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

Grundlage dieser Übung ist die Studie von @wagner2014. Ziel ist es, mithilfe eines
Zweistichproben-T-Tests zu quantifizieren, inwieweit sich die Veränderung der 
Depressionssymptomatik im Verlaufe einer Psychotherapie reliabel zwischen einer 
*Online Studiengruppe* ($n = 25$) und einer *Face-To-Face Studiengruppe* ($n = 28$) 
unterscheidet. Zum Zwecke dieser Übung fokussieren wir auf den *Beck Depression Inventory (BDI)* 
Wert als Ergebnismaß der Studie von @wagner2014.

### Datensatz {-}

Der Datensatz `4-Zweistichproben-T-Test.csv` enthält als Spalten simulierte BDI
Werte zu den Erhebungszeitpunkten *Pre* und *Post* der psychotherapeutischen 
*Online* und *Face-to-Face* Intervention. @tbl-bdi zeigt exemplarisch die Daten
von fünf Patient:innen jeder Studiengruppe.

```{r, echo = F, eval = T}
#| label: Datensimulation
library(MASS)                                                                  # Multivariate Normalverteilung
set.seed(2)                                                                    # Ergebnisreproduzierbarkeit
n_1        <- 25                                                               # Anzahl Patient:innen Online Gruppe
n_2        <- 28                                                               # Anzahl Patient:innen Face-to-Face Gruppe
n          <- n_1 + n_2                                                        # Gesamtzahl Datenpunkte pro Messpunkt
n_m        <- 2                                                                # Anzahl Messzeitpunkte Pre und Post
betas      <- matrix(c(22,23, 12, 12), nrow = 2)                               # Erwartungswertparameter Online/Face-to-Face und Pre/Post
sigsqr     <- 8                                                                # VarianzparameterOnline/Face-to-Face und Pre/Post
Y           <- matrix(rep(NaN,n*2), nrow = n)                                  # Datensatzarrayinitialisierung
X           <- matrix(c(rep(1,n_1), rep(0,n_2),                                # Zweistichproben-T-Test Designmatrix 
                       rep(0,n_1), rep(1,n_2)),
                       nrow = n)
for(i in 1:n_m){                                                               # Pre-Post Iteration                
    Y[,i] <- mvrnorm(1, X %*% betas[,i], sigsqr*diag(n))                       # Datenrealisierung
}
Y           <- round(Y)                                                        # diskrete BDI Werte
Y[Y < 0]    <- 0                                                               # natürliche BDI Werte
D           <- data.frame(Condition  = c(rep("Online"      , n_1), 
                                        rep("Face-to-Face", n_2)),
                         Pre        = Y[,1],                                   # Pre BDI Werte
                         Post       = Y[,2])                                   # Post BDI Werte 
fname       <- file.path("4-Zweistichproben-T-Test.csv")                       # Dateiname
write.csv(D, file = fname, row.names = FALSE)                                  # Speichern 
```

```{r echo = F, eval = T, warning = F}
#| label: tbl-bdi
#| tbl-cap : "Exemplarische Pre- und Post-Intervention BDI Werte der Studiengruppen."
library(knitr)                                                                 # knitr für Tabellen
fname       <- "4-Zweistichproben-T-Test.csv"                                  # Dateiname
D           <- read.table(file.path(fname), sep = ",", header = TRUE)          # Laden des Datensatzes
kable(head(D[c(1:5,31:35),], n = 10L), digits = 2, align = "c")                # Markdowntabellenoutput für head(D)
``` 

\newpage
## Programmieraufgaben {-}

