Datensatz

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

Code-Book:

  • pain.red: Schmerzreduktion, gemessen auf einer Skala von 0 bis 10. Positive Werte bedeuten eine Reduktion.

  • 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")
d.pain <- import("~/Documents/Git/stats-msc-bfh/Data/back_pain.csv")
summary(d.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 daher in einen Faktor um:

d.pain$treatment <- factor(d.pain$treatment)
summary(d.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(d.pain$pain.red[d.pain$treatment == "Trt1"]) # Mittelwert für die erste Gruppe
## [1] 1.9

Analog erfolgt die Berechnung für alle anderen Statistiken. Etwas eleganter geht das mit dem tidyverse-Paket (für Interessierte, wird für die Prüfung nicht vorausgesetzt).

library(tidyverse)
d.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, d.pain)

Wer es auf die Spitzen treiben möchte, stellt zusätzlich noch die einzelnen Werte dar:

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

Die Abbildung ist kongruent zur Interpretation der oben berechneten Statistiken.

Modell

Aufgabe

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

  1. Wie lautet \(H_0\)?

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

  3. Erkläre, warum eine ANOVA/multiple lineare Regression besser ist, als einzelne t-Tests durchzuführen.

  4. Passe ein Modell mit der Funktion aov() 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 Mittelwerte gleich sind, also \(\mu_{Trt1} = \mu_{Trt2} = \mu_{Trt3}\). Die Alternativhypothese besagt folglich, dass sich mindestens ein Mittelwert unterscheidet. Weil wir drei Gruppen haben, können insgesamt 3 Vergleiche gemacht werden (1 vs 2, 1 vs 3 und 2 vs3). Generell gilt, dass es bei \(g\) Gruppen \(g * (g-1) / 2\) paarweise Vergleiche gibt.

Da wir bei jedem Test ein gewisses Risiko für einen \(\alpha\)-Fehler in Kauf nehmen (oft 5%), kumuliert sich das Risiko für einen \(\alpha\)-Fehler mit zunehmender Anzahl Tests. Wir erwarten also bei 20 Tests 1 signifikantes Ergebnis, wenn in Wirklichkeit immer \(H_0\) zutrifft. Aus diesem Grund ist es besser, alle Informationen in einem Modell zu berechnen. Dies wird für dieses Beispiel wie folgt gemacht:

aov.pain <- aov(pain.red ~ treatment, d.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 angezeigte p-Wert? (alle anderen Informationen zum F-Test werden in einer anderen Übung besprochen)

  2. Lass dir mit der Funktion summary.lm() 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(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

Wir sehen im Output oben, dass der p-Wert sehr klein ist, was bedeutet, dass \(H_0\) nicht plausibel ist. Inhaltlich bedeutet das, dass sich mindestens ein Mittelwert von den anderen unterscheidet. Weil der Informationsgehalt des Outputs sonst sehr mager ist, schauen wir uns die Koeffizienten an:

summary.lm(aov.pain)
## 
## Call:
## aov(formula = pain.red ~ treatment, data = d.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 Output erinnert stark an jenen einer linearen Regression. Tatsächlich sind eine ANOVA und eine lineare Regression mathematisch äquivalent!

lm.pain <- lm(pain.red ~ treatment, d.pain)
summary(lm.pain)
## 
## Call:
## lm(formula = pain.red ~ treatment, data = d.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

  • \(\hat{\beta_0}\) entspricht dem Mittelwert der Gruppe Trt1. In diesem Fall hat das Intercept also sehr wohl eine inhaltliche Bedeutung.

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

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

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)*0 + (-1.84)*0 = 1.9\)

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

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

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

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

Wir entnehmen dem Output, dass sich Trt3, jedoch nicht Trt2 statistisch signifikant von Trt1 unterscheidet. 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, in dem wir die Levels von Treatment neu definieren. Wie wir oben gesehen haben, spielt es keine Rolle, ob wir die aov() oder die lm() Funktion brauchen.

d.pain$treatment <- relevel(d.pain$treatment, ref = "Trt3")
lm.pain.2 <- lm(pain.red ~ treatment, d.pain)
summary(lm.pain.2)
## 
## Call:
## lm(formula = pain.red ~ treatment, data = d.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

Möchten wir alle Gruppen mit einander vergleichen, müssen wir also mindestens einmal das Referenzlevel von treatment ändern. Somit haben wir wiederum zwei Modelle, was wir eigentlich nicht wollen (siehe oben).

Aufgabe

Mit dem Tukey-Test kann man alle Gruppen miteinander vergleichen, wobei P-Werte und 95% CIs für multiples Testen korrigiert werden. 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 den Test durch lassen eine Zusammenfassung davon anzeigen.

tukey <- TukeyHSD(aov.pain)
tukey
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = pain.red ~ treatment, data = d.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

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

plot(tukey)

Trt3 unterscheidet sich stat. sign. von Trt1 und Trt2 (die 95% CI’s schneiden 0 nicht) Trt2 und Trt1 unterscheiden sich hingegen nicht stat. sign. (das 95% CI schneidet 0).

Modelldiagnostik

Aufgabe

  1. Welche zwei Voraussetzungen sollten ü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

Die Frage nach einem linearen Zusammenhang ist bei diesem Beispiel weit hergeholt (ein Scatterplot macht in diesem Fall wenig Sinn). Wir fokussieren uns daher auf die Residuenanalyse:

  1. Liegt Heteroskedastizität vor?

Um dies zu überprüfen, erstellen wir einen Residuen-Plot (auch Tukey Anscombe Plot genannt):

plot(aov.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. Sind die Residuen ungefähr normalverteilt?

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

plot(aov.pain, which = 2)

Auch hier besteht keine Evidenz, dass die Verteilung der Residuen von einer Normalverteilung abweicht.

Anmerkung: Dass die Annahmen praktisch perfekt erfüllt sind liegt daran, dass ich die Daten aus einer Normalverteilung mit gleicher Varianz simuliert habe.