[R-sig-ME] Puzzled by simple example.
Rolf Turner
r.turner at auckland.ac.nz
Wed Apr 1 01:31:11 CEST 2009
I continue to struggle to understand the syntax of lmer() from the
lme4 package.
In an attempt to enlighten myself I simulated a rather simple data
set. The
imagined scenario is children's heights (in cm.) at ages 4, 5, 6, and
7, although
no attempt has been made to make the data realistic.
Since the data are balanced, I thought I could check/compare the lmer
() results
with results produced by aov().
The code I used is as follows:
#
# Script scr.01
#
# Generate the data. No attempt is made to make these
# data realistic.
set.seed(42)
Beta <- c(40,50,60,70)
B <- rnorm(80,0,2)
E <- rnorm(80,0,0.25)
Y <- Beta + B + E
ht.dat <- data.frame(y=Y,age=factor(rep(4:7,20)),child=factor(rep
(1:20,each=4)))
# Fit with lmer():
f1 <- lmer(y ~ 0 + age + (1|child),data=ht.dat)
# Fit with aov():
f2 <- aov(y ~ 0 + age + Error(child),data=ht.dat)
# Clean up:
rm(Beta,B,E,Y)
When one does summary(f1) one gets values of the estimates of the
mean heights
which are very close to the ``true'' values used to generate the
data. I am however
puzzled by the random effects estimates. These are:
Random effects:
Groups Name Variance Std.Dev.
child (Intercept) 1.0362e-07 0.0003219
Residual 4.6663e+00 2.1601657
Number of obs: 80, groups: child, 20
I thought at the ``child'' variance would be around 4 (the variance
used to simulate B) and that the residual variance would be around
0.0625
(the variance used to simulate E). This is clearly not so, so there is
clearly something I am misunderstanding here. I am not interpreting the
random effects results correctly. Can someone please point me in the
right direction?
In particular, how can I extract estimates of the variance of B and E?
I am assuming that the model I am fitting in both instances is
y = X*Beta + Z*B + E
where X = kronecker(One_20,I_4) and Z = kronecker(I_20,One_4), where
I_k is the
k-x-k identity and One_k is a k-x-1 matrix all of whose entries are 1.
Have I got this much right? (In the forgoing Beta is a 4-vector of
fixed effects
and B is a 20-vector of random effects.)
I am even more puzzled by the results of aov() --- I cannot connect
anything
with anything from these results. Where/how do I obtain estimates of
the
variance of B and the variance of E? Where/how do I obtain estimates
of the
fixed effects?
If I do coeff(f2) I get
child :
age4
220.0565
Within :
age4 age5 age6
-29.86516 -20.91547 -11.09382
The ``Within'' coefficients are the differences age4 - age7, age5 -
age7, age6 - age7,
respectively. I don't understand the ``child'' coefficient (age4 =
220.0565) at all.
Can anyone explain this to me?
One other puzzlement: If one just types ``f2'' or ``print(f2)'' one
gets messages,
for stratum 1 (child) ``3 out of 4 effects not estimable / Estimable
effects are
balanced'', and for stratum 2 (Within) ``1 out of 4 effects not
estimable / Estimated
effects may be unbalanced'' Huh? What is going on? How could a
design possibly
be more balanced?
I would be grateful for some enlightenment.
cheers,
Rolf Turner
######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
More information about the R-sig-mixed-models
mailing list