如何使用logit函数编写JAGS二项式模型文件 [英] How to write model file for JAGS binomial using logit function

问题描述

我正在使用JAGS进行分配，以建模二项式分布，其 p 参数是另一个变量 d 的函数.

这就是我想要做的:

1. 从后验中为两个参数alpha/beta生成10000个样本
2. 在dist = 25的情况下进行100次尝试，从后验预测成功次数产生样本
3. 计算25英尺距离处的成功率的95个可信区间

我已经编写了模型，但是出现了错误.

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

 #R-code距离= seq(从= 2，到= 20，乘= 1)Ntrys = c(1443,694,455,353,272,256,240,217,200,237,202,192,174,167,201,195,191,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/Ntrysglm1.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文件模型{为(1:N中的i){Nsucc [i]〜dbern(psucc [i])log((psucc [i])/(1-psucc [i]))<-alpha + beta *(距离[i])}alpha〜dunif(-10,10)Beta〜邓尼夫(-10,10)} 

我收到错误

jags.model("glm1.model"，glm1.data，n.chains = 2)中的错误:
运行时错误:
第4行出现编译错误.
pmiss [1]是一个逻辑节点，无法观察到

我认为模型文件甚至都没有设置为执行我想做的事情.

解决方案

您无需计算 rjags 之外的概率，但可以使用二项式分布函数 dbin(p，N)，其中使用参数 p (成功概率)和 N (尝试次数).另外， logit 函数可以用作链接函数.

更新后的模型功能便是

  mod<-模型{#可能性为(1:N中的i){Nsucc [i]〜dbin(p [i]，Ntrys [i])logit(p [i])<-alpha + beta * distance [i]}#先验阿尔法〜邓尼夫(-10,10)Beta〜邓尼夫(-10,10)}" 

通过将预测变量的值添加到数据，并将相关数量的 NA 附加到结果向量，可以在给定某些预测变量值的情况下生成

预测.因此，传递给 rjags 的数据变为

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

然后编译并运行模型

 #set.seed，因此可重现采样图书馆(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))更新(glm1.model，n.iter = 1000，progress.bar ="none")#样本:监视未知预测Nsucc [20]，p [20]glm1.samps<-coda.samples(glm1.model，variable.names = c("alpha"，"beta"，"Nsucc [20]"，"p [20]")，1e5) 

然后您可以根据分位数生成间隔

  s<-摘要(glm1.samps)s $分位数 

或最高密度间隔

  library(HDInterval)hdi(glm1.samps) 

(只是为了好玩，比较 glm 中的系数: summary(glm(cbind(Nsucc，Ntrys-Nsucc)〜距离，family = binomial)))

I am working on an assignment using JAGS to model a binomial distribution who's p parameter is a function of another variable d.

This is what I am trying to do:

1. generate 10000 samples from the posterior for the two parameters alpha/beta
2. produce samples to from the posterior predicted number of success when dist = 25 for 100 attempts
3. calculate 95 credible interval for success rate at 25 feet distance

I have written the model but it is giving an error.

Below is the code I have already tried

#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)
}

I get an error

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

I don't think the model file is even setup to do what I'm trying to do.

解决方案

You do not need to calculate the probabilities outside of rjags but can use the binomial distribution function, dbin(p,N) which takes the arguments, p, the probability of success, and N, the number of tries. Additionally, the logit function can be used as the link function.

The updated model function is then

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)

}"

Predictions can be generated given some value of the predictors by adding the values of the predictors to the data, and appending the relevant number of NA's to the outcome vector. So the data passed to rjags becomes

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

Then compile and run the model

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

You can then generate intervals from the quantiles

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

or the highest density interval

library(HDInterval)
hdi(glm1.samps)

(just for fun, compare coefficients from glm: summary(glm(cbind(Nsucc, Ntrys-Nsucc) ~ distance, family=binomial)))

这篇关于如何使用logit函数编写JAGS二项式模型文件的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！