[R] LME-glmmPQL formulation

[Andrew Beckerman]

Hi all -
R2.0.1 on OSX;MASS library;nlme library

I am trying to emulate the solution to a problem set that has normally  
been run in Genstat, using R.  The problem that I am having at the  
moment is with the following glmm question (using glmmPQL from the MASS  
library):

"We have two different forest habitats (first rotation thicket, and  
high forest) which we want to survey for the presence of our study  
animal. We survey both habitats on each of 10 days, and within each  
habitat we have five transects. The sampling unit is the number of  
animals counted per transect.  We therefore have two sources of random  
variation:
·   counts will vary between days due, say, to variation in weather
·   counts will vary between transects, within sites, for any number of  
(known and unknown) reasons
There is no relationship between transects at the two sites: transect 1  
in site 1 has no link to transect 1 in site 2, etc. The random term for  
transectis therefore nested within site, while the main effect of site,  
which is what we are interested in, is a fixed effect."

 > summary(dat)
     site         day           trans       count
  Here :50   Min.   : 1.0   Min.   :1   Min.   : 48.00
  There:50   1st Qu.: 3.0   1st Qu.:2   1st Qu.: 79.00
             Median : 5.5   Median :3   Median : 95.00
             Mean   : 5.5   Mean   :3   Mean   : 95.85
             3rd Qu.: 8.0   3rd Qu.:4   3rd Qu.:112.25
             Max.   :10.0   Max.   :5   Max.   :165.00

In Genstat, the (supposed) procedure is to fit a model with site as a  
fixed effect and then a random effects model of day+transect.site,  
where the transect.site indicates that there are 5 transects nested  
within each site.

My first thought was the following:

glmmPQL(count~site,data=dat,random=~day|site/transect, family="poisson")

however, the random effects are not separated into day and  
site/transect.  Instead, there is day|site and day|site %in% transect,  
which I realize makes sense in light of the model formulation.

my second guess was

glmmPQL(count~site,random=list(~day|site,~1|trans),family="poisson",data 
=dat2)

which estimates a random effect on ~day|site and on  
~1|trans%in%site..... which seems more appropriate, but does not give  
the same answers as I have for the genstat; nor does it estimate the  
p-value for site.  I guess my question is how to separate the two  
random effects so that there is a estimate for day and for  
transect/site.

I would be happy to provide the data if anyone needs/wants IT.

Cheers
andrew

------------------------------------------------------------------------ 
---------
Dr. Andrew Beckerman
Department of Animal and Plant Sciences, University of Sheffield,
Alfred Denny Building, Western Bank, Sheffield S10 2TN, UK
ph +44 (0)114 222 0026; fx +44 (0)114 222 0002
http://www.shef.ac.uk/beckslab
------------------------------------------------------------------------ 
----------