关于JAGS for R

问题描述

我在学习贝叶斯算法的同时尝试了rjags r软件包。我使用rrr包中的烟草数据集拟合了以下多元模型。

其中\ gamma是随机截距，用于测量结果之间的相关性。另外，我故意使响应的某些值丢失，以探索数据结构不平衡时的结果。同时，我添加了AID变量来表示主题ID。

我使用rjags对此模型进行了拟合。

library(rrr)
require(dplyr)
library(rjags)
data("tobacco")
tobacco <- as_data_frame(tobacco)

tobacco$AID=seq(1:25)
tobacco[4,1]=NA
tobacco[14,1]=NA
tobacco[8,1]=NA
tobacco[6,2]=NA
tobacco[1,2]=NA
tobacco[19,2]=NA
tobacco[21,2]=NA


N1=length(tobacco$Y1.BurnRate)
bayes_model_mul="model {
for(i in 1:N1){
Y1.BurnRate[i]~dnorm(mu1[i],tau1)
Y2.PercentSugar[i]~dnorm(mu2[i],tau2)
mu1[i]=beta1[1] + beta1[2]*X2.PercentChlorine[i] + beta1[3]*X3.PercentPotassium[i] + gamma[AID[i]]
mu2[i]=beta2[1] + beta2[3]*X2.PercentChlorine[i] + beta2[2]*X1.PercentNitrogen[i]+
beta2[4]*X3.PercentPotassium[i]+gamma[AID[i]]
gamma[i] ~ dmnorm(0,tau_u)

}
for (l in 1:3) { beta1[l] ~dnorm(0,0.001) }
for (l in 1:4) { beta2[l] ~dnorm(0,0.001) }
tau1 ~ dgamma(.01,.01)
sigma_tau1 = 1/tau1 

tau2 ~ dgamma(.01,.01)
sigma_tau2 = 1/tau2

tau_u ~ dgamma(.01,.01)
sigma_tau_u = 1/tau_u 
}"


model3 <- jags.model(textConnection(bayes_model_mul),data = list(Y1.BurnRate=tobacco$Y1.BurnRate,Y2.PercentSugar=tobacco$Y2.PercentSugar,X1.PercentNitrogen=tobacco$X1.PercentNitrogen,N1=N1,X2.PercentChlorine=tobacco$X2.PercentChlorine,X3.PercentPotassium=tobacco$X3.PercentPotassium,AID=tobacco$AID),n.chains=1)
params <- c('beta1','sigma_tau1','sigma_tau2','beta2','sigma_tau_u','gamma')
samps.1 <- coda.samples(model3,params,n.iter = 10000)
burn.in=1000
summary.model.1=summary(window(samps.1,start = burn.in))

我没有收到任何错误。但是我收到以下警告消息。

Warning message:
In FUN(X[[i]],...) : start value not changed

任何人都可以帮助我弄清楚此错误消息是关于什么的吗？

谢谢。

解决方法

window函数用于选择样本范围。它使用start和end参数。如果明确设置，则start所取的值应大于n.adapt中传递给jags.models和n.iter中update的迭代次数的总和大量的老化样品。

您可以查看由jags.model重新调整的对象，以查看迭代次数如何变化。

# Define the model
model3 <- jags.model(textConnection(bayes_model_mul),n.adapt=1234,data = list(Y1.BurnRate=tobacco$Y1.BurnRate,Y2.PercentSugar=tobacco$Y2.PercentSugar,X1.PercentNitrogen=tobacco$X1.PercentNitrogen,N1=N1,X2.PercentChlorine=tobacco$X2.PercentChlorine,X3.PercentPotassium=tobacco$X3.PercentPotassium,AID=tobacco$AID),n.chains=1)
# Look at the number of iterations (== n.adapt)
model3$iter()
#1234

# Add some burnin iterations: iterations are updated to n.adapt + burnin
update(model3,n.iter=1111)
model3$iter()
# 2345

# Draw more samples -> iterations updated to n.adapt + burnin + n.iter
samps.1 <- coda.samples(model3,n.iter = 10000,variable.names=c('beta1','sigma_tau1','sigma_tau2','beta2','sigma_tau_u','gamma'))
model3$iter()
# 12345

当我们使用window时，起始值应大于n.adapt + burnin，因为这些迭代已被丢弃。如果不是，则从start对象中提取一个coda.samples值（请参阅mcpar(samps.1[[1]])），并且不使用您的手动start值。

因此，使用w=window(samps.1,start=2345)会给出您看到的警告

警告信息：
在FUN（X [[i]]，...）中：起始值未更改

但是以下情况还可以，因为start大于n.adapt + burnin

w=window(samps.1,start=2346)

请记住，这些只是警告，在这种情况下并不重要。但是，如果您在不重新运行coda.samples的情况下重新运行jags.model，这将再次更新迭代，这使得跟踪要传递给window的值（然后会引发错误）变得更加困难