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