[R-sig-ME] Puzzled by simple example.
Sundar Dorai-Raj
sdorairaj at gmail.com
Wed Apr 1 02:24:08 CEST 2009
Hi, Rolf,
Sorry, but I can't absorb all your commentary below, but the
fundamental problem I believe is your simulation is not reflected by
the model you are fitting with lmer (also, I don't use aov so I won't
comment there). Here's the simulation you want, which I hope will show
you what lmer is fitting:
set.seed(42)
Beta <- c(40,50,60,70)
B <- rnorm(20,0,2) ### note that length(B) = 20
E <- rnorm(80,0,0.25)
ht.dat <- data.frame(age=factor(rep(4:7,20)),child=factor(rep(1:20,each=4)))
ht.dat$y <- with(ht.dat, Beta[age] + B[child] + E)
library(lme4)
# Fit with lmer():
f1 <- lmer(y ~ 0 + age + (1|child),data=ht.dat)
summary(f1)
<snip>
Random effects:
Groups Name Variance Std.Dev.
child (Intercept) 6.888065 2.62451
Residual 0.063592 0.25218
</snip>
HTH,
--sundar
On Tue, Mar 31, 2009 at 4:31 PM, Rolf Turner <r.turner at auckland.ac.nz> wrote:
>
> 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}}
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list