Afsnit 6.6: Tosidet variansanalyse

Den globale opvarmning fører mange steder til mere nedbør og til kraftigere regnbyer (skybrud). Især i byområder, hvor der ikke er naturlige nedsivningsområder, kan dette give anledning til oversvømmelser. Et bidrag til løsning af dette kan være "grønne tage", det vil sige kasser, hvor planter vokser i passende vækstmedier. Der er mange faktorer, der spiller ind i designet af et grønt tag. For at få bedre viden om nogle af disse faktorer er der i artiklen Influence of water storage and plant crop factor on green roof retention and plant drought stress lavet et eksperiment, hvor evapotranspirationen studeres, når højden af vækstkasserne ændres og når plantetætheden ændres. Plantetætheden kan enten være 0, 1, 2 eller 4 planter per potte (faktoren taethed), og potterne kan enten være 15 cm høje eller 30 cm høje (faktoren hoejde). Der foretages en vægtmåling om morgenen og om aftenen, hvorefter plantekasserne vandes grundigt. Forskellen mellem de to vægtmålinger angiver evapotranspirationen for denne dag i gram, der efterfølgende omregnes til et volumen og endelig en højde ved at dividere med pottens overflade. Det endelige fordampningstal, ET, er summen af de daglige evapotranspirationer hen over eksperimentets varighed.

Image by Alexa from Pixabay

For hver af de otte kombinationer af plantetæthed og pottehøjde er der målt på 5 potter. Data er ikke angivet i artiklen, men kan aflæses fra figur i denne. I analysen nedenfor vil jeg som respons bruge logaritmen til den målte ET-værdi ud fra en formodning om, at den procentvise ændring fra for eksempel 1 til 2 planter kan være den samme for de to forskellige pottehøjder.

I kodevinduet nedenfor laves en figur med boxplot for de otte grupper. At lave et boxplot baseret på blot 5 observationer er ækvivalent med blot at vise de 5 dataværdier. Vi kan se, at gruppen med 1 plante og den store pottehøjde skiller sig ud, ved at der er to observationer, der er meget lave i relation til det generelle billede for alle data. Vi kan ikke lave et qqplot for hver gruppe, da der kun er 5 observationer i grupperne. I stedet laves der i kodevinduet et qqplot for hver pottehøjde baseret på residualerne fra en ensidet variansanalyse for hver højde, hvor der deles ind i de fire plantetæthedsgrupper. De to qqplots viser, at det er rimeligt at bruge en normalfordeling til beskrivelse af data. Der er en lille tendens til større varians i den del med den store pottehøjde, hvilket især skyldes gruppen med 1 plante og den store pottehøjde som omtalt ovenfor.

6.6.1 Boxplot og qqplot opdelt efter to faktorer

Se opstartskoden (til/fra)

Jeg vil nu beskrive den tosidede variansanalysemodel generelt. Data består af målinger fra $n$ uafhængige stokastiske variable $X_1,\ldots,X_n.$ Disse inddeles i grupper ved hjælp af to faktorer $G=(G_1,\ldots,G_n)$ og $H=(H_1,\ldots,H_n).$ Vi starter med en model, hvor både middelværdi og varians afhænger af, hvilken af de $k\cdot m$ grupper observationen tilhører. For nemheds skyld betegnes niveauerne i de to faktorer blot med tal.

Statistisk Model 6.6.2. (Grundlæggende tofaktor gruppemodel)

$X_{i}\sim N(\mu_{G_i,H_i},\sigma_{G_i,H_i}^2),\enspace i=1,\ldots,n,\enspace (\mu_{11},\ldots,\mu_{km},\sigma_{11},\ldots,\sigma_{km}) \in \mathbf{R}^{k\cdot m}\times \mathbf{R}^{k\cdot m}_+.$

Når vi reducerer til modellen, hvor alle grupperne har den samme varians, kaldes dette den tosidede variansanalysemodel.

Statistisk Model 6.6.3. (Tosidet variansanalysemodel (twoway anova))

Vi betragter uafhængige stokastiske variable $X_1,\ldots,X_n,$ der deles ind i grupper efter en faktor $G$ med faktorniveauerne $1,2,\ldots,k$ og en faktor $H$ med niveauerne $1,2,\ldots,m.$ Hver gruppe har sin egen middelværdi og alle grupperne har samme varians,

