7. Analysis of Variance

Author

Lieven Clement

Published

Invalid Date

1 Prostacycline dataset

Effect van arachidonzuur op het prostacycline niveau in het bloedplasma.
3 verschillende concentraties van arachidonzuur:
- laag (L, 10 eenheden)
- gemiddeld (M, 25 eenheden)
- hoge dosis (H, 50 eenheden)
Prostacycline concentratie in bloed plasma via gecalibreerde elisa fluorescentie meting
12 ratten worden at random toegekend aan elke behandelingsgroep.
Factoriële proef, volledige gerandomiseerde proefopzet, “completely randomized design” CRD.

prostacyclin <- read_tsv("https://raw.githubusercontent.com/statOmics/sbc/master/data/prostacyclin.txt")

prostacyclin <- prostacyclin |>
  mutate(dose = as.factor(prostacyclin$dose))

head(prostacyclin)

# A tibble: 6 × 2
  prostac dose 
    <dbl> <fct>
1    19.2 10   
2    10.8 10   
3    33.6 10   
4    11.9 10   
5    15.9 10   
6    33.3 10

1.1 Data exploration

prostacyclin |>
  ggplot(aes(x = dose, y = prostac, fill = dose)) +
  geom_boxplot() +
  geom_point(position = "jitter") +
  ylab("prostacyclin (ng/ml)")

prostacyclin |>
  ggplot(aes(sample = prostac)) +
  geom_qq() +
  geom_qq_line() +
  facet_grid(~ dose)

Data in de drie groepen lijkt normaal verdeeld en de variantie is ongeveer gelijk: [Y_i N(_j,^2),] met $j= \text{1, 2, 3}$

1.2 Onderzoeksvraag

Vraagstelling kan vertaald worden in volgende hypotheses

$H_0$: De arachidonzuurconcentratie heeft geen effect op het gemiddelde prostacycline niveau bij ratten [ H_0:_1=_2 = _3 ]
$H_1$: De arachidonzuurconcentratie heeft een effect op het gemiddelde prostacycline niveau bij ratten. Wat betekent dat minstens twee gemiddelden verschillend zijn: [H_1: j,k {1,,g} : _j_k]

naïeve benadering: nulhypothese op splitsen in partiële hypotheses [ H_{0jk}: _j=k H{1jk}: _j _k ]

Elk van deze partiële hypotheses testen met two-sample $t$-testen

$\rightarrow$ Probleem van meervoudig toetsen + verlies van power.

$\rightarrow$ $H_0:\mu_1=\mu_2=\mu_3$ testen met 1 enkele test.

2 Analyse van variantie

Correcte oplossing voor het testprobleem: variantie-analyse, afgekort door ANOVA (ANalysis Of VAriance)
We leiden de methode af voor de meest eenvoudige uitbreiding met 3 groepen (prostacycline voorbeeld)
Data modelleren a.d.h.v een lineair model door gebruik te maken van dummy variabelen.
1 dummy variable minder nodig hebben dan het aantal groepen. Hier dus 2 dummy variabelen.
De veralgemening naar g groepen $g>3$ is triviaal (extra dummy variabelen)

2.1 Model

\[\begin{eqnarray} Y_i &=& g(x_{i1},x_{i2}) + \epsilon_i\\ Y_i &=& \beta_0+\beta_1 x_{i1} +\beta_2 x_{i2} +\epsilon_i \end{eqnarray}\]

