---
title: "Kovarianzanalyse"
author: "Übungsaufgabe zu Analyse und Dokumentation SoSe 2024"
bibliography: 8-Referenzen.bib
lang: de
format: 
  pdf:
    include-in-header:
      text: |
       \usepackage[font=small,format=plain,labelfont=bf, labelsep=period,justification=justified,singlelinecheck=false]{caption}
---

Inspiration für diese Übungsaufgabe ist die Studie von @blackburn1981. Ziel ist es, 
mithilfe einer Kovarianzanalyse von Pre-Post-Therapie-BDI-Differenzwerten
die Effekte von Kognitiver Verhaltenstherapie und Pharmakotherapie vor dem 
Hintergrund der Dauer einer bestehenden Depressionssymptomatik bei Therapiebeginn 
zu quantizifizieren. 

## Datensatz {-}

Der Datensatz `8-Kovarianzanalyse.csv` enthält simulierte Daten von insgesamt 
$n = 40$ Patient:innen als Zeilen. Die erste Spalte des Datensatzes enthält eine 
patient:innenspezifische Therapievariante (`CBT`: Kognitive Verhaltenstherapie, 
`PHC`: Pharmakotherapie), die zweite Spalte die patient:innenspezifische Dauer 
der Depressionsymptomatik zu Therapiebeginn in Monaten und die dritte Spalte 
schließlich die Pre-Post-Therapie-BDI-Differenzwerten (`BDI`). @tbl-bdi zeigt 
exemplarisch die Daten der ersten acht Patient:innen jeder Therapiegruppe.

```{r, echo = F, eval = F}
#| label: Datensimulation
library(MASS)                                                                  # Multivariate Normalverteilung
set.seed(1)                                                                    # Ergebnisreproduzierbarkeit
n_i       <- c(20,20)                                                          # Anzahl Patient:innen pro Gruppe
p         <- 3                                                                 # Anzahl Betaparameter
n         <- sum(n_i)                                                          # Gesamtzahl Datenpunkte
beta      <- matrix(c(16,0,-4), nrow = p)                                      # wahrer, aber unbekannter, Betaparameter
sigsqr    <- 10                                                                # Varianzparameter 
x         <- c(rnorm(n_i[1],2,1), rnorm(n_i[2],4,1))                           # Depressionsdauerkovariate
X         <- matrix(c(rep(1,n_i[1]), rep(1,n_i[2]),                            # \mu_0    Regressor
                    rep(0,n_i[1]), rep(1,n_i[2]),                              # \alpha_2 Regressor
                    x),                                                        # \beta_1  Regressor
               nrow = n)                                                       # Anzahl Zeilen
y        <- mvrnorm(1, X %*% beta, sigsqr*diag(n))                             # Datenrealisierung
y        <- round(y)                                                           # diskrete BDI Werte
D        <- data.frame(THP  = c(rep("CBT", n_i[1]),                            # CBT Bedingung
                                rep("PHC", n_i[2])),                           # PHC Bedingung 
                       DUR  = x,                                               # Depressionsdauer vor Therapie
                       BDI  = y)                                               # BDI Werte                   
fname    <- file.path("8-Kovarianzanalyse.csv")                                # Dateiname
write.csv(D, file = fname, row.names = FALSE)                                  # Speichern 
```

\footnotesize
```{r echo = F}
#| label: tbl-bdi
#| tbl-cap : "Auszug des Beispieldatensatzes"
fname    <- file.path("8-Kovarianzanalyse.csv")           # Dateiname
D        <- read.table(file.path(fname), sep = ",", header = TRUE)             # Laden des Datensatzes
knitr::kable(D[c(1:8,21:28),], digits = 2, align = "c")                        # erste acht Patient:innen jeder Therapiebedingung
``` 
\normalsize


\newpage
## Programmieraufgaben {-}
\noindent 1. Bestimmen Sie zunächst die Stichprobenmittel der `BDI` und `DUR`
Werte für jede Therapievariante. Bestimmen Sie dann für ein lineares Modell aus 
einem Offset-Regressor und einem kategorialen Therapieeffektregressor die Beta- und 
Varianzparameterschätzer sowie die T-Statistk für den Therapieeffektregressor (`ALM 1`).
Bestimmen Sie dann für ein lineares Modell aus einem Offset-Regressor, einem kategorialen
Therapieeffektregressor und dem kontinuierlichen Regressor der `DUR` Werte die Beta- 
und Varianzparameterschätzer sowie die T-Statistk für den Therapieeffektregressor
(`ALM 2`). Sie sollten folgende Ergebnisse erhalten.

