Datensatz

Für diese Übung brauchst du den Datensatz “back-pain.csv”. Der Datensatz beinhaltet simulierte Daten einer randomisierten Studie, in der der Effekt von Kinesio-Tape auf die kurzfristige Schmerzreduktion untersucht wurde.

Codebuch:

  • pain.red: Schmerzreduktion in Skalenpunkten (positive Werte) bzw. Schmerzzunahme (negative Werte). Die Variable ist eine Veränderungsgrösse und daher nicht auf Werte zwischen 0 und 10 beschränkt.

  • treatment: Trt1 = Tape, Trt2 = Placebo-Tape, Trt3 = Kontrollgruppe (kein Tape).

Aufgaben

  1. Lade und inspiziere den Datensatz.
  2. Wandle die Variable treatment in einen Faktor um.
  3. Berechne den Mittelwert und die Standardabweichung pro Gruppe und interpretiere.
  4. Stelle die Daten pro Gruppe geeignet dar.

Lösungen

Wir laden den Datensatz wie gewohnt mit import() aus dem Paket rio und schauen diesen kurz an (hier mit summary()).

library("rio")
df.pain <- import("../data/back-pain.csv")
summary(df.pain)
##   treatment            pain.red     
##  Length:90          Min.   :-4.600  
##  Class :character   1st Qu.:-0.075  
##  Mode  :character   Median : 1.100  
##                     Mean   : 1.270  
##                     3rd Qu.: 2.550  
##                     Max.   : 5.800

Wir sehen, dass die Variable treatment momentan noch als character hinterlegt ist und wandeln sie daher in einen Faktor um:

df.pain$treatment <- factor(df.pain$treatment)
summary(df.pain)
##  treatment    pain.red     
##  Trt1:30   Min.   :-4.600  
##  Trt2:30   1st Qu.:-0.075  
##  Trt3:30   Median : 1.100  
##            Mean   : 1.270  
##            3rd Qu.: 2.550  
##            Max.   : 5.800

Wir sehen nun, dass es in jeder Gruppe 30 Beobachtungen gibt. Als Nächstes berechnen wir die Mittelwerte und Standardabweichungen für alle Gruppen. Es gibt verschiedene Möglichkeiten, das zu tun. Die Indexierung kannst du wie gewohnt mit den eckigen Klammern vornehmen, z.B.

mean(df.pain$pain.red[df.pain$treatment == "Trt1"]) # Mittelwert für die erste Gruppe
## [1] 1.9

Mit der Funktion tapply() kannst du die Werte gleichzeitig für alle Gruppen berechnen:

tapply(df.pain$pain.red, df.pain$treatment, mean)
##       Trt1       Trt2       Trt3 
## 1.90000000 1.85333333 0.05666667
tapply(df.pain$pain.red, df.pain$treatment, sd)
##     Trt1     Trt2     Trt3 
## 1.959416 1.663136 1.736362

Etwas eleganter geht das mit dem tidyverse-Paket (für Interessierte).

library(tidyverse)
df.pain %>% 
  group_by(treatment) %>% 
  summarise(M = mean(pain.red), S = sd(pain.red))
## # A tibble: 3 × 3
##   treatment      M     S
##   <fct>      <dbl> <dbl>
## 1 Trt1      1.9     1.96
## 2 Trt2      1.85    1.66
## 3 Trt3      0.0567  1.74

Die Mittelwerte der ersten beiden Gruppen sind sehr ähnlich, derjenige der dritten Gruppe ist deutlich tiefer. Die Standardabweichungen sind relativ ähnlich.

Zum Schluss stellen wir die Daten graphisch dar. Weil es sich bei der Variable treatment um eine kategoriale Variable handelt, eignen sich Boxplots:

boxplot(pain.red ~ treatment, df.pain)

Man kann zusätzlich noch die einzelnen Werte darstellen (für Interessierte):

library(ggplot2)
ggplot(df.pain, aes(x = treatment, y = pain.red)) +
  geom_boxplot() +
  geom_jitter(width = 0.1)

Die Abbildung stimmt mit der Interpretation der oben berechneten deskriptiven Statistiken überein.

Modell

Aufgabe