$Y_i$ de uitkomst voor observatie $i$ ($i=1,\ldots, n$)
$\epsilon_i\text{i.i.d.} N(0,\sigma^2)$
en dummyvariabelen \[x_{i1} = \left\{ \begin{array}{ll} 1 & \text{ als observatie $i$ tot de middelste dosisgroep behoort (M)} \\ 0 & \text{ als observatie $i$ behoort tot een andere dosisgroep} \end{array}\right.\] \[x_{i2} = \left\{ \begin{array}{ll} 1 & \text{ als observatie $i$ behoort tot de groep met hoge doses (H)} \\ 0 & \text{als observatie $i$ behoort tot een andere dosisgroep} \end{array}\right. .\]
Lage dosis groep (L) with $x_{i1}=x_{i2}=0$ is referentie groep

Regressiemodel kan worden herschreven als een model voor elke groep: \[\begin{eqnarray*} Y_{i\vert \text{dose=L}} &=& \beta_0+\epsilon_i \\ Y_{i\vert \text{dose=M}} &=& \beta_0+\beta_1+ \epsilon_i \\ Y_{i\vert \text{dose=H}} &=& \beta_0+\beta_2 + \epsilon_i \end{eqnarray*}\] met $\epsilon_i \sim N(0,\sigma^2)$

Interpretatie van model parameters: \[\begin{eqnarray*} \beta_0 &=& \text{E}\left[Y_i \mid \text{Behandeling met lage dosis groep L}\right] \\ \beta_1 &=& (\beta_0+\beta_1)-\beta_0 = \text{E}\left[Y_i \mid \text{behandeling M}\right] - \text{E}\left[Y_i \mid \text{behandeling L}\right] \\ \beta_2 &=& (\beta_0+\beta_2)-\beta_0 = \text{E}\left[Y_i \mid \text{behandeling H}\right]-\text{E}\left[Y_i \mid \text{behandeling L}\right]. \end{eqnarray*}\]

$\beta_0$ is gemiddelde uitkomst in groep L
$\beta_1$ is effect (verschil in gemiddelde concentratie) van groep M t.o.v. groep L
$\beta_2$ is effect van groep H t.o.v. groep L

We herformuleren het model door $\mu$-notaties te gebruiken: \[\begin{eqnarray*} Y_{i\vert \text{dose=L}} &=& \beta_0+\epsilon_i = \mu_1+\epsilon_i \\ Y_{i\vert \text{dose=M}} &=& \beta_0+\beta_1+ \epsilon_i = \mu_2+\epsilon_i \\ Y_{i\vert \text{dose=H}} &=& \beta_0+\beta_2 + \epsilon_i = \mu_3+\epsilon_i . \end{eqnarray*}\] met $\epsilon_i \sim N(0,\sigma^2)$ en \[ \mu_j = \text{E}\left[Y_i \mid \text{treatment group } j\right].\]

De oorspronkelijk nulhypothese $H_0:\mu_1=\mu_2=\mu_3$ kan equivalent geformuleerd worden als

[H_0: _1=_2=0]

Model laat toe om methoden van lineaire regressie te gebruiken voor meervoudig vergelijken van gemiddelden.

Parameterschatting van parameters, varianties en standard errors uit theorie van lineaire regressie
Inferentie: Betrouwbaarheidsintervallen, hypothesetests
Test $H_0: \beta_1=\beta_2=0$ met $F$-test.

2.2 Prostacyclin voorbeeld

model1 <- lm(prostac ~ dose, data = prostacyclin)
summary(model1)


Call:
lm(formula = prostac ~ dose, data = prostacyclin)

Residuals:
    Min      1Q  Median      3Q     Max 
-35.167 -17.117  -4.958  17.927  41.133 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   40.108      6.150   6.521 2.10e-07 ***
dose25         8.258      8.698   0.949    0.349    
dose50        43.258      8.698   4.974 1.99e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.3 on 33 degrees of freedom
Multiple R-squared:  0.458, Adjusted R-squared:  0.4252 
F-statistic: 13.94 on 2 and 33 DF,  p-value: 4.081e-05

3 Kwadratensommen en Anova

Zoals bij enkelvoudige regressie kwadratensom van regressie gebruiken bij het opstellen van de F-test. \[\begin{eqnarray*} \text{SSR}&=&\sum\limits_{i=1}^n (\hat Y_i -\bar Y)^2\\ &=& \sum\limits_{i=1}^n (\hat{g} (x_{i1},x_{i2}) - \bar Y)^2\\ &=& \sum\limits_{i=1}^n (\hat\beta_0+\hat\beta_1x_{i1}+\hat\beta_2x_{i2}) - \bar Y)^2\\ &=& \sum\limits_{i=1}^{n_1} (\hat\beta_0 - \bar Y)^2 +\sum\limits_{i=1}^{n_2} (\hat\beta_0 + \hat\beta_1 - \bar Y)^2+\sum\limits_{i=1}^{n_3} (\hat\beta_0 + \hat\beta_2 - \bar Y)^2\\ &=& \sum\limits_{i=1}^{n_1} (\bar Y_1- \bar Y)^2 +\sum\limits_{i=1}^{n_2} (\bar Y_2- \bar Y)^2+\sum\limits_{i=1}^{n_3} (\bar Y_3 - \bar Y)^2\\ \end{eqnarray*}\] met $n_1$, $n_2$ en $n_3$ het aantal waarnemingen in elke groep (hier $n-1=n_2=n_3=12$).

\[\begin{eqnarray*} \text{SSR}&=&\sum\limits_{i=1}^n (\hat Y_i -\bar Y)^2 \end{eqnarray*}\]

Kwadratensom opnieuw equivalent aan vergelijken van model (1) en een gereduceerd model met enkel een intercept.
Voor gereduceerd model zal intercept worden geschat door steekproefgemiddelde.
Deze kwadratensom heeft dus g-1=2 vrijheidsgraden:
- g=3 model parameters - 1 parameter voor steekproefgemiddelde of
- g=3 par. van complexe model - 1 par. van gereduceerde model.

3.1 Ontbinding van de Totale Kwadratensom

De conventie in een Anova setting is om de kwadratensom te noteren als SST, de kwadratensom van de behandeling (treatment) of als SSBetween.
De kwadratensom van de regressie geeft voor model (1) inderdaad de variabiliteit weer tussen de groepen.
De overeenkomstige gemiddelde kwadratensom wordt dan $\text{MST}=\text{SST}/2$.

De decompositie van SSTot wordt dan geschreven als [ = + ]

3.2 SSTot

3.3 SST

3.4 SSE

3.5 Anova test

Test $H_0: \beta_1=\beta_2=0$ met $F$-test. [ F = ]

met

$\text{MST}=\text{SST}/(g-1)$
$\text{MSE}=\text{SSE}/(n-p)$
Teststatistiek vergelijkt de variabiliteit verklaard door model (MST) met de residuele variabiliteit (MSE)

Variabiliteit tussen groepen (MST) tot variabiliteit binnen groepen (MSE)
onder $H_0$: $F \sim F_{g-1,n-g}$, met g=3.

3.6 Anova Tabel

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Treatment	d.f. SST	SST	MST	F-statistiek	p-waarde
Error	d.f. SSE	SSE	MSE

anova(model1)

Analysis of Variance Table

Response: prostac
          Df Sum Sq Mean Sq F value    Pr(>F)    
dose       2  12658  6329.0  13.944 4.081e-05 ***
Residuals 33  14979   453.9                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.6.1 F-verdeling met kritieke waarde ($\alpha$=5%) en geobserveerde F-statistiek voor prostacycline voorbeeld

3.6.2 F-verdeling met verschillend aantal vrijheidsgraden in de noemer en teller

3.6.3 Prostacyclin voorbeeld: welke groepen zijn verschillend?

summary(model1)


Call:
lm(formula = prostac ~ dose, data = prostacyclin)

Residuals:
    Min      1Q  Median      3Q     Max 
-35.167 -17.117  -4.958  17.927  41.133 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   40.108      6.150   6.521 2.10e-07 ***
dose25         8.258      8.698   0.949    0.349    
dose50        43.258      8.698   4.974 1.99e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.3 on 33 degrees of freedom
Multiple R-squared:  0.458, Adjusted R-squared:  0.4252 
F-statistic: 13.94 on 2 and 33 DF,  p-value: 4.081e-05

Met output van model kunnen we ook testen of de gemiddelde prostacycline concentratie verschillend is tussen de matige en lage dosis groep ($\beta_1$: dose25) en tussen de hoge en lage dosis groep ($\beta_2$: dose50).

De p-waarden houden geen rekening met het feit dat we meervoudig toetsen.

4 Post hoc analysis: Meerdere vergelijkingen van gemiddelden

4.1 Naïeve methode

In eerste deel van dit hoofdstuk hebben we $F$-test besproken voor het testen van

\[ H_0: \mu_1=\cdots = \mu_g \text{ versus } H_1: H_1: \exists\ j,k \in \{1,\ldots,g\} : \mu_j\neq\mu_k\]

Als we $H_0$ verwerpen besluiten we dat er minstens twee gemiddelden verschillen van elkaar.
De methode niet in staat om te identificeren welke gemiddelden van elkaar verschillen.

Een eerste, maar naïeve benadering: $H_0$ opsplitsen in partiële hypotheses \[H_{0jk}: \mu_j=\mu_k \text{ versus } H_{1jk}: \mu_j \neq \mu_k\]

partiële hypotheses testen met two-sample $t$-testen
Voor vergelijken van groep $j$ met groep $k$ wordt klassieke two-sample $t$-test onder gelijkheid van variantie gegeven door: \[T_{jk} = \frac{\bar{Y}_j-\bar{Y}_k}{S_p\sqrt{\frac{1}{n_j}+\frac{1}{n_k}}} \sim t_{n-2}\]

Met

$S_p^2$ de gepoolde variantieschatter is, \[S_p^2 = \frac{(n_j-1)S_j^2 + (n_k-1)S_k^2}{n_j+n_k-2}\]
met $S_j^2$ en $S_k^2$ de steekproefvarianties van respectievelijk de uitkomsten uit groep $j$ en $k$.

In ANOVA context veronderstellen we dat variantie in alle $g$ groepen dezelfde is: de residuele variantie $\sigma^2$.

Gebruik van $S_p^2$ is niet efficiënt omdat die niet van alle data gebruik maakt
Aan efficiëntie winnen door door MSE te gebruiken \[\text{MSE}= \sum_{j=1}^g \frac{(n_j-1)S_j^2}{n-g}\]
De $t$-testen worden dus best gebaseerd op \[T_{jk} = \frac{\bar{Y}_j-\bar{Y}_k}{\text{MSE}\sqrt{\frac{1}{n_j}+\frac{1}{n_k}}} \sim t_{n-g}.\]

with(
  prostacyclin,
  pairwise.t.test(prostac, dose, "none")
  )


    Pairwise comparisons using t tests with pooled SD 

data:  prostac and dose 

   10      25     
25 0.34927 -      
50 2e-05   0.00031

P value adjustment method: none

Het werken met $m$-testen op het $\alpha$ significantieniveau is echter een foute aanpak die de kans op een type I fout niet onder controle kan houden.

4.2 We tonen aan dat naïve methode niet werkt via simulatie

We simuleren uit een ANOVA model met $g=3$ groepen.
De gemiddelden in het ANOVA model zijn gelijk aan elkaar, zodat de nulhypothese \[H_0: \mu_1=\mu_2=\mu_3\] opgaat.
Voor iedere gesimuleerde dataset zijn er $m=3$ paarsgewijze two-sample $t$-testen
Zodra minstens één van de $p$-waarden kleiner is dan het significantieniveau $\alpha=5\%$, wordt de nulhypothese $H_0: \mu_1=\mu_2=\mu_3$ verworpen omdat er minstens twee gemiddelden verschillend zijn volgens de $t$-testen.
We rapporteren de relatieve frequentie van het verwerpen van de globale nulhypothese, meer bepaald de kans op een type I fout van de test voor $H_0: \mu_1=\mu_2=\mu_3$.

g <- 3 # aantal behandelingen (g=3)
ni <- 12 # aantal herhalingen in iedere groep
n <- g*ni # totaal aantal observaties
alpha <- 0.05 # significantieniveau van een individuele test
N <- 10000 #aantal simulaties
set.seed(302) #seed zodat resultaten exact geproduceerd kunnen worden
trt <- factor(rep(1:g, ni)) #factor
cnt <- 0 #teller voor aantal foutieve verwerpingen

for(i in 1:N) {
if (i %% 1000 == 0) cat(i, "/", N, "\n")
y <- rnorm(n)
tests <- pairwise.t.test(y, trt, "none")
verwerp <- min(tests$p.value, na.rm = T) < alpha
if(verwerp) cnt <- cnt+1
}

1000 / 10000 
2000 / 10000 
3000 / 10000 
4000 / 10000 
5000 / 10000 
6000 / 10000 
7000 / 10000 
8000 / 10000 
9000 / 10000 
10000 / 10000

cnt/N

[1] 0.1209

Kans op een type I fout gelijk is aan 12.1%
Is meer dan dubbel zo groot is dan vooropgestelde $\alpha=5$%.
Als we simulatiestudie herhalen met g = 5 groepen (i.e. 10 paarsgewijze t-testen) dan vinden we 28.0% in plaats van de gewenste 5%.
Deze simulaties illustreren het probleem van multipliciteit (Engels: multiplicity)
- klassieke $p$-waarden mogen enkel met het significantieniveau $\alpha$ vergeleken worden, indien het besluit op exact één $p$-waarde gebaseerd is.
- Finale besluit (al dan niet verwerpen van $H_0: \mu_1=\cdots =\mu_g$) gebaseerd op $m=g\times(g-1)/2$ $p$-waarden.
We bespreken eerst een uitbreiding van het begrip van type I fout en vervolgens introduceren we enkele oplossingen.

4.3 Family-wise error rate

Wanneer $m>1$ toetsen worden aangewend om 1 beslissing te vormen, is het noodzakelijk te corrigeren voor het risico op vals positieve resultaten (type I fout).
Meeste procedures voor meervoudig toetsen gaan ervan uit dat alle $m$ nulhypotheses waar zijn.
Er wordt dan geprobeerd om het risico op minstens 1 vals positief resultaat te controleren op experimentgewijs significantieniveau $\alpha_E$, typisch $\alpha_E=0.05$.
In de Engelstalige literatuur wordt het experimentgewijs significantieniveau family-wise error rate (FWER) genoemd.

4.4 Bonferroni correction

Bij het uitvoeren van $m$ onafhankelijke toetsen met elk significantieniveau $\alpha$, is \[\begin{eqnarray*} \alpha_E&=&\text{P}[\text{minstens 1 Type I fout}]\\ &=&1-(1-\alpha)^m \leq m\alpha \end{eqnarray*}\]

Als we 3 toetsen uitvoeren op het 5% significantieniveau is FWER $\approx 15\%$.
Door ze op het 1% significantieniveau uit te voeren, bekomen we FWER $\approx 5\%$.
De Bonferroni correctie houdt de FWER begrensd op $\alpha_E$ door \[\alpha=\alpha_E/m\] te kiezen voor het uitvoeren van de $m$ paarsgewijze vergelijkingen.

Als alternatieve methode kunnen we ook

aangepaste p-waarden rapporteren zodat we deze met het experimentgewijze $\alpha_E$ niveau kunnen vergelijken: \[\tilde{p}=min(m\times p,1)\]
$(1-\alpha_E/m)100\%$ betrouwbaarheidsintervallen rapporteren.

4.4.1 Prostacyclin voorbeeld

with(
  prostacyclin,
  pairwise.t.test(
    prostac,
    dose,
    p.adjust.method="bonferroni")
  )


    Pairwise comparisons using t tests with pooled SD 

data:  prostac and dose 

   10      25     
25 1.00000 -      
50 6e-05   0.00094

P value adjustment method: bonferroni

Conclusies blijven gelijk behalve dat FWER nu gecontroleerd is $\alpha_E=5\%$ en $\tilde{p}$-waarden zijn factor 3 groter

Zelfde analyse kan via multcomp R package ontwikkeld voor multipliciteit in lineaire modellen.

library(multcomp)

library(multcomp)
model1.mcp <- glht(model1, linfct = mcp(dose = "Tukey"))
summary(model1.mcp, test = adjusted("bonferroni"))


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = prostac ~ dose, data = prostacyclin)

