---
fontsize: 8pt
bibliography: 11_Referenzen.bib
citation_package: natbib
output:
  beamer_presentation:
    keep_tex: true
    includes:
      in_header: 11_header.tex
---


```{r, include = F}
source("11_R_common.R")
```

#  {.plain}
\center
```{r, echo = FALSE, out.width = "20%"}
knitr::include_graphics("11_Abbildungen/mvda_11_otto.png")
```

\vspace{2mm}

\Huge
Multivariate Datenanalyse
\vspace{6mm}

\Large
MSc Psychologie WiSe 2021/22

\vspace{6mm}
\large
Prof. Dr. Dirk Ostwald

# {.plain}
\vfill
\center
\huge
\textcolor{black}{(11) Einfaktorielle Varianzanalyse}
\vfill

#
\setstretch{2}
\vfill
\large
Anwendungsszenario und Datendeskription

Modellformulierung und Modellschätzung

Modellevaluation

* Wilks' $\Lambda$
* Pillai's $V$  

Selbstkontrollfragen
\vfill

#
\setstretch{2}
\vfill
\large
**Anwendungsszenario und Datendeskription**

Modellformulierung und Modellschätzung

Modellevaluation mit

* Wilks' $\Lambda$
* Pillai's $V$  

Selbstkontrollfragen
\vfill


# Anwendungsszenario und Datendeskription
\setstretch{1.2}
Anwendungsszenario

\small
* Zwei oder mehr Stichproben (oft Gruppen genannt) experimenteller Einheiten.

* Annahme der unabhängigen und identischen Normalverteilung $N(\mu_i,\Sigma)$ der Daten.

* $\mu_i \in \mathbb{R}^m, i = 1,...,k$ und $\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d. }$ unbekannt.

* Absicht des inferentiellen Testens der Nullhypothese identischer Gruppenerwartungswerte.

\center
Generalisierung des Zweistichproben-T$^2$-Tests bei unabhängigen Stichproben mit
einfacher Nullhypothese für mehr als zwei Stichproben

\normalsize
\flushleft
Anwendungsbeispiele

\small
* BDI/Glukokortikoid Analyse von drei Gruppen psychiatrischer Patient:innen nach PA, CBT,ST

$\quad\quad\quad\circ$ Inferentielle Evidenz für Gruppenerwartungswertunterschiede?

$\quad\quad\quad\circ$ Evidenz für unterschiedliche Therapiewirksamkeit?

* Gruppenvergleich von Testdaten bei gutem, befriedigendem, ungenügenden Studienabschluss

$\quad\quad\quad\circ$ Inferentielle Evidenz für Gruppenerwartungswertunterschiede?

$\quad\quad\quad\circ$ Evidenz für Studienerfolgsprädiktivität der Testdaten?


# Anwendungsszenario und Datendeskription
\textcolor{darkblue}{Anwendungsbeispiel}

\small
Nach @rudolf_2020 Kapitel 4, vgl. Einheiten zur Prädiktiven Modellierung

Datensatz zum Verhältnis psychologischer Testdiagnostik und Studienerfolg

Faktor Studienerfolg mit $k = 3$ Leveln (*Ungenügend, Befriedigend, Gut*)

Datendimension $m = 2$ mit $y_{ij} \in \mathbb{R}^2$ ($y_{ij_1}$ *Intelligenztestscore*, $y_{ij_2}$ *Mathematiktestscore*)
\vspace{4mm}

\setstretch{1}
\tiny
```{r, message = F}
# Einlesen des Datensatzes aus Rudolf & Buse (2019) Multivariate Verfahren
library("foreign")
S        = read.spss(file.path(getwd(), "11_Daten", "studienerfolg.sav"), to.data.frame = T)

# Datenpräprozessierung
m        = 2                                       # Datendimension von Interesse
k        = 3                                       # Anzahl Gruppen
l        = 15                                      # Anzahl Datenpunkte pro Gruppe
Y        = array(dim = c(m,l,k))                   # Datenarrayinitialisierung
Y[,,1]   = rbind(S$X1[S$Gruppe == "ungenügend"],   # Y_{1j_1} Intelligenztestscore
                 S$X2[S$Gruppe == "ungenügend"])   # Y_{1j_2} Mathematiktestscore
Y[,,2]   = rbind(S$X1[S$Gruppe == "befriedigend"], # Y_{2j_1} Intelligenztestscore
                 S$X2[S$Gruppe == "befriedigend"]) # Y_{2j_2} Mathematiktestscore
Y[,,3]   = rbind(S$X1[S$Gruppe == "gut"],          # Y_{3j_1} Intelligenztestscore
                 S$X2[S$Gruppe == "gut"])          # Y_{3j_2} Mathematiktestscore
```


# Anwendungsszenario und Datendeskription
\textcolor{darkblue}{Anwendungsbeispiel}

\small
Datendeskription
\vspace{1mm}
\setstretch{.9}
\tiny
```{r}
Y_bar_i = array(dim = c(m,k))                                            # Stichprobenmittelarray
C_i     = array(dim = c(m,m,k))                                          # Stichprobenkovarianzmatrizenarray
j_l     = matrix(rep(1,l), nrow = l)                                     # 1_{l}
I_l     = diag(l)                                                        # Einheitsmatrix I_l
J_l     = matrix(rep(1,l^2), nrow = l)                                   # 1_{ll}
for (i in 1:k){                                                          # Gruppeniterationen
    Y_bar_i[,i] = (1/l)*(Y[,,i] %*% j_l)                                 # Stichprobenmittel \bar{Y}_i
    C_i[,,i]    = (1/(l-1))*(Y[,,i] %*% (I_l-(1/l)*J_l) %*% t(Y[,,i]))}  # Stichprobenkovarianzmatrix C_i
```


```{r, echo = F, eval = F}
library(ellipse)
graphics.off()
dev.new()
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.2,
cex.main    = 1.2)
cols        = c("Gray20", "Gray50", "Gray80")

plot(
NaN,
NaN,
xlim = c(30,80),
ylim = c(0,100),
xlab = TeX("$Y_{ij_1}$ Intelligenztestscore"),
ylab = TeX("$Y_{ij_2}$ Mathematiktestscore"))
for (i in 1:k){
    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 = 0.5, centre = Y_bar_i[,i])
    lines(
    iso[,1],
    iso[,2],
    col = cols[i])
}
legend(
"topleft",
c("Ungenügend", "Befriedigend", "Gut"),
pch         = c(21,21,21),
bg          = cols,
col         = cols,
bty         = "n",
cex         = 1,
x.intersp   = 1)
fdir        =  file.path(getwd(), "11_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_11_anwendungsbeispiel.pdf"),
width       = 6,
height      = 6)
```


```{r, echo = FALSE, out.width = "45%"}
knitr::include_graphics("11_Abbildungen/mvda_11_anwendungsbeispiel.pdf")
```


#
\setstretch{2}
\vfill
\large
Anwendungsszenario und Datendeskription

**Modellformulierung und Modellschätzung**

Modellevaluation

* Wilks' $\Lambda$
* Pillai's $V$  

\large
Selbstkontrollfragen
\vfill


# Modellformulierung und Modellschätzung
\setstretch{1.2}
\footnotesize
\begin{definition}[Modell der einfaktoriellen Varianzanalyse]
\justifying
Für $i = 1,...,k$ sei
\begin{equation}\label{eq:manova_1}
Y_{i1},...,Y_{il} \sim N(\mu_i,\Sigma)
\end{equation}
eine Stichprobe eines multivariaten Normalverteilungsmodells von Größe $l$ mit
unbekanntem Erwartungswertparameter $\mu_i \in \mathbb{R}^m$ und unbekannten
Kovarianzmatrixparameter $\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d.}$ Dann
heißt \eqref{eq:manova_1} \textit{Klassisches Modell der einfaktoriellen
Varianzanalyse}. Die Gesamtstichprobengröße
bezeichnen wir mit $n := kl$. Äquivalent sei
\begin{equation}\label{eq:manova_2}
Y_{ij}  = \mu_i + \varepsilon_{ij}
\mbox{ mit }
\varepsilon_{ij} \sim N(0_m,\Sigma)
\mbox{ für }
i = 1,...,k \mbox{ und } j = 1,...,l,
\end{equation}
wobei
\begin{itemize}
\itemsep-0.5mm
\item $i$ die Stichproben und $j$ die experimentellen Einheiten indizieren,
\item $l$ die Stichprobengrößen und $n := lk$ die Gesamtstichprobengröße sind,
\item $Y_{ij}$ beobachtbare Zufallsvektoren sind,
\item $\mu_i \in \mathbb{R}^m$ feste Erwartungswertparameter der Stichprobenvariablen sind, und
\item $\varepsilon_{ij}$ unabhängige normalverteilte nicht-beobachtbare Zufallsvariablen mit $\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d. }$ sind.
\end{itemize}
Dann heißt \eqref{eq:manova_2} \textit{Generatives Modell der einfaktoriellen
Varianzanalyse.}
\end{definition}
\vspace{-2mm}
Bemerkungen

* Der Einfachheit halber setzen wir hier identische Stichprobengrößen voraus.
* Die Kovarianzmatrixparameter aller Stichproben werden als identisch vorausgesetzt.
* Wir verzichten auf einen Beweis der Äquivalenz von \eqref{eq:manova_1} und \eqref{eq:manova_2}.
* Die Gesamtheit aller Stichprobenzufallsvektoren bezeichen wir mit $Y := \{Y_{ij}\}_{i = 1,...,k, j = 1,...,l}$.

