Hello mixed-effects community,
First, thank you in advance for any assistance you can provide us. I am glad to see there is a r-sig mixed effects forum! I am working with a colleague on a mixed-effects design, both of us are relatively new to the mixed effects world. We were wondering if the forum can provide some advice on how to proceed/code a mixed effect model in R. My colleague previously posted to the general R-sig-eco but this forum was just suggested as a better venue.
We have a repeated measures design, 8 regions, 8 plots (sampling plates) nested within region, sampled every 2 months for 8 months (4 times). We have 2 covariates that have only 1 measure for each region per sampling period (COV1, COV2); these are estimates of gamma diversity. We also have 3 variables that are measured at the plot scale (two measures of alpha diversity and one measure of habitat, I have 8 individual measures for each of these per region per time period), one of these is the response (RESP) and the other two are additional covariates (exp1, exp2). Each regional and plot variable is measured at each time. So we have a repeated measures model with fixed effects including COV1, COV2, exp1, exp2, time. Random effects are plot nested within region. All variables are counts except one covariate. I sqrt transformed the response and ran the following model using lme (following the examples in Zuur et al. 2009 and Crawley 2007):
MODEL.LME <- LME(RESP ~ COV1 + COV2 + EXP1 + EXP2 + TIME, RANDOM = ~1|REGION/PLOT, METHOD="REML")
I also tried some interactions in the model, e.g.:
MODEL.LME2 <- LME(RESP ~ COV1 + COV2 + EXP1 + EXP2 + TIME + TIME*EXP2, RANDOM = ~1|REGION/PLOT, METHOD="REML")
Both of these models ran, but being new to the world of mixed models I want to make sure we set up the models correctly (prior to moving to model selection and validation). Our ultimate question is how does both measures of gamma diversity (COV1, COV2) and the two plot level measures (EXP1, EXP2) affect our diversity measure of interest? And do these relationships change with time - the easiest option would be to remove time and run an individual lme for reach time period but believe this to be the worst option. Is it possible to perform a post-hoc test with an interaction in lme (like the TIME*EXP2 effect) and are post-hoc tests in a mixed effects models valid? This may be sufficient for our needs. I have found examples in the R-help for main effects using multicomp and the glht function but have not been able to code the summary function for an interaction (was successful for a main effect time). I also am thinking we may need to add a within-group correlation structure using the "correlation = " but not sure if it is necessary in our case (possibly corAR1?).
I am also interested in using the lmer function in lme4 (reading through the archives of this forum provided many examples!). I am interested in using this package as my response is a count (richness) and would like to test if a model using a Poisson will fit my data better (I would like to compare output from all these models if possible). The model I coded was:
MODEL.LMER <- LMER(RESP ~ COV1 + COV2 + EXP1 + EXP2 + TIME + (1|Region/Plot), data = dataset, family = poisson)
But I get an error that states:
Warning in Plot:Region :
numerical expression has 256 elements: only the first used
Warning in Plot:Region :
numerical expression has 256 elements: only the first used
Error: length(f1) == length(f2) is not TRUE
First time coding an lmer equation, not sure I have successfully followed the example in Zuur and Crawleys' books nor the reference manual.
Thank you again,
Kelly
Kelly O. Maloney
Post-Doctoral Research Ecologist
Smithsonian Environmental Research Center
647 Contees Wharf Rd, PO Box 28
Edgewater MD 21037
443.482.2297
[[alternative HTML version deleted]]