[R-sig-ME] glm module in

[Martyn Plummer]

Hi Ben,

There are two reasons that your model is not recognized by the the glm
module. I'll repeat the essential parts of the model:

SibNeg[i] ~ dpois(mu[i])
mu[i] <- lambda[i] * z[i] + 0.00001
log(lambda[i]) <- offset[i] + alpha + inprod(X[i,],beta) + a[nest[i]]

Firstly, this is not a log-link GLM due to the perturbation you
have added to mu[i].  Secondly, and this is the real problem, even
without the perturbation JAGS still would not recognize it. 

Normally you can work around this problem by adding log(z) as an offset,
but in your case this will not work since z may be zero.

I am aware of this issue: it is a commonly occurring set-up in
epidemiological applications, where lambda would a be in incidence rate
and z would be the person-time at risk. I'll put this on the TODO list.

Martyn


On Mon, 2011-09-19 at 16:53 -0400, Ben Bolker wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
> 
> 
>   I'm curious whether anyone knows whether, or how, the following model
> can be made to use the 'glm' module in JAGS ... I followed the advice in
> <https://stat.ethz.ch/pipermail/r-sig-mixed-models/2011q1/005780.html>
> and looked at list.samplers() on the result, with the following outcome:
> 
> > unique(names(list.samplers(ooo$model)))
> [1] "FiniteMethod"  "RealSlicer"    "ConjugateBeta"
> 
>   i.e. no glm module.
> 
>   I'm using
> [1] R2jags_0.02-15   rjags_2.2.0-4    R2WinBUGS_2.1-18 coda_0.14-4
> 
>  which uses JAGS 2.2.0 (I haven't spent the time yet to figure out how
> to upgrade), but a glance at the manual suggests that the glm module
> hasn't changed radically since then ...
> 
>   For background, the model is a mixed zero-inflated Poisson, with the
> model matrix constructed in R and passed into JAGS in the data.
> 
>   If I have time I will work to break this down to a simpler example,
> but in the meantime if anyone happens to have any bright ideas ...
> 
>   sincerely
>     Ben Bolker
> ============
> 
> model {
> 
>   ## PRIORS
>   for (m in 1:5){
>     beta[m] ~ dnorm(0, 0.01)           # Linear effects
>   }
>   alpha ~ dnorm(0, .01)
>   sigma ~ dunif(0, 5)
>   tau <- 1/(sigma*sigma)
>   psi ~ dunif(0, 1)
>   for(j in 1:nnests){
>      a[j] ~ dnorm(0, tau)
>    }
> 
>    for(i in 1:N){
>        SibNeg[i] ~ dpois(mu[i])
>        mu[i] <- lambda[i]*z[i]+0.00001
>      ## hack required for Rjags -- otherwise 'incompatible'
>        z[i] ~ dbern(psi)
>        log(lambda[i]) <- offset[i] + alpha +
>               inprod(X[i,],beta) + a[nest[i]]
>    }
> }



-----------------------------------------------------------------------
This message and its attachments are strictly confidenti...{{dropped:8}}