---
title: "Einstichproben-T-Test"
author: "Übungsaufgabe zu Analyse und Dokumentation SoSe 2024"
bibliography: 3-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
Einstichproben-T-Tests zu quantifizieren, inwieweit sich die Depressionssymptomatik
einer Gruppe von Patient:innen zwischen Beginn und Ende einer Psychotherapie verändert 
hat. Diese Anwendung eines Einstichproben-T-Tests wird häufig auch auch als 
*Zweistichproben-T-Test bei abhängigen Stichproben* bezeichnet. Zum Zwecke dieser 
Übung fokussieren wir auf die *Online Studiengruppe* ($n = 32$) und den 
*Beck Depression Inventory (BDI)* Wert als Ergebnismaß der Studie von @wagner2014.

### Datensatz {-}

Der Datensatz `3-Einstichproben-T-Test.csv` enthält als Spalten simulierte BDI
Werte zu den Erhebungszeitpunkten *Pre* und *Post* der psychotherapeutischen 
*Online* Intervention. @tbl-bdi zeigt exemplarisch die Daten der ersten zwölf Patient:innen.

```{r, echo = F, eval = F}
#| label: Datensim
library(MASS)                                                                  # Multivariate Normalverteilung
set.seed(10)                                                                   # Ergebnisreproduzierbarkeit
n           <- 32                                                              # Anzahl Patient:innen
k           <- 2                                                               # Anzahl Gruppen (Pre/Post)
betas       <- c(23,12)                                                        # Erwartungswertparameter \mu_{pre}/\mu_{post}
sigsqrs     <- c(6 ** 2, 10 ** 2)                                              # Varianzparameter \sigma^2_{pre}/\sigma^2_{post}
Y           <- matrix(rep(NaN,n*2), nrow = n)                                  # Datensatzarrayinitialisierung
X           <- matrix(rep(1,n), nrow = n)                                      # Designmatrix 
for(i in 1:k){                                                                 # Pre-Post Iteration                
    Y[,i] <- mvrnorm(1, X %*% betas[i], sigsqrs[i]*diag(n))                    # Datenrealisierung
}
Y           <- round(Y)                                                        # diskrete BDI Werte
Y[Y < 0]    <- 0                                                               # natürliche BDI Werte
D           <- data.frame(Pre  = Y[,1],                                        # Pre BDI Werte
                         Post = Y[,2])                                         # Post BDI Werte 
fname       <- file.path("3-Einstichproben-T-Test.csv")                        # Dateiname
write.csv(D, file = fname, row.names = FALSE)                                  # Speichern 
```

```{r echo = F, warning = F}
#| label: tbl-bdi
#| tbl-cap : "Pre- und Post-Intervention BDI Werte"
library(knitr)                                                                 # knitr für Tabellen
fname       <- "3-Einstichproben-T-Test.csv"                                   # Dateiname
D           <- read.table(file.path(fname), sep = ",", header = TRUE)          # Laden des Datensatzes
kable(head(D, n = 12L), digits = 2, align = "c")                               # Markdowntabellenoutput für head(D)
``` 

\newpage
## Programmieraufgaben {-}

\noindent 1. Bestimmen Sie die Differenzen der Pre und Post BDI Werte. Führen Sie dann
basierend auf diesen Differenzwerten einen zweiseitigen Einstichproben-T-Test 
mit Nullhypothesenparameter $\mu_0 = 0$ durch. Bestimmen Sie dabei insbesondere die Beta- und 
Varianzparameterschätzer des Einstichproben-T-Testmodells, den Wert der
Einstichproben-T-Teststsatitik, den korrespondierenden p-Wert und als alternative 
Effektgröße Cohens' $d$. Geben Sie weiterhin das 95\%-Konfidenzintervall für den 
Erwartungswert der Pre-Post-Differenzen 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 Einstichproben-T-Test bei einer Stichprobengröße von $n = 32$ 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, warning = F}
#| label: t-Test (Theorem)
# Datenanalyse
fname       <- "3-Einstichproben-T-Test.csv"                                   # Dateiname
D           <- read.table(file.path(fname), sep = ",", header = TRUE)          # Laden des Datensatzes
y           <- D$Post - D$Pre                                                  # Post-Pre Differenzwerte
n           <- length(y)                                                       # Anzahl Datenpunkte
p           <- 1                                                               # Anzahl Betaparameter
c           <- 1                                                               # Kontrastgewichtsvektor
mu_0        <- 0                                                               # Nullhypothesenparameter       
delta       <- 0.95                                                            # Konfidenzlevel
alpha_0     <- 0.05                                                            # Signifikanzlevel  
X           <- matrix(rep(1,n), nrow <- n)                                     # Einstichproben-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-1)                                             # \Psi^{-1}(1+\delta)/2, n-1)
lambda      <- diag(solve(t(X) %*% X))                                         # \lambda_j Werte
kappa_u     <- beta_hat - sqrt(sigsqr_hat*lambda)*t_delta                      # untere Konfidenzintervallgrenze
kappa_o     <- beta_hat + sqrt(sigsqr_hat*lambda)*t_delta                      # obere Konfidenzintervallgrenze
t_num       <- t(c) %*% beta_hat - mu_0                                        # Zähler der Einstichproben-T-Teststatistik
t_den       <- sqrt(sigsqr_hat %*% t(c) %*% solve(t(X) %*% X) %*% c)           # Nenner der Einstichproben-T-Teststatistik
t           <- t_num/t_den                                                     # Wert der Einstichproben-T-Teststatistik
pval        <- 2*(1 - pt(abs(t), n-1))                                         # p-Wert bei zweiseitigem Einstichproben-T-Test
d           <- t/sqrt(n)                                                       # Cohen's d
k_alpha_0   <- qt(1-alpha_0/2, n-1)                                            # kritischer Wert 
p_delta_hat <- 1-pt(k_alpha_0, n-1, t)+pt(-k_alpha_0, n-1, t)                  # "Post-hoc power"

