Anova: 2-factor model.
Voorbeeld 12.9 - Buijs p. 346
2-factor model met replicaties.
> bron <- "http://www.mzandee.net/~zandee/statistiek/data/inkomen.txt"
> inkom <- read.table(bron, header=T)
> attach(inkom)
> names(inkom)
[1] "inkomen" "sexe" "opleiding"
> model <- aov(inkomen ~ sexe * opleiding, data=inkom)
# is identiek met de volgende aanroep van de functie aov:
> model <- aov(inkomen ~ sexe + opleiding + sexe:opleiding, data=inkom)
# dwz: onderzoek, op basis van een additief model, de effecten van sexe en opleiding,
# plus dat van hun interactie
> summary(model)
Df Sum Sq Mean Sq F value Pr(>F)
sexe 1 968.00 968.00 17.2600 0.0003563 ***
opleiding 3 1035.00 345.00 6.1516 0.0029656 **
sexe:opleiding 3 827.00 275.67 4.9153 0.0084018 **
Residuals 24 1346.00 56.08
De interactie van sexe en opleiding blijkt dus significant. De hoofdeffecten zijn daardoor minder interessant geworden.
De kwadraatsommen uit de tabel kunnen we ook 'op de hand' uitrekenen, mbv R als rekenmachine:
We beginnen met de Correctiefactor, de gekwadrateerde som van alle waarnemingen, gedeeld door het aantal waarnemingen. We hebben deze term nodig als we op de alternatieve manier (zie Buijs, p. 85) kwadraatsommen uit willen rekenen.
CF <- sum(inkomen)^2/length(inkomen)
De totale kwadraatsom is makkelijk uit te rekenen:
KStot <- sum(inkomen^2) - CF
Dat zelfde is het geval met de kwadraatsom die overeenkomt met alle variatie die wordt veroorzaakt door beide factoren sexe en opleiding plus hun interactie:
KS<-sum(as.vector(tapply(inkomen,list(sexe,opleiding),sum))^2)/4
KSso <- KS - CF
We delen door 4 omdat elk van de subtotalen de som is van 4 termen.
Let op het gebruik van as.vector om de tabel die door tapply wordt gemaakt om te zetten in een serie getallen [zonder tekst] .
De residuele kwadraatsom, die de binnenvariatie van de 8 steekproeven met elk 4 waarnemingen beschrijft, is nu uit te rekenen als:
KSres <- KStot - KSso
ofwel:
KSres <- sum(inkomen^2) - sum(KSso)
We moeten nu nog de kwadraatsommen van de hoofdeffecten sexe en opleiding uitrekenen en die van hun interactie. Daarvoor hebben we weer de functie tapply nodig om de gegevens uit de datamatrix op een slimme manier samen te vatten, en op te tellen. Eerst het effect van sexe:
KSs<-sum(((as.vector(tapply(inkomen,list(sexe),sum)))^2)/16)-CF
We delen door 16 omdat elk van de subtotalen (via tapply) de som is van 16 termen.
Nu het effect van opleiding:
KSo<-sum(((as.vector(tapply(inkomen,list(opleiding),sum)))^2)/8)-CF
We delen nu door 8 omdat elk van de subtotalen de som is van 8 termen.
Nu resteert nog de kwadraatsom van de interactie. Die kunnen we uitrekenen als de rest die overblijft als we alle nu bekende termen in mindering brengen op de totale kwadraatsom:
KSinter <- KStot - (KSo + KSs + KSres)
Voorbeeld Crawley: 2-factor model zonder replicaties.
> bron <- "http://www.mzandee.net/~zandee/statistiek/data/groei.txt"
> plant <- read.table(bron, header=T)
> attach(plant)
> names(plant)
[1] "Groei" "Licht" "Genotype"
> model2 <- aov(Groei ~ Genotype + Licht, data=plant)
> summary(model2)
Df Sum Sq Mean Sq F value Pr(>F)
Genotype 5 27.8750 5.5750 18.0811 7.092e-06 ***
Licht 3 7.1250 2.3750 7.7027 0.002404 **
Residuals 15 4.6250 0.3083
Voorbeeld Crawley: 2-factor model met replicaties.
> bron <- "http://www.mzandee.net/~zandee/statistiek/data/dieet.txt"
> dieet <- read.table(bron, header=T)
> attach(dieet)
> names(dieet)
[1] "growth" "diet" "coat"
> dieet
growth diet coat
1 6.6 A light
2 7.2 A light
3 6.9 B light
4 8.3 B light
5 7.9 C light
6 9.2 C light
7 8.3 A dark
8 8.7 A dark
9 8.1 B dark
10 8.5 B dark
11 9.1 C dark
12 9.0 C dark
> model <- aov(growth~diet*coat)
> summary(model)
Df Sum Sq Mean Sq F value Pr(>F)
diet 2 2.66000 1.33000 3.6774 0.09069 .
coat 1 2.61333 2.61333 7.2258 0.03614 *
diet:coat 2 0.68667 0.34333 0.9493 0.43833
Residuals 6 2.17000 0.36167
Hoe komen we aan de kwadraatsommen in deze tabel, als we ze 'op de hand' uit moeten rekenen, met R als rekenmachine?
We beginnen met de Correctiefactor, de gekwadrateerde som van alle waarnemingen, gedeeld door het aantal waarnemingen. We hebben deze term nodig als we op de alternatieve manier (zie Buijs, p. 85) kwadraatsommen uit willen rekenen.
> CF <- sum(growth)^2/length(growth)
De totale kwadraatsom is makkelijk uit te rekenen:
> KStot <- sum(growth^2) - CF
Dat zelfde is het geval met de kwadraatsom die overeenkomt met alle variatie die wordt veroorzaakt door beide factoren diet en coat plus hun interactie:
> KS <- sum(as.vector(tapply(growth,list(coat,diet),sum))^2)/2
> KScd <- KS - CF
We delen door 2 omdat elk van de subtotalen de som is van 2 termen.
Let op het gebruik van as.vector om de tabel die door tapply wordt gemaakt om te zetten in een serie getallen [zonder tekst] . Dit is wat tapply doet:
> tapply(growth,list(coat,diet),sum)
A B C
dark 17.0 16.6 18.1
light 13.8 15.2 17.1
De residuele kwadraatsom, die de variatie in de 6 steekproeven met elk 2 waarnemingen beschrijft, is nu uit te rekenen als:
> KSres <- KStot - KScd
> KSres
[1] 2.17
ofwel:
> KSres <- sum(growth^2) - sum(KScd)
Van de 5 kwadraatsom termen uit de variantieanalysetabel hebben we er nu 2. We moeten nu nog de kwadraatsommen van de hoofdeffecten diet en coat uitrekenen en die van hun interactie. Daarvoor hebben we weer de functie tapply nodig om de gegevens uit de datamatrix op een slimme manier samen te vatten, en op te tellen. Eerst het effect van diet:
> KSd<-sum(((as.vector(tapply(growth,list(diet),sum)))^2)/4)-CF
> KSd
[1] 2.66
We delen door 4 omdat elk van de subtotalen de som is van 4 termen.
Dit is wat tapply hier doet:
> tapply(growth,list(diet),sum)
A B C
30.8 31.8 35.2
Nu het effect van coat:
> KSc<-sum(((as.vector(tapply(growth,list(coat),sum)))^2)/6)-CF
> KSc
[1] 2.613333
We delen nu door 6 omdat elk van de subtotalen de som is van 6 termen.
Dit is wat tapply hier doet:
> tapply(growth,list(coat),sum)
dark light
51.7 46.1
Nu resteert nog de kwadraatsom van de interactie. Die kunnen we uitrekenen als de rest die overblijft als we alle nu bekende termen in mindering brengen op de totale kwadraatsom:
> KSinter <- KStot - (KSd + KSc + KSres)
> KSinter
[1] 0.6866667
Bron:
Buijs, A. - Statistiek om mee te werken. Stenfert Kroese, Groningen (2003)
Crawley, M.J. - Statistics. An introduction using R. Wiley, Hoboken, NJ, USA (2005)
