Datensatz

In einer randomisierten, kontrollierten Studie wurde untersucht, ob das «Beat the Blues» Programm zur Behandlung von Menschen mit Depression effektiver ist als die Standardbehandlung. Bei den vorliegenden Daten handelt es sich um einen Teildatensatz. Mehr Informationen zur Methodik dieser Studie finden Sie hier.

Der Datensatz umfasst folgende Variablen:

  • drug: Ob der Proband Antidepressiva nimmt oder nicht
  • length: Länge der aktuellen Depressionsepisode (<6 Monate oder >6 Monate)
  • treatment: Standardtherapie (TAU) oder «Beat the Blues» (BtheB)
  • bdi_pre: Beck Depression Inventory vor der Therapie
  • bdi_2m: Beck Depression Inventory nach zwei Monaten
  • bdi_4m: Beck Depression Inventory nach vier Monaten
  • bdi_6m: Beck Depression Inventory nach sechs Monaten

Eine höherer BDI steht für einen höheren Grad einer Depression.

Aufgaben

  1. Laden Sie den Datensatz in RStudio. In welchem Format liegt der Datensatz aktuell vor: im long- oder im wide-Format?
  2. Stellen Sie die BDI-Werte nach 6 Monaten pro Gruppe zum Vergleich mittels Boxplots dar.
  3. Stellen Sie die BDI-Werte nach Zeitpunkt mittels Boxplots dar.

Lösungen

  1. Laden Sie den Datensatz in RStudio. In welchem Format liegt der Datensatz aktuell vor: im long- oder im wide-Format?
library(rio)
data <- import("../../Data/BtheB.csv")
str(data)
## 'data.frame':    58 obs. of  8 variables:
##  $ Subject  : int  2 4 6 7 8 9 10 11 14 15 ...
##  $ drug     : chr  "Yes" "No" "Yes" "Yes" ...
##  $ length   : chr  ">6m" ">6m" "<6m" "<6m" ...
##  $ treatment: chr  "BtheB" "BtheB" "BtheB" "TAU" ...
##  $ bdi_pre  : int  32 21 7 17 20 18 20 30 30 23 ...
##  $ bdi_2m   : int  16 17 0 7 20 13 5 32 26 13 ...
##  $ bdi_4m   : int  24 16 0 7 21 14 5 24 36 13 ...
##  $ bdi_6m   : int  17 10 0 3 19 20 8 12 27 12 ...

Der Datensatz liegt im wide-Format vor.

  1. Stellen Sie die BDI-Werte nach 6 Monaten pro Gruppe zum Vergleich mittels Boxplots dar.

Die Standardvariante:

boxplot(data$bdi_6m ~ data$treatment)

Für die Motivierten mit ggplot

library(tidyverse)
data %>% 
   ggplot(aes(y = bdi_6m, fill = treatment)) + geom_boxplot() + 
   theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

  1. Stellen Sie die BDI-Werte nach Zeitpunkt mittels Boxplots dar.

Die Standardvariante:

boxplot(data$bdi_pre, data$bdi_2m, data$bdi_4m, data$bdi_6m)


Für die Motivierten mit ggplot:

library(tidyverse)
data %>% 
   pivot_longer(bdi_pre:bdi_6m, values_to = "bdi_scores", names_to = "time") %>% # macht ein long-format
   mutate(time = fct_relevel(time, "bdi_pre","bdi_2m", "bdi_4m", "bdi_6m")) %>%
   ggplot(aes(y = bdi_scores, x = time, fill = treatment)) + geom_boxplot()

Berücksichtigen von einem Zeitpunkt

Da es sich um ein RCT handelt, wäre es möglich die Gruppen nach 6 Monaten zu vergleichen.

Aufgaben

  1. Mit welchem Test können Sie überprüfen, ob sich die Gruppen nach 6 Monaten bzgl. dem BDI-Score unterscheiden?
  2. Führen Sie diesen Test (oben) durch. Was sagt Ihnen das Ergebnis?

Lösungen

  1. Mit welchem Test können Sie überprüfen, ob sich die Gruppen nach 6 Monaten bzgl. dem BDI-Score unterscheiden?

z.B. mit einen ungepaarten t-Test bzw. mit einer einfachen linearen Regression

  1. Führen Sie diesen Test (oben) durch. Was sagt Ihnen das Ergebnis?

t-Test:

t.test(data$bdi_6m ~ data$treatment)
## 
##  Welch Two Sample t-test
## 
## data:  data$bdi_6m by data$treatment
## t = -2.5109, df = 46.969, p-value = 0.01554
## alternative hypothesis: true difference in means between group BtheB and group TAU is not equal to 0
## 95 percent confidence interval:
##  -12.670559  -1.398406
## sample estimates:
## mean in group BtheB   mean in group TAU 
##            9.241379           16.275862


