---
title: "Kovarianzanalyse"
author: "Übungsaufgabe zu Design, Analyse, Dokumentation SoSe 2023"
bibliography: referenzen.bib
lang: en
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}
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}
# 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])), nrow = n)                  # \alpha_2 
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}
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"}  

\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