---
format:
  pdf:
    fontsize: 11pt
    number-sections: true
    include-in-header: ../MV_Header_Seminar.tex
lang: de
---


\setcounter{section}{5}
# Multivariate Varianzanalyse


## Mathematische Grundlagen

Gegeben sei das *Modell der einfaktoriellen multivariaten Varianzanalyse*

\begin{equation} \label{eq:MANOVA}
y_{ij} = \mu_i + \varepsilon_{ij}
\quad \text{mit} \quad
\varepsilon_{ij} \sim N(0_m, \Sigma) \quad \text{u.i.v.} \; ,
\end{equation}

wobei $i = 1,\ldots,p$ ein Index über Gruppen, $j = 1,\ldots,n_i$ ein Index über experimentelle Einheiten innerhalb einer Gruppe und $m$ die Dimensionalität eines Datenvektors $y_{ij} \in \mathbb{R}^m$ ist.

Wir betrachten die *multivariate Quadratsummenzerlegung* für dieses Modell mit

\vspace{-0.5em}
\begin{equation} \label{eq:TBW}
\begin{split}
T &= \sum_{i=1}^p \sum_{j=1}^{n_i} \left(y_{ij}-\bar{y}\right) \left(y_{ij}-\bar{y}\right)^\mathrm{T}     \\
B &= \sum_{i=1}^p n_i \left(\bar{y}_{i}-\bar{y}\right) \left(\bar{y}_{i}-\bar{y}\right)^\mathrm{T}        \\
W &= \sum_{i=1}^p \sum_{j=1}^{n_i} \left(y_{ij}-\bar{y}_i\right) \left(y_{ij}-\bar{y}_i\right)^\mathrm{T}
\end{split}
\end{equation}

und die darauf basierende *Wilks'-Lambda-Statistik*

\begin{equation} \label{eq:Lambda}
\Lambda = \frac{|W|}{|B+W|} \; ,
\end{equation}

wobei $B$ die *between sum-of-squares Matrix*, $W$ die *within sum-of-squares Matrix* und $T = B + W$ die *total sum of-squares Matrix* ist.

Das Theorem zu speziellen Verteilungen von Wilks'-$\Lambda$-Transformationen gibt Spezialfälle von Datendimension $m$ und Gruppenanzahl $p$ an, in denen die Verteilung einer Transformation von $\Lambda$ unter der Nullhypothese $H_0 : \mu_1 = \ldots = \mu_p$ exakt bekannt ist. Es besagt, dass für den in den ersten beiden Tabellenspalten aufgeführten Spezialfall die in der dritten Tabellenspalte genannte Teststatistik bei Vorliegen der Nullhypothese einer $f$-Verteilung mit den Freiheitsgradparametern in der vierten Tabellenspalte folgt:

\begin{center}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{cccc}
Datendimension $m$  & Gruppenanzahl $p$  & Statistik                                   & $f$-Verteilungsparameter \\ \hline
beliebig  & 2         & $\frac{1 - \Lambda}{\Lambda} \frac{n-p-m+1}{m}$                & $m, \; n-p-m+1$          \\
beliebig  & 3         & $\frac{1 - \sqrt{\Lambda}}{\sqrt{\Lambda}} \frac{n-p-m+1}{m}$  & $2m, \; 2(n-p-m+1)$      \\
1         & beliebig  & $\frac{1 - \Lambda}{\Lambda} \frac{n-p}{p-1}$                  & $p-1, \; n-p$            \\
2         & beliebig  & $\frac{1 - \sqrt{\Lambda}}{\sqrt{\Lambda}} \frac{n-p-1}{p-1}$  & $2(p-1), \; 2(n-p-1)$
\end{tabular}
\end{center}

Damit gilt z.B. für $m>2$ und $p=3$ (zweite Tabellenzeile): Überschreitet $\frac{1 - \sqrt{\Lambda}}{\sqrt{\Lambda}} \frac{n-p-m+1}{m}$ den kritischen Wert $k_{\alpha_0}$, der mit einem Signifikanzniveau $\alpha_0$ aus der inversen kumulativen Verteilungsfunktion der $f$-Verteilung mit Freiheitsgradparametern $2m$ und $2(n-p-m+1)$ errechnet werden kann, wird die Nullhypothese identischer Erwartungswerte über Gruppen abgelehnt. Anderenfalls wird sie nicht abgelehnt.


## Analyse in R

Die Datei `FADE_SAME.csv` enthält den in der ersten Seminarsitzung vorgestellten Datensatz. Erklären Sie die Funktion des folgenden R-Codes:

