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

【问题标题】：How to write model file for JAGS binomial using logit function如何使用 logit 函数为 JAGS 二项式编写模型文件
【发布时间】：2019-04-25 07:06:53
【问题描述】：

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

这就是我想要做的：

从后验中为两个参数 alpha/beta 生成 10000 个样本
当 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)
}

我收到一个错误

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

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

【问题讨论】：

感谢您的帮助，是的，事实上我传递了导致错误的预先计算的 p 值，它被过度指定了。
不客气。

标签： r bayesian montecarlo jags

【解决方案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))）

【讨论】：