\noindent 1. Bestimmen Sie die Differenzen der Pre und Post BDI Werte für beide Studiengruppen. 
Führen Sie dann basierend auf diesen Differenzwerten einen zweiseitigen 
Zweistichproben-T-Test mit Nullhypothesenparameter $\mu_0 = 0$ durch. Bestimmen 
sie dabei insbesondere die Beta- und Varianzparameterschätzer des Zweistichproben-T-Testmodells, 
den Wert der Zweistichproben-T-Teststsatitik, sowie den korrespondierenden p-Wert. 
Geben Sie weiterhin das 95\%-Konfidenzintervall für den Erwartungswert der 
Pre-Post-Testdifferenzen an. Bestimmen Sie schließlich unter der Annahme, dass die 
Werte der Erwartungswert- und Varianzparameterschätzer den wahren, aber unbekannten, 
Parametern gleichen, die Wahrscheinlichkeit dafür, dass der Zweistichproben-T-Test bei den 
Stichprobengröße von $n_1 = 25$ und $n_2 = 29$ und einem kritischen Wert, der 
einem Signifikanzlevel von $\alpha_0 := 0.05$ entspricht, den Wert 1 annimmt. 
Diese geschätzte Wahrscheinlichkeit wird manchmal als *Post-hoc power* bezeichnet. 
Sie sollten folgende Ergebnisse erhalten:

\footnotesize
```{r echo = F, eval  = T}
#| label: t-Test (Theorem)
# Datenanalyse
fname       <- "4-Zweistichproben-T-Test.csv"                                  # Dateiname
D           <- read.table(file.path(fname), sep = ",", header = TRUE)          # Laden des Datensatzes
n_1         <- 25                                                              # Anzahl Patient:innen Online Gruppe
n_2         <- 28                                                              # Anzahl Patient:innen Face-to-Face Gruppe
n           <- n_1 + n_2                                                       # Gesamtzahl Datenpunkte pro Messpunkt
y           <- D$Post - D$Pre                                                  # Post-Pre Differenzwerte
p           <- 2                                                               # Anzahl Betaparameter
c           <- matrix(c(1,-1), nrow = 2)                                       # Kontrastgewichtsvektor
delta       <- 0.95                                                            # Konfidenzlevel
mu_0        <- 0                                                               # Nullhypothesenparameter       
alpha_0     <- 0.05                                                            # Signifikanzlevel  
X           <- matrix(c(rep(1,n_1),rep(0,n_2),rep(0,n_1),rep(1,n_2)),nrow=n)   # Zweistichproben-T-Test Designmatrix           
beta_hat    <- solve(t(X)%*%X)%*%t(X)%*%y                                      # Betaparameterschätzer
eps_hat     <- y - X %*% beta_hat                                              # Residuenvektor 
sigsqr_hat  <- (t(eps_hat) %*% eps_hat)/(n-p)                                  # Varianzparameterschätzer
t_delta     <- qt((1+delta)/2,n-2)                                             # Psi^{-1}((1+delta)/2,n-2)
lambda      <- diag(solve(t(X) %*% X))                                         # \lambda_j Werte
kappa       <- matrix(rep(NaN, p*2), nrow = p)                                 # Konfidenzintervallarray
for(j in 1:p){                                                                 # \beta_j Iteration
  kappa[j,1] <- beta_hat[j]-sqrt(sigsqr_hat*lambda[j])*t_delta                 # untere Konfidenzintervallgrenze
  kappa[j,2] <- beta_hat[j]+sqrt(sigsqr_hat*lambda[j])*t_delta                 # obere Konfidenzintervallgrenze
}
t_num       <- t(c) %*% beta_hat - mu_0                                        # 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
pval        <- 2*(1 - pt(abs(t), n-2))                                         # p-Wert bei zweiseitigem Zweistichproben-T-Test
k_alpha_0   <- qt(1-alpha_0/2, n-2)                                            # kritischer Wert 
p_delta_hat <- 1-pt(k_alpha_0, n-2, t)+pt(-k_alpha_0, n-2, t)                  # "Post-hoc power"

# Ausgabe
cat("                 Betaparameterschätzer           : ", beta_hat,
    "\n                 95%-Konfidenzintervall beta_1   : ", kappa[1,],
    "\n                 95%-Konfidenzintervall beta_2   : ", kappa[2,],
    "\n                 Varianzparameterschätzer        : ", sigsqr_hat,
    "\n                 Zweistichproben-T-Teststatistik : ", t,
    "\n                 p-Wert                          : ", pval,
    "\n                 Post-hoc power                  : ", p_delta_hat)
```
\normalsize

