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

[Steven McKinney]

Hi Rolf, 

Quick observation in-line below

> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org
[mailto:r-sig-mixed-models-
> bounces at r-project.org] On Behalf Of Rolf Turner
> Sent: Tuesday, March 31, 2009 4:31 PM
> To: R-Sig Mixed-models
> Subject: [R-sig-ME] Puzzled by simple example.
> 
> 
> 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

This does not appear to contain a random effects component.
This just looks like average height data with two independent
gaussian rv's added as noise (so in essence it's just an rnorm(80, 0,
sqrt(4.0625)) random error).  This model should only need 20
rnorm errors at the child level.  So it's a matter of getting
this model data constructed appropriately according to what
you outline below.  I'll try to get this done and post it,
but you or someone else may get there before I have time.

As I see it, this is why aov() is choking - there's no different
error structure at the whole child level.

Also note the following from the aov() help page
" The default 'contrasts' in R are not orthogonal contrasts, and aov and
its helper functions will work better with such contrasts: see the
examples for how to select these. "

so you might need to set orthog contrasts first
## Set orthogonal contrasts.
op <- options(contrasts=c("contr.helmert", "contr.poly"))

HTH

Steve McKinney


> 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