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

Question

我正在使用 JAGS 对二项分布建模，其中 p 参数是另一个变量 d.

的函数

这就是我想要做的：

为两个参数从后验生成 10000 个样本alpha/beta
当 dist = 25 进行 100 次尝试时，根据后验预测的成功次数生成样本
计算 25 英尺距离成功率的 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)
}

我收到一个错误

Error in jags.model("glm1.model", glm1.data, n.chains = 2) :
RUNTIME ERROR:
Compilation error on line 4.
pmiss[1] is a logical node and cannot be observed

我认为模型文件甚至没有设置为执行我正在尝试执行的操作。

Answer 1

您不需要计算 rjags 之外的概率，但可以使用二项分布函数 dbin(p,N)，它采用参数 p、成功概率，以及N，尝试次数。此外，logit 函数可以用作 link 函数。

更新后的模型函数为

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))）

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

How to write model file for JAGS binomial using logit function

r

bayesian

montecarlo

jags