[R-sig-ME] Puzzled by simple example.
Rolf Turner
r.turner at auckland.ac.nz
Wed Apr 1 02:28:47 CEST 2009
On 1/04/2009, at 12:40 PM, Robert Kushler wrote:
> You need
>
> B <- rep(rnorm(20,0,2),each=4)
>
> (the "child effect" should be the same for all four ages).
Right. That was a bit of duhhhhh, wasn't it?
Anyhow --- I fixed that, and now the lmer() results
make sense. I get a standard deviation of 2.62 for
the child effect (close enuff!) and 0.25 for the residual
effect (bang on!).
I still cannot make head or tail of the aov results however.
Is my call to aov() correct? I'm still getting the ``Estimated
effects may be unbalanced'', usw, messages.
cheers,
Rolf Turner
> I'm generally skeptical of forcing the intercept to be zero,
> but it's not illegal. :-)
>
> Regards, Rob Kushler
>
> Rolf Turner 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
>>
######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
More information about the R-sig-mixed-models
mailing list