---
fontsize: 11pt
format: pdf
bibliography: 6_Referenzen_S.bib
include-in-header: 6_Header_S.tex
lang: de
---

# (6) Einfaktorielle Varianzanalyse

Ziel dieses Seminar ist es, die Durchführung einer einfaktoriellen multivariaten
Varianzanalyse an einem Anwendungsbeispiel zu demonstrieren. Dazu betrachten wir 
die Anlayse von Prä-Post-Interventions-BDI-Score und Prä-Post-Interventions-
Glukokortikoidplasmalevel Differenzenwerten von drei Gruppen von jeweils 15 
Patienti:innen, die unterschiedliche Psychotherapiesettings (Face-to-face und Online) 
bzw. eine Wartelistenkontrollbedingung durchlaufen haben. Wir stellen dazu in 
@tbl-bdi-glu einen simulierten Beispieldatensatz dar. Die erste Spalte von 
@tbl-bdi-glu (`COND`) listet das spezifische Therapiesetting (`F2F`: Face-to-face, 
`ONL`: online, `WLC`: waitlist control) der Patient:innen auf. Die zweite Spalte
(`dBDI`) listet die entsprechenden BDI-Score-Differenzwerte  und die dritte 
Spalte (`dGLU`) die entsprechenden Glukokortikoidplasmalevel-Differenzwerte auf. 
In beiden Fällen sollen positive Werte eine Abnahme der Depressionssymptomatik, negative 
Werte dagegen einer Zunahme der Depressionssymptomatik anzeigen. 


\tiny
```{r echo = F, warning = F}
#| label: tbl-bdi-glu
#| tbl-cap : "Prä-Post-Interventions-BDI-II-Score und -Glukokortikoidplasmalevel Differenzenwerte von drei Studiengruppen (F2F: Face-to-face, ONL: online, WLC: waitlist control) jeweils 15 Patienti:innen"
library(knitr)                                                                  # knitr für Tabellen
D = read.csv("./6_Daten_S/6_Einfaktorielle_Varianzanalyse.csv")                 # Dateneinlesen
kable(D, digits = 1, align = "c")                                               # Datentabelle
``` 

\normalsize
Folgender **R** Code demonstriert die Auswertung gruppenspezifischer 
Deskriptivstatistiken für diesen Datensatz, die wir in  @fig-anwendungsbeispiel 
visualisieren.

\tiny
```{r, warning = F}
# Studiengruppenspezifische Deskriptivstatistiken
D       = read.csv("./6_Daten_S/6_Einfaktorielle_Varianzanalyse.csv")           # Dateneinlesen
m       = 2                                                                     # Datendimension von Interesse
p       = 3                                                                     # Anzahl Gruppen
k       = 15                                                                    # Anzahl Datenpunkte pro Gruppe
Y       = array(dim = c(m,k,p))                                                 # Datenarrayinitialisierung
Y[,,1]  = rbind(D$dBDI[D$COND == "F2F"],                                        # F2F dBDI Werte
                D$dGLU[D$COND == "F2F"])                                        # F2F dGLU Werte
Y[,,2]  = rbind(D$dBDI[D$COND == "ONL"],                                        # ONL dBDI Werte
                D$dGLU[D$COND == "ONL"])                                        # ONL dGLU Werte
Y[,,3]  = rbind(D$dBDI[D$COND == "WLC"],                                        # WLC dBDI Werte
                D$dGLU[D$COND == "WLC"])                                        # WLC dGLU Werte 
y_bar_i = array(dim = c(m,p))                                                   # Stichprobenmittelarray
C_i     = array(dim = c(m,m,p))                                                 # Stichprobenkovarianzmatrizenarray
j_k     = matrix(rep(1,k), nrow = k)                                            # 1_{l}
I_k     = diag(k)                                                               # Einheitsmatrix I_l
J_k     = matrix(rep(1,k^2), nrow = k)                                          # 1_{ll}
for (i in 1:p){                                                                 # Gruppeniterationen
    y_bar_i[,i] = (1/k)*(Y[,,i] %*% j_k)                                        # Stichprobenmittel \bar{\ups}_i
    C_i[,,i]    = (1/(k-1))*(Y[,,i] %*% (I_k-(1/k)*J_k) %*% t(Y[,,i]))}         # Stichprobenkovarianzmatrix C_i
```
\normalsize