Linear Hypotheses:
             Estimate Std. Error t value Pr(>|t|)    
25 - 10 == 0    8.258      8.698   0.949 1.000000    
50 - 10 == 0   43.258      8.698   4.974 5.98e-05 ***
50 - 25 == 0   35.000      8.698   4.024 0.000943 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- bonferroni method)

Om Bonferroni aangepaste betrouwbaarheidsintervallen te verkrijgen moeten we eerst zelf functie definiëren in R om bonferroni kritische waarde te bepalen.

Bonferonni-betrouwbaarheidsintervallen worden niet geïmplementeerd omdat er betere methoden bestaan voor meervoudige tests.
De onderstaande functie is hier alleen voor de volledigheid toegevoegd, maar we zullen over het algemeen de standaardmethode gebruiken voor meervoudige tests in multcomp.

calpha_bon_t <- function(object, level)
{
  abs(
    qt(
      (1-level)/2/nrow(object$linfct),
      object$df)
    )
}

confint(model1.mcp, calpha = calpha_bon_t)


     Simultaneous Confidence Intervals

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = prostac ~ dose, data = prostacyclin)

Quantile = 2.5222
95% confidence level
 

Linear Hypotheses:
             Estimate lwr      upr     
25 - 10 == 0   8.2583 -13.6790  30.1957
50 - 10 == 0  43.2583  21.3210  65.1957
50 - 25 == 0  35.0000  13.0626  56.9374

4.4.2 Simulatie om Bonferroni methode te evalueren

g <- 3 # aantal behandelingen (g=3)
ni <- 12 # aantal herhalingen in iedere groep
n <- g*ni # totaal aantal observaties
alpha <- 0.05 # significantieniveau van een individuele test
N <- 10000 #aantal simulaties
set.seed(302) #seed zodat resultaten exact geproduceerd kunnen worden
trt <- factor(rep(1:g, ni)) #factor
cnt <- 0 #teller voor aantal foutieve verwerpingen

for(i in 1:N) {
if (i %% 1000 == 0) cat(i, "/", N, "\n")
y <- rnorm(n)
tests <- pairwise.t.test(y, trt, "bonferroni")
verwerp <- min(tests$p.value, na.rm = T) < alpha
if(verwerp) cnt <- cnt+1
}

1000 / 10000 
2000 / 10000 
3000 / 10000 
4000 / 10000 
5000 / 10000 
6000 / 10000 
7000 / 10000 
8000 / 10000 
9000 / 10000 
10000 / 10000

cnt/N

[1] 0.0457

We vinden een FWER van 4.6% (een beetje conservatief)
Bij simulaties voor $g=5$ groepen, vinden we een FWER van $4.1\%$ (conservatiever).
Door Bonferroni correctie is kans op minstens één vals positief resultaat $< \alpha_E$.
Power verlies aangezien werkelijke niveau lager is dan het vooropgestelde 5% experimentsgewijs significantieniveau.

