[R-sig-ME] glm module in

[Ben Bolker]

On 09/20/2011 05:34 AM, Martyn Plummer wrote:
> 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

  Thanks, Martyn, for the quick response.  Nice to know I wasn't missing
something too obvious.  It's not critical for this application (I can
afford to let the model run for 10 minutes), but I wonder if a
workaround would be to use something like

   logzzvals <- c(-10,0)
   lzcat[i] ~ dbern(prob)
   log(lambda[i]) <- offset[i] + logzvals[lzcat[i]+1] + ...

?



> 
> 
> 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 confiden...{{dropped:8}}