# Re: [R] Odp: Anova

From: Peter Dalgaard <P.Dalgaard_at_biostat.ku.dk>
Date: Fri, 11 May 2007 11:54:43 +0200

Petr PIKAL wrote:
> Hi
>
> r-help-bounces@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

______________________________________________
R-help_at_stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help