\footnotesize
```{r, echo = F}
#| label: Kovarianzanalyse (Theorem)

# Laden des Datensatzes
fname       <- "8-Kovarianzanalyse.csv"                                        # Dateiname
D           <- read.table(file.path(fname), sep = ",", header = TRUE)          # Laden des Datensatzes
n_i         <- c(20,20)                                                        # Anzahl Datenpunkte pro Gruppe
y           <- D$BDI                                                           # Daten
n           <- length(y)                                                       # Gesamtanzahldatenpunkte
XS          <- list()                                                          # Designmatrixliste
CS          <- list()                                                          # Kontrastgewichtsvektorenliste
XS[[1]]     <- matrix(c(rep(1,n_i[1]),rep(1,n_i[2]),                           # \mu_0
                        rep(0,n_i[1]),rep(1,n_i[2])),                          # \alpha_2
                      nrow = n)
XS[[2]]     <- matrix(c(rep(1,n_i[1]),rep(1,n_i[2]),                           # \mu_0
                        rep(0,n_i[1]),rep(1,n_i[2]),                           # \alpha_2
                        D$DUR),
                      nrow = n)                                                # x
CS[[1]]     <- matrix(c(0,1)  , nrow = 2)                                      # Kontrastgewichtsvektor \alpha_2
CS[[2]]     <- matrix(c(0,1,0), nrow = 3)                                      # Kontrastgewichtsvektor \alpha_2
B           <- list()                                                          # Betaparameterschätzerliste
S           <- rep(NaN,2)                                                      # Varianzparameterschätzervektor
T           <- rep(NaN,2)                                                      # T-Statistikvektor
for(i in 1:2){                                                                 # ALM Iterationen
  X           <- XS[[i]]                                                       # Designmatrix
  p           <- ncol(X)                                                       # Anzahl Betaparameter
  beta_hat    <- solve(t(X) %*% X) %*% t(X) %*% y                              # Betaparameterschätzer vollständiges Modell
  eps_hat     <- y - X %*% beta_hat                                            # Residuenvektor 
  sigsqr_hat  <- (t(eps_hat) %*% eps_hat)/(n-p)                                # Varianzparameterschätzer
  c           <- CS[[i]]                                                       # Kontrastgewichtsvektor
  t_num       <- t(c) %*% beta_hat                                             # Zähler der Zweistichproben-T-Teststatistik
  t_den       <- sqrt(sigsqr_hat %*% t(c) %*% solve(t(X) %*% X) %*% c)         # Nenner der Zweistichproben-T-Teststatistik
  t           <- t_num/t_den                                                   # Wert der Zweistichproben-T-Teststatistik
  B[[i]]      <- beta_hat                                                      # Betaparameterschätzer
  S[i]        <- sigsqr_hat                                                    # Varianzparameterschätzer
  T[i]        <- t                                                             # T-Statistik
}
cat("Stichprobenmittel BDI (CBT, PHC) : ",   round(c(mean(D$BDI[D$THP == "CBT"]), mean(D$BDI[D$THP == "PHC"])),2),
    "\nStichprobenmittel DUR (CBT, PHC) : ", round(c(mean(D$DUR[D$THP == "CBT"]), mean(D$DUR[D$THP == "PHC"])),2),  
    "\nBetaparameterschätzer      ALM 1 : ", round(B[[1]],2),
    "\nBetaparameterschätzer      ALM 2 : ", round(B[[2]],2),
    "\nVarianzparameterschätzer   ALM 1 : ", round(S[[1]],2),
    "\nVarianzparameterschätzer   ALM 2 : ", round(S[[2]],2),
    "\nT-Statistik                ALM 1 : ", round(T[[1]],2),
    "\nT-Statistik                ALM 2 : ", round(T[[2]],2))
```

\normalsize
\noindent 2. Visualisieren Sie die Stichprobenmittel- und Standardabweichungen 
sowie die einzelnen Datenpunkte der `BDI `Werte für jede Therapievariante ("Rohdaten"). 
Subtrahieren Sie dann von den `BDI` Werten die aufgrund des Betaparameterschätzers
für den `DUR` Regressor durch `ALM 2` prädizierten `BDI` Anteile und visualisieren
Sie die so adjustierten `BDI` Werte ("Adjustierte Daten"). Ihre Abbildungen sollten
in etwa aussehen wie @fig-abbildung.

