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
------------------------------------------------------------------------
----------
