[R] Odp: Anova

Peter Dalgaard P.Dalgaard at biostat.ku.dk
Fri May 11 11:54:43 CEST 2007


Petr PIKAL wrote:
> Hi
>
> r-help-bounces at stat.math.ethz.ch napsal dne 11.05.2007 07:26:00:
>
>   
>> Doea anyone know how to compute components of variance analysis? For 
>>     
> example, 
>   
>> I have the score of pupils on a test. Each pupil belongs to a school and 
>>     
>
>   
>> within each school I may have several classes? How can I estimate the 
>>     
> variance
>   
>> of the pupils, schools, classes and the errror variance?
>>
>>
>> Any examples?
>>     
>
> What is wrong with examples in ?aov or ?lm pages?
>
> e.g.
> ## From Venables and Ripley (2002) p.165.
> data(npk, package="MASS")
>
> ## Set orthogonal contrasts.
> op <- options(contrasts=c("contr.helmert", "contr.poly"))
> ( npk.aov <- aov(yield ~ block + N*P*K, npk) )
> summary(npk.aov)
> coefficients(npk.aov)
>
> Regards
> Petr
>   
This does not give out the variance components. Actually, the model has
none.
What we do need (and have not, to the best of my knowledge) is
variations of the following:

> npk.aov <- aov(yield ~ Error(block) + N*P*K, npk)
> summary(npk.aov)

Error: block
          Df  Sum Sq Mean Sq F value Pr(>F)
N:P:K      1  37.002  37.002  0.4832 0.5252
Residuals  4 306.293  76.573

Error: Within
          Df  Sum Sq Mean Sq F value   Pr(>F)
N          1 189.282 189.282 12.2587 0.004372 **
P          1   8.402   8.402  0.5441 0.474904
K          1  95.202  95.202  6.1657 0.028795 *
N:P        1  21.282  21.282  1.3783 0.263165
N:K        1  33.135  33.135  2.1460 0.168648
P:K        1   0.482   0.482  0.0312 0.862752
Residuals 12 185.287  15.441

In this model, the estimated interblock variance can be computed as
(76.573-15.441)/4 = 15.283 qsince there are 4 observations per block.
This is easy to work out when you've learned that the residual mean
square for the blocks stratum is based on the variation of the block
means,  and that it is scaled so as to  be comparable to the Within 
stratum residuals. It does get trickier if there are crossed random effects.

The lmer() function from lme4 allows you to obtain variance components
(but it has issues with small-sample cases.)

-- 
   O__  ---- Peter Dalgaard             Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907



More information about the R-help mailing list