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

[Rolf Turner]

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}}