Lineares Modell (abgesehen von der Welch-Korrektur identisch):

summary(lm(bdi_6m ~ treatment, data = data))
## 
## Call:
## lm(formula = bdi_6m ~ treatment, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.2759  -8.0000   0.2414   6.4828  30.7241 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     9.241      1.981   4.665 1.96e-05 ***
## treatmentTAU    7.034      2.802   2.511    0.015 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.67 on 56 degrees of freedom
## Multiple R-squared:  0.1012, Adjusted R-squared:  0.08514 
## F-statistic: 6.305 on 1 and 56 DF,  p-value: 0.01495

Die Mittelwertsdifferenz beträgt -7 (BtheB - TAU) und das 95% CI geht von -12.2 bis -1.4. Wir vertrauen also zu 95% darauf, dass die wahre Mittelwertsdifferenz in diesem Intervall liegt. Weil dieses die Nullhypothese nicht einschliesst, ist diese Mittelwertsdifferenz nach 6 Monaten statistisch signifikant (p = 0.01554).


Berücksichtigen von zwei Zeitpunkten

Der Ansatz oben ist nur bei grossen randomisierten Studien vertretbar, wenn die Gruppen zu Beginn perfekt vergleichbar sind. Diese Abbildung zeigt, dass dies in dieser Studie nicht ganz der Fall ist:

boxplot(data$bdi_pre ~ data$treatment, range = 0)
stripchart(data$bdi_pre ~ data$treatment, vertical = TRUE,
           method = "jitter", add = TRUE, pch = 20, cex = 0.3)

Aufgaben

  1. Berechnen Sie eine neue Variable bdi_diff, welche die Differenz zwischen bdi_pre und bdi_6m darstellt. Stellen Sie die Differenzen der beiden Gruppen auch wieder mittels Boxplots dar.
  2. Führen Sie einen t-Test durch, um zu untersuchen, ob sich die Mittelwerte der Differenz zwischen den beiden Gruppen unterscheiden. Interpretieren Sie das Resultat.
  3. Rechnen Sie ein lineares Modell, in welchem die BDI-Differenz die abhängige Variable ist und die Gruppe die unabhängige Variable. Korrigieren Sie zusätzlich für die Baseline. Ändert sich die Aussage im Vergleich zum Test ohne Baseline-Adjustment?

Lösungen

  1. Berechnen Sie eine neue Variable bdi_diff, welche die Differenz zwischen bdi_pre und bdi_6m darstellt. Stellen Sie die Differenzen der beiden Gruppen auch wieder mittels Boxplots dar.
data$bdi_diff <- data$bdi_6m - data$bdi_pre
boxplot(data$bdi_diff ~ data$treatment)

  1. Führen Sie einen t-Test durch, um zu untersuchen, ob sich die Mittelwerte der Differenz zwischen den beiden Gruppen unterscheiden. Interpretieren Sie das Resultat.
t.test(data$bdi_diff ~ data$treatment)
## 
##  Welch Two Sample t-test
## 
## data:  data$bdi_diff by data$treatment
## t = -1.8522, df = 53.063, p-value = 0.06956
## alternative hypothesis: true difference in means between group BtheB and group TAU is not equal to 0
## 95 percent confidence interval:
##  -10.5579177   0.4199867
## sample estimates:
## mean in group BtheB   mean in group TAU 
##          -12.241379           -7.172414

Die Differenz ist in der BtheB Gruppe etwa 5 grösser. Jedoch schneidet das 95% CI die 0 und der p-Wert beträgt 7%. Nach der klassischen 5%-Regel würde man hier nicht von einem statistisch signifikanten Unterschied sprechen.

  1. Rechnen Sie ein lineares Modell, in welchem die BDI-Differenz die abhängige Variable ist und die Gruppe die unabhängige Variable. Korrigieren Sie zusätzlich für die Baseline. Ändert sich die Aussage im Vergleich zum Test ohne Baseline-Adjustment?
my_lm2 <- lm(bdi_diff ~ treatment + bdi_pre, data)
summary(my_lm2)
## 
## Call:
## lm(formula = bdi_diff ~ treatment + bdi_pre, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.3365  -6.3056  -0.3402   5.3605  24.1388 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -2.0297     3.1156  -0.651 0.517465    
## treatmentTAU   6.0033     2.4508   2.450 0.017517 *  
## bdi_pre       -0.4753     0.1208  -3.936 0.000235 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.288 on 55 degrees of freedom
## Multiple R-squared:  0.2648, Adjusted R-squared:  0.2381 
## F-statistic: 9.904 on 2 and 55 DF,  p-value: 0.000212


Und die Vertrauensintervalle der Koeffizienten:

