2017-08-14 128 views
1

我想用rjagsrjags错误无效向量参数ilogit

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同样回归到一个比较回归betareg主场迎战同回归

#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)如何设置先验的测试回归

谢谢。

回答

1

发生此错误是因为您已将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,bphi)。