\footnotesize
```{r, echo = T, eval = T, results = 'hide'}
# Daten einlesen
fname = 'FADE_SAME.csv'                             # Dateiname
D     = read.csv(fname)                             # Dataframe

# Datenmatrix extrahieren
rows  = startsWith(D$subject, 'subA')               # Personen aus Studie A
cols  = c('novelty.FADE','novelty.SAME',            # Definition Variablen
          'memory.FADE', 'memory.SAME')
n     = sum(rows)                                   # Anzahl Datenpunkte
m     = length(cols)                                # Anzahl Variablen
Y     = t(as.matrix(D[rows,cols]))                  # m x n Datenmatrix
print(dim(Y))                                       # Überprüfung Datenmatrix
```
\normalsize


## Erste Programmieraufgabe

Wir wollen systematische Unterschiede der FADE-Scores und SAME-Scores zwischen den drei Altersgruppen "junge Personen" (bis 35 Jahre), "mittelalte Personen" (50-59 Jahre) und "ältere Personen" (ab 60 Jahre) untersuchen.

Wir werden dazu zunächst eine Variable $x$ definieren, die die unabhängige Variable "Altersgruppe" repräsentiert und die Zugehörigkeit der Datenpunkte zu den verschiedenen Studiengruppen kodiert. Gehen Sie dazu wie folgt vor:

- Extrahieren Sie die Variable "age" aus dem Dataframe, aber nur für diejenigen Datenpunkte, die durch die Variable \verb|rows| beschrieben werden.

- Bilden Sie drei logische Variablen \verb|g1|, \verb|g2|, \verb|g3|, d.h. Vektoren von derselben Länge wie die Datenmatrix, die jeweils \verb|TRUE| sind, wenn der entsprechende Datenpunkt zu Gruppe 1, 2, 3 gehört. Beispielsweise ist \verb|g1 = (x <= 35)|.

- Bilden Sie einen Vektor, der für jeden Datenpunkt den entsprechenden Gruppenindex enthält, d.h. einen Vektor, der 1 ist, wenn Altersgruppe \verb|g1| vorliegt, 2 ist, wenn Altersgruppe \verb|g2| vorliegt und 3 ist, wenn Altersgruppe \verb|g1| vorliegt. Sie erreichen dies z.B., indem Sie die logischen Vektoren mit den Gruppenindizes multiplizieren und aufsummieren. Speichern Sie das Ergebnis in die Variable \verb|x|.

- Speichern Sie das Maximum von \verb|x|, d.h. die Anzahl der Gruppen, in die Variable \verb|p|.

- Geben Sie die ersten 100 Einträge des Vektors \verb|x| aus. Sie sollten folgende Ergebnisse erhalten:

\footnotesize
```{r, echo = F, eval = T}
# Gruppen definieren
x     = D$age[rows]                                 # gruppendefinierende Variable
g1    = (x <= 35)                                   # Gruppe 1: Alter <= 35
g2    = (x >  50) & (x < 60)                        # Gruppe 2: 50 < Alter < 60
g3    = (x >= 60)                                   # Gruppe 3: Alter >= 60
x     = g1*1 + g2*2 + g3*3                          # diskrete Gruppenvariable
p     = max(x)                                      # Anzahl Gruppen
print(x[1:100])                                     # Überprüfung UV
```
\normalsize


## Analyse in R

Im zweiten Teil wollen wir nun mit der Datenmatrix $Y$ und dem Gruppenvektor $x$ eine einfaktorielle multivariate Varianzanalyse durchführen. Erklären Sie dazu zunächst den folgenden R-Code und die Funktionen, die er definiert:

