---
fontsize: 8pt
format:
    beamer:
        include-in-header: "6_Header.tex"
bibliography: 6_Referenzen.bib
---

#  {.plain}

\center
```{r, echo = FALSE, out.width = "20%", fig.align = "center"}
knitr::include_graphics("6_Abbildungen/mv_6_otto.png")
```
\vspace{2mm}

\huge
Multivariate Verfahren
\vspace{4mm}

\large
MSc Psychologie | MSc Klinische Psychologie und Psychotherapie   

WiSe 2023/24

\vspace{4mm}
\normalsize
Prof. Dr. Dirk Ostwald

#  {.plain}

\vfill
\center
\huge
\textcolor{black}{(6) Einfaktorielle Varianzanalyse}
\vfill

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

Modellformulierung und Modellschätzung

Modellevaluation mit der Wilks'-$\Lambda$-Statistik

Selbstkontrollfragen
\vfill

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

Modellformulierung und Modellschätzung

Modellevaluation mit der Wilks'-$\Lambda$-Statistik

Selbstkontrollfragen
\vfill


# Anwendungsszenario und Datendeskription
\setstretch{1.2}
Anwendungsszenario

\small
* Zwei oder mehr Gruppen experimenteller Einheiten mit Datendimension $m>1$ 

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

* $\mu_i \in \mathbb{R}^m, i = 1,...,p$ und $\Sigma \in \mathbb{R}^{m \times m} \mbox{ pd }$ 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
Anwendungsbeispiel

\small
* Analyse von Daten dreier Therapiegruppen (F2F, ONL, WLC) von `dBDI` und `dGLU` Daten

* Testen der Nullhypothese $\mu_1 = \mu_2 = \mu_3$

# Anwendungsszenario und Datendeskription
\footnotesize

```{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 Patient:innen. Die Tabelle zeigt exemplarisch die ersten fünf Datenpunkte jeder Gruppe."
library(knitr)                                                                  # knitr für Tabellen
D = read.csv("./6_Daten/6_Einfaktorielle_Varianzanalyse.csv")                   # Dateneinlesen
kable(D[c(1:5, 16:20, 31:35),], digits = 1, align = "c")                                               # Datentabelle
``` 

# Anwendungsszenario und Datendeskription
Anwendungsbeispiel

```{r, echo = F, warning = F, eval = F}
# Studiengruppenspezifische Deskriptivstatistiken
D       = read.csv("./6_Daten/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

# Visualisierung
library(ellipse)
library(latex2exp)
pdf(
file        = "6_Abbildungen/mv_6_anwendungsbeispiel.pdf",
width       = 5.5,
height      = 5.5)
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(16,16,16),
bg          = cols,
col         = cols,
bty         = "n",
cex         = 1,
x.intersp   = 1)
dev.off()
```

\center
```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics("6_Abbildungen/mv_6_anwendungsbeispiel.pdf")
```

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

**Modellformulierung und Modellschätzung**

Modellevaluation mit der Wilks'-$\Lambda$-Statistik

\large
Selbstkontrollfragen
\vfill

# Modellformulierung und Modellschätzung
\setstretch{1.2}
\footnotesize
\justifying
\begin{definition}[Modell der einfaktoriellen multivariaten Varianzanalyse]
Für $i = 1,...,p$ und $j = 1,...,n_i$ seien $\upsilon_{ij}$ $m$-dimensionale
Zufallsvektoren, die die $n := \sum_{i=1}^p n_i$ $m$-dimensionalen Datenpunkte
eines einfaktoriellen multivariaten Varianzanalyseszenarios modellieren. Dann hat
das Modell der einfaktoriellen multivariaten Varianzanalyse die strukturelle Form
\begin{equation}
\ups_{ij} = \mu_i + \varepsilon_{ij} 
\mbox{ mit } \varepsilon_{ij} \sim N(0_m,\Sigma) \mbox{ u.i.v. } 
\mbox{ mit } \mu_i \in \mathbb{R}^m \mbox{ und } \Sigma \in \mathbb{R}^{m \times m} \mbox{ pd }
\end{equation}
und die Datenverteilungsform
\begin{equation}
\ups_{ij} \sim N(\mu_i,\Sigma) \mbox{ u.v.  mit }  \mu_i \in \mathbb{R}^m  \mbox{ und } \Sigma \in \mathbb{R}^{m \times m} \mbox{ pd}.
\end{equation}
\end{definition}
\vspace{-2mm}
Bemerkungen