4.5 Tukey Methode

Minder conservatieve methode
Implementatie benadert de nuldistributie van de posthoc test d.m.v. simulaties.
Resultaten kunnen lichtjes verschillen wanneer je posthoc analyse opnieuw uitvoert.
Details van de methode vallen buiten het bestek van deze cursus.
Implementatie in multcomp package:
- aangepaste p-waarden
- aangepaste betrouwbaarheidsintervallen
- Niet nodig om functies te definiëren voor Tukey gecorrigeerde BIs

4.5.1 Captopril voorbeeld

model1.mcp <- glht(model1, linfct = mcp(dose = "Tukey"))
summary(model1.mcp)


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = prostac ~ dose, data = prostacyclin)

Linear Hypotheses:
             Estimate Std. Error t value Pr(>|t|)    
25 - 10 == 0    8.258      8.698   0.949 0.613390    
50 - 10 == 0   43.258      8.698   4.974  < 1e-04 ***
50 - 25 == 0   35.000      8.698   4.024 0.000835 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

confint(model1.mcp)


     Simultaneous Confidence Intervals

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = prostac ~ dose, data = prostacyclin)

Quantile = 2.4539
95% family-wise confidence level
 

Linear Hypotheses:
             Estimate lwr      upr     
25 - 10 == 0   8.2583 -13.0849  29.6016
50 - 10 == 0  43.2583  21.9151  64.6016
50 - 25 == 0  35.0000  13.6567  56.3433

model1.mcp |>
  confint() |>
  plot()

4.5.2 Evalueer Tukey methode in een simulatiestudie

g <- 3 # aantal behandelingen (g=3)
ni <- 12 # aantal herhalingen in iedere groep
n <- g*ni # totaal aantal observaties
alpha <- 0.05 # significantieniveau van een individuele test
N <- 10000 #aantal simulaties
set.seed(302) #seed zodat resultaten exact geproduceerd kunnen worden
trt <- factor(rep(1:g, ni)) #factor
cnt <- 0 #teller voor aantal foutieve verwerpingen

for(i in 1:N) {
if (i %% 1000 == 0) cat(i, "/", N, "\n")
y <- rnorm(n)
m <- lm(y ~ trt)
m.mcp <- glht(m, linfct = mcp(trt = "Tukey"))
tests <- summary(m.mcp)$test
verwerp <- min(
  as.numeric(tests$pvalues),
  na.rm=T) < alpha
if(verwerp) cnt <- cnt+1
}

1000 / 10000 
2000 / 10000 
3000 / 10000 
4000 / 10000 
5000 / 10000 
6000 / 10000 
7000 / 10000 
8000 / 10000 
9000 / 10000 
10000 / 10000

cnt/N

[1] 0.0503

De methode wordt mooi op het 5% experimentsgewijs significantieniveau gecontroleerd.

5 Conclusies: Prostacycline voorbeeld

Volledige analyse voor voorbeeld prostacycline

Anova vóór posthoc-tests: F-test heeft een hogere power dan paarsgewijze t-test
- F-test gebruikt alle gegevens
- Voor F-test hoeven we niet te corrigeren voor meervoudige testen: er wordt één test uitgevoerd voor de algemene omnibushypothese

model1 <- lm(prostac ~ dose, data = prostacyclin)
anova(model1)

Analysis of Variance Table

Response: prostac
          Df Sum Sq Mean Sq F value    Pr(>F)    
dose       2  12658  6329.0  13.944 4.081e-05 ***
Residuals 33  14979   453.9                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model1.mcp <- glht(model1, linfct = mcp(dose = "Tukey"))
summary(model1.mcp)


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = prostac ~ dose, data = prostacyclin)

Linear Hypotheses:
             Estimate Std. Error t value Pr(>|t|)    
25 - 10 == 0    8.258      8.698   0.949 0.613433    
50 - 10 == 0   43.258      8.698   4.974  < 1e-04 ***
50 - 25 == 0   35.000      8.698   4.024 0.000922 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

confint(model1.mcp)


     Simultaneous Confidence Intervals

Multiple Comparisons of Means: Tukey Contrasts


Fit: lm(formula = prostac ~ dose, data = prostacyclin)

Quantile = 2.4526
95% family-wise confidence level
 

Linear Hypotheses:
             Estimate lwr      upr     
25 - 10 == 0   8.2583 -13.0736  29.5902
50 - 10 == 0  43.2583  21.9264  64.5902
50 - 25 == 0  35.0000  13.6681  56.3319

Er is een extreem significant effect van arachidonzuur op de gemiddelde bloedconcentratie van prostacycline bij ratten ($p<0.001$). De gemiddelde prostacyclineconcentratie is hoger in de hoge dosisgroep dan in de lage en matige dosisgroep (beide p-waardes zijn kleiner dan $p<0.001$).
De gemiddelde concentratie in de hoge dosisgroep is 43.3ng/ml (95% CI [21.9,64.6]ng/ml) en 35ng/ml (95% BI [13.6,56.4]ng/ml) hoger dan in de lage en matige dosisgroep, respectievelijk.
Het verschil in gemiddelde prostacyclineconcentratie tussen de matige en lage dosisgroep is niet significant (p=0.61). (Alle p-waardes en betrouwbaarheidsintervallen voor post-hoc-tests worden gecorrigeerd voor meervoudige tests met behulp van de Tukey-methode).

Merk op dat het belangrijk is om ook de niet significante resultaten te vermelden!