\footnotesize
```{r, echo = T, eval = T, results = 'hide'}
# Parameter schätzen
estimate = function(Y,x) {

    # Datendimensionalität
    n = ncol(Y)                                     # Anzahl Datenpunkte
    m = nrow(Y)                                     # Anzahl Variablen
    p = max(x)                                      # Anzahl Gruppen

    # Parameterschätzer
    mu_hat    = matrix(rep(0,m*p), nrow = m)        # Erwartungswertparameter
    Sigma_hat = matrix(rep(0,m*m), nrow = m)        # Kovarianzmatrixparameter
    for (i in 1:p) {
        Y_i        = Y[,x==i]                       # Datenmatrix Gruppe i
        n_i        = ncol(Y_i)                      # Anzahl Datenpunkte Gruppe i
        mu_hat[,i] = as.matrix(rowMeans(Y_i))       # Erwartungswerte Gruppe i
        for (j in 1:n_i) {
            Sigma_hat = Sigma_hat + (Y_i[,j] - mu_hat[i]) %*% t(Y_i[,j] - mu_hat[i])
        }
    }
    Sigma_hat = (1/(n-p)) * Sigma_hat

    # Funktionswerte
    return(list(mu_hat = mu_hat, Sigma_hat = Sigma_hat))
}

# Multivariate Quadratsummenzerlegung
sumofsqr = function(Y,x) {

    # Datendimensionalität
    n = ncol(Y)                                     # Anzahl Datenpunkte
    m = nrow(Y)                                     # Anzahl Variablen
    p = max(x)                                      # Anzahl Gruppen

    # sum-of-squares Matrizen
    y_bar   = as.matrix(rowMeans(Y))                # Gesamtstichprobenmittel
    y_i_bar = matrix(rep(0,m*p), nrow = m)          # Gruppenstrichprobenmittel
    B       = matrix(rep(0,m*m), nrow = m)          # between-group sum-of-squares Matrix
    W       = matrix(rep(0,m*m), nrow = m)          # within-group sum-of-squares Matrix
    for (i in 1:p) {
        Y_i         = Y[,x==i]                      # Datenmatrix Gruppe i
        n_i         = ncol(Y_i)                     # Anzahl Datenpunkte Gruppe i
        y_i_bar[,i] = as.matrix(rowMeans(Y_i))      # Erwartungswerte Gruppe i
        for (j in 1:n_i) {
            B = B + (y_i_bar[,i] - y_bar) %*% t(y_i_bar[,i] - y_bar)
            W = W + (Y_i[,j] - y_i_bar[,i]) %*% t(Y_i[,j] - y_i_bar[,i])
        }
    }
    T         = B + W                               # total sum-of-squares Matrix

    # Funktionswerte
    return(list(T = T, B = B, W = W))
}
```
\normalsize


## Zweite Programmieraufgabe

Nutzen Sie die Funktionen der Parameterschätzung und der multivariaten Quadratsummenzerlegung, um die MANOVA-Modellparameter zu schätzen. Gehen Sie dazu wie folgt vor:

- Wenden Sie die Funktion \verb|estimate| auf $Y$ und $x$ an und extrahieren Sie die Felder \verb|mu_hat| sowie \verb|Sigma_hat|.

- Wenden Sie die Funktion \verb|sumofsqr| auf $Y$ und $x$ an und extrahieren Sie die Felder \verb|T| sowie \verb|B| und \verb|W|.

- Geben Sie die Schätzer für die Gruppenerwartungswertparameter und den Kovarianzmatrixparameter aus. Sie sollten folgende Ergebnisse erhalten:

\footnotesize
```{r, echo = F, eval = T}
# Modellschätzung und Modellevaluation
mu_hat    = estimate(Y,x)$mu_hat                    # Erwartungswertparameterschätzer
Sigma_hat = estimate(Y,x)$Sigma_hat                 # Kovarianzmatrixparameterschätzer
T         = sumofsqr(Y,x)$T                         # total sum-of-squares Matrix
B         = sumofsqr(Y,x)$B                         # between sum-of-squares Matrix
W         = sumofsqr(Y,x)$W                         # within sum-of-squares Matrix
colnames(mu_hat) = c("young","middle.aged","older") # Spaltennamen \hat{\mu}
print(mu_hat)
print(Sigma_hat)
```
\normalsize

\pagebreak
Führen Sie nun eine einfaktorielle multivariate Varianzanalyse (MANOVA) durch. Gehen Sie dazu wie folgt vor:

- Legen Sie das Signifikanzniveau $\alpha_0 = 0.05$ fest.

- Berechnen Sie die Wilks'-Lambda-Statistik $\Lambda$ anhand der obigen Formel. Die Determinante ist in R als die Funktion \verb|det()| implementiert.

- Suchen Sie aus der Tabelle in Abschnitt 6.1 den für diese Datenanalyse passenden Spezialfall heraus, berechnen Sie aus $\Lambda$ die Teststatistik und speichern Sie das Ergebnis in die Variable \verb|eff|.

- Berechnen Sie die zum Spezialfall gehörenden Freiheitsgradparameter und speichern Sie diese in die Variablen \verb|df_1| und \verb|df_2|.

- Ermitteln Sie den kritischen Wert $k_{\alpha_0}$ des Tests. Die inverse kumulative Verteilungsfunktion der F-Verteilung ist in R als \verb|qf()| implementiert. Informieren Sie sich über die Eingabeparameter dieser Funktion, um den kritischen Wert korrekt zu bestimmen.

- Ermitteln Sie den p-Wert für den Test. Die kumulative Verteilungsfunktion der F-Verteilung ist in R als \verb|pf()| implementiert. Informieren Sie sich über die Eingabeparameter dieser Funktion, um den p-Wert korrekt zu bestimmen.

