Modèles de régression Bayésiens avec OpenBUGS

Régression linéaire avec des erreurs normales
Diagnostoque de la convergence package coda
Exercice: A comparer avec un autre modèle Bayésien

Régression linéaire avec des erreurs normales

Dans cette section on considère les données contenant
- \(y_i,\;i=1,\ldots,n\) une variable réponse
- \(x_{i1},\ldots,x_{ip}\), \(i=1,\ldots,n\) \(p\) variables explicatives.
On considère le modèle suivant : pour tout \(i=1,\ldots,n\) \[y_i\sim\mbox{Normal}(\mu_i,\sigma^2),\;\;\mu_i=\beta_0+\displaystyle\sum_{k=1}^p \beta_k x_{ip}\] où \(\beta_0\), \(\beta_1,\ldots,\beta_p\) et \(\sigma\) sont les paramètres inconnus.
Si on veut une solution bayésienne au probléme de régression avec une loi a priori non informative on peut considérer les lois a priori suivantes \[\beta_k\sim\mbox{Normal}(0,10^5)\;\mbox{ et }\;\log \sigma \sim\mbox{Unif}(-100,100)\]
Dans une regression bayésienne il est preférable de centrer toutes les variables expliquatives. Ceci pourrait réduire les corrélation entre les coefficients \((\beta_1,\ldots,\beta_p)\) et \(\beta_0\)

Exemple : Courbe de croissance

Dans Glefand et al. (1990) on examine la croissance de 5 jeunes rats où on a mesuréd’une facon hebdomadaire leurs poids pendant 5 semaines.
Dans cet exemple ajuste un modèle de régression linéaire pour le neuvième rat.

Le modèle Bayésien