* Der Einfachheit halber setzen wir meist identische Gruppengrößen $k := n_i$ voraus.
* Die Gesamtheit aller Datenzufallsvektoren bezeichnen wir mit $\Upsilon$.

# Modellformulierung und Modellschätzung
\setstretch{1.2}
\footnotesize
\begin{theorem}[Parameterschätzer]
Gegeben sei das Modell der einfaktoriellen multivariaten Varianzanalyse. Dann
ist für $i = 1,...,p$
\begin{equation}
\hat{\mu}_i := \frac{1}{n_i}\sum_{j = 1}^{n_i} \ups_{ij}
\end{equation}
ein unverzerrte Schätzer des gruppenspezifischen Erwartungswertparameters $\mu_i$ und
\begin{equation}
\hat{\Sigma} := \frac{1}{n-p}\sum_{i = 1}^p \sum_{j = 1}^{n_i} \left(\ups_{ij} - \hat{\mu}_i\right)\left(\ups_{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 $1/(n-p)$ skalierte Within-Group Sum-of-Squares Matrix (siehe unten).
* Anstelle eines Beweis validieren wir die Aussage des Theorems mithilfe einer Simulation.

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

Modellformulierung und Modellschätzung

**Modellevaluation mit der Wilks'-$\Lambda$-Statistik**

Selbstkontrollfragen
\vfill

# Modellevaluation mit der Wilks'-$\Lambda$-Statistik
Überblick

\small
* Ziel einer einfaktoriellen multivariaten Varianzanalyse ist meist das Testen der Nullhypothese
\begin{equation}
H_0 : \mu_1 = \cdots = \mu_p.
\end{equation}

* Es ergibt sich damit die Alternativhypothese 
\begin{align}
\begin{split}
H_1 : \mu_{i_l} \neq \mu_{j_l}
& \mbox{ für mindestens ein Paar } i,j
\mbox{ mit } i \neq j, 1 \le i,j \le p
\\
& \mbox{ und mindestens ein } l
\mbox{ mit } 1 \le l \le m.
\end{split}
\end{align}

* Zur Evaluation von $H_0$ wurden eine Reihe von Teststatistiken vorgeschlagen.

* Alle Teststatistiken beruhen auf der Kreuzproduktsummenmatrixzerlegung.

* Wir betrachten hier exemplarisch die Wilks'-$\Lambda$-Statistik.

* Die Verteilungen der Wilks'-$\Lambda$-Statistik unter der Nullhypothese sind zum Teil analytisch angebbar,
zum Teil müssen sie approximiert werden und haben die Form von $f$-Verteilungen.



# Modellevaluation mit der Wilks'-$\Lambda$-Statistik
\small
\begin{theorem}[Kreuzproduktsummenmatrizenzerlegung]
\justifying
\normalfont
Für $i = 1,...,p$ und $j = 1,...,n_i$ bezeichne $\ups_{ij}$ den $j$ten Stichprobenvektor
der $i$ten Stichprobengruppe eines einfaktoriellen multivariaten Varianzanalysemodells. 
Weiterhin seien
\begin{equation}
\bar{\ups} := \frac{1}{n}\sum_{i=1}^p \sum_{j=1}^{n_i} \ups_{ij}
\mbox{ und }
\bar{\upsilon}_i := \frac{1}{n_i} \sum_{j=1}^{n_i} \ups_{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}^p \sum_{j=1}^{n_i}\left(\ups_{ij}-\bar{\ups}\right)\left(\ups_{ij}-\bar{\ups}\right)^T$
&
die \textit{Totale Sum-of-Squares Matrix}
\\
\\
$B := \sum_{i=1}^p n_i\left(\bar{\ups}_{i}-\bar{\ups}\right)\left(\bar{\ups}_{i}-\bar{\ups}\right)^T$
&
die \textit{Between-Group Sum-of-Squares Matrix}
\\
\\
$W := \sum_{i=1}^p \sum_{j=1}^{n_i}\left(\ups_{ij}-\bar{\upsilon}_i\right)\left(\ups_{ij}-\bar{\upsilon}_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 mit der Wilks'-$\Lambda$-Statistik
\small
\setstretch{1.8}
Bemerkungen

$T \in \mathbb{R}^{m \times m}$ misst die totale Variabilität der Datenvektoren um das Gesamtstichprobenmittel.
\vspace{1mm}

$B \in \mathbb{R}^{m \times m}$ misst die Variabilität der Gruppenstichprobenmittel um das Gesamtstichprobenmittel.
\vspace{1mm}

$W \in \mathbb{R}^{m \times m}$ misst die Variabilität der Datenvektoren um ihre jeweiligen Gruppenstichprobenmittel.
\vspace{1mm}

Die totale Variabilität wird hier also in zwei unabhängige Beiträge von Variabilität zerlegt.
\vspace{1mm}

$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 = (n-p)\hat{\Sigma}.
\end{equation}

Ein Beweis gibt sich durch algebraische Umformung.

# Modellevaluation mit der Wilks'-$\Lambda$-Statistik
\tiny
\setstretch{.8}
\begin{align}
\begin{split}
\mbox{\underline{Beweis} } \quad\quad
T
& = \sum_{i=1}^p \sum_{j=1}^{n_i}(\ups_{ij}-\bar{\ups})(\ups_{ij}-\bar{\ups})^T
\\
& = \sum_{i=1}^p \sum_{j=1}^{n_i} (\ups_{ij}-\bar{\ups}_i+\bar{\ups}_i-\bar{\ups})(\ups_{ij}-\bar{\ups}_i+\bar{\ups}_i-\bar{\ups})^T
\\
& = \sum_{i=1}^p \sum_{j=1}^{n_i} \left((\ups_{ij}-\bar{\ups}_i)+(\bar{\ups}_i-\bar{\ups})\right)\left((\ups_{ij}-\bar{\ups}_i)+(\bar{\ups}_i-\bar{\ups})\right)^T
\\
& = \sum_{i=1}^p \sum_{j=1}^{n_i} \left(
                                      (\ups_{ij}-\bar{\ups}_i)(\ups_{ij}-\bar{\ups}_i)^T
								    +2(\ups_{ij}-\bar{\ups}_i)(\bar{\ups}_i-\bar{\ups})^T
				    		        + (\bar{\ups}_i-\bar{\ups})(\bar{\ups}_i-\bar{\ups})^T
				  		          \right)
\\
& = \sum_{i=1}^p \left(
                       \sum_{j=1}^{n_i} (\ups_{ij}-\bar{\ups}_i)(\ups_{ij}-\bar{\ups}_i)^T
   				      +\sum_{j=1}^{n_i}2(\ups_{ij}-\bar{\ups}_i)(\bar{\ups}_i-\bar{\ups})^T
				      +\sum_{j=1}^{n_i}(\bar{\ups}_i-\bar{\ups})(\bar{\ups}_i-\bar{\ups})^T
				 \right)
\\
& = \sum_{i=1}^p \left(
                       \sum_{j=1}^{n_i} (\ups_{ij}-\bar{\ups}_i)(\ups_{ij}-\bar{\ups}_i)^T
   				      +2\left(\sum_{j=1}^{n_i}(\ups_{ij}-\bar{\ups}_i)\right)(\bar{\ups}_i-\bar{\ups})^T
				      +n_i(\bar{\ups}_i-\bar{\ups})(\bar{\ups}_i-\bar{\ups})^T
				 \right)
\\
& = \sum_{i=1}^p \left(
                       \sum_{j=1}^{n_i} (\ups_{ij}-\bar{\ups}_i)(\ups_{ij}-\bar{\ups}_i)^T
   				      +2\left(\sum_{j=1}^{n_i}\left(\ups_{ij}-\frac{1}{n_i}\sum_{j=1}^{n_i} \ups_{ij}\right)\right)(\bar{\ups}_i-\bar{\ups})^T
				      +n_i(\bar{\ups}_i-\bar{\ups})(\bar{\ups}_i-\bar{\ups})^T
				 \right)
\\
& = \sum_{i=1}^p \left(
                       \sum_{j=1}^{n_i} (\ups_{ij}-\bar{\ups}_i)(\ups_{ij}-\bar{\ups}_i)^T
   				      +2\left(\sum_{j=1}^{n_i} \ups_{ij}-\sum_{j=1}^{n_i} \ups_{ij}\right) (\bar{\ups}_i-\bar{\ups})^T
				      +n_i(\bar{\ups}_i-\bar{\ups})(\bar{\ups}_i-\bar{\ups})^T
				 \right)
\\
& = \sum_{i=1}^p \left(
                       \sum_{j=1}^{n_i} (\ups_{ij}-\bar{\ups}_i)(\ups_{ij}-\bar{\ups}_i)^T
   				       +n_i(\bar{\ups}_i-\bar{\ups})(\bar{\ups}_i-\bar{\ups})^T
				 \right)
\\
& = \sum_{i=1}^p n_i(\bar{\ups}_i-\bar{\ups})(\bar{\ups}_i-\bar{\ups})^T
   +\sum_{i=1}^p \sum_{j=1}^{n_i} (\ups_{ij}-\bar{\ups}_i)(\ups_{ij}-\bar{\ups}_i)^T
\\
& = B + W.
\end{split}
\end{align}

# Modellevaluation mit der Wilks'-$\Lambda$-Statistik 
\small
\begin{definition}[Wilks'-$\Lambda$-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 Wilks'-$\Lambda$-Statistik Teststatistik definiert als
\begin{equation}
\Lambda := \frac{|W|}{|B+W|},
\end{equation}
wobei $|\cdot|$ die Determinante bezeichnet.
\end{definition}

\footnotesize
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{\ups}_1 = \cdots = \bar{\ups}_p = \bar{\ups}$ gilt $B = 0_{mm}$ und damit $\Lambda = 1$.
* Für steigende Unterschiede zwischen den $\bar{\upsilon}_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
\footnotesize
\begin{theorem}[Spezielle $H_0$ Verteilungen von $\Lambda$ Transformationen]
\justifying
\normalfont
Es seien das Modell der einfaktoriellen Varianzanalyse und Wilks'-$\Lambda$-Statistik
definiert wie oben und es gelte außerdem
\begin{equation}
H_0 : \mu_1 = \cdots = \mu_p.
\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 $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{theorem}

Bemerkungen

* Die Verteilungen gehen zurück auf @wilks1932.


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

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

# Szenariensimulationen
library(MASS)                                      # R Paket für multivariate Normalverteilungen
nsc = length(M)                                    # Szenarienanzahl
T   = matrix(rep(NaN,nsm*nsc), ncol = nsc)         # Teststatistik Array
for(sc in 1:nsc){                                  # Szenarioiterationen

    # Modellparameter
    m     = M[sc]                                  # Datendimension
    p     = P[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,k,p))              # Datenanarrayinitialisierung
        for(i in 1:p){                             # Gruppeniterationen
            Y[,,i] =  t(mvrnorm(k,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 Teststatistik ("Prüfgröße")
        if     (sc == 1){t = ((1-L)/L)*((n-p-m+1)/m)}
        else if(sc == 2){t = ((1-sqrt(L))/sqrt(L))*((n-p-m+1)/m)}
        else if(sc == 3){t = ((1-L)/L)*((n-p)/(p-1))}
        else if(sc == 4){t = ((1-sqrt(L))/sqrt(L))*((n-p-1)/(p-1))}

        # Statistikrealisation
        T[sm,sc] = t}}
```

```{r, echo = F, eval = F}
# figure parameters
pdf(
file        = "./6_Abbildungen/mv_6_spezielle_lambda_verteilungen.pdf",
width       = 6,
height      = 6)
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]
    p     = P[sc]
    n     = N[sc]
    if(sc == 1){
        df = c(m,n-p-m+1)}
    else if(sc == 2){
        df = c(2*m, 2*(n-p-m+1))}
    else if(sc == 3){
        df = c(p-1, n-p)}
    else if(sc == 4){
        df = c(2*(p-1), 2*(n-p-1))}
    pdf         = df(f,df[1],df[2])
    hist(
    T[,sc],
    breaks      = 50,
    col         = "gray90",
    prob        = TRUE,
    xlim        = c(f_min, f_max),
    ylim        = ylimits[[sc]],
    xlab        = TeX("Statistik"),
    ylab        = "",
    main        = sprintf("m = %d, p = %d", M[sc], P[sc]))
    lines(
    f,
    pdf,
    lwd         = 2,
    col         = "darkorange")
}
dev.off()
```

\center
```{r, echo = FALSE, out.width = "62%"}
knitr::include_graphics("6_Abbildungen/mv_6_spezielle_lambda_verteilungen.pdf")
```

# Modellevaluation mit der Wilks'-$\Lambda$-Statistik
\footnotesize
\begin{theorem}[Approximative $H_0$ Verteilungen von $\Lambda$ Transformationen]
\normalfont
\justifying
Es seien das Modell der einfaktoriellen Varianzanalyse und Wilks'-$\Lambda$-Statistik
definiert wie oben und es gelte außerdem
\begin{equation}
H_0 : \mu_1 = \cdots = \mu_p.
\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(p-1)
\mbox{ und }
\nu_2 := wt-\frac{1}{2}(m(p-1)-2)
\end{equation}
sowie
\begin{equation}
w := n-1-\frac{1}{2}(m+k)
\mbox{ und }
t   := \sqrt{\frac{m^2(p-1)^2 - 4}{m^2 + (p-1)^2 - 5}}
\end{equation}
approximativ $f$-verteilt mit Freiheitsgradparametern $\nu_1$ und $\nu_2$.
\end{theorem}

Bemerkungen

* Die Approximation geht zurück auf @rao1951.

# 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
library(latex2exp)
pdf(
file        = "./6_Abbildungen/mv_6_tau_lambda.pdf",
width       = 6,
height      = 4)
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, p = 4, k = 15",
  "m = 3, p = 9, k = 15",
  "m = 4, p = 4, k = 15",
  "m = 9, p = 4, k = 15"),
lty = 1,
col = cols,
pch = NA,
bty = "n",
y.intersp = 1.2,
cex = 0.9)
dev.off()
```

\center
```{r, echo = FALSE, out.width = "80%"}
knitr::include_graphics("6_Abbildungen/mv_6_tau_lambda.pdf")
```

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

```{r, eval = F, echo = F}
# Szenarioparameter
library(MASS)                                          # R Paket für multivariate Normalverteilungen
nsm = 1e4                                              # Datensimulationsanzahl
M   = c(3,3,4,9)                                       # Datendimension
P   = c(4,9,4,4)                                       # Gruppenanzahl
k   = 15                                               # Datenpunkte pro Gruppe
N   = k*P                                              # 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
    p       = P[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+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
    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,k,p))                  # Datenanarrayinitialisierung
        for(i in 1:p){                                 # Gruppeniterationen
            Y[,,i] =  t(mvrnorm(k,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
```

```{r, echo = F, eval = F}
library(latex2exp)
pdf(
file        = "./6_Abbildungen/mv_6_approximative_lambda_verteilungen.pdf",
width       = 6,
height      = 6)
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)
    pdf   = 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, p = %d", M[sc], P[sc]))
    lines(
    f,
    pdf,
    lwd         = 2,
    col         = "darkorange")
}
dev.off()
```

\center
```{r, echo = FALSE, out.width = "62%"}
knitr::include_graphics("6_Abbildungen/mv_6_approximative_lambda_verteilungen.pdf")
```

# Modellevaluation mit der Wilks'-$\Lambda$-Statistik
\small
\begin{theorem}[Wilks'-$\Lambda$-Statistik-basierter Test]
\justifying
\normalfont
Es seien das Modell der einfaktoriellen Varianzanalyse und die Wilks'-$\Lambda$-Statistik
basierte Teststatistik $\tau$ mit Verteilungsparametern $\nu_1,\nu_2$ wie oben definiert.
Weiterhin sei der kritische Wert-basierte Test
\begin{equation}
\phi(\Upsilon) := 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{theorem}


\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, p = 4, k = 15 \Rightarrow \nu_1 = 9, \nu_2 = 132$ und $\alpha_0 = 0.05$.
\vspace{2mm}

```{r,echo = F, eval = F}
# Visualisierungsparameter
library(latex2exp)
pdf(
file        = "./6_Abbildungen/mv_6_kritischer_wert_lambda.pdf",
width       = 8,
height      = 4)
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")
dev.off()
```

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

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

Praktisches Vorgehen
\footnotesize

* \justifying Man nimmt an, dass ein vorliegender Datensatz von $p$ Gruppen von 
$m$-dimensionalen Datenvektoren für jeweils $j = 1,...,n_i$ Realisationen von 
$\upsilon_{ij} \sim N(\mu_i,\Sigma)$ mit unbekannten Parametern $\mu_i \in 
\mathbb{R}^m, i = 1,...,p$ und $\Sigma \in \mathbb{R}^{m \times m} \mbox{ pd }$ sind.

* Man möchte entscheiden ob $H_0 : \mu_1 = \cdots = \mu_p$ 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, p = 4, k = 15$ für $i = 1,...p$ 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,p$ und $k$ 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}

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

Modellformulierung und Modellschätzung

Modellevaluation mit der Wilks'-$\Lambda$-Statistik

**Selbstkontrollfragen**
\vfill

# Selbstkontrolfragen
\footnotesize
\setstretch{1.6}

1. Erläutern Sie das Anwendungsszenario einer einfaktoriellen multivariaten Varianzanalyse.
1. Geben Sie die Definition des Modells der einfaktoriellen multivariaten Varianzanalyse wieder.
1. Geben Sie das Theorem zu den Parameterschätzern der einfaktoriellen multivariaten Varianzanalyse wieder.
1. Erläutern Sie die Null- und Alternativhypothesen einer einfaktoriellen multivariaten Varianzanalyse.
1. Geben Sie das Theorem zur Kreuzproduktsummenmatrizenzerlegung wieder.
1. Was messen die Totale, Between-Group und die Within-Group Sum-of-Squares Matrizen, respektive?
1. Geben Sie die Definition der Wilks'-$\Lambda$-Statistik wieder.
1. Erläutern Sie Gemeinsamkeiten und Unterschiede zwischen speziellen und approximativen $H_0$ Verteilungen
von Wilks-$\Lambda$-Transformationen bei der einfaktoriellen multivariaten Varianzanalyse.
1. Geben Sie das Theorem zum Wilks-$\Lambda$-Statistik-basierten Test im Rahmen der einfaktoriellen multivariaten Varianzanalyse wieder.
1. Erläutern Sie das praktische Vorgehen zur Durchführung eines Wilks-$\Lambda$-Statistik-basierten Tes tim Rahmen der einfaktoriellen multivariaten Varianzanalyse.

# Referenzen
\footnotesize