Dear all,
I'm not very sure on the use of repeated measures in R, so some advice
would be very appreciate.
Here is a simple example similar to my real problem (R 2.6.0 for
windows): Lets supose I have annual tree production measured in 9
trees during 3 years; the 9 trees are located in 3 different mountains
(sites), and each tree receive different annual rainfall (different
locations). I would like to known the parameters that explain the
variability in production. The data would be something like:
set.seed(111)
mydat <- data.frame(tree= factor(rep(1:9,3)), yeargl(3,9,lab=2001:2003,
orde=T), site= gl(3,3,27,lab=c("A","B","C")),
rain= c(rnorm(9, 100), rnorm(9, 200),rnorm(9, 300)), prod51:77+rnorm(27, 1),
pr01= rbinom(27, 1, 0.5))
mydat
# see for instance
interaction.plot(mydat$year, mydat$site, mydat$prod)
#My first attempt was to use aov:
summary(aov(prod ~ rain + year + site + Error(tree), data=mydat))
#
# Error: tree
# Df Sum Sq Mean Sq F value Pr(>F)
# rain 1 36.814 36.814 6.4423 0.05201 .
# site 2 112.588 56.294 9.8513 0.01843 *
# Residuals 5 28.572 5.714
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
#
# Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
# rain 1 1688.50 1688.50 1142.6437 1.422e-15 ***
# year 2 12.74 6.37 4.3103 0.03319 *
# Residuals 15 22.17 1.48
#
# The results seems OK to me,
# Is there a way to get predictions from the model? (as in predict.lm)
# And to get the explained variance?
# Using the lme (which allows to predict), I suppose the same model
would be (is it right?):
library(nlme)
res <- lme(prod ~ rain + year + site, data=mydat, random= ~ 1 |
tree); anova(res)
# numDF denDF F-value p-value
# (Intercept) 1 15 23027.700 <.0001
# rain 1 15 1144.461 <.0001
# year 2 15 2.288 0.1358
# site 2 6 17.267 0.0032
# The rain, varies within tree and between trees (in time), thus the
aov give me the significance of each part (less significant Between
than Within)
# Does lme give me Within only? or it include both Between and Within?
Year was significant in the aov and not in the lme
# I also want to test a binary (binomial) variable (no production vs
production), so I guess I should use lmer
# First the same model as above but with lmer:
library(lme4)
res2 <- lmer(prod ~ rain + year + site + (1|tree), data=mydat); anova(res2)
res2
# Is this the correct way to fit the model above? I'm unsure; results
are slightly different.
# I understand that there is no predict for lmer models. anova does
not give me the significance for lmer modes, so I guess I should enter
the variables by steps and compare models with anova(m1, m2, ...)
# and now the binary data. Is this correct?
res3 <- lmer(pr01 ~ rain + year + site + (1|tree), data=mydat,
family=binomial); anova(res3)
res3
Many thanks for any comments
Juli
--
http://www.ceam.es/pausas