$X_{i}\sim N(\mu_{G_i,H_i},\sigma^2),\enspace i=1,\ldots,n,\enspace (\mu_{11},\ldots,\mu_{km},\sigma) \in \mathbf{R}^{k\cdot m}\times \mathbf{R}_+.$

Den tosidede variansanalysemodel har en meget vigtig undermodel kaldet den additive model. Modellen er vigtig, da den giver mulighed for simple fortolkninger af de parametre, der indgår, hvilket beskrives detaljeret i næste afsnit.

Statistisk Model 6.6.4. (Additive model)

Vi betragter uafhængige stokastiske variable $X_1,\ldots,X_n,$ der deles ind i grupper efter en faktor $G$ med faktorniveauerne $1,2,\ldots,k$ og en faktor $H$ med niveauerne $1,2,\ldots,m.$ Middelværdien i den additive model kan skrives som et bidrag fra faktor $G$ $(\zeta_1,\ldots,\zeta_k)$ plus et bidrag fra faktor $H$ $(\eta_1,\ldots,\eta_m),$

$X_{i}\sim N(\zeta_{G_i}+\eta_{H_i},\sigma^2),\enspace i=1,\ldots,n,\enspace (\zeta_{1},\ldots,\zeta_{k},\eta_1,\ldots,\eta_m,\sigma) \in \mathbf{R}^{k+m}\times \mathbf{R}_+.$ Modellen har $d(M)=k+m-1$ frie parametre i middelværdispecifikatioen.

Fra den additive model kan vi prøve at reducere modellen til en model, hvor middelværdien kun har et bidrag fra faktoren $G$ (eller kun har et bidrag fra faktoren $H$ ). Dette fører os tilbage til den ensidede variansanalysemodel 6.2.3.

For at få en fornemmelse af om data kan beskrives med den additive model, kan man lave et interaktionsplot. Den indbyggede funktion i R er lidt mangelfuld på dette punkt, så i stedet anbefaler jeg en funktion additivitetsPlot, som findes i filen Rfunktioner.txt, jævnfør punktet Egne funktioner i R i afsnit 1.2.

I et interaktionsplot beregner man gennemsnit i alle grupperne givet ved opdeling efter $G*H.$ Gennemsnit afsættes mod niveauerne for den ene faktor, og alle gennemsnit, der ligger på det samme niveau af den anden faktor, forbindes. Hvis data kan beskrives med den additive model, Statistisk Model 6.6.4 ovenfor, afspejler gennemsnittene i figuren altså $\zeta_u+\eta_v$ afsat mod for eksempel $u=1,\ldots,k,$ og punkterne med samme værdi af $v$ forbindes. De kurver, der fremkommer, svarer altså til kurven $(u,\zeta_u),$ der parallelforskydes med værdierne fra $\eta_v.$ I et interaktionsplot prøver vi derfor at vurdere, om kurverne ser ud til at være parallelle.

6.6.5 Interaktionsplot

I det følgende kodevindue vises interaktionsplots baseret på funktionen additivitetsPlot der er indskrevet i opstartsdelen. Input til denne funktion er de to faktorer, der bruges til at dele data op i undergrupper, og vektoren med responsværdierne. Efter at I har kørt koden så prøv at skifte additivitetsPlot ud med den indbyggede R-funktion interaction.plot.

Se opstartskoden (til/fra)

Begge figurer viser approksimative parallelle kurver, hvilket tyder på, at data kan beskrives med den additive model.

6.6.1 Analyse i R og parametrisering

Den tosidede variansanalysemodel 6.6.3 i foregående afsnit, hvor hver gruppe har sin egen middelværdi, analyseres med kaldet x $\sim$ G*H. Den additive model 6.6.4 analyseres med kaldet x $\sim$ G+H. Lad os starte med at forstå output fra summary for den sidste model. Vi kan forstå output ved, for $u=1,\ldots,k$ og $v=1,\ldots,m,$ at skrive

