Andrew Dolman
2007-Mar-30 03:16 UTC
[R] problem using mcmcsamp() with glmer models containing interaction terms in fixed effects
Dear All,
I've been using mcmcsamp() successfully with a few different mixed models
but I can't get it to work with the following. Is there an obvious reason
why it shouldn't work with a model of this structure ?
*brief summary of objective:
I want to test the effect of no-fishing marine reserves on the abundance of
a target species.
I have samples at coral reef sites inside and outside of reserves, from
visits before and after the establishment of the reserves. Samples are
counts of the target species.
I plan to test this by looking at the interaction between reserve zone
(green/blue, green means no-fishing allowed) and the change in abundance
between the before sample (visit 0), and after sample (visit 1). -- please
say if you think there is a better way.
Samples come from several different regions, thus it would be nice to be
able to model differences between regions in the response to no-fishing and
to expand possible inference beyond this sample of regions.
My problem is getting a test for the fixed effects, see
http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests for background
on why I want to use mcmcsamp.
Thank you for your attention,
Andrew Dolman - Australian Institute of Marine Science.
--
andydolman@gmail.com
a.dolman@aims.gov.au
# required packages
library('lme4')
# some dummy data
counts <- c(rpois(75, 2), rpois(25,4)) # number of fish in a sample
[dependent variable - poisson distributed]
visit <- c(rep(0, 50), rep(1,50)) # sample visit number, 0 or 1 to indicate
pre and post no-fishing reserve establishment
zone <- as.factor(rep(c(rep("blue", 25),
rep("green",25)),2)) #
inside/outside no-fishing reserve
region <- as.factor(rep(c("island1", "island2"),50)) #
spatial grouping
factor
data1 <- data.frame(counts, visit, zone, region)
rm (counts, visit, zone, region)
# a quick look
bwplot(I(sqrt(counts))~as.factor(visit)|region, data=data1)
# two lmer models
lmer1 <- lmer (counts~visit+visit:zone + (1|region), data=data1,
family=poisson(link="log"))
lmer2 <- lmer (counts~visit+visit:zone + (visit+visit:zone|region),
data=data1, family=poisson(link="log"))
mcmc1 <- mcmcsamp(lmer1, n=10000) # works
mcmc2 <- mcmcsamp(lmer2, n=10000) # doesn't work
[[alternative HTML version deleted]]