Untersuche, ob sich die Mittelwerte der drei Gruppen statistisch signifikant unterscheiden (\(\alpha = 0.05\)).

  1. Wie lautet \(H_0\)?

  2. Wie viele paarweise Vergleiche können gemacht werden?

  3. Erkläre, warum ein globaler F-Test in einer einfaktoriellen ANOVA sinnvoller ist als mehrere unadjustierte t-Tests.

  4. Passe ein Modell mit der Funktion lm() an. Speichere das Resultat in ein neues R-Objekt.

  5. Lass dir mit summary() eine Zusammenfassung des Modells anzeigen.

Lösung

\(H_0\) besagt, dass alle Erwartungswerte, also alle Mittelwerte in der Population, gleich sind, also

\[ \mu_{Trt1} = \mu_{Trt2} = \mu_{Trt3} \]

bzw.

\[ \beta_1 = \beta_2 = 0. \]

Die Alternativhypothese besagt folglich, dass sich mindestens zwei Gruppenmittelwerte unterscheiden bzw. mindestens einer der beiden Koeffizienten der Gruppenindikatoren nicht 0 ist. Weil wir drei Gruppen haben, können insgesamt 3 Vergleiche gemacht werden (1 vs 2, 1 vs 3 und 2 vs 3). Generell gilt, dass es bei \(k\) Gruppen \(k \times (k-1) / 2\) paarweise Vergleiche gibt.

Bei \(k = 3\) Gruppen gibt es \(k(k-1)/2 = 3\) paarweise Vergleiche. Werden alle drei Tests unadjustiert auf dem Niveau \(\alpha = 0.05\) durchgeführt, ist die Wahrscheinlichkeit, mindestens eine wahre Nullhypothese fälschlich abzulehnen (familienweise Fehlerwahrscheinlichkeit), im Allgemeinen grösser als 5%. Wären die drei Tests unabhängig und alle Nullhypothesen wahr, würde diese Wahrscheinlichkeit \(1 - (1 - \alpha)^3 = 1 - 0.95^3 \approx 0.143\), also 14.3%, betragen. Die paarweisen Tests sind hier jedoch abhängig, da sie jeweils Daten derselben Gruppen verwenden; 14.3% ist deshalb nicht die tatsächliche Fehlerwahrscheinlichkeit dieser drei Vergleiche.

Der globale F-Test prüft stattdessen die Gleichheit aller drei Gruppenmittelwerte mit einem einzigen Test auf dem Niveau 5%. Für anschliessende paarweise Vergleiche ist eine Adjustierung erforderlich, wenn die familienweise Fehlerwahrscheinlichkeit auf 5% kontrolliert werden soll, beispielsweise mit dem weiter unten verwendeten Tukey-Verfahren. Das zugehörige ANOVA-Modell passen wir wie folgt an:

lm.pain <- lm(pain.red ~ treatment, df.pain)
summary(lm.pain)
## 
## Call:
## lm(formula = pain.red ~ treatment, data = df.pain)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.657 -1.057 -0.105  1.087  4.043 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.90000    0.32694   5.811 9.98e-08 ***
## treatmentTrt2 -0.04667    0.46237  -0.101 0.919838    
## treatmentTrt3 -1.84333    0.46237  -3.987 0.000139 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.791 on 87 degrees of freedom
## Multiple R-squared:  0.192,  Adjusted R-squared:  0.1734 
## F-statistic: 10.33 on 2 and 87 DF,  p-value: 9.4e-05

Auch die Variante mit aov() ist korrekt (intern sind beide Varianten äquivalent, nur der Output ist verschieden):

aov.pain <- aov(lm(pain.red ~ treatment, df.pain))
summary(aov.pain)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## treatment    2  66.28   33.14   10.34 9.4e-05 ***
## Residuals   87 278.99    3.21                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretation des Outputs

Nun geht es an die Interpretation des Outputs von oben.

Aufgaben

  1. Was bedeutet der p-Wert des globalen F-Tests? (alle anderen Informationen zum F-Test werden in einer anderen Übung besprochen)

  2. Lass dir die geschätzten Koeffizienten anzeigen.

  3. Interpretiere \(\hat{\beta_0}\), \(\hat{\beta_1}\) und \(\hat{\beta_2}\)

  4. Verwende die Regressionsformel \(\hat{y}=\hat{\beta_0} + \hat{\beta_1}Trt2 + \hat{\beta_2}Trt3\) um den Mittelwert für die drei Gruppen zu schätzen.

  5. Was sind die 95%-Konfidenzintervalle der Koeffizienten? Tipp: confint()

  6. Kannst du das Modell so erstellen, dass die Kontrollgruppe die Referenzgruppe ist (Tipp: relevel()).