\noindent 2. Visualisieren Sie die entsprechenden Gruppenmittelwerte als 
Linienplots mit Fehlerbalken analog zu Figure 2 in @wagner2014. Visualisieren außerdem
die Post-Pre-Differenz Werte als gruppenspezifische *Violinplots* mithilfe des
*R* Pakets `vioplot`.  Die Abbildung sollte in etwa aussehen wie @fig-abbildung.

```{r echo = F, eval = F}
#| label: Visualisierung (low-level)
# Datenanalysevisualisierung
# install.packages("vioplot")
library(latex2exp)
library(vioplot)
fname       <- "4-Zweistichproben-T-Test.csv"
D           <- read.table(file.path(fname), sep = ",", header = TRUE)
D$Diff      = -(D$Post - D$Pre)
pdf(
file        = file.path("4-Zweistichproben-T-Test-Abbildung.pdf"),
width       = 8,
height      = 4)
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = .8)

# Post-Pre Differenz Lineplot
x           = 1:2
groupmeans  = matrix(c(mean(D$Pre[D$Condition   == "Face-to-Face"]),
                       mean(D$Post[D$Condition  == "Face-to-Face"]),
                       mean(D$Pre[D$Condition   == "Online"      ]),
                       mean(D$Post[D$Condition  == "Online"      ])),nrow = 2)
groupstds   = matrix(c(sd(D$Pre[D$Condition   == "Face-to-Face"]),
                       sd(D$Post[D$Condition  == "Face-to-Face"]),
                       sd(D$Pre[D$Condition   == "Online"      ]),
                       sd(D$Post[D$Condition  == "Online"      ])),nrow = 2)
cols        = c("gray30", "gray70")
lwds        = c(3,1)
matplot(
x,
groupmeans,
type        = "b",
pch         = 21,
lty         = 1,
col         = cols,
lwd         = lwds,
xlab        = "",
ylab        = "BDI",
xaxt        = "n",
xlim        = c(0.5,2.5),
ylim        = c(5,30))
for(i in 1:2){
  arrows(                                          
  x0          = x,                                 
  y0          = groupmeans[,i] - groupstds[,i],    
  x1          = x,                                
  y1          = groupmeans[,i] + groupstds[,i],    
  col         = cols[[i]],                         
  lwd         = lwds[[i]],                         
  code        = 3,                                 
  angle       = 90,                                
  length      = 0.05)
}
legend(
  "bottomleft",
  c("Face-to-Face", "Online"),
  lty         = 1,
  col         = c("gray30", "gray90"),
  lwd         = lwds,
  bty         = "n",
  cex         = 1,
  x.intersp  = .3,
  seg.len     = 0.6
)
text(1,1.7, "Pretest" , xpd = T, cex = 1.1)
text(2,1.7, "Posttest", xpd = T, cex = 1.1)

# Post-Pre Differenzwert Violinplot
vioplot(
  D$Diff ~ D$Condition,
  D,
  col         = "gray80",                          
  rectCol     = "black",                           
  lineCol     = "white",                           
  colMed      = "gray80",                         
  border      = "black",                           
  pchMed      = 16,                                
  plotCentre  = "points",
  ylab        = "BDI Differenzwerte",
  ylim        = c(0,30),
  xlab        = ""
)
stripchart(
  D$Post ~ D$Condition,
  D,
  method      = "jitter", 
  xaxt        = "n",
  vertical    = TRUE, 
  pch         = 21,
  bg          = "Black",
  col         = "gray80",
  add         = TRUE
)
dev.off()
```

