How to write model file for JAGS binomial using logit function

Question

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:

generate 10000 samples from the posterior for the two parameters alpha/beta
produce samples to from the posterior predicted number of success when dist = 25 for 100 attempts
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.

Thanks for helping with this yes it was the fact that I was passing in precalulated p values that was causing the error, it was overspecified. — user3159273

user20650 user20650 · Accepted Answer · 2019-04-25T12:40:26

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

How to write model file for JAGS binomial using logit function

1 Answers