```{r, echo = F, eval = F}
#| label: Visualisierung (low-level)

graphics.off()
library(latex2exp) 
pdf(
  file        = file.path("8-Kovarianzanalyse-Abbildung.pdf"),
  width       = 10,
  height      = 5
)
par(
  mfcol       = c(1,2),
  family      = "sans",
  pty         = "s",
  bty         = "l",
  mgp         = c(2,1,0),
  lwd         = 1,
  las         = 1
)
fname       <- file.path("8-Kovarianzanalyse.csv")
D           <- read.table(fname, sep = ",", header = TRUE)
Y           <- data.frame(
              CBT = D$BDI[D$THP == "CBT"],
              PHC = D$BDI[D$THP == "PHC"])
lab         <- c("Rohdaten", "Adjustierte Daten")

# Iteration über unkorrigierte und korrigierte Daten
for(i in 1:2){

  # Datenadjustierung
  if(i == 2){
    Y$CBT <- Y$CBT - X[1:20 ,3] * beta_hat[3]
    Y$PHC <- Y$PHC - X[21:40,3] * beta_hat[3]
  }
 
  # Gruppenmittelwerte und Standardabweichungen
  groupmeans  <- colMeans(Y)
  groupstds   <- apply(Y,2,sd)
  x <- barplot(
    groupmeans,
    ylim        = c(-8,22),
    col         = "gray90",
    ylab        = "BDI",
    xlab        = "Therapievariante",
    main        = lab[i]
  )
  arrows(
    x0          = x,
    y0          = groupmeans - groupstds,
    x1          = x,
    y1          = groupmeans + groupstds,
    code        = 3,
    angle       = 90,
    length      = 0.05
  )
  
  # Einzeldatenpunkte
  for(i in 1:ncol(Y)){
    points(
      jitter(rep(x[i], length(Y[,i])), amount = .3),
      Y[,i],
      pch = 21,
      col = "White",
      bg  = "Black"
    )
  }
}
dev.off()
```

