[R-sig-ME] Mixed model and negative binomial distribution

[Ben Bolker]

Seth W. Bigelow <seth at ...> writes:

> 
> Dear mixed-model brain trust: 
> 
> I am comparing snag (dead tree) densities 1 year and 5 years after
> silvicultural treatment in forest plots to densities prior to treatment. In
> nlme, my model is 
> 
> lme(snagnum~treatment, random=(~1|plot), correlation=corExp(form=~year)). 
> 
> {Treatment is a factor with values of Pre/1-year post/5-year post}. This
> gives reasonable output, but I'm having a niggling doubt that I should be
> using something akin to a negative binomial distribution, since about half
> of the values are zeros (i.e., many plots had no snags prior to treatment,
> and did not gain additional snags as a result of treatment). Can anyone
> suggest an appropriate package and associated syntax for doing this mixed
> model based on an alternative probability density function?

  Negative binomial would be a reasonable distribution; the other
answer gives you some methods for doing this.  *However*, incorporating
both serial correlation and non-Gaussian errors in a model of this form
is a bit of a nuisance.
The model you want might be something like

  snagnum ~ Poisson(lambda)
  lambda ~ MVN(mean=treatment,Sigma=...)

where the variance-covariance matrix gives you both some extra-Poisson
variation (to handle overdispersion) and some correlation between
observations.  I'm hoping Alain Zuur will pop up shortly to point you
to a reference in his new book that will tell you how to do this in
WinBUGS ...

  You *might* be able to do this via glmmPQL in the MASS package ...
but I would check your answers particularly carefully to see that
they make sense.