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...
Geraldine Mabille

From: Mabille, Geraldine
Sent: 11. oktober 2013 15:26
To: r-sig-mixed-models@r-project.org
Subject: glmmadmb- problems with explanatory variable with a lot of zeros?

Dear list members,
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)

AIC: 19722.8

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)
(Intercept)           2.9544     0.1374    21.5  < 2e-16 ***
I(scale(Dens_Lynx))   0.1122     0.0229     4.9  9.5e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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


That was indeed very strange and I tried to inspect the residuals from Mod1 to better understand what could be happening.
DATA$Res<-resid(Mod1)
head(DATA[order(-DATA$Res),c("Killed_Lamb","Dens_Lynx","Res")],10)

     Killed_Lamb    Dens_Lynx        Res

3435         120 0.0000000000 272.519920

4074          66 0.0000000000 149.706007

2494         120 0.0000000000 125.349338

4629          54 0.0000000000 122.414026

990           38 0.0000000000  50.755557

3094         238 0.0006861074   5.626109

80           204 0.0000000000   5.472785

774           10 0.0000000000   4.923185

3780          57 0.0000000000   4.645824

4115          16 0.0000000000   3.991187

As visible from the first and third lines, two datapoints with the same number of Killed_Lamb and the same number of Dens_Lynx get very different residuals. Furthermore, the first 5 lines show residuals much above the range of the other residuals (otherwise between -1.316 and 5.626).

I suspect this happens because the model can't handle so many zeros in the explanatory variable because if I try running a model with slaughtering weight of the lamb as explanatory variable
Mod2 <-glmmadmb(Killed_Lamb~I(scale(Weight)), random=~ 1|KOM, data=DATA,zeroInflation=TRUE,family="nbinom")
I get the following summary:

Call:

glmmadmb(formula = Killed_Lamb ~ I(scale(Weight)), data = DATA,

    family = "nbinom", random = ~1 | KOM, zeroInflation = TRUE)



AIC: 19658.7



Coefficients:

                 Estimate Std. Error z value            Pr(>|z|)

(Intercept)        3.0778     0.1499   20.54 <0.0000000000000002 ***

I(scale(Weight))  -0.0471     0.0338   -1.39                0.16

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Number of observations: total=2151, KOM=257

Random effect variance(s):

Group=KOM

            Variance StdDev

(Intercept)    3.638  1.907



Negative binomial dispersion parameter: 1.7665 (std. err.: 0.07474)

Zero-inflation: 0.10091  (std. err.:  0.0093901 )



Log-likelihood: -9824.37

With a higher log-likelihood for the more complex model (Mod2) compared to the base model (Mod0)
anova(Mod0,Mod2)

Analysis of Deviance Table



Model 1: Killed_Lamb ~ 1

Model 2: Killed_Lamb ~ I(scale(Weight))

  NoPar  LogLik Df Deviance          Pr(>Chi)

1     4 -9849.6

2     5 -9824.4  1    50.46 0.000000000001216 ***

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Has anyone experienced those problems before? Which alternatives do I have to run a model with overdispersion and zero-inflation? I have looked at modeling a log-normal Poisson model using lme4
DATA$Obs<-c(1:nrow(DATA))
Mod3<- glmer(Killed_Lamb~Dens_Lynx+(1|KOM)+(1|Obs), family=poisson,link=log, data=DATA)
But I obtain very ugly validation graphs even though the overdispersion seems to be taken into account
I also tried using MCMCglmm but I'm not so confident with using it since I've never used Bayesian statistics before...

Thanks heaps  if anyone could help with those questions :)
Geraldine Mabille

	[[alternative HTML version deleted]]