--- title: "7. Analysis of Variance" author: "Lieven Clement" date: "statOmics, Ghent University (https://statomics.github.io)" format: html: toc: true number-sections: true theme: cosmo highlight-style: tango code-tools: true embed-resources: true --- ```{r setup, include=FALSE} knitr::opts_chunk$set(include = TRUE, comment = NA, echo = TRUE, message = FALSE, warning = FALSE) library(Rmisc) library(tidyverse) ``` # Prostacycline dataset - Effect van arachidonzuur op het prostacycline niveau in het bloedplasma. - 3 verschillende concentraties van arachidonzuur: - laag (L, 10 eenheden) - gemiddeld (M, 25 eenheden) - hoge dosis (H, 50 eenheden) - Prostacycline concentratie in bloed plasma via gecalibreerde elisa fluorescentie meting - 12 ratten worden at random toegekend aan elke behandelingsgroep. - Factoriële proef, *volledige gerandomiseerde proefopzet*, *"completely randomized design" CRD*. ```{r} prostacyclin <- read_tsv("https://raw.githubusercontent.com/statOmics/sbc/master/data/prostacyclin.txt") prostacyclin <- prostacyclin |> mutate(dose = as.factor(prostacyclin$dose)) head(prostacyclin) ``` *** ## Data exploration ```{r} prostacyclin |> ggplot(aes(x = dose, y = prostac, fill = dose)) + geom_boxplot() + geom_point(position = "jitter") + ylab("prostacyclin (ng/ml)") prostacyclin |> ggplot(aes(sample = prostac)) + geom_qq() + geom_qq_line() + facet_grid(~ dose) ``` Data in de drie groepen lijkt normaal verdeeld en de variantie is ongeveer gelijk: \[Y_i \vert \text{groep j} \sim N(\mu_j,\sigma^2),\] met $j= \text{1, 2, 3}$ ## Onderzoeksvraag Vraagstelling kan vertaald worden in volgende hypotheses - $H_0$: De arachidonzuurconcentratie heeft geen effect op het gemiddelde prostacycline niveau bij ratten \[ H_0:\mu_1=\mu_2 = \mu_3 \] - $H_1$: De arachidonzuurconcentratie heeft een effect op het gemiddelde prostacycline niveau bij ratten. Wat betekent dat minstens twee gemiddelden verschillend zijn: \[H_1: \exists\ j,k \in \{1,\ldots,g\} : \mu_j\neq\mu_k\] **naïeve benadering**: nulhypothese op splitsen in **partiële hypotheses** \[ H_{0jk}: \mu_j=\mu_k \text{ versus } H_{1jk}: \mu_j \neq \mu_k \] Elk van deze partiële hypotheses testen met two-sample $t$-testen $\rightarrow$ Probleem van meervoudig toetsen + verlies van power. $\rightarrow$ $H_0:\mu_1=\mu_2=\mu_3$ testen met **1 enkele test**. # Analyse van variantie - Correcte oplossing voor het testprobleem: **variantie-analyse**, afgekort door ANOVA (ANalysis Of VAriance) - We leiden de methode af voor de meest eenvoudige uitbreiding met 3 groepen (prostacycline voorbeeld) - Data modelleren a.d.h.v een lineair model door gebruik te maken van dummy variabelen. - 1 dummy variable minder nodig hebben dan het aantal groepen. Hier dus 2 dummy variabelen. - De veralgemening naar g groepen $g>3$ is triviaal (extra dummy variabelen) ## Model \begin{eqnarray} Y_i &=& g(x_{i1},x_{i2}) + \epsilon_i\\ Y_i &=& \beta_0+\beta_1 x_{i1} +\beta_2 x_{i2} +\epsilon_i \end{eqnarray} - $Y_i$ de uitkomst voor observatie $i$ ($i=1,\ldots, n$) \vspace{7pt} - $\epsilon_i\text{i.i.d.} N(0,\sigma^2)$ \vspace{7pt} - en dummyvariabelen $$x_{i1} = \left\{ \begin{array}{ll} 1 & \text{ als observatie $i$ tot de middelste dosisgroep behoort (M)} \\ 0 & \text{ als observatie $i$ behoort tot een andere dosisgroep} \end{array}\right.$$ $$x_{i2} = \left\{ \begin{array}{ll} 1 & \text{ als observatie $i$ behoort tot de groep met hoge doses (H)} \\ 0 & \text{als observatie $i$ behoort tot een andere dosisgroep} \end{array}\right. .$$ \vspace{7pt} - Lage dosis groep (L) with $x_{i1}=x_{i2}=0$ is *referentie groep* Regressiemodel kan worden herschreven als een model voor elke groep: \vspace{-20pt} \begin{eqnarray*} Y_{i\vert \text{dose=L}} &=& \beta_0+\epsilon_i \\ Y_{i\vert \text{dose=M}} &=& \beta_0+\beta_1+ \epsilon_i \\ Y_{i\vert \text{dose=H}} &=& \beta_0+\beta_2 + \epsilon_i \end{eqnarray*} met $\epsilon_i \sim N(0,\sigma^2)$ \vspace{10pt} Interpretatie van model parameters: \vspace{-20pt} \begin{eqnarray*} \beta_0 &=& \text{E}\left[Y_i \mid \text{Behandeling met lage dosis groep L}\right] \\ \beta_1 &=& (\beta_0+\beta_1)-\beta_0 = \text{E}\left[Y_i \mid \text{behandeling M}\right] - \text{E}\left[Y_i \mid \text{behandeling L}\right] \\ \beta_2 &=& (\beta_0+\beta_2)-\beta_0 = \text{E}\left[Y_i \mid \text{behandeling H}\right]-\text{E}\left[Y_i \mid \text{behandeling L}\right]. \end{eqnarray*} 1. $\beta_0$ is gemiddelde uitkomst in groep L 2. $\beta_1$ is effect (verschil in gemiddelde concentratie) van groep M t.o.v. groep L 3. $\beta_2$ is effect van groep H t.o.v. groep L We herformuleren het model door $\mu$-notaties te gebruiken: \vspace{-7pt} \begin{eqnarray*} Y_{i\vert \text{dose=L}} &=& \beta_0+\epsilon_i = \mu_1+\epsilon_i \\ Y_{i\vert \text{dose=M}} &=& \beta_0+\beta_1+ \epsilon_i = \mu_2+\epsilon_i \\ Y_{i\vert \text{dose=H}} &=& \beta_0+\beta_2 + \epsilon_i = \mu_3+\epsilon_i . \end{eqnarray*} met $\epsilon_i \sim N(0,\sigma^2)$ en $$ \mu_j = \text{E}\left[Y_i \mid \text{treatment group } j\right].$$ De oorspronkelijk nulhypothese $H_0:\mu_1=\mu_2=\mu_3$ kan equivalent geformuleerd worden als \[H_0: \beta_1=\beta_2=0\] Model laat toe om methoden van lineaire regressie te gebruiken voor meervoudig vergelijken van gemiddelden. - Parameterschatting van parameters, varianties en standard errors uit theorie van lineaire regressie - Inferentie: Betrouwbaarheidsintervallen, hypothesetests - Test $H_0: \beta_1=\beta_2=0$ met $F$-test. ## Prostacyclin voorbeeld ```{r} model1 <- lm(prostac ~ dose, data = prostacyclin) summary(model1) ``` # Kwadratensommen en Anova Zoals bij enkelvoudige regressie kwadratensom van regressie gebruiken bij het opstellen van de F-test. \vspace{-10pt} \begin{eqnarray*} \text{SSR}&=&\sum\limits_{i=1}^n (\hat Y_i -\bar Y)^2\\ &=& \sum\limits_{i=1}^n (\hat{g} (x_{i1},x_{i2}) - \bar Y)^2\\ &=& \sum\limits_{i=1}^n (\hat\beta_0+\hat\beta_1x_{i1}+\hat\beta_2x_{i2}) - \bar Y)^2\\ &=& \sum\limits_{i=1}^{n_1} (\hat\beta_0 - \bar Y)^2 +\sum\limits_{i=1}^{n_2} (\hat\beta_0 + \hat\beta_1 - \bar Y)^2+\sum\limits_{i=1}^{n_3} (\hat\beta_0 + \hat\beta_2 - \bar Y)^2\\ &=& \sum\limits_{i=1}^{n_1} (\bar Y_1- \bar Y)^2 +\sum\limits_{i=1}^{n_2} (\bar Y_2- \bar Y)^2+\sum\limits_{i=1}^{n_3} (\bar Y_3 - \bar Y)^2\\ \end{eqnarray*} met $n_1$, $n_2$ en $n_3$ het aantal waarnemingen in elke groep (hier $n-1=n_2=n_3=12$). \begin{eqnarray*} \text{SSR}&=&\sum\limits_{i=1}^n (\hat Y_i -\bar Y)^2 \end{eqnarray*} - Kwadratensom opnieuw equivalent aan vergelijken van model (1) en een gereduceerd model met enkel een intercept. - Voor gereduceerd model zal intercept worden geschat door steekproefgemiddelde. - Deze kwadratensom heeft dus g-1=2 vrijheidsgraden: - g=3 model parameters - 1 parameter voor steekproefgemiddelde of - g=3 par. van complexe model - 1 par. van gereduceerde model. ## Ontbinding van de Totale Kwadratensom - De conventie in een Anova setting is om de kwadratensom te noteren als SST, de **kwadratensom van de behandeling (treatment)** of als SSBetween. - De kwadratensom van de regressie geeft voor model (1) inderdaad de variabiliteit weer tussen de groepen. - De overeenkomstige gemiddelde kwadratensom wordt dan $\text{MST}=\text{SST}/2$. De decompositie van SSTot wordt dan geschreven als \[ \text{SSTot} = \text{SST} + \text{SSE} \] ## SSTot ```{r echo=FALSE, warning=FALSE} par(mfrow=c(1,2)) jitIk=runif(36,-.2,.2)+rep(1:3,each=12) plot(prostac~dose,data=prostacyclin,xlab="Arachidonic acid dose ",ylab="Prostacyclin (ng/ml)",cex.axis=1.5,cex.lab=1.5,cex.main=1.5) points(jitIk,prostacyclin$prostac,col=col,pch=19) points(jitIk,prostacyclin$prostac,col=4) points(1:3,predict(model1,data.frame(dose=factor(c(10,25,50)))),pch=17,col=c("bisque","coral","darkcyan"),cex=1.5) points(1:3,predict(model1,data.frame(dose=factor(c(10,25,50)))),pch=2,col=1,cex=1.5) abline(h=mean(prostacyclin$prostac),lty=1) for (i in 1:36) lines(rep(jitIk[i],2),c(mean(prostacyclin$prostac),prostacyclin$prostac[i]),col=4,lty=2) jitIk=runif(36,-.2,.2)+rep(1:3,each=12) plot(rep(1,36),prostacyclin$prostac-mean(prostacyclin$prostac),xaxt="none",ylab="Deviations",cex.lab=1.5,cex.main=1.5,cex.axis=1.5,col=as.character(prostacyclin$col),xlim=c(1,3),pch=19,xlab="") points(rep(1,36),prostacyclin$prostac-mean(prostacyclin$prostac),pch=1,col=4) axis(at=1:3,labels=c(expression(paste(y[i]," - ",bar(y))),expression(paste(bar(y)[j]," - ",bar(y))),expression(paste(y[i]," - ",bar(y)[j]))),side=1,cex.axis=1.5) ``` ## SST ```{r echo=FALSE, warning=FALSE} par(mfrow=c(1,2)) plot(prostac~dose,data=prostacyclin,xlab="Arachidonic acid dose ",ylab="Prostacyclin (ng/ml)",cex.axis=1.5,cex.lab=1.5,cex.main=1.5) points(jitIk,prostacyclin$prostac,col=col,pch=19) points(jitIk,prostacyclin$prostac,col=4) points(1:3,predict(model1,data.frame(dose=factor(c(10,25,50)))),pch=17,col=c("bisque","coral","darkcyan"),cex=1.5) points(1:3,predict(model1,data.frame(dose=factor(c(10,25,50)))),pch=2,col=2,cex=1.5) abline(h=mean(prostacyclin$prostac),lty=1) for (i in 1:3) lines(rep(i,2),c(mean(prostacyclin$prostac),predict(model1,data.frame(dose=levels(prostacyclin$dose)[i]))),col=2,lty=2) plot(rep(1,36),prostacyclin$prostac-mean(prostacyclin$prostac),xaxt="none",ylab="Deviations",cex.lab=1.5,cex.main=1.5,cex.axis=1.5,col=as.character(prostacyclin$col),xlim=c(1,3),pch=19,xlab="") points(rep(1,36),prostacyclin$prostac-mean(prostacyclin$prostac),pch=1,col=4) points(rep(2,3),predict(model1,data.frame(dose=factor(c(10,25,50))))-mean(prostacyclin$prostac),pch=17,col=unique(prostacyclin$col),cex=1.5) points(rep(2,3),predict(model1,data.frame(dose=factor(c(10,25,50))))-mean(prostacyclin$prostac),pch=2,cex=1.5,col=2) axis(at=1:3,labels=c(expression(paste(y[i]," - ",bar(y))),expression(paste(bar(y)[j]," - ",bar(y))),expression(paste(y[i]," - ",bar(y)[j]))),side=1,cex.axis=1.5) ``` ## SSE ```{r echo=FALSE, warning=FALSE} par(mfrow=c(1,2)) plot(prostac~dose,data=prostacyclin,xlab="Arachidonic acid dose ",ylab="Prostacyclin (ng/ml)",cex.axis=1.5,cex.lab=1.5,cex.main=1.5) points(jitIk,prostacyclin$prostac,col=col,pch=19) points(jitIk,prostacyclin$prostac,col=1) points(1:3,predict(model1,data.frame(dose=factor(c(10,25,50)))),pch=17,col=c("bisque","coral","darkcyan"),cex=1.5) points(1:3,predict(model1,data.frame(dose=factor(c(10,25,50)))),pch=2,col=2,cex=1.5) for (i in 1:3) lines(c(i-.2,i+.2),rep(predict(model1,data.frame(dose=levels(prostacyclin$dose)[i])),2),col=c("bisque","coral","darkcyan")[i]) abline(h=mean(prostacyclin$prostac),lty=1) for (i in 1:36) lines(rep(jitIk[i],2),c(prostacyclin$prostac[i],model1$fitted[i]),col=1,lty=2) plot(rep(1,36),prostacyclin$prostac-mean(prostacyclin$prostac),xaxt="none",ylab="Deviations",cex.lab=1.5,cex.main=1.5,cex.axis=1.5,col=as.character(prostacyclin$col),xlim=c(1,3),pch=19,xlab="") points(rep(1,36),prostacyclin$prostac-mean(prostacyclin$prostac),pch=1,col=4) points(rep(2,3),predict(model1,data.frame(dose=factor(c(10,25,50))))-mean(prostacyclin$prostac),pch=17,col=unique(prostacyclin$col),cex=1.5) points(rep(2,3),predict(model1,data.frame(dose=factor(c(10,25,50))))-mean(prostacyclin$prostac),pch=2,col=2,cex=1.5) points(rep(3,36),model1$res,pch=19,col=as.character(prostacyclin$col)) points(rep(3,36),model1$res,pch=1) axis(at=1:3,labels=c(expression(paste(y[i]," - ",bar(y))),expression(paste(bar(y)[j]," - ",bar(y))),expression(paste(y[i]," - ",bar(y)[j]))),side=1,cex.axis=1.5) ``` ## Anova test Test $H_0: \beta_1=\beta_2=0$ met $F$-test. \[ F = \frac{\text{MST}}{\text{MSE}} \] met - $\text{MST}=\text{SST}/(g-1)$ \vspace{10pt} - $\text{MSE}=\text{SSE}/(n-p)$ \vspace{10pt} - Teststatistiek vergelijkt de variabiliteit verklaard door model (MST) met de residuele variabiliteit (MSE) of - Variabiliteit tussen groepen (MST) tot variabiliteit binnen groepen (MSE) \vspace{10pt} - onder $H_0$: $F \sim F_{g-1,n-g}$, met g=3. ## Anova Tabel | |Df|Sum Sq|Mean Sq|F value|Pr(>F)| |---|---|---|---|---|---| |Treatment|d.f. SST|SST|MST|F-statistiek|p-waarde| |Error|d.f. SSE|SSE|MSE| | | ```{r} anova(model1) ``` ### F-verdeling met kritieke waarde ($\alpha$=5%) en geobserveerde F-statistiek voor prostacycline voorbeeld ```{r prostacF, out.width='100%', fig.asp=.8, fig.align='center',echo=FALSE} grid <- seq(0,17,.01) df1=anova(model1)[1,1] df2=anova(model1)[2,1] fval=anova(model1)[1,4] crit=qf(0.95,df1,df2) reject=c(crit,grid[which(grid>crit)]) accept=c(grid[which(grid<crit)],crit) plot(grid,df(grid,df1,df2),type="l",ylab="Density",xlab="F-statistic",cex.axis=1.5,cex.lab=1.5) polygon(c(0,accept,crit,0),c(0,df(accept,df1,df2),0,0),col="blue",border="blue") text(crit/2,.97,labels="accept\n95%",col="blue",cex=1.5) polygon(c(crit,reject,15,crit),c(0,df(reject,df1,df2),0,0),col="red",border="red") abline(v=crit,col="red",lwd=2) text(crit+(fval-crit)/2,.97,labels="reject\n5%",col="red",cex=1.5) text(pos=4,crit,df(crit,2,33),labels=paste0("F(0.05,",df1,",",df2,")"),col="red",cex=1.5) text(pos=4,fval,df(crit,df1,df2),labels=paste0("f=",round(fval,1)),col="darkorange",cex=1.5) abline(v=fval,col="darkorange",lwd=2,lty=2) text(15.5,.97,labels=paste0("p-value\n",format(anova(model1)[1,5],digits=2)),col="darkorange",cex=1.5) arrows(x0=17.5,x1=fval,y0=.9,y1=.9,col="darkorange") ``` ### F-verdeling met verschillend aantal vrijheidsgraden in de noemer en teller ```{r ftheo, out.width='100%', fig.asp=.8, fig.align='center',echo=FALSE} plot(grid,df(grid,1,5),type="l",ylab="Density",xlab="F-statistic",xlim=c(0,5),ylim=c(0,1.5),lwd=2,cex.axis=1.5,cex.lab=1.5) lines(grid,df(grid,5,5),type="l",col=2,lwd=2) lines(grid,df(grid,10,30),type="l",col=3,lwd=2) lines(grid,df(grid,20,30),type="l",col=4,lwd=2) lines(grid,df(grid,50,50),type="l",col=5,lwd=2) legend("topright",lty=1,col=c(1,2,3,4,5),legend=c("F(1,5)","F(5,5)","F(10,30)","F(20,30)","F(50,50)"),lwd=2,cex=1.5) ``` ### Prostacyclin voorbeeld: welke groepen zijn verschillend? ```{r} summary(model1) ``` Met output van model kunnen we ook testen of de gemiddelde prostacycline concentratie verschillend is tussen de matige en lage dosis groep ($\beta_1$: dose25) en tussen de hoge en lage dosis groep ($\beta_2$: dose50). De p-waarden houden geen rekening met het feit dat we meervoudig toetsen. # Post hoc analysis: Meerdere vergelijkingen van gemiddelden ## Naïeve methode In eerste deel van dit hoofdstuk hebben we $F$-test besproken voor het testen van $$ H_0: \mu_1=\cdots = \mu_g \text{ versus } H_1: H_1: \exists\ j,k \in \{1,\ldots,g\} : \mu_j\neq\mu_k$$ - Als we $H_0$ verwerpen besluiten we dat er minstens twee gemiddelden verschillen van elkaar. - De methode niet in staat om te identificeren welke gemiddelden van elkaar verschillen. Een eerste, maar naïeve benadering: $H_0$ opsplitsen in partiële hypotheses $$H_{0jk}: \mu_j=\mu_k \text{ versus } H_{1jk}: \mu_j \neq \mu_k$$ - partiële hypotheses testen met two-sample $t$-testen - Voor vergelijken van groep $j$ met groep $k$ wordt klassieke two-sample $t$-test onder gelijkheid van variantie gegeven door: $$T_{jk} = \frac{\bar{Y}_j-\bar{Y}_k}{S_p\sqrt{\frac{1}{n_j}+\frac{1}{n_k}}} \sim t_{n-2}$$ Met - $S_p^2$ de gepoolde variantieschatter is, $$S_p^2 = \frac{(n_j-1)S_j^2 + (n_k-1)S_k^2}{n_j+n_k-2}$$ - met $S_j^2$ en $S_k^2$ de steekproefvarianties van respectievelijk de uitkomsten uit groep $j$ en $k$. *** In ANOVA context veronderstellen we dat variantie in **alle** $g$ groepen dezelfde is: de residuele variantie $\sigma^2$. - Gebruik van $S_p^2$ is niet efficiënt omdat die niet van alle data gebruik maakt - Aan efficiëntie winnen door door MSE te gebruiken $$\text{MSE}= \sum_{j=1}^g \frac{(n_j-1)S_j^2}{n-g}$$ - De $t$-testen worden dus best gebaseerd op $$T_{jk} = \frac{\bar{Y}_j-\bar{Y}_k}{\text{MSE}\sqrt{\frac{1}{n_j}+\frac{1}{n_k}}} \sim t_{n-g}.$$ *** ```{r} with( prostacyclin, pairwise.t.test(prostac, dose, "none") ) ``` Het werken met $m$-testen op het $\alpha$ significantieniveau is echter een foute aanpak die de kans op een type I fout niet onder controle kan houden. ## We tonen aan dat naïve methode niet werkt via simulatie 1. We simuleren uit een ANOVA model met $g=3$ groepen. 2. De gemiddelden in het ANOVA model zijn gelijk aan elkaar, zodat de nulhypothese $$H_0: \mu_1=\mu_2=\mu_3$$ opgaat. 3. Voor iedere gesimuleerde dataset zijn er $m=3$ paarsgewijze two-sample $t$-testen 4. Zodra minstens één van de $p$-waarden kleiner is dan het significantieniveau $\alpha=5\%$, wordt de nulhypothese $H_0: \mu_1=\mu_2=\mu_3$ verworpen omdat er minstens twee gemiddelden verschillend zijn volgens de $t$-testen. 5. We rapporteren de relatieve frequentie van het verwerpen van de globale nulhypothese, meer bepaald de kans op een type I fout van de test voor $H_0: \mu_1=\mu_2=\mu_3$. ```{r} g <- 3 # aantal behandelingen (g=3) ni <- 12 # aantal herhalingen in iedere groep n <- g*ni # totaal aantal observaties alpha <- 0.05 # significantieniveau van een individuele test N <- 10000 #aantal simulaties set.seed(302) #seed zodat resultaten exact geproduceerd kunnen worden trt <- factor(rep(1:g, ni)) #factor cnt <- 0 #teller voor aantal foutieve verwerpingen for(i in 1:N) { if (i %% 1000 == 0) cat(i, "/", N, "\n") y <- rnorm(n) tests <- pairwise.t.test(y, trt, "none") verwerp <- min(tests$p.value, na.rm = T) < alpha if(verwerp) cnt <- cnt+1 } cnt/N ``` *** - Kans op een type I fout gelijk is aan `r round(cnt/N,3)*100`% - Is meer dan dubbel zo groot is dan vooropgestelde $\alpha=5$%. - Als we simulatiestudie herhalen met g = 5 groepen (i.e. 10 paarsgewijze t-testen) dan vinden we 28.0% in plaats van de gewenste 5%. - Deze simulaties illustreren het probleem van **multipliciteit** (Engels: *multiplicity*) - klassieke $p$-waarden mogen enkel met het significantieniveau $\alpha$ vergeleken worden, indien het besluit op exact één $p$-waarde gebaseerd is. - Finale besluit (al dan niet verwerpen van $H_0: \mu_1=\cdots =\mu_g$) gebaseerd op $m=g\times(g-1)/2$ $p$-waarden. - We bespreken eerst een uitbreiding van het begrip van type I fout en vervolgens introduceren we enkele oplossingen. ## Family-wise error rate - Wanneer $m>1$ toetsen worden aangewend om 1 beslissing te vormen, is het noodzakelijk te corrigeren voor het risico op vals positieve resultaten (type I fout). \vspace{10pt} - Meeste procedures voor meervoudig toetsen gaan ervan uit dat *alle $m$ nulhypotheses waar* zijn. \vspace{10pt} - Er wordt dan geprobeerd om het *risico op minstens 1 vals positief resultaat* te controleren op **experimentgewijs significantieniveau $\alpha_E$**, typisch $\alpha_E=0.05$. \vspace{10pt} - In de Engelstalige literatuur wordt het experimentgewijs significantieniveau *family-wise error rate (FWER)* genoemd. ## Bonferroni correction Bij het uitvoeren van $m$ onafhankelijke toetsen met elk significantieniveau $\alpha$, is \begin{eqnarray*} \alpha_E&=&\text{P}[\text{minstens 1 Type I fout}]\\ &=&1-(1-\alpha)^m \leq m\alpha \end{eqnarray*} - Als we 3 toetsen uitvoeren op het 5\% significantieniveau is FWER $\approx 15\%$. \vspace{10pt} - Door ze op het 1\% significantieniveau uit te voeren, bekomen we FWER $\approx 5\%$. \vspace{10pt} - De Bonferroni correctie houdt de FWER begrensd op $\alpha_E$ door $$\alpha=\alpha_E/m$$ te kiezen voor het uitvoeren van de $m$ paarsgewijze vergelijkingen. *** Als alternatieve methode kunnen we ook 1. *aangepaste p-waarden* rapporteren zodat we deze met het experimentgewijze $\alpha_E$ niveau kunnen vergelijken: $$\tilde{p}=min(m\times p,1)$$ 2. $(1-\alpha_E/m)100\%$ betrouwbaarheidsintervallen rapporteren. ### Prostacyclin voorbeeld ```{r} with( prostacyclin, pairwise.t.test( prostac, dose, p.adjust.method="bonferroni") ) ``` - Conclusies blijven gelijk behalve dat FWER nu gecontroleerd is $\alpha_E=5\%$ en $\tilde{p}$-waarden zijn factor 3 groter *** Zelfde analyse kan via `multcomp` R package ontwikkeld voor multipliciteit in lineaire modellen. ```{r} library(multcomp) ``` ```{r} library(multcomp) model1.mcp <- glht(model1, linfct = mcp(dose = "Tukey")) summary(model1.mcp, test = adjusted("bonferroni")) ``` *** Om Bonferroni aangepaste betrouwbaarheidsintervallen te verkrijgen moeten we eerst zelf functie definiëren in R om bonferroni kritische waarde te bepalen. - Bonferonni-betrouwbaarheidsintervallen worden niet geïmplementeerd omdat er betere methoden bestaan voor meervoudige tests. - De onderstaande functie is hier alleen voor de volledigheid toegevoegd, maar we zullen over het algemeen de standaardmethode gebruiken voor meervoudige tests in multcomp. ```{r} calpha_bon_t <- function(object, level) { abs( qt( (1-level)/2/nrow(object$linfct), object$df) ) } ``` ```{r} confint(model1.mcp, calpha = calpha_bon_t) ``` ### Simulatie om Bonferroni methode te evalueren ```{r} g <- 3 # aantal behandelingen (g=3) ni <- 12 # aantal herhalingen in iedere groep n <- g*ni # totaal aantal observaties alpha <- 0.05 # significantieniveau van een individuele test N <- 10000 #aantal simulaties set.seed(302) #seed zodat resultaten exact geproduceerd kunnen worden trt <- factor(rep(1:g, ni)) #factor cnt <- 0 #teller voor aantal foutieve verwerpingen for(i in 1:N) { if (i %% 1000 == 0) cat(i, "/", N, "\n") y <- rnorm(n) tests <- pairwise.t.test(y, trt, "bonferroni") verwerp <- min(tests$p.value, na.rm = T) < alpha if(verwerp) cnt <- cnt+1 } cnt/N ``` - We vinden een FWER van `r round(cnt/N*100,1)`% (een beetje conservatief) - Bij simulaties voor $g=5$ groepen, vinden we een FWER van $4.1\%$ (conservatiever). - Door Bonferroni correctie is kans op minstens één vals positief resultaat $< \alpha_E$. - Power verlies aangezien werkelijke niveau lager is dan het vooropgestelde 5% experimentsgewijs significantieniveau. ## Tukey Methode - Minder conservatieve methode - Implementatie benadert de nuldistributie van de posthoc test d.m.v. simulaties. - Resultaten kunnen lichtjes verschillen wanneer je posthoc analyse opnieuw uitvoert. - Details van de methode vallen buiten het bestek van deze cursus. - Implementatie in multcomp package: - aangepaste p-waarden - aangepaste betrouwbaarheidsintervallen - Niet nodig om functies te definiëren voor Tukey gecorrigeerde BIs ### Captopril voorbeeld ```{r} model1.mcp <- glht(model1, linfct = mcp(dose = "Tukey")) summary(model1.mcp) ``` ```{r} confint(model1.mcp) ``` ```{r out.width='100%', fig.asp=.8, fig.align='center'} model1.mcp |> confint() |> plot() ``` ### Evalueer Tukey methode in een simulatiestudie ```{r} g <- 3 # aantal behandelingen (g=3) ni <- 12 # aantal herhalingen in iedere groep n <- g*ni # totaal aantal observaties alpha <- 0.05 # significantieniveau van een individuele test N <- 10000 #aantal simulaties set.seed(302) #seed zodat resultaten exact geproduceerd kunnen worden trt <- factor(rep(1:g, ni)) #factor cnt <- 0 #teller voor aantal foutieve verwerpingen for(i in 1:N) { if (i %% 1000 == 0) cat(i, "/", N, "\n") y <- rnorm(n) m <- lm(y ~ trt) m.mcp <- glht(m, linfct = mcp(trt = "Tukey")) tests <- summary(m.mcp)$test verwerp <- min( as.numeric(tests$pvalues), na.rm=T) < alpha if(verwerp) cnt <- cnt+1 } cnt/N ``` De methode wordt mooi op het 5% experimentsgewijs significantieniveau gecontroleerd. # Conclusies: Prostacycline voorbeeld Volledige analyse voor voorbeeld prostacycline 1. Anova vóór posthoc-tests: F-test heeft een hogere power dan paarsgewijze t-test - F-test gebruikt alle gegevens - Voor F-test hoeven we niet te corrigeren voor meervoudige testen: er wordt één test uitgevoerd voor de algemene omnibushypothese ```{r} model1 <- lm(prostac ~ dose, data = prostacyclin) anova(model1) ``` ```{r} model1.mcp <- glht(model1, linfct = mcp(dose = "Tukey")) summary(model1.mcp) ``` ```{r} confint(model1.mcp) ``` - Er is een extreem significant effect van arachidonzuur op de gemiddelde bloedconcentratie van prostacycline bij ratten ($p<0.001$). De gemiddelde prostacyclineconcentratie is hoger in de hoge dosisgroep dan in de lage en matige dosisgroep (beide p-waardes zijn kleiner dan $p<0.001$). - De gemiddelde concentratie in de hoge dosisgroep is `r round(confint(model1.mcp)$confint[2,1],1)`ng/ml (95% CI [`r paste(round(confint(model1.mcp)$confint[2,2:3],1),collapse=",")`]ng/ml) en `r round(confint(model1.mcp)$confint[3,1],1)`ng/ml (95% BI [`r paste(round(confint(model1.mcp)$confint[3,2:3],1),collapse=",")`]ng/ml) hoger dan in de lage en matige dosisgroep, respectievelijk. - Het verschil in gemiddelde prostacyclineconcentratie tussen de matige en lage dosisgroep is niet significant (p=`r round(summary(model1.mcp)$test$pvalues[1],2)`). (Alle p-waardes en betrouwbaarheidsintervallen voor post-hoc-tests worden gecorrigeerd voor meervoudige tests met behulp van de Tukey-methode). Merk op dat het belangrijk is om ook de niet significante resultaten te vermelden! ***