使用 rjags 包的贝叶斯多项式回归

Question

我试图拟合一个多项式逻辑回归模型，其rjags结果是一个具有 3 个级别的分类（名义）变量（结果），解释变量是年龄（连续）和组（具有 3 个级别的分类）。在此过程中，我想获得Age和Group的后验均值和基于分位数的 95% 区域。

\n

我不太擅长，for loop我认为这就是我为模型编写的代码无法正常工作的原因。

\n

我的 beta 先验遵循正态分布，\xce\xb2j \xe2\x88\xbc Normal(0,100) for j \xe2\x88\x88 {0, 1, 2}。

\n

可重现的 R 代码

\n

library(rjags)\n\nset.seed(1)\ndata <- data.frame(Age = round(runif(119, min = 1, max = 18)),\n                   Group = c(rep("pink", 20), rep("blue", 18), rep("yellow", 81)), \n                   Outcome = c(rep("A", 45), rep("B", 19), rep("C", 55)))\n\nX <- as.matrix(data[,c("Age", "Group")]) \nJ <- ncol(X)\nN <- nrow(X)\n\n## Step 1: Specify model\ncat("\nmodel {\nfor (i in 1:N){\n\n    ##Sampling model\n    yvec[i] ~ dmulti(p[i,1:J], 1)\n    #yvec[i] ~ dcat(p[i, 1:J])  # alternative\n    for (j in 1:J){\n      log(q[i,j]) <- beta0 + beta1*X[i,1] + beta2*X[i,2] \n      p[i,j] <- q[i,j]/sum(q[i,1:J])  \n    } \n    \n    ##Priors\n    beta0 ~ dnorm(0, 0.001)\n    beta1 ~ dnorm(0, 0.001)\n    beta2 ~ dnorm(0, 0.001)\n}\n}",\nfile="model.txt")\n\n##Step 2: Specify data list \ndat.list <- list(yvec = data$Outcome, X=X, J=J, N=N) \n\n## Step 3: Compile and adapt model in JAGS \njagsModel<-jags.model(file = "model.txt",\n                      data = dat.list,\n                      n.chains = 3,\n                      n.adapt = 3000\n)\n\n

\n

错误信息：

\n

\n

我一直在寻找帮助的来源：

\n

http://people.bu.edu/dietze/Bayes2018/Lesson21_GLM.pdf

\n

JAGS 中带有分类 X 的狄利克雷多项式模型

\n

参考http://www.stats.ox.ac.uk/~nicholls/MScMCMC15/jags_user_manual.pdf ，第31页

\n

\n

我刚刚开始学习如何使用该rjags包，因此任何提示/解释以及相关来源的链接将不胜感激！

\n

Answer 1

我将提供解决您问题的方法。我采用了您为系数定义的相同先验。我只需要提到，由于您有一个因素，Group我将使用其水平之一作为参考（在本例中pink），因此模型中的常量将考虑其影响。接下来是代码：

library(rjags)
#Data
set.seed(1)
data <- data.frame(Age = round(runif(119, min = 1, max = 18)),
                   Group = c(rep("pink", 20), rep("blue", 18), rep("yellow", 81)), 
                   Outcome = c(rep("A", 45), rep("B", 19), rep("C", 55)))

#Input Values we will avoid pink because it is used as reference level
#so constant absorbs the effect of that level
r1 <- as.numeric(data$Group=='pink')
r2 <- as.numeric(data$Group=='blue')
r3 <- as.numeric(data$Group=='yellow')
age <- data$Age
#Output 2 and 3
o1 <- as.numeric(data$Outcome=='A')
o2 <- as.numeric(data$Outcome=='B')
o3 <- as.numeric(data$Outcome=='C')
#Dim, all have the same length
N <- length(r2)

## Step 1: Specify model

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

## outcome levels B, C
o1[i] ~ dbern(pi1[i])
o2[i] ~ dbern(pi2[i]) 
o3[i] ~ dbern(pi3[i]) 

## predictors
logit(pi1[i]) <- b1+b2*age[i]+b3*r2[i]+b4*r3[i]
logit(pi2[i]) <- b1+b2*age[i]+b3*r2[i]+b4*r3[i]
logit(pi3[i]) <- b1+b2*age[i]+b3*r2[i]+b4*r3[i]

} 
## priors
b1 ~ dnorm(0, 0.001)
b2 ~ dnorm(0, 0.001)
b3 ~ dnorm(0, 0.001)
b4 ~ dnorm(0, 0.001)
}
"
#Model
model.spec<-textConnection(model.string)

## fit model w JAGS
jags <- jags.model(model.spec,
                   data = list('r2'=r2,'r3'=r3,
                               'o1'=o1,'o2'=o2,'o3'=o3,
                               'age'=age,'N'=N),
                   n.chains=3,
                   n.adapt=3000)

#Update the model
#Update
update(jags, n.iter=1000,progress.bar = 'none')
#Sampling
results <- coda.samples(jags,variable.names=c("b1","b2","b3","b4"),n.iter=1000,
                        progress.bar = 'none')
#Results
Res <- do.call(rbind.data.frame, results)

使用保存在 中的参数链结果Res，您可以使用以下代码计算后验媒体和可信区间：

#Posterior means
apply(Res,2,mean)

         b1          b2          b3          b4 
-0.79447801  0.00168827  0.07240954  0.08650250

#Lower CI limit
apply(Res,2,quantile,prob=0.05)

         b1          b2          b3          b4 
-1.45918662 -0.03960765 -0.61027923 -0.42674155

#Upper CI limit
apply(Res,2,quantile,prob=0.95)

         b1          b2          b3          b4 
-0.13005617  0.04013478  0.72852243  0.61216838 

参数b属于所考虑的每个变量（age以及 的水平Group）。由于混合链，最终值可能会发生变化！