confint(my_lm2)
##                   2.5 %     97.5 %
## (Intercept)  -8.2735136  4.2141427
## treatmentTAU  1.0918008 10.9147227
## bdi_pre      -0.7173863 -0.2333011

Interpretation

Ja: Die Differenz ist mit einer Baseline-Adjustierung 6 Punkte grösser in der BtheB-Gruppe. Das 95% CI dieser Differenz schneidet 0 nicht, der p-Wert beträgt 1.8%. Zudem sieht man, dass bdi_pre ebenfalls ein wichtiger Prädiktor für die Differenz ist.

Exkurs: mehr als zwei Messzeitpunkte

In unseren Foschungsgebieten kommt es häufig vor, dass mehr als zwei Messzeitpunkte vorliegen. Dies ist beispielsweise der Fall, wenn bei der gleichen Person Ganganalysen mit drei verschiedenen Schuheinlagen gemacht werden.

Auch beim vorliegenden Datensatz “Beat the Blues” gibt es mehr als zwei Messzeitpunkte pro Person. Damit man solche Daten einfacher visualisieren und analysieren kann, muss man sie vorher in ein “langes Format” bringen:

data_long <- data %>% 
   pivot_longer(bdi_pre:bdi_6m, values_to = "bdi", names_to = "time") %>% 
   mutate(time = fct_relevel(time, "bdi_pre","bdi_2m", "bdi_4m", "bdi_6m")) %>% 
  mutate(treatment = factor(treatment, levels = c("TAU", "BtheB")))

ggplot(data_long, aes(x = time, y = bdi, fill = treatment)) + 
  geom_boxplot()

Mixed effects models

Mixed effects models sind eine pupuläre Methode, um solche (korrelierten) Daten zu modellieren. Diese Übung hat nie und nimmer den Anspruch, diese Modelle zu erklären. Es geht einzig darum festzuhalten, dass es in Situationen mit mehr als zwei Messzeitpunkten meistens auf ein “Mixed effects model” hinausläuft.

In der Grafik unten sehen wir, dass beide Gruppen mit der Zeit besser werden. Somit könnte der Haupteffekt Zeit relevant sein. Die BtheB-Gruppe ist grundsätzlich tiefer. Somit könnte auch der Haupteffekt Gruppe relevant sein. Man sieht auch, dass die Verbesserung bei der BtheB-Gruppe über die Zeit stärker ist (die Linien sind nicht parallel). Das spricht für einen relevanten Interaktionseffekt. Also, dass die Verbesserung über die Zeit von der Gruppe abhängig ist.

ggplot(data_long, aes(x = time, y = bdi, group = treatment, color = treatment)) + 
  stat_summary(fun = mean, geom = "line")

Hier wird gezeigt, wie man die zwei Haupteffekte und den Interaktionseffekt mit einem Mixed effets model analysieren kann.

library(lmerTest)
lme.btb <- lmer(bdi ~ time + treatment + time:treatment + (1 | Subject), data = data_long)

Grundätzlich sieht das Modell sehr ähnlich aus wie ein normales lineares Modell (lm()). Mit dem Teil (1 | Subject) können wir der Tatsache Rechnung tragen, dass die Daten von der gleichen Person korreliert sind. Eine Kurzzusammenfassung erhalten wir mit anova():

anova(lme.btb)
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## time           3354.8 1118.26     3   168 28.6500 5.168e-15 ***
## treatment       317.9  317.95     1    56  8.1459   0.00604 ** 
## time:treatment  430.2  143.38     3   168  3.6735   0.01343 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wir entnehmen dem Output, dass die zwei Haupteffekte so wie der Interaktionseffekt statistisch signifikant unterschiedlich von 0 sind. Es liegt demnach Evidnez vor, dass sich die Personen über die Zeit verbessern (Effekt time), dass es einen generellen Unterschied zwischen den Gruppen gibt (Effket treatment) und, dass die Verbesserung über die Zeit bei den beiden Gruppen unterschiedlich ist (Interaktionsffekt treatment:time). Für eine detailliertere Analyse (z.B. zu welchem Zeitpunkt sind die Gruppen unterschiedlich) könnte mit spezifischen Kontrasten gearbeitet werden.

Das Schöne an den Mixed effects models ist, dass wir im Sinne der linearen Regression “weiterdenken” und diese beliebig anpassen können (z.B. für Confounder).

Die konkreten Regressionskoeffizienten erhalten wir wie gewohnt mit summary():

