gpt4 book ai didi

r - 使用 rjags 包进行贝叶斯多项式回归

转载 作者:行者123 更新时间:2023-12-03 23:44:54 24 4
gpt4 key购买 nike

我正在尝试使用 rjags 拟合多项逻辑回归模型因为结果是一个具有 3 个水平的分类(名义)变量( 结果 ),解释变量是 年龄 (连续)和 (分类为 3 个级别)。这样做时,我想获得 的后验均值和 95% 基于分位数的区域。年龄 .
我不是很擅长 for loop我认为这就是我为模型编写的代码无法正常工作的原因。
我的 beta 先验遵循正态分布,βj ∼ Normal(0,100) for j ∈ {0, 1, 2}。
可重现的 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
)

错误信息 :
enter image description here
我一直在寻找帮助的来源 :
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 页
enter image description here
我刚刚开始学习如何使用 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"),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参数属于所考虑的每个变量( ageGroup 的水平)。由于混合链,最终值可能会改变!

关于r - 使用 rjags 包进行贝叶斯多项式回归,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/63444427/

24 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com