[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