# Modellformulierung und Modellschätzung
\setstretch{1.2}
\footnotesize
\begin{theorem}[Parameterschätzer der einfaktoriellen Varianzanalyse]
\normalfont
\justifying
Für $i = 1,...,k$ sei mit
\begin{equation}
Y_{i1},...,Y_{il} \sim N(\mu_i,\Sigma)
\end{equation}
für $\mu_i \in \mathbb{R}^m$ und $\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d.}$
das Modell der einfaktoriellen Varianzanalyse gegeben. Dann ist für $i = 1,...,k$
\begin{equation}
\hat{\mu}_i := \frac{1}{l}\sum_{j = 1}^l Y_{ij}
\end{equation}
ein unverzerrte Schätzer des Erwartungswertparameters $\mu_i$ und
\begin{equation}
\hat{\Sigma} := \frac{1}{lk - k}\sum_{i = 1}^k \sum_{j = 1}^l \left(Y_{ij} - \hat{\mu}_i\right)\left(Y_{ij} - \hat{\mu}_i\right)^T
\end{equation}
ein unverzerrter Schätzer des Kovarianzmatrixparameters $\Sigma$.
\end{theorem}

Bemerkungen

* $\hat{\mu}_i$ ist das Stichprobenmittel der $i$ten Gruppe.
* $\hat{\Sigma}$ ist die mit $(lk - k)$ skalierte Within Group Sum of Squares Matrix (siehe unten).
* Anstelle eines Beweis validieren wir die Aussage des Theorems mithilfe einer Simulation.

# Modellformulierung und Modellschätzung

Parameterschätzer der einfaktoriellen Varianzanalyse
\vspace{2mm}

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

  # Diese Funktion evaluiert die Parameterschätzer einer einfaktoriellen
  # Varianzanalyse basierend auf einen m x l x k Datensatz Y.
  #
  # Input
  #     Y          : m x l x k Datenarray
  #
  # Output
  #     $mu_hat    : m x k \mu_i Parameterschätzer
  #     $Sigma_hat : m x m \Sigma Parametschätzer
  # ----------------------------------------------------------------------------
  # Dimensionsparameter
  d         = dim(Y)                                    # Datensatzdimensionen
  m         = d[1]                                      # Datendimension
  l         = d[2]                                      # Anzahl Datenpunkte pro Gruppe
  k         = 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:k){
      for(j in 1:l){
          Sigma_hat = Sigma_hat + (1/(l*k-k))*(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))}


```

# Modellformulierung und Modellschätzung
Parameterschätzer der einfaktoriellen Varianzanalyse
\vspace{2mm}

\tiny
\setstretch{0.8}
```{r, eval = F}
# R Pakete
library(MASS)                                    # multivariate Normalverteilungen
library(matlib)                                  # Matrixalgebra


# Modellparameter
k          = 3                                   # Anzahl Gruppen
l          = 15                                  # Anzahl Datenpunkte pro Gruppe
m          = 2                                   # Datendimensionalität
mu_i       = matrix(c(1,2,2,1,3,2.5), ncol = k)  # Erwartungswertparameter
Sigma      = matrix(c(1,.5,.5,1)  , ncol = m)    # Kovarianzmatrixparameter

# Simulationsparameter und Arrays
nsm        = 1e2                                 # Anzahl Simulation
mu_hat_is  = array(dim = c(m,k,nsm))             # \hat{\mu}_i Array
Sigma_hats = array(dim = c(m,m,nsm))             # \hat{\Sigma} Array

# Simulationen
for(s in 1:nsm){

    # Datengeneration
    Y               = array(dim = c(m,l,k))      # Datenarray
    for(i in 1:k){
        Y[,,i] = t(mvrnorm(l,mu_i[,i],Sigma))    # Datengeneration
    }
    S               = estimate(Y)                # Parameterschätzung
    mu_hat_is[,,s]  = S$mu_hat_i                 # \hat{\mu}_i
    Sigma_hats[,,s] = S$Sigma_hat                # \hat{\Sigma}
}

# Erwartungswertschätzung
E_hat_mu_i_hat  = apply(mu_hat_is , c(1,2), mean)
E_hat_Sigma_hat = apply(Sigma_hats, c(1,2), mean)
```

# Modellformulierung und Modellschätzung
Parameterschätzer der einfaktoriellen Varianzanalyse
\vspace{2mm}

```{r, echo = F, eval = F}
library(ellipse)
graphics.off()
dev.new()
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.2,
cex.main    = 1.2)
cols        = c("Black", "Gray70")

plot(
NaN,
NaN,
xlim = c(0,4),
ylim = c(0,4),
xlab = TeX("$Y_{ij_1}$"),
ylab = TeX("$Y_{ij_2}$"),
main = "100 Simulationen")
for (i in 1:k){
    points(
    mu_i[1,i],
    mu_i[2,i],
    col = cols[1],
    pch = 19)
    points(
    E_hat_mu_i_hat[1,i],
    E_hat_mu_i_hat[2,i],
    col = cols[2],
    pch = 1)
    iso_1 = ellipse(Sigma, level = 0.2, centre = mu_i[,i])
    lines(
    iso_1[,1],
    iso_1[,2],
    col   = cols[1])
    iso_2 = ellipse(E_hat_Sigma_hat, level = 0.2, centre = E_hat_mu_i_hat[,i])
    lines(
    iso_2[,1],
    iso_2[,2],
    col   = cols[2])
}
legend(
"topleft",
c("Wahre, aber unbekannte, Parameter", "Geschätzte Parameterschätzererwartungswerte"),
lty         = 1,
bg          = cols,
col         = cols,
bty         = "n",
cex         = .7,
x.intersp   = 1)
fdir        =  file.path(getwd(), "11_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_11_parameterschätzer.pdf"),
width       = 6,
height      = 6)
```


```{r, echo = FALSE, out.width = "60%"}
knitr::include_graphics("11_Abbildungen/mvda_11_parameterschätzer.pdf")
```




# Modellformulierung und Modellschätzung
\textcolor{darkblue}{Anwendungsbeispiel}

\small
Parameterschätzung

\footnotesize
```{r}
# Parameterschätzung
S = estimate(Y)