![Rohdaten und Adjustierte Daten](8-Kovarianzanalyse-Abbildung.pdf){#fig-abbildung fig-align="center"}  

\noindent 3. Zeigen Sie in einer kurzen R-Übung, wie Sie die in Aufgabe 1 bestimmten Betaparameterschätzer und T-Teststatistiken mithilfe der **R** Funktionen `stats::lm()` und `afes::aov_ez()` bestimmen können.

```{r echo = F, eval = F}
#| label: tidyverse style + afex
library(tidyverse)                                                             # Für Pipe (%>%), mutate(), filter()
library(ggplot2)                                                               # Für ggplot()
library(httpgd)                                                                # Für hgd() und hgd_browse()
library(afex)                                                                  # Für aov_ez()
library(ggpubr)                                                                # Für ggarrange()


# Daten vorbereiten
fname      <- "6_Übungen/8_Kovarianzanalyse/8-Kovarianzanalyse.csv"            # Dateiname
fname      <- "8-Kovarianzanalyse.csv"                                         # Dateiname
D          <- read.table(fname, sep = ",", header = TRUE)                      # Laden des Datensatzes
n_subs     <- nrow(D)
D_processed <- D %>%
  mutate(ID = seq(n_subs)                                                      # Subject ID hinzufügen
  )

# Mittelwerte bestimmen
means <- D_processed %>%
  group_by(THP) %>%                                                            # Nach THP gruppieren
  summarize(                                                                   # Dataframe mit Spalte für mean erstellen
    mean_BDI = mean(BDI),
    mean_DUR = mean(DUR)
  )

# Beta-parameterschätzer
beta_hats = list()
beta_hats <- D_processed %>%
  summarize(
    alm1_beta = list(lm(BDI ~ THP)$coefficients),
    alm2_beta = list(lm(BDI ~ THP + DUR)$coefficients)
  )
alm1_fit <- lm(BDI ~ THP, data = D_processed)                                  # ALM 1, nur eine kategoriale UV
alm2_fit <- lm(formula = BDI ~ THP + DUR, data = D_processed)                  # ALM 2, eine kategoriale UV plus eine kontinuierliche Kovariate


# Kovarianzanalyse mit aov_ez()
alm1 <- aov_ez(                                                                # ALM 1, nur eine kategoriale UV
  id = "ID", dv = "BDI", data = D_processed, between = "THP")
alm2 <- aov_ez(                                                                # ALM 2, eine kategoriale UV
  id = "ID", dv = "BDI", data = D_processed, between = "THP",
    covariate = "DUR",                                                         # plus eine covariate
    factorize = FALSE                                                          # verhindern, dass var. faktorisiert werden
)
# T-Statistiken für THP
t_stat_alm1 <- sqrt(alm1$anova_table$F)
t_stat_alm2 <- sqrt(alm2$anova_table$F[1])


# Ausgabe der Ergebnisse
cat("Stichprobenmittel BDI (CBT, PHC) : ",   round(c(
                                                filter(means, THP == "CBT")$mean_BDI,
                                                filter(means, THP == "PHC")$mean_BDI),
                                              digits = 2),
    "\nStichprobenmittel DUR (CBT, PHC) : ", round(c(
                                                filter(means, THP == "CBT")$mean_DUR,
                                                filter(means, THP == "PHC")$mean_DUR),
                                              digits = 2),
    "\nBetaparameterschätzer      ALM 1 : ", round(unlist(beta_hats$alm1_beta), 2),
    "\nBetaparameterschätzer      ALM 2 : ", round(unlist(beta_hats$alm2_beta), 2),
    "\nT-Statistik                ALM 1 : ", round(t_stat_alm1, 2),
    "\nT-Statistik                ALM 2 : ", round(t_stat_alm2, 2)
  )


# ------ Visualisierung -------------------------------------------------------
# # Start plot viewer
# hgd()
# hgd_browse()
plot_list <- list()                                                            # leere Liste, um plots zu speichern
# Designmatrix erstellen
n_i <- c(20, 20)                                                               # Anzahl Datenpunkte pro Gruppe
n <- nrow(D_processed)
X <- matrix(c(rep(1, n_i[1]), rep(1, n_i[2]),                                  # \mu_0
              rep(0, n_i[1]), rep(1, n_i[2]),                                  # \alpha_2
              D_processed$DUR),                                                # x
            nrow = n)

# Iteration über unkorrigierte und korrigierte Daten
for (i in c("Rohdaten", "Adjustierte Daten")){

  if (i == "Rohdaten"){
    data_set <- D_processed
  } else {
    data_set <- D_processed %>%
      mutate(BDI = BDI - X[, 3] * beta_hats$alm2_beta[[1]]["DUR"]
      )
  }
  print(i)
  print(data_set)

  # Deskriptive Statistiken
  summary_stats <- data_set %>%
    group_by(THP) %>%                                                          # Nach Treatment gruppieren
    summarize(                                                                 # Kreiert dataframe mit eine Zeile für jede Gruppe
      Mean = mean(BDI),
      Var = var(BDI),
      Std = sd(BDI)
    )

  #plot[i] <- # Erstellen des Plots mit ggp und facet_wrap
  plot_list[[i]] <- ggplot(summary_stats, aes(x = THP, y = Mean)) +
    geom_bar(stat = "identity", fill="lightblue") +
    geom_errorbar(
      aes(ymin = Mean - Std, ymax = Mean + Std),
      width = 0.2, position = position_dodge(width = 0.9)) +
    geom_jitter(data = data_set, aes(y = BDI), width = 0.2, height = 0, shape = 21, fill = "black", color = "white") +
    ylim(-8, 22) +
    labs(x = "Therapievariante", y = "BDI", title = i)
}

gesamt_plot <- ggarrange(                                                      # Plots in eine Grafik zusammenfügen
  plot_list[["Rohdaten"]],                                                     # linkes Panel
  plot_list[["Adjustierte Daten"]] + rremove("y.text") + rremove("ylab"),      # rechtes Panel
  align = "hv",                                                                # Grid aligned
  ncol = 2, nrow = 1)

ggsave(                                                                        # Abbildung speichern
  filename = "8-Kovarianzanalyse-ggplot.pdf",
  width    = 10,
  height   = 5
  )

```


\newpage
## Dokumentation {-}

Bitte beachten Sie bei der Erstellung Ihre Dokumentation folgende Vorgaben und 
orientieren Sie sich in der Darstellung Ihrer datenanalytischer Ergebnisse an
den Empfehlungen des [APA Publication Manuals 7th Edition](https://apastyle.apa.org/products/publication-manual-7th-edition), insbesondere Kapitel 6. 


### Einleitung {-}

Stellen Sie kurz die Ausgangsfrage der Übungsaufgabe dar und erläutern Sie die
Therapievarianten, die `DUR` Variable und die `BDI` Variable. Diskutieren Sie
die Therapie-spezifischen Stichprobenmittel der `BDI` und `DUR` Variablen und 
erläutern Sie, inwiefern es sich bei dem vorliegenden Datensatz nicht um das Resultat
einer randomisierten Studie handelt.

### Methoden {-}

Beschreiben Sie die linearen Modelle, die Sie zur Lösung der ersten Programmieraufgabe
herangezogen haben und stellen Sie die Ajdustierung der Daten zur Lösung der 
zweiten Programmieraufgabe dar. 

### Resultate {-}

Diskutieren Sie die Größen des Betaparameterschätzers und der T-Statistik
für den Therapieeffektregressor im Kontext von `ALM 1` und `ALM 2`. Betrachten 
Sie dabei auch den Varianzparameterschätzer. Erläutern Sie Ihre Ergebnisse im Kontext
des Zusammenspiels der Dauer einer bestehenden Depressionssymptomatik bei 
Therapiebeginn und der in diesem Datensatz abgebildeten Therapieeffektivität.

### Schlußfolgerung {-}

Fassen Sie die von Ihnen erstellte Dokumentation in drei Sätzen zusammen.

### Referenzen {-}
\footnotesize