summary(lme.btb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: bdi ~ time + treatment + time:treatment + (1 | Subject)
##    Data: data_long
## 
## REML criterion at convergence: 1598.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8555 -0.4794 -0.0944  0.4049  3.1077 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 66.05    8.127   
##  Residual             39.03    6.248   
## Number of obs: 232, groups:  Subject, 58
## 
## Fixed effects:
##                           Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                 23.448      1.904 102.506  12.318  < 2e-16 ***
## timebdi_2m                  -4.034      1.641 168.000  -2.459  0.01495 *  
## timebdi_4m                  -5.483      1.641 168.000  -3.342  0.00103 ** 
## timebdi_6m                  -7.172      1.641 168.000  -4.372 2.16e-05 ***
## treatmentBtheB              -1.966      2.692 102.506  -0.730  0.46698    
## timebdi_2m:treatmentBtheB   -7.034      2.320 168.000  -3.032  0.00282 ** 
## timebdi_4m:treatmentBtheB   -6.138      2.320 168.000  -2.645  0.00894 ** 
## timebdi_6m:treatmentBtheB   -5.069      2.320 168.000  -2.185  0.03030 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmbd_2 tmbd_4 tmbd_6 trtmBB t_2:BB t_4:BB
## timebdi_2m  -0.431                                          
## timebdi_4m  -0.431  0.500                                   
## timebdi_6m  -0.431  0.500  0.500                            
## tretmntBthB -0.707  0.305  0.305  0.305                     
## tmbd_2m:tBB  0.305 -0.707 -0.354 -0.354 -0.431              
## tmbd_4m:tBB  0.305 -0.354 -0.707 -0.354 -0.431  0.500       
## tmbd_6m:tBB  0.305 -0.354 -0.354 -0.707 -0.431  0.500  0.500

Die Regressionskoeffizienten sehen wir unter Fixed effects. Sie sind nichts anderes, als der mittlere BDI der TAU-Gruppe zur Baseline (= Intercept = 23.44) und die entsprechenden Mittelwertsdifferenzen. Diese können wir in die gewohnte Regressionsgleichung einsetzen.


Beispiel:

Auf den BDI-Wert der BtheB Gruppe nach 4 Monaten kommen wir, indem wir das Intercept mit dem Koeffizienten für den Zeitpunkt “timebdi_4m”, dem Koeffizienten für “die BtheB Gruppe”treatmentBtheB” und dem Koeffizienten für “timebdi_4m:treatmentBtheB” adieren:

\[23.448276 + (-1.965517) + (-5.482759) + (-6.137931) = 9.862069\]

Dies könnten wir auch ohne kompliziertes Multilevelmodell, also rein deskriptiv ausrechnen:

data_long %>% 
  group_by(time, treatment) %>% 
  summarise(bdi = mean(bdi))
## # A tibble: 8 × 3
## # Groups:   time [4]
##   time    treatment   bdi
##   <fct>   <fct>     <dbl>
## 1 bdi_pre TAU       23.4 
## 2 bdi_pre BtheB     21.5 
## 3 bdi_2m  TAU       19.4 
## 4 bdi_2m  BtheB     10.4 
## 5 bdi_4m  TAU       18.0 
## 6 bdi_4m  BtheB      9.86
## 7 bdi_6m  TAU       16.3 
## 8 bdi_6m  BtheB      9.24

Der einzige Regressionskoeffizient, welcher nicht signifikant ist, ist “treatmentBtheB”. Das ist auch gut so, weil das ist die Differenz der Gruppen bei Baseline!

Was ist jetzt der Unterschied zwischen einem linearen Modell und einem mixed effects model?

Wie man unten sieht, schätzt ein lineares Modell die gleichen Regressionskoeffizienten, kommt aber zu ganz anderen P-Werten und Konfidenzintervallen. Letztere sind für die Situation mit Messwiederholungen nicht korrekt!!!

lm.btb <- lm(bdi ~ time + treatment + time:treatment, data = data_long)
summary(lm.btb)
## 
## Call:
## lm(formula = bdi ~ time + treatment + time:treatment, data = data_long)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.414  -6.888  -1.414   6.060  31.035 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 23.448      1.904  12.318  < 2e-16 ***
## timebdi_2m                  -4.034      2.692  -1.499  0.13536    
## timebdi_4m                  -5.483      2.692  -2.037  0.04286 *  
## timebdi_6m                  -7.172      2.692  -2.664  0.00827 ** 
## treatmentBtheB              -1.966      2.692  -0.730  0.46607    
## timebdi_2m:treatmentBtheB   -7.034      3.807  -1.848  0.06596 .  
## timebdi_4m:treatmentBtheB   -6.138      3.807  -1.612  0.10832    
## timebdi_6m:treatmentBtheB   -5.069      3.807  -1.331  0.18439    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.25 on 224 degrees of freedom
## Multiple R-squared:  0.2099, Adjusted R-squared:  0.1853 
## F-statistic: 8.504 on 7 and 224 DF,  p-value: 3.144e-09