如何解决使用rjags包的贝叶斯多项式回归
我正在尝试使用rjags
来拟合多项式逻辑回归模型,因为结果是具有3个级别的分类(名义)变量(结果),解释变量为年龄(连续)和组(分为3个级别)。这样,我想获得年龄和 Group 的后验均值和95%基于分位数的区域。
我对for loop
并不擅长,这是我为模型编写的代码无法正常工作的原因。
我的beta先验遵循正态分布,对于j∈{0,1,2},βj〜正态(0,100)。
可复制的R代码
library(rjags)
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)))
X <- as.matrix(data[,c("Age","Group")])
J <- ncol(X)
N <- nrow(X)
## Step 1: Specify model
cat("
model {
for (i in 1:N){
##Sampling model
yvec[i] ~ dmulti(p[i,1:J],1)
#yvec[i] ~ dcat(p[i,1:J]) # alternative
for (j in 1:J){
log(q[i,j]) <- beta0 + beta1*X[i,1] + beta2*X[i,2]
p[i,j] <- q[i,j]/sum(q[i,1:J])
}
##Priors
beta0 ~ dnorm(0,0.001)
beta1 ~ dnorm(0,0.001)
beta2 ~ dnorm(0,0.001)
}
}",file="model.txt")
##Step 2: Specify data list
dat.list <- list(yvec = data$Outcome,X=X,J=J,N=N)
## Step 3: Compile and adapt model in JAGS
jagsModel<-jags.model(file = "model.txt",data = dat.list,n.chains = 3,n.adapt = 3000
)
错误消息:
我一直在寻求帮助的来源:
http://people.bu.edu/dietze/Bayes2018/Lesson21_GLM.pdf
Dirichlet Multinomial model in JAGS with categorical X
来自http://www.stats.ox.ac.uk/~nicholls/MScMCMC15/jags_user_manual.pdf,第31页的参考
我刚刚开始学习如何使用rjags
软件包,因此,非常感谢任何提示/解释以及与相关资源的链接!
解决方法
我将介绍解决您问题的方法。我采用了与系数相同的先验条件。我只需要提到,由于您在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"),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,quantile,prob=0.05)
b1 b2 b3 b4
-1.45918662 -0.03960765 -0.61027923 -0.42674155
#Upper CI limit
apply(Res,prob=0.95)
b1 b2 b3 b4
-0.13005617 0.04013478 0.72852243 0.61216838
b
参数属于所考虑的每个变量(age
和Group
的级别)。最终值可能会因为混合链而改变!
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。