如何使用logit函数为JAGS二项式编写模型文件

Question

我正在使用JAGS对二项分布进行建模，该分布的p参数是另一个变量d的函数。

这就是我想要做的：

1. 从后面生成10000个样本，用于两个参数alpha / beta
2. 当dist = 25进行100次尝试时，从后验预测的成功次数产生样本
3. 计算95英尺距离成功率的95可信区间

我写了模型，但它给出了一个错误。

以下是我已经尝试过的代码

#R-code
distance=seq(from=2,to=20,by=1)
Ntrys=c(1443,694,455,353,272,256,240,217,200,237,202,192,174,167,201,195,191,147,152)
Nsucc=c(1346,577,337,208,149,136,111,69,67,75,52,46,54,28,27,31,33,20,24)

psucc=Nsucc/Ntrys

glm1.data=list(N=19, Nsucc=Nsucc,psucc=psucc,distance=distance)

glm1.model=jags.model("glm1.model",glm1.data,n.chains=2)

glm1.samps=coda.samples(glm1.model, variable.names=c("alpha", "beta"), 1e5)

#model file
model{ 
    for (i in 1:N){
            Nsucc[i] ~ dbern(psucc[i])
            log((psucc[i])/(1-psucc[i])) <- alpha + beta*(distance[i])
    }
    alpha ~ dunif(-10,10)
    beta ~ dunif(-10,10)
}

我收到一个错误

jags.model（“glm1.model”，glm1.data，n.chains = 2）出错： 运行时错误： 第4行的编译错误。 pmiss [1]是一个逻辑节点，无法观察到

我不认为模型文件甚至设置为我正在尝试做的事情。

Answer 1

你不需要计算rjags之外的概率，但可以使用二项式分布函数dbin(p,N)，它取参数，p，成功概率，N，尝试次数。另外，logit函数可以用作链接函数。

然后是更新的模型函数

mod <-
"model{ 
    # likelihood
    for (i in 1:N){
            Nsucc[i] ~ dbin(p[i], Ntrys[i])
            logit(p[i]) <- alpha + beta*distance[i]
    }
    # priors
    alpha ~ dunif(-10,10)
    beta ~ dunif(-10,10)

}"

通过将预测变量的值与数据相加，并将相关数量的NA附加到结果向量，可以在给定预测变量值的情况下生成预测。所以传递给rjags的数据变成了

glm1.data <- list(N=20, Nsucc=c(Nsucc, NA), Ntrys=c(Ntrys, 100), distance=c(distance, 25))

然后编译并运行模型

# set.seed so sampling is reproducible
library(rjags)
load.module("glm")

glm1.model <- jags.model(textConnection(mod), glm1.data, 
                         n.chains=2,
                         inits=list(.RNG.name="base::Wichmann-Hill",
                                    .RNG.seed=1))
update(glm1.model, n.iter = 1000, progress.bar="none")

# sample: monitor the unknown predictions, Nsucc[20], p[20]
glm1.samps <- coda.samples(glm1.model, variable.names=c("alpha", "beta", "Nsucc[20]", "p[20]"), 1e5)

然后，您可以从分位数生成间隔

s <- summary(glm1.samps)
s$quantiles 

或最高密度间隔

library(HDInterval)
hdi(glm1.samps)

（只是为了好玩，比较来自glm的系数：summary(glm(cbind(Nsucc, Ntrys-Nsucc) ~ distance, family=binomial))）