7. Analysis of Variance

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/GTPB/PSLS20/master/data/prostacyclin.txt")

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

head(prostacyclin)

# A tibble: 6 x 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 \vert \text{groep j} \sim N(\mu_j,\sigma^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:\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.

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*}\]

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: \[\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.

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 \[ \text{SSTot} = \text{SST} + \text{SSE} \]

##SSTot

##SST

3.2 SSE

3.3 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)\)
  
• \(\text{MSE}=\text{SSE}/(n-p)\)
• Teststatistiek vergelijkt de variabiliteit verklaard door model (MST) met de residuele variabiliteit (MSE)

of

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

3.4 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.4.1 F-verdeling met kritieke waarde (\(\alpha\)=5%) en geobserveerde F-statistiek voor prostacycline voorbeeld

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

3.4.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

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\).

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

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.

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

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

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!

---

Home