[R-sig-ME] Puzzled by simple example.

[Sundar Dorai-Raj]

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
>