model{
for(i in 1:5){
y[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha + beta * (x[i]-mean(x[])
  }
# Prior
  alpha ~ dnorm(0,00001)
  beta ~ dnorm(0,.00001)
  r ~ dunif(-100,100)
  sigma2 <- exp(r)
  tau <- 1/sigma2
}

Le code `R`

Construction du fichier

> library(R2OpenBUGS)
> sink('exempleRegLineaire.txt')
> cat("
+     model{
+ for(i in 1:5){
+ y[i] ~ dnorm(mu[i],tau)
+ mu[i] <- alpha + beta * (x[i]-mean(x[]))
+   }
+ # Prior
+   alpha ~ dnorm(0,.00001)
+   beta ~ dnorm(0,.00001)
+   r ~ dunif(-100,100)
+   sigma2 <- exp(r/2)
+   tau <- 1/sigma2
+ }
+ ",fill=T)
> sink()

> filename<-"exempleRegLineaire.txt"

Les données

> donnees=list(y=c(177,236,285,350,376),x=c(8,15,22,29,36))

Les valeurs initiales

> inits<-function(){
+   inits=list(alpha=c(.5,3,30),beta=c(-2,10,-5),r=c(-90,0,50))
+ }

Les paramètres

> params <- c('alpha','beta','sigma2')

Faire tourner OpenBUGS à partir de R

> outBinom <-bugs(donnees,inits,params,filename,codaPkg=F,n.thin =1,
+                 n.iter=10000,debug=F,n.chains = 3,
+                 working.directory=getwd())

Résumé de la loi a posteriori

> outBinom
Inference for Bugs model at "exempleRegLineaire.txt", 
Current: 3 chains, each with 10000 iterations (first 5000 discarded)
Cumulative: n.sims = 15000 iterations saved
          mean    sd  2.5%   25%   50%   75%  97.5% Rhat n.eff
alpha    284.6   8.1 268.4 280.9 284.7 288.3  300.1    1 15000
beta       7.3   0.8   5.8   6.9   7.3   7.7    8.8    1 15000
sigma2   315.8 923.8  38.4  86.2 149.3 295.5 1577.0    1 13000
deviance  39.9   3.6  35.8  37.2  39.0  41.6   49.2    1  7000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = Dbar-Dhat)
pD = 2.9 and DIC = 42.8
DIC is an estimate of expected predictive error (lower deviance is better).

Estimation des lois a posteriori (avec ggplot2)

> library(ggplot2)
> dt=data.frame(outBinom$sims.matrix)
> dt=dt[-c(1:13000),]
> gr<-ggplot(data=dt,aes(x=alpha))+geom_density(alpha=.5,colour='red',fill='red')+xlab(expression(alpha))+ylab("density")
> gr

> gr<-ggplot(data=dt,aes(x=beta))+geom_density(alpha=.5,colour='lightblue',fill='lightblue')+xlab(expression(beta))+ylab("density")
> gr

> gr<-ggplot(data=dt,aes(x=sigma2))+geom_density(alpha=.5,colour='green',fill='green')+xlab(expression(sigma^2))+ylab("density")
> gr

Diagnostoque de la convergence package `coda`

Diagnostique visuel

Si \((\theta_i)_{i=1,\ldots,M}\) est la chaîne de Markov simulée.
Le traceplot est une représentation graphique de \((i,\theta_i)\).
Nous pouvons voir si notre chaîne se coince dans une certaine région de l’espace des paramètres

> library(coda)
> dim(outBinom$sims.array)
[1] 5000    3    4
> dim(outBinom$sims.array[,,1])
[1] 5000    3
> dimnames(outBinom$sims.array)
[[1]]
NULL

[[2]]
NULL

[[3]]
[1] "alpha"    "beta"     "sigma2"   "deviance"
> alpha.mcmc=mcmc(outBinom$sims.array[,,1])
> matplot(alpha.mcmc,col=c("red","blue","green"),type="l",ylab=expression(alpha),xlab="Iterations")

> beta.mcmc=mcmc(outBinom$sims.array[,,2])
> matplot(beta.mcmc,col=c("red","blue","green"),type="l",ylab=expression(beta),xlab="Iterations")

ou avec traceplot

> traceplot(list(alpha.mcmc[,1],alpha.mcmc[,2],alpha.mcmc[,3]),col=c("red","blue","green"))

Le résumé d’un objet `mcmc`

> alpha.mcmc=mcmc(outBinom$sims.array[,,1])
> summary(alpha.mcmc)

Iterations = 1:5000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

      Mean    SD Naive SE Time-series SE
[1,] 284.6 8.376   0.1185         0.1205
[2,] 284.5 7.969   0.1127         0.1035
[3,] 284.5 8.062   0.1140         0.1191

2. Quantiles for each variable:

      2.5%   25%   50%   75% 97.5%
var1 268.3 281.1 284.7 288.3 299.8
var2 268.1 280.7 284.5 288.3 300.8
var3 269.0 280.9 284.7 288.3 299.8

A comparer avec

La moyenne

> mean(outBinom$sims.array[,1,1])
[1] 284.641

Standard-deviation

> sd(outBinom$sims.array[,1,1])
[1] 8.375994

Naive standard-error.

> sd(outBinom$sims.array[,1,1])/sqrt(nrow(alpha.mcmc))
[1] 0.1184544

L’autocorrélation

C’est une corrélation dans la CM avec le \(k-\)ième décalage \item L’autocorrélation \(\rho_k\) s’exprime \[\rho_k=\displaystyle\frac{\sum_{i=1}^{n-k}(\theta_i-\overline{\theta})(\theta_{i+k}-\overline{\theta})}{\sum_{i=1}^n(\theta_i-\overline{\theta})^2}\]
Sous R

> autocorr(alpha.mcmc)
, , 1

                [,1]         [,2]         [,3]
Lag 0   1.0000000000  0.012698341 -0.007371685
Lag 1  -0.0066858529  0.013116688  0.022311364
Lag 5   0.0059180225  0.033834899 -0.010044573
Lag 10 -0.0135458633 -0.001377333 -0.010339931
Lag 50 -0.0009784652 -0.016211932  0.013690217

, , 2

               [,1]         [,2]         [,3]
Lag 0   0.012698341  1.000000000 -0.007867831
Lag 1   0.034562726 -0.017593199  0.006079726
Lag 5  -0.009497296  0.006777567 -0.013140899
Lag 10 -0.023438346 -0.014230336  0.014795043
Lag 50  0.010218417 -0.008884439  0.023339124

, , 3

               [,1]         [,2]         [,3]
Lag 0  -0.007371685 -0.007867831  1.000000000
Lag 1  -0.012619016  0.004496634 -0.020315248
Lag 5   0.014872212 -0.008786492 -0.011774687
Lag 10  0.001732580 -0.017331261 -0.016412201
Lag 50  0.008282702 -0.001656412 -0.006702233
> autocorr.plot(alpha.mcmc,lag.max = 10)

Si on change le paramètre thin

> alpha.mcmc=mcmc(outBinom$sims.array[,,1],thin=10)
> autocorr(alpha.mcmc)
, , 1

                 [,1]         [,2]         [,3]
Lag 0    1.0000000000  0.012698341 -0.007371685
Lag 10  -0.0066858529  0.013116688  0.022311364
Lag 50   0.0059180225  0.033834899 -0.010044573
Lag 100 -0.0135458633 -0.001377333 -0.010339931
Lag 500 -0.0009784652 -0.016211932  0.013690217

, , 2

                [,1]         [,2]         [,3]
Lag 0    0.012698341  1.000000000 -0.007867831
Lag 10   0.034562726 -0.017593199  0.006079726
Lag 50  -0.009497296  0.006777567 -0.013140899
Lag 100 -0.023438346 -0.014230336  0.014795043
Lag 500  0.010218417 -0.008884439  0.023339124

, , 3

                [,1]         [,2]         [,3]
Lag 0   -0.007371685 -0.007867831  1.000000000
Lag 10  -0.012619016  0.004496634 -0.020315248
Lag 50   0.014872212 -0.008786492 -0.011774687
Lag 100  0.001732580 -0.017331261 -0.016412201
Lag 500  0.008282702 -0.001656412 -0.006702233
> autocorr.plot(alpha.mcmc,lag.max = 10)

Le facteur de réduction `Rhat`

Le facteur de réduction d’échelle s’exprime \[\widehat{R}=\sqrt{\displaystyle\frac{\widehat{\mbox{Var}(\theta)}}{W}}\] où \[\widehat{\mbox{Var}(\theta)}=\left(1-\displaystyle\frac{1}{n}\right) W+\displaystyle\frac{1}{n} B,\]

\[W=\displaystyle\frac{1}{m}\sum_{j=1}^{m} \displaystyle\frac{1}{n-1}\sum_{i=1}^n(\theta_{ij}-\overline{\theta}_j)^2,\] \[B=\displaystyle\frac{n}{m-1}\sum_{j=1}^m (\overline{\theta}_j-\overline{\overline{\theta}})^2\] et \[\overline{\overline{\theta}}=\displaystyle\frac{1}{m}\sum_{j=1}^m\overline{\theta}_j\] \(n\) étant le nombre d’itérations et \(m\) est le nombre chaînes générées.

Interprétation :
- Si \(\widehat{R}\) est grand (plus grand que 1,1) on devrait faire tourner la chaîne un peu plus longtemps jusqu’à atteindre la convergence.
- Si on a plusieurs parametres à estimer, il faut alors calculer un facteur pour chaque paramètre.
- Il faut aussi faire tourner la chaîne un peu plus longtemps pour que \(\widehat{R}\) devient de plus en petit.
- Une valeur de \(\widehat{R}\) proche de 1 indique que la CM a convergé.

> gelman.diag(list(alpha.mcmc[,1],alpha.mcmc[,2],alpha.mcmc[,3]))
Potential scale reduction factors:

     Point est. Upper C.I.
[1,]          1          1

> gelman.plot(list(alpha.mcmc[,1],alpha.mcmc[,2],alpha.mcmc[,3]))

Diagnostique de Geweke

Cette méthode est basée sur un test d’égalité des moyennes entre une première partie de la chaîne (par défaut les premiers 10\(\%\)) et une seconde partie (les 50\(\%\) dernières simulations de la chaînes).
\(Z-\)score est alors calculé à partir de la différence des moyennes empiriques divisées par l’écart-type empirique.
Sous l’hypothèse nulle (les deux parties de la chaînes suivent la même loi de probabilité a posteriori), \(Z\) suit une loi normale \(\mathcal{N}(0,1)\).

> geweke.diag(alpha.mcmc)

Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5 

   var1    var2    var3 
-1.0221 -0.5755 -0.3331 
> geweke.plot(alpha.mcmc)

> geweke.plot(alpha.mcmc)

Et le graphique

Exercice: A comparer avec un autre modèle Bayésien

Avec des lois a priori informatives.

model{
for(i in 1:5){
y[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha + beta * (x[i]-mean(x[])
  }
# Prior
  alpha ~ dnorm(0,001)
  beta ~ dnorm(0,.001)
  r ~ dunif(-10,10)
  sigma2 <- exp(r)
  tau <- 1/sigma2
}

> outBinom2
Inference for Bugs model at "exempleRegLineaire.txt", 
Current: 3 chains, each with 10000 iterations (first 5000 discarded)
Cumulative: n.sims = 15000 iterations saved
          mean   sd  2.5%   25%   50%   75% 97.5% Rhat n.eff
alpha    220.3  4.8 211.0 217.0 220.3 223.5 229.7    1  2600
beta       7.3  0.5   6.2   6.9   7.3   7.7   8.4    1 15000
sigma2   146.3  2.1 140.8 145.5 147.0 147.8 148.4    1  5100
deviance 180.5 21.0 141.3 165.9 179.8 194.5 223.5    1  2700

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = Dbar-Dhat)
pD = 1.7 and DIC = 182.2
DIC is an estimate of expected predictive error (lower deviance is better).

Modèles de régression Bayésiens avec OpenBUGS

Dhafer Malouche

Novembre, 2015

Régression linéaire avec des erreurs normales

Exemple : Courbe de croissance

Le modèle Bayésien

Le code `R`

Diagnostoque de la convergence package `coda`

Diagnostique visuel

Le résumé d’un objet `mcmc`

L’autocorrélation

Le facteur de réduction `Rhat`

Diagnostique de Geweke

Exercice: A comparer avec un autre modèle Bayésien

Modèles de régression Bayésiens avec OpenBUGS

Dhafer Malouche

Novembre, 2015

Régression linéaire avec des erreurs normales

Exemple : Courbe de croissance

Le modèle Bayésien

Le code R

Diagnostoque de la convergence package coda

Diagnostique visuel

Le résumé d’un objet mcmc

L’autocorrélation

Le facteur de réduction Rhat

Diagnostique de Geweke

Exercice: A comparer avec un autre modèle Bayésien

Le code `R`

Diagnostoque de la convergence package `coda`

Le résumé d’un objet `mcmc`

Le facteur de réduction `Rhat`