# Ausgabe
cat("                 Betaparameterschätzer           : ", round(beta_hat, digits = 3),
    "\n                 95%-Konfidenzintervall          : ", round(kappa_u,digits = 3), round(kappa_o,digits = 3),
    "\n                 Varianzparameterschätzer        : ", round(sigsqr_hat, digits = 3),
    "\n                 Einstichproben-T-Teststatistik  : ", round(t, digits = 3),
    "\n                 p-Wert                          : ", round(pval, digits = 3),
    "\n                 Cohen's d                       : ", round(d,digits = 3),
    "\n                 Post-hoc power                  : ", round(p_delta_hat,digits = 3))
```
\normalsize

\noindent 2. Visualisieren Sie die Pre und Post-Interventions BDI Werte aller Patient:innen
als Liniendiagramm. Visualisieren weiterhin die Post-Pre Differenz BDI Werte
als *Violinplot* mithilfe des **R** Pakets `vioplot`. Die Abbildung sollte in 
etwa aussehen wie @fig-abbildung.

```{r echo = F, eval = F, warning = F}
#| label: Visualisierung (low-level)
# install.packages("vioplot")
library(latex2exp)
library(vioplot)
pdf(
  file        = file.path("3-Einstichproben-T-Test-Abbildung.pdf"),
  width       = 7,
  height      = 3.5
)
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
)

# Pre und Post BDI Werte
M <- rbind(matrix(D$Pre, ncol = n), matrix(D$Post, ncol = n))
matplot(
  M,
  type    = "b",
  pch     = 21,
  bg      = "Black",
  col     = "gray",
  bty     = "L",
  xlim    = c(0.8,2.2),
  ylim    = c(0,40),
  ylab    = "BDI",
  lty     = 1,
  xaxt    = "n",
  cex     = .8,
  main    = "Pre Post BDI Werte",
  xpd     = T
)
axis(side = 1, at = c(1, 2), labels = c("Pre", "Post"))

# Pre-Post-Differenzwerte Violinplot
vioplot(
  y,
  ylim        = c(-30,30),
  xlim        = c(0,2),
  col         = "gray80",
  rectCol     = "black",
  lineCol     = "white",
  colMed      = "gray80",
  border      = "black",
  pchMed      = 16,
  plotCentre  = "points",
  xaxt        = "n",
  ylab        = "BDI DIfferenzwerte",
  main        = "Pre Post BDI Differenzwerte"
)
stripchart(
  y, 
  method      = "jitter", 
  xaxt        = "n",
  vertical    = TRUE, 
  pch         = 21,
  bg          = "Black",
  col         = "gray80",
  add         = TRUE
)
dev.off()
```
![Pre-Post-Interventions BDI Werte](3-Einstichproben-T-Test-Abbildung.pdf){#fig-abbildung fig-align="center" width=90%}  

\noindent 3. Zeigen Sie in einer kurzen R-Übung wie Sie die Post-Pre Differenz BDI Werte mit einer \href{https://style.tidyverse.org/pipes.html#pipes}{tydiverse Pipe} berechnen 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/3_Einstichproben-T-Test/3-Einstichproben-T-Test.csv"  # Dateiname
fname       <- "3-Einstichproben-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

# T-Test berechnen
t_test_results <- t.test(D_processed$Diff)                                     # t-test durchführen
t_test_results                                                                 # Ergebnisse ausgeben

# Violinplot erstellen
# # Start plot viewer
# hgd()
# hgd_browse()
ggplot(data = D_processed, aes(x = "", y = Diff)) +
  geom_violin(fill = "gray80", color = "black") +                              # Violinenplit
  geom_point(                                                                  # Streudiagramm
    aes(y = Diff),
    position = position_jitter(width = 0.1, height = 0),
    size = 2,
    fill = "black"
  ) +
  labs(
    x = NULL, y = "BDI DIfferenzwerte", title = "Pre Post BDI Differenzwerte"
  )

ggsave(                                                                        # Abbildung speichern
  filename = "3-Einstichproben-T-Test-ggplot.pdf",
  height   = 5,
  width    = 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 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 Einstichproben-T-Tests im Sinne
eines Zweistichproben-T-Tests bei abhängigen Stichproben. Konsultieren Sie dazu
auch die entsprechenden Abschnitte in den  [Vorlesungsfolien](https://bit.ly/3Uk92c2) 
und das [Vorlesungsvideo](https://youtu.be/ChpvGVE5DSs) aus dem Sommersemester 2021.
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