Regarding Bayesian linear mixed effects models using R

Question

I have fitted following Bayesian linear mixed effects model using R. The analysis is based on sleepstudy dataset in lme4 package and the Bayesian modeling has done using Rjags package.

Where $\beta$ corresponds to fixed effects and b's corresponds to subject specific random effects.

I have fitted this model using two methods and I want to know whether both methods are correct. The main difference in two methods is how the Reaction variable has been used. In the first method , the reaction variable has reshaped so that it is a matrix. But in the second method I have used the reaction variable as it is(as a vector).

Method1

require(lme4)
require(rjags)
sleepstudy1=sleepstudy
M <- 10 
N <- length(unique(sleepstudy1$Subject)) 
Days <-c(0:9)

Reaction<- matrix(as.numeric(sleepstudy1$Reaction),N,10,byrow=TRUE)
model_string.1 <- "model {
for (i in 1:N) {
#
b0[i,1] <- 0
b0[i,2] <- 0
b[i,1:2] ~ dmnorm(b0[i,1:2],ISigma[,])


for (j in 1:M) {
## Response model ###
Reaction[i,j] ~ dnorm(mu[i,j], tau)
mu[i,j] <- beta[1]+b[i,1] + (beta[2]+b[i,2])*Days[j]
}

}

for (l in 1:2) { beta[l] ~dnorm(0, 0.00001) }
#(2) Precision parameters
tau ~dgamma(.01,.01)
sigma.tau <- 1/tau
#(3) Variance-covariance matrix
ISigma[1:2,1:2] ~dwish(R[,],3)
Sigma[1:2,1:2] <- inverse(ISigma[,])
R[1,1] <-1
R[1,2] <-0
R[2,2] <-1
R[2,1] <-0 
}"

model3 <- jags.model(textConnection(model_string.1), 
                     data = list(Reaction=Reaction,N=N,M=M,Days=Days),
                     n.chains=2)
params <- c('b','beta','Sigma','sigma.tau')
samps.2 <- coda.samples(model3, params, n.iter = 10000)


summary.model.2=summary(window(samps.2, start = burn.in))

Stat.model.2=as.data.frame(summary.model.2$statistics)

Method2

N <- length((sleepstudy1$Subject)) 

model_string.1 <- "model {
for (i in 1:N) {

b0[i,1] <- 0
b0[i,2] <- 0
b[i,1:2] ~ dmnorm(b0[i,1:2],ISigma[,])


Reaction[i] ~ dnorm(mu[i], tau)
mu[i] <- beta[1]+b[Subject[i],1] + (beta[2]+b[Subject[i],2])*Days[i]
}
for (l in 1:2) { beta[l] ~dnorm(0, 0.00001) }
#(2) Precision parameters
tau ~dgamma(.01,.01)
sigma.tau <- 1/tau
ISigma[1:2,1:2] ~dwish(R[,],3)
Sigma[1:2,1:2] <- inverse(ISigma[,])
R[1,1] <-1
R[1,2] <-0
R[2,2] <-1
R[2,1] <-0 
}"

model2 <- jags.model(textConnection(model_string.1), 
                     data = list(Reaction=sleepstudy1$Reaction,N=N,Days=sleepstudy$Days,Subject=sleepstudy$Subject),
                     n.chains=2)
params <- c('b','beta','Sigma','sigma.tau')
samps.1 <- coda.samples(model2, params, n.iter = 10000)


summary.model.1=summary(window(samps.1, start = burn.in))

Stat.model.1=as.data.frame(summary.model.1$statistics)

When i compare the results especially the fixed effects, the estimates are more or less same. But the estimates for variances,covariances and for random effects are different. So I want to know whether both these two methods are correct in terms of coding. Especially whether the method 2 is correct. Because method 2 has less number of loops.

Your expert knowledge is highly appreciated.

Thanks.

Regarding Bayesian linear mixed effects models using R

0 Answers