[ADMB Users] help with a negative binomial, zero inflated model formula

[Ben Bolker]

On 12-02-06 05:40 PM, Danson, Bryan wrote:

> I am attempting to run a glmm on data which I have already determined
> display a negative binomial distribution and are zero inflated. My
> question is in the syntax of the formula for the proper analysis.  I
> have a count of the number of legal lobsters within a trap. There are
> five trap designs and were pulled on 12 different occasions. The trap
> pulls are independent, and therefore a random effect.  The
> trap design is my fixed effect.  Of those five trap designs, four are
> alternatives to a standard trap (control).
> 


> I am trying to determine if the legal lobster catch rate differs
> between the alternative trap designs and the standard wood trap.  I
> have run the model as follows:
> 
> mod <- glmmadmb(legal~trap_type+(1|pull/trap_type), data=l.lob,
> family="nbinom", zeroInflation=T); summary(mod)

  This seems reasonable.  In general it is slightly safer to use TRUE
rather than abbreviating it to T ...

> 
> Is this correct? This is the output:
> 
> summary(mod)
> 
> Call: glmmadmb(formula = legal ~ trap_type + (1 | pull/trap_type), 
> data = l.lob, family = "nbinom", zeroInflation = T)
> 
> 
> Coefficients: Estimate   Std. Error   z value   Pr(>|z|) (Intercept)
> 0.2450        0.2612        0.94          0.348 trap_typewire basket
> -0.0949       0.1923       -0.49         0.622 trap_typewire w/ wood
> frame   0.3045       0.1761        1.73          0.084       . 
> trap_typewire w/ wood slats     0.1134        0.1804        0.63
> 0.530 trap_typewood                                 0.1037
> 0.1797        0.58          0.564 --- Signif. codes:  0 '***' 0.001
> '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> 
> Number of observations: total=1004, date=12, pull:trap_type=60 Random
> effect variance(s): Group=date Variance   StdDev (Intercept)
> 0.5152        0.7178 Group=pull:trap_type Variance   StdDev 
> (Intercept)                          0.04389      0.2095 Negative
> binomial dispersion parameter: 2.8077 (std. err.: 0.84525) 
> Zero-inflation: 0.40333  (std. err.:  0.036329 )
> 
> Log-likelihood: -1316.88
> 
> In the output, it appears that there is no difference between the
> alternative trap types and the standard wood trap (control).  But
> which p-value do I refer to?  

  I don't remember whether drop1() works for glmmadmb or not (I don't
think so), but you should be able to test the overall effect of trap
type by fitting the model with legal ~ 1 + (1|pull/trap_type) [i.e. no
fixed effect of trap] and using anova(full_model, reduced_model) to test
the difference.

  Be aware that computation of correct degrees of freedom, significance,
etc., is quite a can of worms for generalized linear mixed models -- see
<http://glmm.wikidot.com/faq> for some starting points.  The bottom line
is that the p-values you get from summary(glmmADMBfit) assume (1)
arbitrarily large data sets (i.e. large denominator degrees of freedom
in the classical ANOVA sense) and (2) quadratic likelihood surfaces
(i.e., they are Wald tests).  #2 is relatively easy to relax by using
anova() instead but #1 is difficult.

> Also, the last alternative trap design
> (vertical wood) is not on the coefficients list, is this because all
> of the other trap types are being compared to it? 

  Yes.

 If so, how do I
> change the syntax in the formula to test the alternatives against the
> standard wood trap?

  l.lob$trap_type <- relevel(l.lob$trap_type,"vertical wood")

  Questions about glmmADMB are probably slightly more appropriate on the
r-sig-mixed-models at r-project.org mailing list.

  cheers
    Ben Bolker