![Post-Pre BDI Differenz Gruppenanalyse.](4-Zweistichproben-T-Test-Abbildung.pdf){#fig-abbildung fig-align="center" width=70%}  

\noindent 3. Zeigen Sie in einer kurzen R-Übung wie Sie mit einer \href{https://style.tidyverse.org/pipes.html#pipes}{tydiverse Pipe} die Post-Pre-Differenz der BDI-Werte berechnen und diese direkt als neue Spalte in das Dataframe integrieren können. Konsultieren Sie hierfür auch die Einführung zu Data transformation in \href{https://r4ds.hadley.nz/data-transform#sec-the-pipe}{R for Data Science}. Zeigen Sie weiterhin, wie Sie die Post-Pre Differenz BDI Werte als *Violinplot* mithilfe des **R** Pakets `ggplot2` visualisieren können.

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

# Daten vorbereiten
fname     <- "6_Übungen/4_Zweistichproben-T-Test/4-Zweistichproben-T-Test.csv" # Dateiname
fname     <- "4-Zweistichproben-T-Test.csv"                                    # Dateiname
D           <- read.table(fname, sep = ",", header = TRUE)                     # Laden des Datensatzes
D_processed <- D %>%
  mutate(Diff = as.numeric(Post - Pre)) %>%                                    # Neue Variable berechnen
  mutate(Diff_betrag = as.numeric(-(Post - Pre)))

# T-Test durchführen
alpha_0 <- 0.05
t_test_results <- t.test(
        Diff_betrag ~ Condition, data = D_processed,
        var.equal = TRUE,
        alternative = "two.sided",
        conf.level = 1 - alpha_0
        )
print(t_test_results)                                                          # Ergebnisse ausgeben

# Violinplots erstellen
# # Start plot viewer
# hgd()
# hgd_browse()
# Violinplot erstellen
ggplot(data = D_processed, aes(x = "", y = Diff_betrag, fill = Condition)) +
  geom_violin(trim = T) +                                                      # Violinplots
  geom_point(                                                                  # Datenpunkte
    aes(y = Diff_betrag),
    position = position_jitter(width = 0.15, height = 0),
    size = 2,
    fill = "black"
    ) +
  labs(x = "", y = "Differenz der Werte", title = "Face-to-Face vs. Online") + # Label hinzufügen
  facet_wrap(~ Condition, ncol = 2) +                                          # Panel nach Condition in zwei Spalten soirtern
  theme(legend.position = "none")                                              # Legende ausblenden, da die Füllfarbe schon die Gruppen darstellt

ggsave(                                                                        # Abbildung speichern
  filename = "4-Zweistichproben-T-Test-ggplot.pdf",
  width    = 5,
  height   = 5
  )
```

\newpage
## Dokumentation

Bitte beachten Sie bei der Erstellung Ihre Dokumentation folgende Vorgaben.

### Einleitung {-}

Stellen Sie die Ausgangsfrage von  @wagner2014 dar und erläutern Sie kurz die
Therapieprinzipien der *Online* und der *Face-to-Face* Studiengruppen. 

### Methoden {-}

Beschreiben Sie die Patient:innen- und Therapiebedingungsgruppen. Erläutern Sie 
kurz die Logik der Anwendung eines Zweistichproben-T-Tests bei unabhängigen Stichproben.
Dokumentieren Sie Ihre Datenanalyse in Form kommentierten **R** Codes zur Lösung
von Programmieraufgabe 1.

### Resultate {-}

Reportieren Sie die von Ihnen bestimmten Statistiken aus Programmieraufgabe 1
und beziehen Sie zur Validität der Nullhypothese $\mu_0 = 0$ Stellung. Kommentieren
Sie weiterhin vor diesem Hintergrund den resultierenden Wert der Post-hoc Power.
Beschreiben Sie die in Programmieraufgabe 2 erstellte Abbildung.

### Schlußfolgerung {-}

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

## Referenzen
\footnotesize