Lösungen

summary(lm.pain)
## 
## Call:
## lm(formula = pain.red ~ treatment, data = df.pain)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.657 -1.057 -0.105  1.087  4.043 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1.90000    0.32694   5.811 9.98e-08 ***
## treatmentTrt2 -0.04667    0.46237  -0.101 0.919838    
## treatmentTrt3 -1.84333    0.46237  -3.987 0.000139 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.791 on 87 degrees of freedom
## Multiple R-squared:  0.192,  Adjusted R-squared:  0.1734 
## F-statistic: 10.33 on 2 and 87 DF,  p-value: 9.4e-05

Der p-Wert des globalen F-Tests beträgt ungefähr 0.000094. Unter \(H_0\) und den Modellannahmen ist dies die Wahrscheinlichkeit, einen mindestens so grossen F-Wert wie den beobachteten zu erhalten. Wir verwerfen \(H_0\) auf dem Niveau 5%: Die Daten liefern Evidenz dafür, dass sich mindestens zwei Gruppenmittelwerte unterscheiden. Der Globaltest sagt nicht, welche Gruppen sich unterscheiden.

  • \(\hat{\beta_0}\) entspricht dem Mittelwert der Gruppe Trt1 (Referenzkodierung mit Trt1 als Referenzgruppe). Der Achsenabschnitt (Intercept) hat hier also eine inhaltliche Bedeutung.

  • \(\hat{\beta_1}\) entspricht der Mittelwertsdifferenz Trt2 − Trt1.

  • \(\hat{\beta_2}\) entspricht der Mittelwertsdifferenz Trt3 - Trt1.

Dabei sind \(Trt2\) und \(Trt3\) Indikatorvariablen: Sie sind für die jeweilige Gruppe 1 und sonst 0. Wir setzen nun in die Gleichung \(\hat{y}=\hat{\beta_0} + \hat{\beta_1}Trt2 + \hat{\beta_2}Trt3\) ein, um die Mittelwerte für die Gruppen zu schätzen.

  • \(\hat{y_{Trt1}} = 1.9 + (-0.047) \times 0 + (-1.84) \times 0 = 1.9\)

  • \(\hat{y_{Trt2}} = 1.9 + (-0.047) \times 1 + (-1.84) \times 0 = 1.853\)

  • \(\hat{y_{Trt3}} = 1.9 + (-0.047) \times 0 + (-1.84) \times 1 = 0.06\)

Bis auf Rundungsdifferenzen sind das genau die Werte, welche oben zur Deskription berechnet wurden:

df.pain %>% 
  group_by(treatment) %>% 
  summarise(M = mean(pain.red), S = sd(pain.red))
## # A tibble: 3 × 3
##   treatment      M     S
##   <fct>      <dbl> <dbl>
## 1 Trt1      1.9     1.96
## 2 Trt2      1.85    1.66
## 3 Trt3      0.0567  1.74

Die unadjustierten t-Tests im Output zeigen einen statistisch signifikanten Unterschied zwischen Trt3 und Trt1, jedoch nicht zwischen Trt2 und Trt1. Ein nicht signifikanter Unterschied belegt keine Gleichheit der Gruppenmittelwerte. Die 95%-Konfidenzintervalle erhalten wir so:

confint(aov.pain)
##                    2.5 %     97.5 %
## (Intercept)    1.2501643  2.5498357
## treatmentTrt2 -0.9656731  0.8723398
## treatmentTrt3 -2.7623398 -0.9243269

In diesem Fall erscheint es sinnvoll, die Gruppe Trt3 (=Kontrollgruppe) als Referenzlevel zu wählen. Dies erreichen wir, indem wir die Referenzgruppe von treatment ändern.