# Ausgabe
print(S$mu_hat_i)
print(S$Sigma_hat)
```

#
\setstretch{2}
\vfill
\large
Anwendungsszenario und Datendeskription

Modellformulierung und Modellschätzung

**Modellevaluation**

* Wilks' $\Lambda$
* Pillai's $V$  

Selbstkontrollfragen
\vfill


# Modellevaluation
Überblick

\footnotesize
Primäres Ziel der einfaktoriellen Varianzanalyse ist zumeist das Testen der Nullhypothese
\begin{equation}
H_0 : \mu_1 = \mu_2 = \cdots = \mu_k.
\end{equation}

Die Nullhypothese gesagt also, dass zwischen den Gruppen keine 
Erwartungswertparameterunterschiede bestehen. Die Alternativhypothese lautet somit
\begin{equation}
H_1 : \mu_{i_c} \neq \mu_{j_c}
\mbox{ für mindestens ein Paar } (i,j)
\mbox{ mit } i \neq j
\mbox{ und mindestens ein } c
\mbox{ mit } 1 \le c \le m.
\end{equation}

Die Alternativhypothese besagt also, dass sich mindestens zwei Erwartungswertparameter
in mindestens einer ihrer Komponenten unterscheiden. 

Kritische Wert-basierte Tests der Nullhypothesen können mit verschiedenen
Teststatistiken (\textit{Wilks' $\Lambda$, Pillai Statisik}) konstruiert werden.
Diesen Teststatistiken ist gemein, dass sie auf eine Generalisierung der vom univariaten
Fall bekannten  Quadratsummenzerlegung zurück gehen. Wir führen also zunächst
diese sogenannte *Kreuzproduktsummenmatrizenzerlegung* ein und betrachten dann
die durch die Wilks' $\Lambda$ und die Bartlett-Pillai-Spur induzierten Tests
anhand der Gliederung (1) Teststatistik und Test, (2) Analyse der Teststatistik
und (3) Testumfangkontrolle.

Einen Überblick über die Modellevaluation bei der univariaten einfaktoriellen Varianzanalyse
gibt [\textcolor{red}{(12) Varianzanalysen}](https://bit.ly/3rdQsUE).

# Modellevaluation

\footnotesize
\begin{theorem}[Kreuzproduktsummenmatrizenzerlegung]
\justifying
\normalfont
Für $i = 1,...,k$ und $j = 1,...,l$ bezeichne $Y_{ij}$ den $j$ten Stichprobenvektor
der $i$ten Stichprobengruppe eines einfaktoriellen Varianzanalysemodells. Weiterhin
seien
\begin{equation}
\bar{Y} := \frac{1}{n}\sum_{i = 1}^k \sum_{j = 1}^l Y_{ij}
\mbox{ und }
\bar{Y}_i := \frac{1}{l} \sum_{j = 1}^l Y_{ij}
\end{equation}
das \textit{Gesamtstichprobenmittel} und das \textit{$i$te Gruppenstichprobenmittel}, respektive.
Schließlich seien
\vspace{1mm}

\center
\begin{tabular}{ll}
$T := \sum_{i = 1}^k \sum_{j = 1}^l \left(Y_{ij} - \bar{Y}\right)\left(Y_{ij} - \bar{Y}\right)^T$
&
die \textit{Totale Sum of Squares Matrix}
\\
\\
$B := \sum_{i = 1}^k l \left(\bar{Y}_{i} - \bar{Y}\right)\left(\bar{Y}_{i} - \bar{Y}\right)^T$
&
die \textit{Between Group Sum of Squares Matrix}
\\
\\
$W := \sum_{i = 1}^k \sum_{j = 1}^l \left(Y_{ij} - \bar{Y}_i\right)\left(Y_{ij} - \bar{Y}_i\right)^T$
&
die \textit{Within Group Sum of Squares Matrix}.
\end{tabular}
\vspace{1mm}

\flushleft
Dann gilt
\begin{equation}
T = B + W.
\end{equation}
\end{theorem}

# Modellevaluation
\footnotesize
\setstretch{1.8}
Bemerkungen

* \justifying $T \in \mathbb{R}^{m \times m}$ repräsentiert die totale Variabilität der Datenvektoren um das Gesamtstichprobenmittel.

* $B \in \mathbb{R}^{m \times m}$ repräsentiert die Variabilität der Gruppenstichprobenmittel um das Gesamtstichprobenmittel.

* $W \in \mathbb{R}^{m \times m}$ repräsentiert die Variabilität der Datenvektoren um ihre jeweiligen Gruppenstichprobenmittel.

* Die totale Variabilität wird hier also in zwei unabhängige Beiträge von Variabilität zerlegt.

* $W$ heißt auch \textit{Residualvariabiliät}, weil sie die verbleibende Variabilität nach
Schätzung der Gruppenerwartungswertparameter quantifiziert und gilt das
\begin{equation}
W = (lk -k)\hat{\Sigma}.
\end{equation}

* Die Kreuzproduktsummenmatrixzerlegung ergibt sich anhand von algebraischen Identitäten wie unten gezeigt.

#  Modellevaluation
\vspace{2mm}
\tiny
\underline{Beweis}
\vspace{-2mm}
\begin{align*}
\begin{split}
T
& = \sum_{i=1}^k \sum_{j=1}^{l}(Y_{ij}-\bar{Y})(Y_{ij}-\bar{Y})^T
\\
& = \sum_{i=1}^k \sum_{j=1}^{l} (Y_{ij}-\bar{Y}_i+\bar{Y}_i-\bar{Y})(Y_{ij}-\bar{Y}_i+\bar{Y}_i-\bar{Y})^T
\\
& = \sum_{i=1}^k \sum_{j=1}^{l} \left((Y_{ij}-\bar{Y}_i)+(\bar{Y}_i-\bar{Y})\right)\left((Y_{ij}-\bar{Y}_i)+(\bar{Y}_i-\bar{Y})\right)^T
\\
& = \sum_{i=1}^k \sum_{j=1}^{l} \left(
                                      (Y_{ij}-\bar{Y}_i)(Y_{ij}-\bar{Y}_i)^T
								    +2(Y_{ij}-\bar{Y}_i)(\bar{Y}_i-\bar{Y})^T
				    		        + (\bar{Y}_i-\bar{Y})(\bar{Y}_i-\bar{Y})^T
				  		        \right)
\\
& = \sum_{i=1}^k \left(
                       \sum_{j=1}^{l} (Y_{ij}-\bar{Y}_i)(Y_{ij}-\bar{Y}_i)^T
   				      +\sum_{j=1}^{l}2(Y_{ij}-\bar{Y}_i)(\bar{Y}_i-\bar{Y})^T
				      +\sum_{j=1}^{l}(\bar{Y}_i-\bar{Y})(\bar{Y}_i-\bar{Y})^T
				 \right)
\end{split}
\end{align*}


#  Modellevaluation
\vspace{2mm}
\tiny
\underline{Beweis (fortgeführt)}
\vspace{-2mm}
\begin{align*}
\begin{split}
& = \sum_{i=1}^k \left(
                       \sum_{j=1}^{l} (Y_{ij}-\bar{Y}_i)(Y_{ij}-\bar{Y}_i)^T
   				      +2\left(\sum_{j=1}^{l}(Y_{ij}-\bar{Y}_i)\right)(\bar{Y}_i-\bar{Y})^T
				      +l(\bar{Y}_i-\bar{Y})(\bar{Y}_i-\bar{Y})^T
				 \right)
\\
& = \sum_{i=1}^k \left(
                       \sum_{j=1}^{l} (Y_{ij}-\bar{Y}_i)(Y_{ij}-\bar{Y}_i)^T
   				      +2\left(\sum_{j=1}^{l}\left(Y_{ij}-\frac{1}{l}\sum_{j=1}^l Y_{ij}\right)\right)(\bar{Y}_i-\bar{Y})^T
				      +l(\bar{Y}_i-\bar{Y})(\bar{Y}_i-\bar{Y})^T
				 \right)
\\
& = \sum_{i=1}^k \left(
                       \sum_{j=1}^{l} (Y_{ij}-\bar{Y}_i)(Y_{ij}-\bar{Y}_i)^T
   				      +2\left(\sum_{j=1}^{l} Y_{ij}-\sum_{j=1}^{l} Y_{ij}\right) (\bar{Y}_i-\bar{Y})^T
				      +l(\bar{Y}_i-\bar{Y})(\bar{Y}_i-\bar{Y})^T
				 \right)
\\
& = \sum_{i=1}^k \left(
                       \sum_{j=1}^{l} (Y_{ij}-\bar{Y}_i)(Y_{ij}-\bar{Y}_i)^T
   				       +l(\bar{Y}_i-\bar{Y})(\bar{Y}_i-\bar{Y})^T
				 \right)
\\
& = \sum_{i=1}^kl(\bar{Y}_i-\bar{Y})(\bar{Y}_i-\bar{Y})^T
   +\sum_{i=1}^k \sum_{j=1}^{l} (Y_{ij}-\bar{Y}_i)(Y_{ij}-\bar{Y}_i)^T
\\
& = B + W.
\end{split}
\end{align*}


# Modellevaluation
\vspace{3mm}

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

  # Diese Funktion evaluiert die Kreuzproduktsummenmatrizen T,B,W einer
  # einfaktoriellen Varianzanalyse basierend auf einen m x l x k Datensatz Y.
  #
  # Input
  #     Y        : m x l x k Datenarray
  #
  # Output
  #     $Y_bar   : m x 1 Gesamtmittelwert
  #     $Y_bar_i : m x k 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
  l        = d[2]                                      # Anzahl Datenpunkte pro Gruppe
  k        = 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:k){
      for(j in 1:l){
          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:k){
    B = B + l*(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:k){
      for(j in 1:l){
          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))}
```

# Modellevaluation


\textcolor{darkblue}{Überblick zu Teststatistiken und ihren Verteilungen}
\footnotesize

Basierend auf der Kreuzproduktsummenmatrixzerlegung $T = B + W$ wurden eine 
Reihe von Teststatisken für Hypothesentests der Nullhypothese $H_0 : \mu_1 = \cdots \mu_k$
vorgeschlagen. Wir betrachten hier (nur)
\begin{equation}
\mbox{Wilks' } \Lambda = \frac{|W|}{|B + W|} \mbox{ und  Pillai's } V := \mbox{tr}\left((B+W)^{-1}B \right)
\end{equation}
Im Gegensatz zur $T^2$-Teststatistik und zur F-Teststatistik der univariaten 
Varianzanalyse sind die Verteilungen von $\Lambda$ und $V$ nur 
für bestimmte Anwendungsszenarien, d.h. kleine Werte der Datendimensionalität $m$
und die Anzahl der Gruppen $k$, analytisch exakt beschreibbar und führen, wie im Fall
der T$^2$-Teststatistik auf $f$-Verteilungen. 

Für Anwendungsszenarien mit größeren Werten von $m$ und oder $k$ existieren 
lediglich Approximationen der Verteilungen von $\Lambda$ und $V$, die für 
unendlich große Stichprobenumfänge $n \to \infty$  exakt sind und wiederum 
durch $f$-Verteilungen gegeben sind.

@anderson_2003, Kapitel 8 gibt eine Einführung in die Approximationstheorie für 
multivariate Modelle, das Wissen um die exakten Verteilungen der Teststatistiken 
um 1970 wird von @rao_1972 zusammengefasst. Im Sinne der Anwendung unterscheiden 
sich Testentscheidungen basierend auf exakten oder approximativen Verteilungen 
von Wilk's $\Lambda$ und Pillai's $V$ nicht. Zur Absicherung dieser Aussage mögen
im konkreten Fall von $m$ und $k$ Simulationen wie im Folgenden diskutiert helfen, 
Unterschiede zwischen $\Lambda$ und $V$ und der Approximation ihrer Verteilungen abzuschätzen.

#
\setstretch{2}
\vfill
\large
Anwendungsszenario und Datendeskription

Modellformulierung und Modellschätzung

**Modellevaluation**

* **Wilks' $\Lambda$**
* Pillai's $V$

\large
Selbstkontrollfragen
\vfill


# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\small
\begin{definition}[Wilks' $\Lambda$]
Es seien das Modell der einfaktoriellen Varianzanalyse sowie die Between Sum of
Squares Matrix $B$ und die Within Sum of Squares Matrix $W$ definiert wie oben.
Dann ist die Wilks' $\Lambda$ Teststatistik definiert als
\begin{equation}
\Lambda := \frac{|W|}{|B+W|},
\end{equation}
wobei $|\cdot|$ die Determinante bezeichnet.
\end{definition}

Bemerkungen

* Intuitiv misst $\Lambda$ das Verhältnis von Residualvariabilität und Gesamtvariabilität.
* Ohne Beweis halten wir fest, dass $\Lambda \in [0,1]$
* Für $\bar{Y}_1 = \cdots = \bar{Y}_p = \bar{Y}$ gilt $B = 0_{mm}$ und damit $\Lambda = 1$.
* Für steigende Unterschiede zwischen den $\bar{Y}_i$ nimmt $|B+W|$ gegenüber $|W|$ zu, $\Lambda$ also ab.
* Kleine Werte von $\Lambda$ sprechen also für eine Abweichung von der Nullhypothese.

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\small
\begin{theorem}[Eigenwertform von Wilks' $\Lambda$]
\justifying
\normalfont
Es seien das Modell der einfaktoriellen Varianzanalyse, die Between Sum of
Squares Matrix $B$, die Within Sum of Squares Matrix $W$ und Wilks' $\Lambda$ definiert
wie oben. Weiterhin seien $\lambda_1,...,\lambda_s$ die Eigenwerte von $W^{-1}B$.
Dann gilt
\begin{equation}
\Lambda = \prod_{i=1}^s \frac{1}{1 + \lambda_i}
\end{equation}
\end{theorem}

Bemerkungen

* Wir verzichten auf einen Beweis.
* Die Matrix $W^{-1}B$ ist das multivariate Analogon zu $\frac{\mbox{SQB}}{\mbox{SQB}}$.

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\footnotesize
\begin{theorem}[Spezielle $H_0$ Verteilungen von Wilks' $\Lambda$ Transformationen]
\justifying
\normalfont
Es seien das Modell der einfaktoriellen Varianzanalyse und Wilks' $\Lambda$
definiert wie oben und es gelte außerdem
\begin{equation}
H_0 : \mu_1 = \mu_2 = \cdots = \mu_p
\mbox{ bzw. }
H_0 : \alpha_1 = \alpha_2 = \cdots = \alpha_p = 0_m
\end{equation}
Dann sind für die in den ersten beiden Tabellenspalten aufgeführten Spezialfällen
die in der dritten Tabellenspalte aufgeführten Statistiken $f$-Zufallsvariablen
und zwar mit den in der vierten Tabellenspalte aufgeführten Parametern.
\vspace{1mm}
\center
\renewcommand{\arraystretch}{1.5}
\begin{tabular}{cccc}
Datendimension $m$
& Gruppenanzahl $k$
& Statistik
& $f$-Verteilungsparameter
\\\hline
Beliebig
& 2
& $\frac{1 - \Lambda}{\Lambda} \frac{n-k-m+1}{m}$
& $m, n-k-m+1$
\\
Beliebig
& 3
& $\frac{1 - \sqrt{\Lambda}}{\sqrt{\Lambda}} \frac{n-k-m+1}{m}$
& $2m, 2(n-k-m+1)$
\\
1
& Beliebig
& $\frac{1 - \Lambda}{\Lambda} \frac{n-k}{k-1}$
& $k-1,n-k$
\\
2
& Beliebig
& $\frac{1 - \sqrt{\Lambda}}{\sqrt{\Lambda}} \frac{n-k-1}{k-1}$
& $2(k-1),2(n-k-1)$
\end{tabular}
\end{theorem}

Bemerkungen

* Die Verteilungen gehen zurück auf @wilks_1932.

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\vspace{1mm}
\small
Simulation spezieller $H_0$ Verteilungen von Wilks' $\Lambda$ Transformationen 
\vspace{1mm}
\setstretch{.8}
\tiny

```{r, eval = F}
# Szenarioparameter
nsm = 1e4                                          # Datensimulationsanzahl
M   = c(3,3,1,2)                                   # Datendimension
K   = c(2,3,4,4)                                   # Gruppenanzahl
l   = 15                                           # Datenpunkte pro Gruppe
N   = l*K                                          # Gesamtanzahl Datenpunkte

# Szenariensimulationen
library(MASS)                                      # R Paket für multivariate Normalverteilungen
nsc = length(M)                                    # Szenarienanzahl
P   = matrix(rep(NaN,nsm*nsc), ncol = nsc)         # Statistik ("Prüfgröße") Array
for(sc in 1:nsc){                                  # Szenarioiterationen

    # Modellparameter
    m     = M[sc]                                  # Datendimension
    k     = K[sc]                                  # Gruppenanzahl
    n     = N[sc]                                  # Gesamtanzahl Datenpunkte
    mu_i  = matrix(rep(1,m), nrow = m)             # Identische Gruppenerwartungswertparameter bei H_0
    Sigma = diag(m)                                # Identische Gruppenkovarianzmatrixparameter

    # Datensimulationen
    for(sm in 1:nsm){

        # Datengeneration
        Y     = array(dim = c(m,l,k))              # Datenanarrayinitialisierung
        for(i in 1:k){                             # Gruppeniterationen
            Y[,,i] =  t(mvrnorm(l,mu_i,Sigma))}    # Datensimulation

        # Analyse
        S        = sos(Y)                          # Stichprobenmittel und Sum of Squares Matrizen
        L        = det(S$W)/det(S$W + S$B)         # Wilks' Lambda

        # Szenarioabhängige Statistik ("Prüfgröße")
        if     (sc == 1){p = ((1-L)/L)*((n-k-m+1)/m)}
        else if(sc == 2){p = ((1-sqrt(L))/sqrt(L))*((n-k-m+1)/m)}
        else if(sc == 3){p = ((1-L)/L)*((n-k)/(k-1))}
        else if(sc == 4){p = ((1-sqrt(L))/sqrt(L))*((n-k-1)/(k-1))}

        # Statistikrealisation
        P[sm,sc] = p }}
```

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\vspace{1mm}
\small
Simulation spezieller $H_0$ Verteilungen von Wilks' $\Lambda$ Transformationen


```{r, echo = F, eval = F}
# figure parameters
dev.new()
graphics.off()
par(
family      = "sans",
mfcol       = c(2,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex.main    = 1)
ylimits     = list(c(0,.8),c(0,1),c(0,1),c(0,1))

# Szenarioiterationen
for(sc in 1:nsc){

    # Statistik Histogramm und WDF
    f_min = 0
    f_max = 5
    f_res = 1e3
    f     = seq(f_min,f_max,len =f_res)
    m     = M[sc]
    k     = K[sc]
    n     = N[sc]
    if(sc == 1){
        df = c(m,n-k-m+1)}
    else if(sc == 2){
        df = c(2*m, 2*(n-k-m+1))}
    else if(sc == 3){
        df = c(k-1, n-k)}
    else if(sc == 4){
        df = c(2*(k-1), 2*(n-k-1))}
    p           = df(f,df[1],df[2])
    hist(
    P[,sc],
    breaks      = 50,
    col         = "gray90",
    prob        = TRUE,
    xlim        = c(f_min, f_max),
    ylim        = ylimits[[sc]],
    xlab        = TeX("Statistik"),
    ylab        = "",
    main        = sprintf("m = %d, k = %d", M[sc], K[sc]))
    lines(
    f,
    p,
    lwd         = 2,
    col         = "darkorange")
}
fdir        =  file.path(getwd(), "11_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_11_spezielle_lambda_verteilungen.pdf"),
width       = 6,
height      = 6)
```


```{r, echo = FALSE, out.width = "62%"}
knitr::include_graphics("11_Abbildungen/mvda_11_spezielle_lambda_verteilungen.pdf")
```

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\footnotesize
\begin{theorem}[Approximative $H_0$ Verteilungen von Wilks' $\Lambda$ Transformationen]
\normalfont
\justifying
Es seien das Modell der einfaktoriellen Varianzanalyse und Wilks' $\Lambda$
definiert wie oben und es gelte außerdem die Nullhypothese
\begin{equation}
H_0 : \mu_1 = \mu_2 = \cdots = \mu_p
\mbox{ bzw. }
H_0 : \alpha_1 = \alpha_2 = \cdots = \alpha_p = 0_m
\end{equation}
Dann ist die Statistik
\begin{equation}
\tau := \frac{1 - \Lambda^{1/t}}{\Lambda^{1/t}} \frac{\nu_2}{\nu_1}
\end{equation}
mit
\begin{equation}
\nu_1 := m(k-1)
\mbox{ und }
\nu_2 := wt-\frac{1}{2}(m(k-1)-2)
\end{equation}
sowie
\begin{equation}
w := n-1-\frac{1}{2}(m+k)
\mbox{ und }
t   := \sqrt{\frac{m^2(k-1)^2 - 4}{m^2 + (k-1)^2 - 5}}
\end{equation}
approximativ $f$-verteilt mit Freiheitsgradparametern $\nu_1$ und $\nu_2$.
\end{theorem}

Bemerkungen

* Die Approximation geht zurück auf @rao_1951.

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\vspace{2mm}
\small
Simulation approximativer $H_0$ Verteilungen von Wilks' $\Lambda$ Transformationen 
\vspace{1mm}
\setstretch{.8}
\tiny

```{r, eval = F}
# Szenarioparameter
library(MASS)                                          # # R Paket für multivariate Normalverteilungen
nsm = 1e4                                              # Datensimulationsanzahl
M   = c(3,3,4,9)                                       # Datendimension
K   = c(4,9,4,4)                                       # Gruppenanzahl
l   = 15                                               # Datenpunkte pro Gruppe
N   = l*K                                              # Gesamtanzahl Datenpunkte
nsc = length(M)                                        # Szenarienanzahl
TAU = matrix(rep(NaN,nsm*nsc), ncol = nsc)             # Statistik Array
NU  = matrix(rep(NaN,2*nsc)  , ncol = nsc)             # Parameter Array
WL  = seq(0,1,len = 1e3)                               # Wilk's Lambda Values
TL  = matrix(rep(NaN,length(WL)*nsc), nrow = nsc)      # \tau(\Lambda) Array
for(sc in 1:nsc){                                      # Szenarioiterationen

    # Modellparameter
    m       = M[sc]                                    # Datendimension
    k       = K[sc]                                    # Gruppenanzahl
    n       = N[sc]                                    # Gesamtanzahl Datenpunkte
    mu_i    = matrix(rep(1,m), nrow = m)               # Identische Gruppenerwartungswertparameter bei H_0
    Sigma   = diag(m)                                  # Identische Gruppenkovarianzmatrixparameter

    # Varianzanalyse Parameter
    w       = n-1-(1/2)*(m+k)                          # w
    t       = sqrt((m^2*(k-1)^2-4)/(m^2+(k-1)^2-5))    # t
    nu_1    = m*(k-1)                                  # \nu_1
    nu_2    = w*t-(1/2)*(m*(k-1)-2)                    # \nu_2
    TL[sc,] = ((1-WL^(1/t))/WL^(1/t))*(nu_2/nu_1)      # \tau(\Lambda)

    # Datensimulationen
    for(sm in 1:nsm){

        # Datengeneration
        Y     = array(dim = c(m,l,k))                  # Datenanarrayinitialisierung
        for(i in 1:k){                                 # Gruppeniterationen
            Y[,,i] =  t(mvrnorm(l,mu_i,Sigma))}        # Datensimulation

        # Varianzanalyse
        S          = sos(Y)                            # Stichprobenmittel und Sum of Squares Matrizen
        L          = det(S$W)/det(S$W + S$B)           # Wilks' Lambda
        tau        = ((1-L^(1/t))/L^(1/t))*(nu_2/nu_1) # Statistik
        TAU[sm,sc] = tau                               # Statistik
        NU[1,sc]   = nu_1                              # \nu_1
        NU[2,sc]   = nu_2}}                            # \nu_2
```

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\vspace{1mm}
$\tau$ als Funktion von $\Lambda$ : $\Lambda \downarrow \Rightarrow \tau \uparrow$
\vspace{2mm}
```{r, echo = F, eval = F}
# figure parameters
dev.new()
graphics.off()
par(
family      = "sans",
mfcol       = c(1,1),
pty         = "m",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex.main    = 1)

cols        = c("gray20","gray40","gray60","gray80")
matplot(
as.matrix(WL),
t(TL),
ylim = c(0,20),
type = "l",
lty  = 1,
col  = cols,
xlab = TeX("$\\Lambda$"),
ylab = TeX("$\\tau$"))

legend(
0.6,
20,
c("m = 3, k = 4, n = 15",
  "m = 3, k = 9, n = 15",
  "m = 4, k = 4, n = 15",
  "m = 9, k = 4, n = 15"),
lty = 1,
col = cols,
pch = NA,
bty = "n",
y.intersp = 1.2,
cex = 0.9)

fdir        =  file.path(getwd(), "11_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_11_tau_lambda.pdf"),
width       = 6,
height      = 4)
```

```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("11_Abbildungen/mvda_11_tau_lambda.pdf")
```


# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\vspace{1mm}
\small
Simulation approximativer $H_0$ Verteilungen von Wilks' $\Lambda$ Transformationen


```{r, echo = F, eval = F}
# figure parameters
dev.new()
graphics.off()
par(
family      = "sans",
mfcol       = c(2,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex.main    = 1)
ylimits     = list(c(0,1.6),c(0,1.6),c(0,1.6),c(0,1.6))

# Szenarioiterationen
for(sc in 1:nsc){

    # Statistik Histogramm und WDF
    f_min = 0
    f_max = 3
    f_res = 1e3
    f     = seq(f_min,f_max,len =f_res)
    p     = df(f,NU[1,sc], NU[2,sc])
    hist(
    TAU[,sc],
    breaks      = 50,
    col         = "gray90",
    prob        = TRUE,
    xlim        = c(f_min, f_max),
    ylim        = ylimits[[sc]],
    xlab        = TeX("Statistik"),
    ylab        = "",
    main        = sprintf("m = %d, k = %d", M[sc], K[sc]))
    lines(
    f,
    p,
    lwd         = 2,
    col         = "darkorange")
}
fdir        =  file.path(getwd(), "11_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_11_approximative_lambda_verteilungen.pdf"),
width       = 6,
height      = 6)
```


```{r, echo = FALSE, out.width = "62%"}
knitr::include_graphics("11_Abbildungen/mvda_11_approximative_lambda_verteilungen.pdf")
```


# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\small
\begin{definition}[Wilks' $\Lambda$-basierter Test, Testumfangkontrolle, p-Wert]
\justifying
Es seien das Modell der einfaktoriellen Varianzanalyse und die Wilks' $\Lambda$
basierte Teststatistik $\tau$ mit Verteilungsparametern $\nu_1,\nu_2$ wie oben definiert.
Weiterhin sei der kritische Wert-basierte Test
\begin{equation}
\phi(Y) := 1_{\{\tau > k\}} :=
\begin{cases}
1 & \tau >   k \\
0 & \tau \le k
\end{cases}
\end{equation}
definiert. $\phi$ ist genau dann ein Level-$\alpha_0$-Test mit Testumfang $\alpha$, wenn
\begin{equation}
k := k_{\alpha_0} := F^{-1}(1-\alpha_0; \nu_1,\nu_2)
\end{equation}
ist und der p-Wert einer realisiertern $\tau$-Teststatistik $\tilde{\tau}$ ergibt sich zu
\begin{equation}
\mbox{p-Wert} = \mathbb{P}(\tau \ge \tilde{\tau}) = 1 - F(\tilde{\tau}; \nu_1,\nu_2)
\end{equation}
\end{definition}


\footnotesize
Bemerkungen

* Ein Beweis kann in Analogie zum Einstichproben-T$^2$-Test Fall geführt werden.
* Wir validieren die Testumfangkontrolle mithilfe von $k_{\alpha_0}$ in untenstehender Simulation.

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
Testumfangkontrolle

\vspace{2mm}
\small
\center
Wahl von $k_{\alpha_0}$ bei $m = 3,k = 4, n = 15 \Rightarrow \nu_1 = 9, \nu_2 = 132$ und $\alpha_0 = 0.05$.
\vspace{2mm}

```{r,echo = F, eval = F}

par(
family      = "sans",
mfcol       = c(1,2),
pty         = "m",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 1.2)

# Parameter
alpha_0     = 0.05
k_alpha_0   = qf(1 - alpha_0, NU[1,1],NU[2,1])
tau         = seq(0,3,length=1e4 )
Ptau        = pf(tau,NU[1,1],NU[2,1])
ptau        = df(tau,NU[1,1],NU[2,1])

# KVF Perspektive
plot(
tau,
Ptau,
type        = "l",
xlab        = TeX("$\\tau$"),
ylab        = " ",
ylim        = c(0,1),
main        = TeX("$F(\\tau;9,132)$"))

lines(
k_alpha_0,
0,
type        = "p",
pch         = 16,
xpd         = TRUE)

lines(
min(tau),
1 - alpha_0,
type        = "p",
pch         = 16,
xpd         = TRUE)

arrows(
x0          = min(tau),
y0          = 1 - alpha_0,
x1          = k_alpha_0,
y1          = 1 - alpha_0,
col         = "darkorange",
angle       = 45,
length      = .1)

arrows(
x0          = k_alpha_0,
y0          = 1-alpha_0,
x1          = k_alpha_0,
y1          = 0,
col         = "darkorange",
angle       = 45,
length      = .1)

text(k_alpha_0-.2, .05 , TeX("$\\k_{\\alpha_0}$"), xpd = TRUE)
text(1      , 1 , TeX("$1 - \\alpha_0$"), xpd = TRUE)

# WDF Perspektive
plot(
tau,
ptau,
type        = "l",
ylab        = " ",
xlab        = TeX("$\\tau$"),
ylim        = c(0,1),
main        = TeX("$f(\\tau;9,132)$"))

polygon(
c(tau[tau  >= k_alpha_0] , max(tau), k_alpha_0),
c(ptau[tau >= k_alpha_0],       0, 0),
col = "gray90",
border = NA)


lines(
seq(k_alpha_0, max(tau), len = 1e2),
rep(0,1e2),
type        = "l",
lwd         = 5,
col         = "darkorange")

lines(
k_alpha_0,
0,
type        = "p",
pch         = 16,
xpd         = TRUE)

text( k_alpha_0-.2,.05, TeX("$\\k_{\\alpha_0}$") , xpd = TRUE)
text( 2.35, .2, TeX("$P(\\tau > k_{\\alpha_0}) = \\alpha_0$"), xpd = TRUE, cex = 1, col = "gray50")

fdir        =  file.path(getwd(), "11_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_11_kritischer_wert_lambda.pdf"),
width       = 8,
height      = 4)
```

```{r, echo = FALSE, out.width = "90%"}
knitr::include_graphics("11_Abbildungen/mvda_11_kritischer_wert_lambda.pdf")
```

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
Testumfangkontrolle
\vspace{1mm}

\setstretch{0.8}
\tiny
```{r}
# Szenarioparameter
nsm       = 1e4                                        # Datensimulationsanzahl
M         = c(3,3,4,9)                                 # Datendimension
K         = c(4,9,4,4)                                 # Gruppenanzahl
l         = 15                                         # Datenpunkte pro Gruppe
N         = l*K                                        # Gesamtanzahl Datenpunkte
alpha_0   = 0.05                                       # \alpha_0
nsc       = length(M)                                   # Szenarienanzahl
TAU       = matrix(rep(NaN,nsm*nsc), ncol = nsc)       # Statistik Array
NU        = matrix(rep(NaN,2*nsc)  , ncol = nsc)       # Parameter Array
KA        = rep(NaN, nsc)                              # Kritische Werte
PHI       = matrix(rep(0,nsm*nsc)  , ncol = nsc)       # Testarray

# Simulationen
for(sc in 1:nsc){                                      # Szenarioiterationen

    # Modellparameter
    m         = M[sc]                                  # Datendimension
    k         = K[sc]                                  # Gruppenanzahl
    n         = N[sc]                                  # Gesamtanzahl Datenpunkte
    mu_i      = matrix(rep(1,m), nrow = m)             # Identische Gruppenerwartungswertparameter bei H_0
    Sigma     = diag(m)                                # Identische Gruppenkovarianzmatrixparameter

    # Varianzanalyse Parameter
    w         = n-1-(1/2)*(m+k)                        # w
    t         = sqrt((m^2*(k-1)^2-4)/(m^2+(k-1)^2-5))  # t
    nu_1      = m*(k-1)                                # \nu_1
    nu_2      = w*t-(1/2)*(m*(k-1)-2)                  # \nu_2
    KA[sc]    = qf(1-alpha_0,nu_1,nu_2)                # kritischer Wert

    # Datensimulationen
    for(sm in 1:nsm){
        Y     = array(dim = c(m,l,k))                  # Datenanarrayinitialisierung
        for(i in 1:k){                                 # Gruppeniterationen
            Y[,,i] =  t(mvrnorm(l,mu_i,Sigma))}        # Datensimulation
        S          = sos(Y)                            # Stichprobenmittel und Sum of Squares Matrizen
        L          = det(S$W)/det(S$W + S$B)           # Wilks' Lambda
        tau        = ((1-L^(1/t))/L^(1/t))*(nu_2/nu_1) # Statistik
        PHI[sm,sc] = tau > KA[sc]}}                    # Test

```

```{r, echo = F}
cat("Kritische Werte       : ", KA,
    "\nGeschätzte Testumfänge: ", apply(PHI, 2, mean))
```


# Modellevaluation mit der Wilks' $\Lambda$ Statistik

Praktisches Vorgehen
\footnotesize

* \justifying Man nimmt an, dass ein vorliegender Datensatz von $k$ Gruppendatensätzen
$y_{11}, ...,y_{1l}$, $y_{21}, ...,y_{2l}$, ..., $y_{k1}, ...,y_{kl}$
Realisationen von $Y_{ij} \sim N(\mu_i,\Sigma)$ mit unbekannten Parametern
$\mu_i \in \mathbb{R}^m, i = 1,...,k$ und $\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d.}$ sind.

* Man möchte entscheiden ob $H_0 : \mu_1 = \cdots = \mu_k$ eher zutrifft oder eher nicht.

* Man wählt ein Signifikanzniveau $\alpha_0$ und bestimmt den zugehörigen
Freiheitsgradparameter-abhängigen kritischen Wert $k_{\alpha_0}$. Zum Beispiel gilt
bei Wahl von $\alpha_0 := 0.05, m = 3, k = 4, l = 15$ und somit $n = 60$ sowie
$\nu_1 = 9, \nu_2 = 132$, dass $k_{\alpha_0} = F^{-1}(1 - 0.05; 9, 132) \approx 1.95$ ist.

* Anhand der Wilk's $\Lambda$ Statistik sowie $m,k$ und $n$ berechnet man  man den realisierten
Wert der $\tau$-Teststatistik, den wir hier mit $\tilde{\tau}$ bezeichnen.

* Wenn $\tilde{\tau}$ größer als $k_{\alpha_0}$ ist, lehnt man die Nullhypothese ab,
andernfalls nicht.

* Die oben entwickelte Theorie garantiert dann, dass man im Mittel in höchstens
$\alpha_0 \cdot 100$ von 100 Fällen die Nullhypothese fälschlicherweise ablehnt.

* Schließlich ergibt sich der assoziierte p-Wert der realisiertern $\tau$-Teststatistik
$\tilde{\tau}$ zu
\begin{equation}
\mbox{p-Wert} = \mathbb{P}(\tau \ge \tilde{\tau}) = 1 -  F(\tilde{\tau}; \nu_1,\nu_2)
\end{equation}

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\textcolor{darkblue}{Anwendungsbeispiel}
\vspace{1mm}

Einfaktorielle Varianzanalyse mit R's `lm()` und `Manova()`
\vspace{2mm}

\footnotesize
```{r, message = F}
# Einlesen des Datensatzes aus Rudolf & Buse (2019) Multivariate Verfahren
library(foreign)
library(car)
D        = read.spss(file.path(getwd(), "11_Daten", "studienerfolg.sav"), to.data.frame = T)
model    = lm(cbind(D$X1,D$X2) ~ D$Gruppe, D)      # Modellspezifikation
Manova(model, test.statistic = "Wilks")            # Einfaktorielle Varianzanalyse
```

# Modellevaluation mit der Wilks' $\Lambda$ Statistik
\textcolor{darkblue}{Anwendungsbeispiel}

\small
Einfaktorielle Varianzanalyse mit `sos()`
\vspace{1mm}
\tiny
\setstretch{.8}
```{r, message = F}
# Einlesen des Datensatzes aus Rudolf & Buse (2019) Multivariate Verfahren
library("foreign")
D      = read.spss(file.path(getwd(), "11_Daten", "studienerfolg.sav"), to.data.frame = T)

# Datepräprozessierung
m      = 2                                       # Datendimension von Interesse
k      = 3                                       # Anzahl Gruppen
l      = 15                                      # Anzahl Datenpunkte pro Gruppe
n      = k*l                                     # Gesamtdatenpunktanzahl
Y      = array(dim = c(m,l,k))                   # Datenarrayinitialisierung
Y[,,1] = rbind(D$X1[D$Gruppe == "ungenügend"],   # Y_{1j_1} Intelligenztestscore
               D$X2[D$Gruppe == "ungenügend"])   # Y_{1j_2} Mathematiktestscore
Y[,,2] = rbind(D$X1[D$Gruppe == "befriedigend"], # Y_{2j_1} Intelligenztestscore
               D$X2[D$Gruppe == "befriedigend"]) # Y_{2j_2} Mathematiktestscore
Y[,,3] = rbind(D$X1[D$Gruppe == "gut"],          # Y_{3j_1} Intelligenztestscore
               D$X2[D$Gruppe == "gut"])          # Y_{3j_2} Mathematiktestscore

# Einfaktorielle Varianzanalyse
S      = sos(Y)                                  # Sum of Squares Matrizen
L      = det(S$W)/det(S$W + S$B)                 # Wilks' Lambda
w      = n-1-(1/2)*(m+k)                         # w
t      = sqrt((m^2*(k-1)^2-4)/(m^2+(k-1)^2-5))   # t
nu_1   = m*(k-1)                                 # \nu_1
nu_2   = w*t-(1/2)*(m*(k-1)-2)                   # \nu_2
tau    = ((1-L^(1/t))/L^(1/t))*(nu_2/nu_1)       # Teststatistik
P      = 1-pf(tau,nu_1,nu_2)                     # Überschreitungswahrscheinlichkeit

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

#
\setstretch{2}
\vfill
\large
Anwendungsszenario und Datendeskription

Modellformulierung und Modellschätzung

**Modellevaluation**

* Wilks' $\Lambda$
* **Pillai's $V$**

Selbstkontrollfragen
\vfill

# Modellevaluation mit der Pillai's $V$ Statistik
\small
\begin{definition}[Pillai's $V$ Statistik]
Es seien das Modell der einfaktoriellen Varianzanalyse sowie die Between Sum of
Squares Matrix $B$ und die Within Sum of Squares Matrix $W$ definiert wie oben.
Dann ist die \textit{Pillai's $V$ Statistik} definiert als
\begin{equation}
V := \mbox{tr}\left((B + W)^{-1}B \right)
\end{equation}
wobei $\mbox{tr}(\cdot)$ die Spur, also die Summe der Diagonalelemente, einer Matrix
bezeichnet.
\end{definition}

\footnotesize
Bemerkungen

* Die Pillai's $V$ Statistik betrachtet die Diagonalelemente von $T^{-1}B$
* Ein hoher Wert von $V$ spricht also für einen großen Anteil der Between-Varianz an der Gesamtvarianz.
* Für einen hohen Wert von $V$ würde man also die Nullhypothese $B = 0_{mm}$ verwerfen.

# Modellevaluation mit der Pillai's $V$ Statistik
\small
\begin{theorem}[Eigenwertform der Pillai's $V$ Statistik]
\justifying
\normalfont
Es seien das Modell der einfaktoriellen Varianzanalyse, die Between Sum of
Squares Matrix $B$, die Within Sum of Squares Matrix $W$ und die Pillai's $V$ Statistik
definiert wie oben. Weiterhin seien $\lambda_1,...,\lambda_s$ die Eigenwerte
von $W^{-1}B$. Dann gilt
\begin{equation}
V = \sum_{i=1}^s \frac{\lambda_i}{1 + \lambda_i}
\end{equation}
\end{theorem}

footnotesize
Bemerkungen

* Wir verzichten auf einen Beweis.
* Die Matrix $W^{-1}B$ ist das multivariate Analogon zu $\frac{\mbox{SQB}}{\mbox{SQW}}$

# Modellevaluation mit der Pillai's $V$ Statistik
\small
\begin{theorem}[Approximative $H_0$ Verteilungen von Pillai's $V$ Transformationen]
\normalfont
\justifying
Es seien das Modell der einfaktoriellen Varianzanalyse und Pillai's $V$ definiert
wie oben und es gelte außerdem die Nullhypothese
\begin{equation}
H_0 : \mu_1 = \mu_2 = \cdots = \mu_p
\mbox{ bzw. }
H_0 : \alpha_1 = \alpha_2 = \cdots = \alpha_p = 0_m.
\end{equation}
Weiterhin seien die Parameter
\begin{equation}
s := \min(k-1,m),
t := \frac{1}{2}(|k-1-m|-1)
\mbox{ und }
w := \frac{1}{2}(n-k-m-1)
\end{equation}
definiert. Dann ist
\begin{equation}
\tau = \frac{(2w+s+1)V}{(2t+s+1)(s-V)}
\end{equation}
approximativ $f$-verteilt mit Freiheitsgradparametern
\begin{equation}
\nu_1 := s(2t+s+1)
\mbox{ und }
\nu_2:=  s(2w+s+1)
\end{equation}
\end{theorem}

\footnotesize
Bemerkungen

* Die Approximation geht zurück auf @pillai_1955.

# Modellevaluation mit der Pillai's $V$ Statistik
\vspace{2mm}
\small
Simulation approximativer $H_0$ Verteilungen von Pillai's $V$ Transformatione
\vspace{1mm}
\setstretch{.8}
\tiny

```{r, eval = F}
# Szenarioparameter
library(MASS)                                          # R Paket für multivariate Normalverteilungen
library(matlib)                                        # R Paket für Matrixalgebra
nsm = 1e1                                              # Datensimulationsanzahl
M   = c(3,3,4,9)                                       # Datendimension
K   = c(4,9,4,4)                                       # Gruppenanzahl
l   = 15                                               # Datenpunkte pro Gruppe
N   = l*K                                              # Gesamtanzahl Datenpunkte
nsc = length(M)                                        # Szenarienanzahl
TAU = matrix(rep(NaN,nsm*nsc), ncol = nsc)             # Statistik Array
NU  = matrix(rep(NaN,2*nsc)  , ncol = nsc)             # Parameter Array
for(sc in 1:nsc){                                      # Szenarioiterationen

    
    # Modellparameter
    m     = M[sc]                                      # Datendimension
    k     = K[sc]                                      # Gruppenanzahl
    n     = N[sc]                                      # Gesamtanzahl Datenpunkte
    mu_i  = matrix(rep(1,m), nrow = m)                 # Identische Gruppenerwartungswertparameter bei H_0
    Sigma = diag(m)                                    # Identische Gruppenkovarianzmatrixparameter

    # Varianzanalyseparameter
    s     = min(k-1,m)                                 # s
    t     = (1/2)*(abs(k-1-m)-1)                       # t
    w     = (1/2)*(n-k-m-1)                            # w
    nu_1  = s*(2*t+s+1)                                # \nu_1
    nu_2  = s*(2*w+s+1)                                # \nu_2

    # Datensimulationen
    for(sm in 1:nsm){

        # Datengeneration
        Y     = array(dim = c(m,l,k))                  # Datenanarrayinitialisierung
        for(i in 1:k){                                 # Gruppeniterationen
            Y[,,i] =  t(mvrnorm(l,mu_i,Sigma))}        # Datensimulation

        # Varianzanalyse
        S          = sos(Y)                            # Stichprobenmittel und Sum of Squares Matrizen
        V          = sum(diag(inv(S$B+S$W) %*% S$B))   # Pillai's V
        tau        = ((2*w+s+1)*V)/((2*t+s+1)*(s-V))   # Statistik
        TAU[sm,sc] = tau                               # Statistik
        NU[1,sc]   = nu_1                              # \nu_1
        NU[2,sc]   = nu_2}}                            # \nu_2
```

# Modellevaluation mit der Pillai's $V$ Statistik
\vspace{1mm}
\small
Simulation approximativer $H_0$ Verteilungen von Pillai's $V$ Transformationen
\vspace{1mm}

```{r, echo = F, eval = F}
# figure parameters
dev.new()
graphics.off()
par(
family      = "sans",
mfcol       = c(2,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex.main    = 1)
ylimits     = list(c(0,1.6),c(0,1.6),c(0,1.6),c(0,1.6))

# Szenarioiterationen
for(sc in 1:nsc){

    # Statistik Histogramm und WDF
    f_min = 0
    f_max = 3
    f_res = 1e3
    f     = seq(f_min,f_max,len =f_res)
    p     = df(f,NU[1,sc], NU[2,sc])
    hist(
    TAU[,sc],
    breaks      = 50,
    col         = "gray90",
    prob        = TRUE,
    xlim        = c(f_min, f_max),
    ylim        = ylimits[[sc]],
    xlab        = TeX("Statistik"),
    ylab        = "",
    main        = sprintf("m = %d, k = %d", M[sc], K[sc]))
    lines(
    f,
    p,
    lwd         = 2,
    col         = "darkorange")
}
fdir        =  file.path(getwd(), "11_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_11_approximative_v_verteilungen.pdf"),
width       = 6,
height      = 6)
```

```{r, echo = FALSE, out.width = "62%"}
knitr::include_graphics("11_Abbildungen/mvda_11_approximative_v_verteilungen.pdf")
```


# Modellevaluation mit der Pillai's $V$ Statistik
\small
\begin{definition}[Pillai's $V$ -basierter Test, Testumfangkontrolle, p-Wert]
\justifying
Es seien das Modell der einfaktoriellen Varianzanalyse und die Pillai's $V$
basierte Teststatistik $\tau$ mit Verteilungsparametern $\nu_1,\nu_2$ wie oben definiert.
Weiterhin sei der kritische Wert-basierte Test
\begin{equation}
\phi(Y) := 1_{\{\tau > k\}} :=
\begin{cases}
1 & \tau >   k \\
0 & \tau \le k
\end{cases}
\end{equation}
definiert. $\phi$ ist genau dann ein Level-$\alpha_0$-Test mit Testumfang $\alpha$, wenn
\begin{equation}
k := k_{\alpha_0} := F^{-1}(1-\alpha_0; \nu_1,\nu_2)
\end{equation}
ist und der p-Wert einer realisiertern $\tau$-Teststatistik $\tilde{\tau}$ ergibt sich zu
\begin{equation}
\mbox{p-Wert} = \mathbb{P}(\tau \ge \tilde{\tau}) = 1 - F(\tilde{\tau}; \nu_1,\nu_2)
\end{equation}
\end{definition}

# Modellevaluation mit der Pillai's $V$ Statistik

Testumfangkontrolle
\vspace{2mm}

\small
\center
Wahl von $k_{\alpha_0}$ bei $m = 3,k = 4, n = 15 \Rightarrow \nu_1 = 9, \nu_2 = 168$ und $\alpha_0 = 0.05$.
\vspace{2mm}

```{r,echo = F, eval = F}
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "m",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 1.2)

# Parameter
alpha_0     = 0.05
k_alpha_0   = qf(1 - alpha_0, NU[1,1],NU[2,1])
tau         = seq(0,3,length=1e4 )
Ptau        = pf(tau,NU[1,1],NU[2,1])
ptau        = df(tau,NU[1,1],NU[2,1])

# KVF Perspektive
plot(
tau,
Ptau,
type        = "l",
xlab        = TeX("$\\tau$"),
ylab        = " ",
ylim        = c(0,1),
main        = TeX("$F(\\tau;9,168)$"))

lines(
k_alpha_0,
0,
type        = "p",
pch         = 16,
xpd         = TRUE)

lines(
min(tau),
1 - alpha_0,
type        = "p",
pch         = 16,
xpd         = TRUE)

arrows(
x0          = min(tau),
y0          = 1 - alpha_0,
x1          = k_alpha_0,
y1          = 1 - alpha_0,
col         = "darkorange",
angle       = 45,
length      = .1)

arrows(
x0          = k_alpha_0,
y0          = 1-alpha_0,
x1          = k_alpha_0,
y1          = 0,
col         = "darkorange",
angle       = 45,
length      = .1)

text(k_alpha_0-.2, .05 , TeX("$\\k_{\\alpha_0}$"), xpd = TRUE)
text(1      , 1 , TeX("$1 - \\alpha_0$"), xpd = TRUE)

# WDF Perspektive
plot(
tau,
ptau,
type        = "l",
ylab        = " ",
xlab        = TeX("$\\tau$"),
ylim        = c(0,1),
main        = TeX("$f(\\tau;9,168)$"))

polygon(
c(tau[tau  >= k_alpha_0] , max(tau), k_alpha_0),
c(ptau[tau >= k_alpha_0],       0, 0),
col = "gray90",
border = NA)


lines(
seq(k_alpha_0, max(tau), len = 1e2),
rep(0,1e2),
type        = "l",
lwd         = 5,
col         = "darkorange")

lines(
k_alpha_0,
0,
type        = "p",
pch         = 16,
xpd         = TRUE)

text( k_alpha_0-.2,.05, TeX("$\\k_{\\alpha_0}$") , xpd = TRUE)
text( 2.35, .2, TeX("$P(\\tau > k_{\\alpha_0}) = \\alpha_0$"), xpd = TRUE, cex = 1, col = "gray50")

fdir        =  file.path(getwd(), "11_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_11_kritischer_wert_v.pdf"),
width       = 8,
height      = 4)
```

```{r, echo = FALSE, out.width = "90%"}
knitr::include_graphics("11_Abbildungen/mvda_11_kritischer_wert_v.pdf")
```

# Modellevaluation mit der Pillai's $V$ Statistik
Testumfangkontrolle
\vspace{1mm}

\setstretch{0.8}
\tiny
```{r, eval = F}
# Szenarioparameter
library(MASS)                                          # R Paket für multivariate Normalverteilungen
library(matlib)                                        # R Paket für Matrixalgebra
nsm     = 1e3                                          # Datensimulationsanzahl
M       = c(3,3,4,9)                                   # Datendimension
K       = c(4,9,4,4)                                   # Gruppenanzahl
l       = 15                                           # Datenpunkte pro Gruppe
N       = l*K                                          # Gesamtanzahl Datenpunkte
alpha_0 = 0.05                                         # Signifikanzlevel
nsc     = length(M)                                    # Szenarienanzahl
KA      = rep(NaN, nsc)                                # kritische Werte
PHI     = matrix(rep(0,nsm*nsc), ncol = nsc)           # Testarray

# Szenariosimulationen
for(sc in 1:nsc){                                      # Szenarioiterationen

    # Modell- und Varianzanalyseparameter
    m      = M[sc]                                     # Datendimension
    k      = K[sc]                                     # Gruppenanzahl
    n      = N[sc]                                     # Gesamtanzahl Datenpunkte
    mu_i   = matrix(rep(1,m), nrow = m)                # Identische Gruppenerwartungswertparameter bei H_0
    Sigma  = diag(m)                                   # Identische Gruppenkovarianzmatrixparameter
    s      = min(k-1,m)                                # s
    t      = (1/2)*(abs(k-1-m)-1)                      # t
    w      = (1/2)*(n-k-m-1)                           # w
    nu_1   = s*(2*t+s+1)                               # \nu_1
    nu_2   = s*(2*w+s+1)                               # \nu_2
    KA[sc] = qf(1-alpha_0,nu_1,nu_2)                   # kritischer Wert 
    
    # Datensimulationen
    for(sm in 1:nsm){

        # Datengeneration und Varianzanalyse
        Y     = array(dim = c(m,l,k))                  # Datenanarrayinitialisierung
        for(i in 1:k){                                 # Gruppeniterationen
            Y[,,i] =  t(mvrnorm(l,mu_i,Sigma))}        # Datensimulation
        S          = sos(Y)                            # Stichprobenmittel und Sum of Squares Matrizen
        V          = sum(diag(inv(S$B+S$W) %*% S$B))   # Pillai's V
        tau        = ((2*w+s+1)*V)/((2*t+s+1)*(s-V))   # Statistik
        PHI[sm,sc] = tau > KA[sc]}}                    # Test
```

```{r, echo = F}
cat("Kritische Werte       : ", KA,
    "\nGeschätzte Testumfänge: ", apply(PHI, 2, mean))
```

# Modellevaluation mit der Pillai's $V$ Statistik 
Praktisches Vorgehen
\footnotesize

* \justifying Man nimmt an, dass ein vorliegender Datensatz von $k$ Gruppendatensätzen
$y_{11}, ...,y_{1l}$, $y_{21}, ...,y_{2l}$, ..., $y_{k1}, ...,y_{kl}$
Realisationen von $Y_{ij} \sim N(\mu_i,\Sigma)$ mit unbekannten Parametern
$\mu_i \in \mathbb{R}^m, i = 1,...,k$ und $\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d.}$ sind.

* Man möchte entscheiden ob $H_0 : \mu_1 = \cdots = \mu_k$ eher zutrifft oder eher nicht.

* Man wählt ein Signifikanzniveau $\alpha_0$ und bestimmt den zugehörigen
Freiheitsgradparameter-abhängigen kritischen Wert $k_{\alpha_0}$. Zum Beispiel gilt
bei Wahl von $\alpha_0 := 0.05, m = 3, k = 4, l = 15$ und somit $n = 60$ sowie
$\nu_1 = 9, \nu_2 = 168$, dass $k_{\alpha_0} = F^{-1}(1 - 0.05; 9, 168) \approx 1.94$ ist.

* Anhand der Pillai's $V$ Statistik sowie $m,k$ und $n$ berechnet man  man den realisierten
Wert der $\tau$-Teststatistik, den wir hier mit $\tilde{\tau}$ bezeichnen.

* Wenn $\tilde{\tau}$ größer als $k_{\alpha_0}$ ist, lehnt man die Nullhypothese ab,
andernfalls nicht.

* Die oben entwickelte Theorie garantiert dann, dass man im Mittel in höchstens
$\alpha_0 \cdot 100$ von 100 Fällen die Nullhypothese fälschlicherweise ablehnt.

* Schließlich ergibt sich der assoziierte p-Wert der realisiertern $\tau$-Teststatistik
$\tilde{\tau}$ zu
\begin{equation}
\mbox{p-Wert} = \mathbb{P}(\tau \ge \tilde{\tau}) = 1 -  F(\tilde{\tau}; \nu_1,\nu_2)
\end{equation}


# Modellevaluation mit der Pillai's $V$ Statistik
\textcolor{darkblue}{Anwendungsbeispiel}
\vspace{1mm}

Einfaktorielle Varianzanalyse mit R's `lm()` und `Manova()`
\vspace{2mm}

\footnotesize
```{r, message = F}
# Einlesen des Datensatzes aus Rudolf & Buse (2019) Multivariate Verfahren
library(foreign)
library(car)
D        = read.spss(file.path(getwd(), "11_Daten", "studienerfolg.sav"), to.data.frame = T)
model    = lm(cbind(D$X1,D$X2) ~ D$Gruppe, D)      # Modellspezifikation
Manova(model, test.statistic = "Pillai")           # Einfaktorielle Varianzanalyse
```

# Modellevaluation mit der Pillai's $V$ Statistik
\textcolor{darkblue}{Anwendungsbeispiel}

\small
Einfaktorielle Varianzanalyse mit `sos()`
\vspace{1mm}
\tiny
\setstretch{.8}
```{r, message = F}
# Einlesen des Datensatzes aus Rudolf & Buse (2019) Multivariate Verfahren
library(foreign)
library(matlib)
D      = read.spss(file.path(getwd(), "11_Daten", "studienerfolg.sav"), to.data.frame = T)

# Datepräprozessierung
m      = 2                                       # Datendimension von Interesse
k      = 3                                       # Anzahl Gruppen
l      = 15                                      # Anzahl Datenpunkte pro Gruppe
n      = k*l                                     # Gesamtdatenpunktanzahl
Y      = array(dim = c(m,l,k))                   # Datenarrayinitialisierung
Y[,,1] = rbind(D$X1[D$Gruppe == "ungenügend"],   # Y_{1j_1} Intelligenztestscore
               D$X2[D$Gruppe == "ungenügend"])   # Y_{1j_2} Mathematiktestscore
Y[,,2] = rbind(D$X1[D$Gruppe == "befriedigend"], # Y_{2j_1} Intelligenztestscore
               D$X2[D$Gruppe == "befriedigend"]) # Y_{2j_2} Mathematiktestscore
Y[,,3] = rbind(D$X1[D$Gruppe == "gut"],          # Y_{3j_1} Intelligenztestscore
               D$X2[D$Gruppe == "gut"])          # Y_{3j_2} Mathematiktestscore

# Einfaktorielle Varianzanalyse
S      = sos(Y)                                  # Sum of Squares Matrizen
V      = sum(diag(inv(S$B+S$W) %*% S$B))         # Pillai's V
s      = min(k-1,m)                              # s
t      = (1/2)*(abs(k-1-m)-1)                    # t
w      = (1/2)*(n-k-m-1)                         # w
nu_1   = s*(2*t+s+1)                             # \nu_1
nu_2   = s*(2*w+s+1)                             # \nu_2
tau    = ((2*w+s+1)*V)/((2*t+s+1)*(s-V))         # Statistik
P      = 1-pf(tau,nu_1,nu_2)                     # Überschreitungswahrscheinlichkeit

# Ausgabe
cat("Pillai's  V        "  , V,
    "\ntau                ", tau,
    "\nnu_1               ", nu_1,
    "\nnu_2               ", nu_2,
    "\nP(tau > tau_tilde) ", P)
```

#
\setstretch{2}
\vfill
\large
Anwendungsszenario und Datendeskription

Modellformulierung und Modellschätzung

Modellevaluation

* Wilks' $\Lambda$
* Pillai's $V$

**Selbstkontrollfragen**
\vfill

# Selbskontrolfragen
\setstretch{1.8}
\footnotesize
\begin{enumerate}
\item Erläutern Sie das Anwendungsszenario einer multivariaten einfaktoriellen Varianzanalyse.
\item Definieren Sie das Klassische Modell der einfaktoriellen Varianzanalyse.
\item Geben Sie das Theorem zur Kreuzproduktsummenmatrixzerlegung wieder.
\item Erläutern Sie die intuitive Bedeutung der $T,B$ und $W$ Matrizen der Kreuzproduktsummenmatrixzerlegung. 
\item Geben Sie einen Überblick zu den Teststatistiken der einfaktoriellen Varianzanalyse und ihren Verteilungen.
\item Definieren Sie die Wilk's $\Lambda$ Teststatistik.
\item Erläutern Sie die intuitive Bedeutung der Wilk's $\Lambda$ Teststatistik.
\item Definieren Sie die Pillai's $V$ Teststatistik.
\item Erläutern Sie die intuitive Bedeutung der Pillai's $V$ Teststatistik
\end{enumerate}

# References
\footnotesize