[R] Mixed Model nested ANOVA (lme help)
Stephen Cole
swbcole at gmail.com
Fri Mar 14 14:34:30 CET 2008
Hello Again R help readers - I have posted in the past and want to
thank all those who replied with helpful suggestions.
This regards mixed model anova, and the nlme package.
My purpose is to make a comparison of barnacle recruit density from 3
different regions. In each of these regions i have sampled 6 randomly
chosen sites. I sampled 20 quadrats at each site giving me a total
sample size of 360. Thus i have as a design
Response: density
Factors:
regions (3 levels, fixed)
sites (18 levels, random)
I started with creating a dataframe that included density, region
(coded 1, 2, 3), and site (coded 1,2...18)
I then created a grouped data object with the code;
datag <- groupedData(density_recruit ~ region | site, data=data,
labels = list(x="Region", y="Recruit Density"), units =
list(y="dm^2"))
1: my first question is this code correctly creating a grouped data
object for a nested design with one nested random factor?
I than used lme to build my full model with the code; (i had to
transform my response)
lmm1 <- lme(asinh(density_recruit) ~ region, data=datag, random=list(site=~1))
2: Is this code correct, i have tried other variations such as
lmm1 <- lme(asinh(density_recruit) ~ region, data=datag, random=~1| site)
but i essentially get the same results.
I can assess significance of my fixed factor using anova(), however i
would also like to know if the random factor is significant or not.
>From the help archives and the P/B book on mixed models i think the
correct procedure is to build a model without the random effect and
compare it to the model with the random effect using anova()
Thus i created a model without the random effect using the code and
compared lmm1 to lmod
lmod <- lm(asinh(density_recruit) ~ region, data=datag)
3. my printout from anova() was
Model df AIC BIC logLik Test L.Ratio p-value
lmm1 1 5 1458.319 1477.708 -724.1596
lm1 2 4 1677.033 1692.544 - 834.5165 1 vs 2 220.7138 <.0001
Does this mean that lmm1 is a superior model given the data? I would
say yes based on the lower AIC and BIC however i do not know what the
L.ratio is and what the p-value of <.0001 means.
Thanks again for any help on this issue.
Stephen Cole
M.Sc candidate
Marine Ecology Lab
Saint Francis Xavier University
More information about the R-help
mailing list