关于使用R的贝叶斯线性混合效应模型

如何解决关于使用R的贝叶斯线性混合效应模型

我已经使用R拟合了以下贝叶斯线性混合效果模型。该分析基于lme4包中的sleepstudy数据集，而贝叶斯建模已使用Rjags包来完成。

其中$ \ beta $对应于固定效果，b对应于受试对象特定的随机效果。

我已经使用两种方法拟合了该模型，我想知道这两种方法是否正确。 两种方法的主要区别是如何使用反应变量。在第一种方法中，反应变量已重塑为矩阵。但是在第二种方法中，我按原样使用了反应变量（作为向量）。

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)

方法2

N <- length((sleepstudy1$Subject)) 

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

b0[i,])


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,.01)
sigma.tau <- 1/tau
ISigma[1:2,1] <-0 
}"

model2 <- jags.model(textConnection(model_string.1),data = list(Reaction=sleepstudy1$Reaction,Days=sleepstudy$Days,Subject=sleepstudy$Subject),'sigma.tau')
samps.1 <- coda.samples(model2,n.iter = 10000)


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

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

当我比较结果，尤其是固定效果时，估计值大致相同。但是方差，协方差和随机效应的估计是不同的。所以我想知道这两种方法在编码方面是否正确。特别是方法2是否正确。因为方法2的循环次数较少。

我们非常感谢您的专业知识。

谢谢。