df.pain$treatment <- relevel(df.pain$treatment, ref = "Trt3")
lm.pain.2 <- lm(pain.red ~ treatment, df.pain)
summary(lm.pain.2)
## 
## Call:
## lm(formula = pain.red ~ treatment, data = df.pain)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.657 -1.057 -0.105  1.087  4.043 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.05667    0.32694   0.173 0.862801    
## treatmentTrt1  1.84333    0.46237   3.987 0.000139 ***
## treatmentTrt2  1.79667    0.46237   3.886 0.000199 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.791 on 87 degrees of freedom
## Multiple R-squared:  0.192,  Adjusted R-squared:  0.1734 
## F-statistic: 10.33 on 2 and 87 DF,  p-value: 9.4e-05

Wir sehen, dass sich Trt1 und Trt2 signifikant von Trt3 unterscheiden.

Paarweise Vergleiche

Durch einen Wechsel der Referenzgruppe können wir weitere paarweise Vergleiche ausgeben lassen. Dabei wird dasselbe Modell anders parametrisiert; die vorhergesagten Werte bleiben gleich. Die unadjustierten Tests kontrollieren jedoch nicht die familienweise Fehlerwahrscheinlichkeit für alle drei Vergleiche.

Aufgabe

Mit dem Tukey-Verfahren kann man alle Gruppen paarweise vergleichen. Es liefert adjustierte p-Werte und simultane 95%-Konfidenzintervalle und kontrolliert die familienweise Fehlerwahrscheinlichkeit unter den ANOVA-Annahmen. Führe diesen Test durch und stelle das Resultat grafisch dar. Speichere dazu das Resultat des Tests wiederum in einem neuen R-Objekt ab. Beachte, dass dieser Test ein aov-Objekt braucht. Tipp: TukeyHSD(), plot()

Lösung

Zuerst führen wir das Verfahren durch und lassen eine Zusammenfassung davon anzeigen.

tukey <- TukeyHSD(aov.pain)
tukey
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = lm(pain.red ~ treatment, df.pain))
## 
## $treatment
##                  diff       lwr        upr     p adj
## Trt2-Trt1 -0.04666667 -1.149174  1.0558404 0.9943999
## Trt3-Trt1 -1.84333333 -2.945840 -0.7408263 0.0004054
## Trt3-Trt2 -1.79666667 -2.899174 -0.6941596 0.0005772

Diese Ergebnisse können mit der plot()-Funktion grafisch dargestellt werden:

plot(tukey)

Trt3 unterscheidet sich statistisch signifikant von Trt1 und Trt2 (die simultanen 95%-Konfidenzintervalle enthalten 0 nicht). Für Trt2 und Trt1 ist der Unterschied nicht statistisch signifikant (das Intervall enthält 0). Dies belegt keine Gleichheit der Gruppenmittelwerte.

Modelldiagnostik

Aufgabe

  1. Welche zwei Verteilungsannahmen sollten anhand der Residuen überprüft werden, damit das Modell (insbesondere die inferenzstatistischen Angaben) korrekt interpretiert werden kann?

  2. Überprüfe die Voraussetzungen durch visuelle Hilfsmittel und interpretiere diese. Tipp: plot(…, which = …)

Lösung

Mit einem kategorialen Prädiktor werden separate Gruppenmittelwerte modelliert; eine lineare Beziehung zu numerischen Gruppenlabels ist nicht erforderlich. Zusätzlich setzt die klassische ANOVA unabhängige Fehler voraus. Diese Annahme muss anhand des Studiendesigns beurteilt werden und lässt sich nicht mit den folgenden Residuenplots prüfen. Wir untersuchen:

  1. Liegt Heteroskedastizität vor?

Um dies zu überprüfen, erstellen wir einen Residuen-Plot (TA-Plot):

plot(lm.pain, which = 1)

Die Varianzen der Gruppen sind sehr ähnlich. Das überrascht uns nicht, da wir bei der deskriptiven Statistik schon gesehen haben, dass die Standardabweichungen ähnlich sind.

  1. Ist die Annahme normalverteilter Fehler plausibel?

Um dies zu überprüfen, eignet sich ein QQ-Plot:

plot(lm.pain, which = 2)

Der QQ-Plot zeigt keine auffälligen Abweichungen von der Normalverteilung. Dies stützt die Annahme normalverteilter Fehler, beweist sie jedoch nicht.

Anmerkung: Die Daten wurden aus Normalverteilungen mit gleicher Varianz simuliert. Die Verteilungsannahmen gelten somit für den datengenerierenden Prozess; die Stichprobe kann dennoch zufällige Abweichungen zeigen.