10

I am trying to fit a logistic regression model in JAGS, but I have data in the form of (# success y, # attempts n), rather than a binary variable. In R, one can fit a model to data such as these by using glm(y/n ~ ) with the "weights" argument, but I am not sure how to fit this in JAGS.

Here is a simple example that I hope addresses what I am trying to ask. Note that I am using the rjags package. Thanks for any help!

y <- rbinom(10, 500, 0.2)
n <- sample(500:600, 10)
p <- y/n
x <- sample(0:100, 10) # some covariate

data <- data.frame(y, n, p, x)

model <- "model{
# Specify likelihood
for(i in 1:10){
    y[i] ~ dbin(p[i], n[i])
    logit(p[i]) <- b0 + b1*x
}

# Specify priors
b0 ~ dnorm(0, 0.0001)
b1 ~ dnorm(0, 0.0001)
}"
Kirk Fogg
  • 521
  • 5
  • 14
  • Your model is wrapped in quotation marks. I'm not familiar with RJags, but this looks incorrect to me. – Phil Apr 30 '15 at 21:00
  • @Phil, BUGS/JAGS models are sometimes specified that way (they would then need to be written to a temporary file) – Ben Bolker Apr 30 '15 at 21:10
  • That's precisely why I thought I'd flag it rather than dive in the deep end and edit it! Glad you got a solution. – Phil Apr 30 '15 at 21:30

1 Answers1

10

You don't need to compute p in your data set at all. Just let it be a logical node in your model. I prefer the R2jags interface, which allows you to specify a BUGS model in the form of an R function ...

jagsdata <- data.frame(y=rbinom(10, 500, 0.2),
                   n=sample(500:600, 10),
                   x=sample(0:100, 10))
model <- function() {
    ## Specify likelihood
    for(i in 1:10){
        y[i] ~ dbin(p[i], n[i])
        logit(p[i]) <- b0 + b1*x[i]
    }
    ## Specify priors
    b0 ~ dnorm(0, 0.0001)
    b1 ~ dnorm(0, 0.0001)
}

Now run it:

library("R2jags") 
jags(model.file=model,data=jagsdata,
     parameters.to.save=c("b0","b1"))
jbaums
  • 27,115
  • 5
  • 79
  • 119
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453