我想比较 betareg 回归与使用 rjags 的相同回归
library(betareg)
d = data.frame(p= sample(c(.1,.2,.3,.4),100, replace= TRUE),
id = seq(1,100,1))
# I am looking to reproduce this regression with jags
b=betareg(p ~ id, data= d,
link = c("logit"), link.phi = NULL, type = c("ML"))
summary(b)
下面我正在尝试使用 rjags 进行相同的回归
#install.packages("rjags")
library(rjags)
jags_str = "
model {
#model
y ~ dbeta(alpha, beta)
alpha <- mu * phi
beta <- (1-mu) * phi
logit(mu) <- a + b*id
#priors
a ~ dnorm(0, .5)
b ~ dnorm(0, .5)
t0 ~ dnorm(0, .5)
phi <- exp(t0)
}"
id = d$id
y = d$p
model <- jags.model(textConnection(jags_str),
data = list(y=y,id=id)
)
update(model, 10000, progress.bar="none"); # Burnin for 10000 samples
samp <- coda.samples(model,
variable.names=c("mu"),
n.iter=20000, progress.bar="none")
summary(samp)
plot(samp)
我在这一行上收到错误
model <- jags.model(textConnection(jags_str),
data = list(y=y,id=id)
)
Error in jags.model(textConnection(jags_str), data = list(y = y, id = id)) :
RUNTIME ERROR:
Invalid vector argument to ilogit
你能建议一下吗
(1(如何修复错误
(2( 如何设置贝塔回归的先验
谢谢。
发生此错误是因为您已将id
向量提供给标量函数logit
。在 Jags 中,逆链接函数无法矢量化。要解决此问题,您需要使用for
循环来遍历id
的每个元素。为此,我可能会在数据列表中添加一个额外的元素,表示id
有多长。
d = data.frame(p= sample(c(.1,.2,.3,.4),100, replace= TRUE),
id = seq(1,100,1), len_id = length(seq(1,100,1)))
从那里,您只需对jags代码进行少量编辑即可。
for(i in 1:(len_id)){
y[i] ~ dbeta(alpha[i], beta[i])
alpha[i] <- mu[i] * phi
beta[i] <- (1-mu[i]) * phi
logit(mu[i]) <- a + b*id[i]
}
但是,如果您跟踪mu
它将是一个矩阵,该矩阵为 20000(# 次迭代(乘以 100(id
的长度(。您可能对实际参数(a
、b
和phi
(更感兴趣。