$\zeta_u+\eta_v=(\zeta_1+\eta_1)+(\zeta_u-\zeta_1)+(\eta_v-\eta_1) =\text{Intercept+Gu+Hv},$ hvor højresiden viser de parametre, der bruges ved kaldet x $\sim$ G+H. Vi ser her at modellen kan parametriseres med $k+m-1$ parametre, nemlig $\zeta_1+\eta_1$ , $\zeta_2-\zeta_1,\ldots,\zeta_k-\zeta_1$ og $\eta_2-\eta_1,\ldots,\eta_m-\eta_1$ . Går vi nu tilbage til den tosidede variansanalysemodel, skriver vi i stedet

$\begin{aligned} \mu_{u,v} &=\mu_{1,1}+(\mu_{u,1}-\mu_{1,1})+(\mu_{1,v}-\mu_{1,1})+ (\mu_{u,v}-\mu_{u,1}-\mu_{1,v}+\mu_{1,1}) \\ & = \text{Intercept+Gu+Hv+Gu:Hv}, \end{aligned}$ hvor den anden linje viser de parametre, der bruges ved kaldet x $\sim$ G*H. Det sidste led kaldes interaktionen mellem de to faktorer. I den tosidede variansanalysemodel er Gu således forskel mellem niveau $u$ og niveau 1 for faktoren $G,$ når faktoren $H$ ligger på niveau 1, og vice versa for Hv. Det nyttige ved den additive model er, at Gu nu er forskel mellem niveau $u$ og niveau 1 for faktoren $G,$ uanset hvilket niveau faktoren $H$ befinder sig på, og Hv er forskel mellem niveau $v$ og niveau 1 for faktoren $H,$ uanset hvilket niveau faktoren $G$ befinder sig på.

Parametrisering i R

Som for den ensidede variansanalyse i afsnit 6.4 betragtes simulerede data med spredning $\sigma=0,$ således at vi direkte kan se parametrene, der bruges i R. Vi betragter den additive model inden for den tosidede variansanalysemodel.

Hvad er middelværdien for en observation med gFak på niveau "2" og hFak på niveau "B" ? Hvad er værdien af Intercept+gFak2+hFakB ?
Udtryk R-parametrene Intercept, gFak2, hFakB og hFakC ud fra eta og zeta.
Opskriv de statistiske modeller svarende til henholdsvis kaldet lm(x~gFak*hFak) og til kaldet lm(x~gFak+hFak).

Svar: Forstå output

Fra gFak kommer bidraget 3 til middelværdien og fra hFak bidraget 0, hvorfor middelværdien er $3+0=3.$ $\text{Intercept+gFak2+hFakB}=3+2-2=3,$ som er middelværdien, vi lige har udregnet.
Intercept=eta[1]+zeta[1], gFak2=eta[2]-eta[1], hFakB=zeta[2]-zeta[1] og hFakC=zeta[3]-zeta[1].
Det $i$ 'te respons er $X_i,$ og den $i$ 'te værdi af de to faktorer er $\text{gFak}_i$ og $\text{hFak}_i.$ Model svarende til det første kald siger, at hver undergruppe givet ved de to faktorer har sin egen middelværdi, $X_i\sim N(\mu_{\text{gFak}_i,\text{hFak}_i},\sigma^2).$ Model svarende til det andet kald siger, at middelværdien består af et bidrag fra gruppen bestemt af faktoren gFak plus et bidrag fra gruppen bestemt af faktoren hFak, $X_i\sim N(\eta_{\text{gFak}_i}+\zeta_{\text{hFak}_i},\sigma^2).$

Eksempel 6.6.6. (Grønne tage)

Vi betragter data omtalt i starten af dette afsnit omkring evapotranspirationen fra grønne tage. Lad $\text{logET}_i$ være den stokastiske variabel, der angiver respons (logaritmen til ET) for den $i$ 'te måling, og lad $\text{taethed}_i$ og $\text{hoejde}_i$ være de tilhørende værdier for de to faktorer taethed og hoejde, hvor taethed har niveauerne "t0", $\ldots$ , "t4", og hoejde har niveauerne "h1" og "h2". Lad os starte med den grundlæggende tofaktormodel 6.6.2, her skrevet som

$\text{LogET}_i\sim N\big(\mu_{\text{taethed}_i,\text{hoejde}_i}, \sigma^2_{\text{taethed}_i,\text{hoejde}_i}\big),\enspace i=1,\ldots,40,$ hvor hver gruppe bestemt af taethed*hoejde har sin egen middelværdi og sin egen varians, og de otte middelværdier og varianser kan variere frit. Først undersøges hypotesen om fælles varians:

