# Experimental design & questions on use of generalized linear models

I have an ecological experiment for which I need to analyze bird count data. Here is the set up:

2 treatments (open/control), 3 regions. Not quite a full 3x2 factorial because in 2 regions there are 3 plots (250m x 250m) for both open & control (6 plots). For “reference” region there are just 3 plots (more similar to open plots in the upstream region).

Birds were counted weekly (can be considered independent) in each plot over 4 years. Bird type & densities are highly dependent on water depth (have depth for each survey date) – would like to use as covariate (note: the response isn’t linear, it’s unimodal, and it differs by species).

The goal is to examine treatment vs control differences in numbers of types of birds, & examine differences/interactions of region & year (as open plots mature & vegetation grows in). The closed plots contain many zeros (I may need zero-inflated model).

The literature is pointing me towards negative binomial or zero-inflated Poisson GLM. If I choose to do this, though, I have following questions:

• Can I nest plots within region for this type of model? I wanted unit of replication to be the week for each treatment, without sacrificing replication for each treatment x region plot.
• Also, any suggestions on what to do with a covariate like depth? It could really dominate the terms of a GLM.

I'm still learning R, but assuming based on this forum & speaking with others that I will need to do this in R or SAS. Advice is welcome on the software front, including what R packages I'll need to load.

• By "plot", do you mean area of land within which you will count the birds? – jrennie Aug 15 '12 at 13:03
• What is a "closed" plot? – jrennie Aug 15 '12 at 13:03
• Sorry -- sometimes I used "closed" synonymously w/ "control": – kayaker Aug 15 '12 at 16:59
• So control/closed plots are those in which vegetation wasn't removed. Open (treated) plots are those in which vegetation was removed (via herbicide/fire regimen). Yes a plot is an area of land 250m x 250m in which birds were counted by helicopter ea wk. – kayaker Aug 15 '12 at 17:00
• 2 regions had 3 sets of paired open/closed plots (side-by-side), for total of 6 pairs or 12 plots. There were also 3 "reference" plots in another region which are naturally more open, thus we hope to compare bird densities in the open plots only (from the 2 higher nutrient regions, n = 6 plots) w/ those 3 reference plots. – kayaker Aug 15 '12 at 17:03

Software: R is certainly a good choice. I use python for this sort of thing; I write my own objective/gradient function(s) and use one of the scipy optimizers, like L-BFGS. But, R is better if you aren't a strong programmer.

Caveat: I'm a machine learning guy, not a statistician, so please consider my answer to be one opinion, not the "right answer".

It sounds like your model should have at least coefficients for (1) is treatment?, (2) is control?, (3) each plot, (4) each region, (5) week-of-year, (6) week-of-year-and-region, (7) water depth, (8) bird species, (9) water-depth-and-bird-species. After including all of these, I'd look at residuals to try to determine any obvious ones I missed. Though, it sounds like you have a pretty good idea of all of the major covariates.

I would try different models (Poisson, negative binomial, zero-inflated Poisson) and use a hold-out set to determine which is more appropriate. I would use L2 regularization and seriously consider L2 normalizing the covariates (at least approximately).

If there is a strong covariate like depth, I'd definitely want to include it in my model, else other covariates may appear stronger than they really are.

If you choose a negative bionomial, then the package glmmADMB in R is the package for you. It is fantastic! It works just like a lmer model (where you can have fixed, random effects and nested components if you choose, along with the negative bionomial distribution/family. Best of luck.

Here is an example from the package description:

om <- glmmadmb(SiblingNegotiation ~ FoodTreatment * SexParent(1|Nest) +
offset(log(BroodSize)),
zeroInflation=TRUE, family="nbinom", data=Owls)

• (+1) Does this package also implement regularised regression (shrinkage) for negative binomial family? Any idea? Any other packages you're aware of? – Arun Dec 26 '12 at 14:32