CODEWITHSUNDEEP
rbayesianmontecarlojags
user3159273
user3159273

Reputation: 23

How to write model file for JAGS binomial using logit function

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.

Upvotes: 2

Views: 4242

Answers (1)

user20650
user20650

Reputation: 25844

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

Upvotes: 2

Related Questions