[R-sig-ME] glmmadmb- problems with explanatory variable with a lot of zeros?

[Ben Bolker]

Mabille, Geraldine <Geraldine.Mabille at ...> writes:


> Edit...I have just understood the "problem" with the residuals that
> were different even though they had the same number of Lamb Killed
> and the same Dens_Lynx. It is only that they come from different
> municipalities!  The rest of the problem, with log likelihood
> behaving weirdly is still their though...

[snip]

> I'm just starting to use glmmadmb to try to model the number of
> lambs killed by lynx (Killed_Lamb) in each Norwegian
> municipality. We have repeated data over 11 years so I use
> municipality as a random effect in the analysis. The data contains a
> lot of zeros (403 over a total of 2151 lines) and evidence for
> overdispersion (tested using the overdisp_fun() after fitting of a
> poisson model in lme4, ratio obtained=23.8). I therefore decided to
> try using a negative binomial distribution with zero inflation.


> I first modeled the number of killed lambs as a function of the
> density of lynx (Dens_Lynx) in the municipality. We also have a lot
> of zero in Dens_Lynx (1518 over a total of 2151 lines).

> Mod1 <-glmmadmb(Killed_Lamb ~I(scale(Dens_Lynx)), random=~ 1|KOM, 
> data=DATA,zeroInflation=TRUE,family="nbinom")
> I obtain the following summary for Mod 1:
> Call:
> glmmadmb(formula = Killed_Lamb ~ I(scale(Dens_Lynx)), data = DATA,
>     family = "nbinom", random = ~1 | KOM, zeroInflation = TRUE)
 

[snip]
 
> Number of observations: total=2151, KOM=257
> Random effect variance(s):
> Group=KOM
>             Variance StdDev
> (Intercept)    4.053  2.013
> 
> Negative binomial dispersion parameter: 1.7665 (std. err.: 0.075737)
> Zero-inflation: 0.090062  (std. err.:  0.0083659 )
> 
> Log-likelihood: -9856.38
 
> I then tried to compare AIC for Mod1, with AIC for a base model
>  containing only an intercept as fixed effect:

> Mod0 <-glmmadmb(Killed_Lamb ~1, random=~ 1|KOM, 
> data=DATA,zeroInflation=TRUE,family="nbinom")
> anova(Mod0,Mod1).

> I get the following result and a warning because the more complex model 
> (Mod1) has a lower log-likelihood
> than the base model (Mod0)
> 
> Analysis of Deviance Table
> 
> Model 1: Killed_Lamb ~ 1
> 
> Model 2: Killed_Lamb ~ I(scale(Dens_Lynx))
> 
>   NoPar  LogLik Df Deviance Pr(>Chi)
> 
> 1     4 -9849.6
> 
> 2     5 -9856.4  1   -13.56        1
> 


  I don't know exactly what's going on here, but I agree it's weird.
Some suggestions:

 * try fitting the more and less complex models with starting
values set to the parameters of the other model (i.e. fit the 
more complex model starting at the reduced model, and the reduced
model starting at the estimate of the full model [forcing the
slope variable to zero])

* see whether adding a binary lynx present/absent variable helps
(this would add flexibility in the lynx density variable -- you
could also add a spline based on the lynx density; e.g. library(splines),
response ~ ns(Dens_Lynx,3); when
you plot the data, do the fitted parameters seem reasonable?)

* you can try a zero-inflated model via expectation-maximization 
with lme4:

https://groups.nceas.ucsb.edu/non-linear-modeling/projects/owls/
    WRITEUP/owls.Rnw
https://groups.nceas.ucsb.edu/non-linear-modeling/projects/owls/R/
  owls_R_funs.R