```{r, echo = F, eval = F}
library(ellipse)
library(latex2exp)
pdf(
file        = "6_Abbildungen_S/mv_6_anwendungsbeispiel_s.pdf",
width       = 6,
height      = 6)
par(
family      = "sans",
mfcol       = c(1,1),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 1)
cols        = c("Gray20", "Gray50", "Gray80")
plot(
NaN,
NaN,
xlim = c(-5,15),
ylim = c(-1,6),
xlab = TeX("dBDI"),
ylab = TeX("dGLU"))
for (i in 1:p){
    points(
    Y[1,,i],
    Y[2,,i],
    col = "White",
    bg  = cols[i],
    pch = 21)
    points(
    y_bar_i[1,i],
    y_bar_i[2,i],
    col = cols[i],
    pch = 1)
    iso = ellipse(C_i[,,i], level = .6, centre = y_bar_i[,i])
    lines(
    iso[,1],
    iso[,2],
    col = cols[i])
}
legend(
"topleft",
c("F2F", "ONL", "WLC"),
pch         = c(21,21,21),
bg          = cols,
col         = cols,
bty         = "n",
cex         = 1,
x.intersp   = 1)
dev.off()
```

![Deskriptivstatistiken der `dBDI`, `dGLU` Daten des Beispieldatensatzes. Jeder Punkt visualisiert die Daten einer Patient:in.](./6_Abbildungen_S/mv_6_anwendungsbeispiel_s.pdf){#fig-anwendungsbeispiel fig-align="center" width=70%}

\newpage
Zur Parameterschätzung nutzen wir folgende **R** Funktion.

\tiny
```{r}
estimate = function(Y){

  # Diese Funktion evaluiert die Parameterschätzer einer einfaktoriellen
  # multivariaten Varianzanalyse basierend auf einen m x l x p Datensatz Y.
  #
  # Input
  #     Y          : m x k x p Datenarray
  #
  # Output
  #     $mu_hat    : m x p \mu_i Parameterschätzer
  #     $Sigma_hat : m x m \Sigma Parametschätzer
  # ---------------------------------------------------------------------------
  # Dimensionsparameter
  d         = dim(Y)                                    # Datensatzdimensionen
  m         = d[1]                                      # Datendimension
  k         = d[2]                                      # Anzahl Datenpunkte pro Gruppe
  p         = d[3]                                      # Anzahl Gruppen

  # Erwartungswertparameterschätzer
  mu_hat_i  = matrix(apply(Y,3,rowMeans), nrow = m)

  # Kovarianzmatrixparameterschätzer
  Sigma_hat = matrix(rep(0,m*m), nrow = m)
  for(i in 1:p){
      for(j in 1:k){
          Sigma_hat = Sigma_hat + (1/(k*p-p))*(Y[,j,i] - mu_hat_i[,i]) %*% t(Y[,j,i] - mu_hat_i[,i])}}


  # Outputspezifikation
  return(list(mu_hat_i = mu_hat_i, Sigma_hat = Sigma_hat))}
```
\normalsize

Zur Kreuzproduktsummenmatrixzerlegung nutzen wir folgende **R** Funktion.

\tiny
```{r}
sos = function(Y){

  # Diese Funktion evaluiert die Kreuzproduktsummenmatrizen T,B,W einer
  # einfaktoriellen Varianzanalyse basierend auf einen m x k x p Datensatz Y.
  #
  # Input
  #     Y        : m x k x p Datenarray
  #
  # Output
  #     $y_bar   : m x 1 Gesamtmittelwert
  #     $y_bar_i : m x p Gruppenmittelwerte
  #     $T       : m x m Total Sum of Squares Matrix
  #     $B       : m x m Between-Group Sum-of-Squares Matrix
  #     $W       : m x m Within  Group  Sum of Squares Matrix
  # ---------------------------------------------------------------------------
  d        = dim(Y)                                    # Datensatzdimensionen
  m        = d[1]                                      # Datendimension
  k        = d[2]                                      # Anzahl Datenpunkte pro Gruppe
  p        = d[3]                                      # Anzahl Gruppen

  # Mittelwerte
  y_bar_i  = matrix(apply(Y,3,rowMeans), nrow = m)     # Gruppenstichprobenmittel
  y_bar    = matrix(rowMeans(y_bar_i)  , nrow = m)     # Gesamtstichprobenmittel

  # Totale Sum-of-Squares Matrix
  T = matrix(rep(0,m*m), nrow = m)
  for(i in 1:p){
      for(j in 1:k){
          T = T + (Y[,j,i] - y_bar) %*% t(Y[,j,i] - y_bar)}}

  # Between Sum of Squares Matrix
  B = matrix(rep(0,m*m), nrow = m)
  for(i in 1:p){
    B = B + k*(y_bar_i[,i] - y_bar) %*% t(y_bar_i[,i] - y_bar)}

  # Within Sum of Squares Matrix
  W = matrix(rep(0,m*m), nrow = m)
  for(i in 1:p){
      for(j in 1:k){
          W = W + (Y[,j,i] - y_bar_i[,i]) %*% t(Y[,j,i] - y_bar_i[,i])}}

  # Outputspezifikation
  return(list(y_bar_i = y_bar_i, y_bar = y_bar, T = T, B = B, W = W))}
```

\normalsize
Zur Durchführung der einfaktoriellen Varianzanalyse nutzen wir folgenden **R** Code.

\tiny
```{r, message = F}
# Einlesen und Präprozessierung des Datensatzes
D           = read.csv("./6_Daten_S/6_Einfaktorielle_Varianzanalyse.csv")       # Dateneinlesen
m           = 2                                                                 # Datendimension von Interesse
p           = 3                                                                 # Anzahl Gruppen
k           = 15                                                                # Anzahl Datenpunkte pro Gruppe
n           = p*k                                                               # Gesamtanzahl Datenpunkte
Y           = array(dim = c(m,k,p))                                             # Datenarrayinitialisierung
Y[,,1]      = rbind(D$dBDI[D$COND == "F2F"],                                    # F2F dBDI Werte
                    D$dGLU[D$COND == "F2F"])                                    # F2F dGLU Werte
Y[,,2]      = rbind(D$dBDI[D$COND == "ONL"],                                    # ONL dBDI Werte
                    D$dGLU[D$COND == "ONL"])                                    # ONL dGLU Werte
Y[,,3]      = rbind(D$dBDI[D$COND == "WLC"],                                    # WLC dBDI Werte
                    D$dGLU[D$COND == "WLC"])                                    # WLC dGLU Werte         

# Einfaktorielle Varianzanalyse
alpha_0     = 0.05                                                              # Signifikanzlevel
S           = sos(Y)                                                            # Sum of Squares Matrizen
L           = det(S$W)/det(S$W + S$B)                                           # Wilks' Lambda
w           = n-1-(1/2)*(m+p)                                                   # w
t           = sqrt((m^2*(p-1)^2-4)/(m^2+(p-1)^2-5))                             # t
nu_1        = m*(p-1)                                                           # \nu_1
nu_2        = w*t-(1/2)*(m*(p-1)-2)                                             # \nu_2
tau         = ((1-L^(1/t))/L^(1/t))*(nu_2/nu_1)                                 # Teststatistik
k_alpha_0   = qf(1-alpha_0,nu_1,nu_2)                                           # kritischer Wert
phi         = as.numeric(tau > k_alpha_0)                                       # Nullhypothesentest 
P           = 1-pf(tau,nu_1,nu_2)                                               # Überschreitungswahrscheinlichkeit

# Ausgabe
cat("Wilks' Lambda     "  , L,
    "\ntau               ", tau,
    "\nnu_1              ", nu_1,
    "\nnu_2              ", nu_2,
    "\nphi               ", phi,
    "\nP(tau > tau_tilde)", P)
```
\normalsize

Im vorliegenden Fall wird die Nullhypothese identischer Gruppenerwartungswertparameter
also verworfen. Schließlich validieren wir obige Analyse im Sinne eines Black-Box-Verfahrens
mithilfe der **R** `lm()` und `Manova()` Funktionen.

\tiny
```{r, message = F}
library(car)
D        = read.csv("./6_Daten_S/6_Einfaktorielle_Varianzanalyse.csv")          # Dateneinlesen
model    = lm(cbind(D$dBDI,D$dGLU) ~ D$COND, D)                                 # Modellspezifikation
Manova(model, test.statistic = "Wilks")                                         # Einfaktorielle Varianzanalyse
```