- Ermitteln Sie das Testergebnis. Es ist 1, wenn die Teststatistik den kritischen Wert überschrittet. Andernfalls ist es 0.

- Geben Sie die Resultate ihrer Analyse aus. Sie sollten folgende Ergebnisse erhalten:

\footnotesize
```{r, echo = F}
# Einfaktorielle multivariate Varianzanalyse (manuell)
alpha_0   = 0.05                                    # Signifikanzlevel
L         = det(W)/det(W + B)                       # Wilks' Lambda
eff       = (1-sqrt(L))/sqrt(L) * ((n-p-m+1)/m)     # Teststatistik
df_1      = 2*m                                     # Freiheitsgradparameter 1
df_2      = 2*(n-p-m+1)                             # Freiheitsgradparameter 2
k_alpha_0 = qf(1-alpha_0, df_1, df_2)               # kritischer Wert
phi       = as.numeric(eff > k_alpha_0)             # Nullhypothesentest
p_val     = 1-pf(eff, df_1, df_2)                   # p-Wert
cat(  "Wilks' Lambda          : ", L,
    "\nTeststatistik          : ", eff,
    "\nFreiheitsgradparameter : ", c(df_1, df_2),
    "\nSignifikanzniveau      : ", alpha_0,
    "\nkritischer Wert        : ", k_alpha_0,
    "\nTestergebnis           : ", phi,
    "\np-Wert                 : ", p_val, "\n")
```
\normalsize

- Überprüfen Sie Ihre multivariate Varianzanalyse, indem Sie folgenden Code mit der in R eingebauten Funktion \verb|Manova()| ausführen:

\footnotesize
```{r, eval = F, warning = F, results = 'hide'}
# Einfaktorielle multivariate Varianzanalyse (automatisch)
library(car)                                        # R-Paket "car"
group  = as.factor(x)                               # faktorielle Variable
model  = lm(t(Y) ~ group)                           # Modellformulierung
result = Manova(model, test.statistic = "Wilks")    # Modellevaluation
print(result)                                       # Ergebnisausgabe
```
\normalsize


\pagebreak
## Lückentext

Füllen Sie mit den in der Übung gewonnenen Erkenntnissen den folgenden Lückentext aus und präsentieren Sie die Ergebnisse im Seminar:

\vspace{1em}
**Lückentext:** Die Anzahl der beobachteten experimentellen Einheiten ist \underline{ 259 }, die Anzahl der abhängigen Variablen pro Einheit ist \underline{ 4 } und die Anzahl der Level der unabhängigen Variable ist \underline{ 3 }. Die Wilks'-Lambda-Statistik hat den Wert \underline{ 0.596 }, die daraus sich ergebende Teststatistik ist \underline{ 18.67 }. Der kritische Wert für die multivariate ANOVA ergibt sich aus \underline{ der inversen kumulativen Verteilungsfunktion der F-Verteilung mit den Freiheitsgradparametern } \underline{ 8 und 506 } und beträgt \underline{ 1.96 }. Die Nullhypothese identischer Gruppenerwartungswerte wird daher bei einem Signifikanznivea von 5 Prozent \underline{ abgelehnt } (abgelehnt/nicht abgelehnt). Der p-Wert für die multivariate ANOVA ergibt sich aus \underline{ der kumulativen Verteilungsfunktion der } \underline{ F-Verteilung mit den Freiheitsgradparametern 8 und 506 } und beträgt \underline{ praktisch 0 ($< 2.2 \cdot 10^{-16}$) }.


## Mögliche Klausurfrage

Präsentieren Sie im Seminar folgende Klausurfrage und erklären Sie die richtige Antwort:

\vspace{1em}
**Frage:** Welche Aussage über die Parameterschätzer $\hat{\mu}_i$ und $\hat{\Sigma}$ im Modell der einfaktoriellen multivariaten Varianzanalyse ist **nicht** korrekt?
\begin{enumerate}[a)]
\item $\hat{\mu}_i$ ist das multivariate Stichprobenmittel der $i$-ten Gruppe.
\item $\hat{\mu}_i$ entspricht dem Gesamtstichprobenmittel $\bar{y}$. $\checkmark$
\item $\hat{\Sigma}$ entspricht einer skalierten Form der der within-group sum-of-squares Matrix $W$.
\item $\hat{\Sigma}$ ist ein unverzerrter Schätzer des Kovarianzmatrixparameters $\Sigma$.
\end{enumerate}


## Kinderwitz

Gestern haben wir zwei Biber beim Essen beobachtet:

\begin{turn}{180}
Gab Steg.
\end{turn}