$H:\enspace \sigma^2_{\text{t0},\text{h1}}= \sigma^2_{\text{t0},\text{h2}}= \sigma^2_{\text{t1},\text{h1}}= \sigma^2_{\text{t1},\text{h2}}= \sigma^2_{\text{t2},\text{h1}}= \sigma^2_{\text{t2},\text{h2}}= \sigma^2_{\text{t4},\text{h1}}= \sigma^2_{\text{t4},\text{h2}}.$ Beregningerne nedenfor i R viser, at Bartlett teststørrelsen er 20.5, og den tilhørende $p$ -værdi fra en $\chi^2(7)$ -fordeling er 0.005. Dette er en noget lav $p$ -værdi. Som vi har set i figur fra den skjulte kode ovenfor, er det især situationen med 1 plante og den store pottehøjde, der har en afvigende stor varians. Vi kan også se dette direkte ved at lave et test for ens varianser i de syv andre grupper. Fra beregningerne nedenfor ses, at Bartlett teststørrelsen nu er 12.2, og den tihørende $p$ -værdi er 0.058. Jeg vil nedenfor betragte modellen med samme varians i alle 8 grupper og til sidst i næste afsnit i eksempel 6.7.2 se på, at de konklusioner, vi kommer frem til, er robuste over for den afvigende varians i den ene gruppe. Den tosidede variansanalysemodel 6.6.3, med sammme varians i de 8 grupper, formuleres som

$\text{LogET}_i\sim N\big(\mu_{\text{taethed}_i,\text{hoejde}_i}, \sigma^2\big),\enspace i=1,\ldots,40,$ hvor de 8 middelværdier og den fælles varians kan variere frit. Vi ønsker nu at teste reduktionen til den additive model 6.6.4 skrevet som

$\text{LogET}_i\sim N\big(\zeta_{\text{taethed}_i}+\eta_{\text{hoejde}_i}, \sigma^2\big),\enspace i=1,\ldots,40.$ Vi har ovenfor lavet interaktionsplots, der viser overensstemmelse med den additive model. Parametertabellen hørende til model $M_1$ viser, at der er tre parametre, der vedrører interaktionen mellem taethed og hoejde, nemlig taethedt1:hoejdeh2, taethedt2:hoejdeh2 og taethedt4:hoejdeh2. Test for, at hver enkel af disse tre parametre kan sættes lig med nul, aflæses under $t$ -testet i parametertabellen, og giver alle en meget høj $p$ -værdi, hvilket understøtter observationen fra interaktionsplottet. I næste afsnit indføres et $F$ -test for at teste hypotesen om additivitet, og for data her fås en $p$ -værdi på 0.80. Konklusionen er derfor, at data ikke strider mod hypotesen om additivitet.

Vi laver nu en parametertabel under den additive model. Fra rækken hoejdenh2 ses, at et test for ingen effekt af pottehøjden giver en $p$ -værdi på $1.1\cdot 10^{-6}$ , hvorfor vi konkluderer, at pottehøjden har en stor effekt. Samtidig viser tabellen en kraftig forskel mellem nul planter og enten 1 2 eller 4 planter, hvorfor vi konkluderer, at plantetaetheden har en effekt. Et formelt $F$ -test for denne hypotese laves i eksempel 6.7.2 i næste afsnit. Et 95%-konfidensinterval for forskel mellem middelværdierne af logaritmen til ET for den store og den lille pottehøjde er $[0.20,\, 0.41].$ Konfidensinterval for forskel mellem 1 og ingen plante er $[1.41,\, 1.71],$ og mellem 4 planter og ingen plante er $[1.73,\, 2.03].$ Der er som forventeligt en stor forskel mellem nul planter og 1 plante, hvorimod stigningen fra 1 til 2 planter er cirka 0.16, og igen en stigning på cirka 0.16 fra 2 til 4 planter. Da vi kigger på logaritmen til ET, vil en fordobling svare til en stigning på cirka 0.7. Vi er således langt fra en fordobling, når vi går fra 1 plante til 2 planter.

6.6.7 Bartletts test og parametertabeller

Se opstartskoden (til/fra)

Foregående Næste