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


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

#  {.plain}
\center
```{r, echo = FALSE, out.width = "20%"}
knitr::include_graphics("10_Abbildungen/mvda_10_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}{(10) T$^2$-Tests}
\vfill

#
\vfill
\setstretch{2.4}
\large
**Vorbemerkungen**

Einstichproben-T$^2$-Tests

Zweistichproben-T$^2$-Tests

Univariates vs. multivariates Testen

Selbstkontrollfragen
\vfill

# Vorbemerkungen
Ausgewählte multivariate Methoden der Frequentistischen Inferenz

\small
Multivariate Generalisierungen bekannter Frequentistischer Verfahren

**T$^2$-Tests** als Generalisierung von T-Tests

* Inferentieller Vergleich von ein bis zwei Gruppen multivariater Daten

Multivariate **Einfaktorielle Varianzanalyse** als Generalisierung der einfaktoriellen Varianzanalyse

* Inferentieller Vergleich von drei oder mehr Gruppen multivariater Daten

**Kanonische Korrelationsanalyse** als Generalisierung (multipler) Korrelation

* Korrelative Zusammenhänge von Zufallsvektoren

Zur Revision univariater Frequentistischer Verfahren

[\textcolor{red}{$\quad\bullet$ Inferenzstatistik SoSe 21}](https://bit.ly/3sTMP8U)

[\textcolor{red}{$\quad\bullet$ Wahrscheinlichkeitstheorie und Frequentistische Inferenz WiSe 21/22}](https://bit.ly/3qJKlan)

# Vorbemerkungen
```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("10_Abbildungen/mvda_10_frequentistische_inferenz.pdf")
```

# Vorbemerkungen
Standardprobleme Frequentistischer Inferenz

\small
\noindent (1) **Parameterschätzung**

Ziel der Parameterschätzung ist es, einen möglichst guten Tipp für den wahren,
aber unbekannten, Parameterwert (oder eine Funktion dessen) abzugeben,
typischerweise basierend auf der Beobachtung einer Realisierung von
$Y_1,...,Y_n \sim \mathbb{P}_\theta$.
\vspace{2mm}

\noindent (2) **Konfidenzintervalle**

Das Ziel der Bestimmung von Konfidenzintervallen ist es, basierend auf der
Verteilung möglicher Parameterschätzwerte eine quantitative Aussage über die
mit dem Schätzwert assoziierte Unsicherheit zu treffen.
\vspace{2mm}

\noindent (3) **Hypothesentests**

Das Ziel der Auswertung von Hypothesentests ist es, basierend auf der angenommenen
Verteilung der Beobachtungen $Y_1,...,Y_n$ in einer möglichst sinnvollen Form
zu entscheiden, ob der wahre, aber unbekannte Parameterwert, in einer von zwei
sich gegenseitig ausschließenden Untermengen des Parameterraumes, welche man als
Hypothesen bezeichnet, liegt.

# Vorbemerkungen
Standardannahmen Frequentistischer Inferenz

\footnotesize
$\mathcal{M}$ sei ein statistisches Modell mit unabhängig und identisch verteilten
Zufallsvektoren $Y_1,...,Y_n \sim \mathbb{P}_\theta$. Es wird angenommen, dass ein
vorliegender Datensatz eine der möglichen Realisierungen von $Y_1,...,Y_n \sim
\mathbb{P}_\theta$ ist. Aus frequentistischer Sicht kann man die Erhebung
von Datensätzen unendlich oft wiederholen und zu jedem Datensatz Statistiken
auswerten.

\footnotesize
\begin{itemize}
\item[] Datensatz (1) : $y^{(1)} = \left(y_1^{(1)}, y_2^{(1)}, ...,y_n^{(1)}\right)$, Statistik (1): $S : \mathbb{R}^{m\times n} \to \Sigma, y^{(1)} \mapsto S\left(y^{(1)}\right)$
\item[] Datensatz (2) : $y^{(2)} = \left(y_1^{(2)}, y_2^{(2)}, ...,y_n^{(2)}\right)$, Statistik (2): $S : \mathbb{R}^{m\times n} \to \Sigma, y^{(2)} \mapsto S\left(y^{(2)}\right)$
\item[] Datensatz (3) : $y^{(3)} = \left(y_1^{(3)}, y_2^{(3)}, ...,y_n^{(3)}\right)$, Statistik (3): $S : \mathbb{R}^{m\times n} \to \Sigma, y^{(3)} \mapsto S\left(y^{(3)}\right)$
\item[] Datensatz (4) : $y^{(4)} = \left(y_1^{(4)}, y_2^{(4)}, ...,y_n^{(4)}\right)$, Statistik (4): $S : \mathbb{R}^{m\times n} \to \Sigma, y^{(4)} \mapsto S\left(y^{(4)}\right)$
\item[] ...
\end{itemize}

Um die Qualität statistischer Methoden zu beurteilen betrachtet die frequentistische
Statistik deshalb die Wahrscheinlichkeitsverteilungen von Statistiken und Schätzern
unter der Annahme von $Y_1,...,Y_n \sim \mathbb{P}_\theta$.

Wenn eine statistische Methode im Sinne der frequentitischen Standardannahmen "gut" ist,
dann heißt das also, dass sie bei häufiger Anwendung "im Mittel gut" ist. Im Einzelfall,
also im realen Normalfall nur eines vorliegenden Datensatzes, kann sie auch "schlecht" sein.

# Vorbemerkungen
```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("10_Abbildungen/mvda_10_hypothesentests.pdf")
```

# Vorbemerkungen
\small
\begin{definition}[Mahalanobis Distanz]
\justifying
$\xi_1$ sei ein Zufallsvektor, eine Realisiation eines Zufallsvektors,
ein multivariater Erwartungswert oder ein multivariates Stichprobenmittel,
$\xi_2$ sei ein Zufallsvektor, eine Realisiation eines Zufallsvektors,
ein multivariater Erwartungswert oder ein multivariates Stichprobenmittel
und $\Xi$ sei eine Kovarianzmatrix oder eine Stichprobenkovarianzmatrix. Dann heißt
\begin{equation}
D = \left(\xi_1 - \xi_2 \right)^T\Xi^{-1}\left(\xi_1 - \xi_2\right)
\end{equation}
\textit{Mahalanobis Distanz} von $\xi_1$ und $\xi_2$ hinsichtlich $\Xi$.
\end{definition}

\footnotesize
Bemerkungen

* Eine Mahalanobis Distanz ist eine Kovarianzmatrix-normalisierte quadrierte Euklidische Distanz.
* Ähnliche Maße in der univariaten Statistik sind die $z$-Transformation $z = \frac{y - \mu}{\sigma}$ und Cohen's $d = \frac{\bar{y}_1 - \bar{y}_2}{s_{12}}$.
* Ähnlich wie bei $z$-Werten wird bei der Mahalanobis Distanz in "Einheiten von Kovarianzen" gemessen.
* Stark variante Komponenten von $\xi_1$ und $\xi_2$ tragen weniger zur Distanz bei.
* Stark kovariante Komponenten von $\xi_1$ und $\xi_2$ tragen weniger zur Distanz bei.

# Vorbemerkungen

Mahalanobis Distanzen als Funktion von Komponentenvarianzen
\vspace{2mm}

\footnotesize
\begin{equation*}
\quad
\Sigma := \begin{pmatrix} 1.0 & 0.0 \\ 0.0 & 1.0\end{pmatrix}
\quad\quad\quad\quad\quad\quad\quad
\Sigma := \begin{pmatrix} 0.5 & 0.0 \\ 0.0 & 0.5\end{pmatrix}
\quad\quad\quad\quad\quad\quad\quad
\Sigma := \begin{pmatrix} 1.5 & 0.0 \\ 0.0 & 1.5\end{pmatrix}
\end{equation*}
\vspace{2mm}

```{r, echo = F, eval = F}
library(ellipse)
library(matlib)
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,3),
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)

# Mahalanobis Distanzszenarien
xi          = matrix(c(-1,-1,1,1), nrow = 2)
Xi          = array(dim = c(2,2,3))
Xi[,,1]     = matrix(c(1.0, 0.0, 0.0,1.0), nrow = 2)
Xi[,,2]     = matrix(c(0.5, 0.0, 0.0,0.5), nrow = 2)
Xi[,,3]     = matrix(c(1.5, 0.0, 0.0,1.5), nrow = 2)
for(i in 1:3){
    plot(
    ellipse(as.matrix(Xi[,,i]), level = 0.40, centre = c(0,0)),
    type = "l",
    col = "Black",
    xlim = c(-2,2),
    ylim = c(-2,2),
    xlab = TeX("$\\xi_{i_1}$"),
    ylab = TeX("$\\xi_{i_2}$"),
    main = sprintf("D = %0.2f ", t(xi[,1] - xi[,2]) %*% inv(Xi[,,i]) %*% (xi[,1] - xi[,2])))
    points(
    xi[1,],
    xi[2,],
    col  = "Black",
    bg   = "Black",
    pch  = 21,
    cex  = 1.2)
}
dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_mahalanobis_distanzen_1.pdf"),
width  = 9,
height = 3)
```

```{r, echo = FALSE, out.width="100%"}
knitr::include_graphics("10_Abbildungen/mvda_10_mahalanobis_distanzen_1.pdf")
```


# Vorbemerkungen

Mahalanobis Distanzen als Funktion von Komponentenkovarianzen
\vspace{2mm}

\footnotesize
\begin{equation*}
\quad\quad\quad
\Sigma := \begin{pmatrix} 1.0 & 0.0 \\ 0.0 & 1.0\end{pmatrix}
\quad\quad\quad\quad\quad\quad
\Sigma := \begin{pmatrix} 1.0 & 0.9 \\ 0.9 & 1.0\end{pmatrix}
\quad\quad\quad\quad\quad\quad
\Sigma := \begin{pmatrix*}[r] 1.0 & -0.9 \\ -0.9 & 1.0\end{pmatrix*}
\end{equation*}
\vspace{2mm}

```{r, echo = F, eval = F}
library(ellipse)
library(matlib)
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,3),
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)

# Mahalanobis Distanzszenarien
xi          = matrix(c(-1,-1,1,1), nrow = 2)
Xi          = array(dim = c(2,2,3))
Xi[,,1]     = matrix(c(1.0, 0.0, 0.0,1.0), nrow = 2)
Xi[,,2]     = matrix(c(1.0, 0.9, 0.9,1.0), nrow = 2)
Xi[,,3]     = matrix(c(1.0,-0.9,-0.9,1.0), nrow = 2)
for(i in 1:3){
    plot(
    ellipse(as.matrix(Xi[,,i]), level = 0.40, centre = c(0,0)),
    type = "l",
    col = "Black",
    xlim = c(-2,2),
    ylim = c(-2,2),
    xlab = TeX("$\\xi_{i_1}$"),
    ylab = TeX("$\\xi_{i_2}$"),
    main = sprintf("D = %0.2f ", t(xi[,1] - xi[,2]) %*% inv(Xi[,,i]) %*% (xi[,1] - xi[,2])))
    points(
    xi[1,],
    xi[2,],
    col  = "Black",
    bg   = "Black",
    pch  = 21,
    cex  = 1.2)
}
dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_mahalanobis_distanzen_2.pdf"),
width  = 9,
height = 3)
```

```{r, echo = FALSE, out.width="100%"}
knitr::include_graphics("10_Abbildungen/mvda_10_mahalanobis_distanzen_2.pdf")
```

# Vorbemerkungen
\footnotesize
\begin{definition}[$f$-Zufallsvariable]
\justifying
$X$ sei eine Zufallsvariable mit Ergebnisraum $\mathbb{R}_{>0}$ und Wahrscheinlichkeitsdichtefunktion
\begin{equation}
p_X : \mathbb{R} \to \mathbb{R}_{>0}, x \mapsto p_X(x)
:= \nu_1^{\frac{\nu_1}{2}}\nu_2^{\frac{\nu_2}{2}}
   \frac{\Gamma\left(\frac{\nu_1+\nu_2}{2}\right)}{\Gamma\left(\frac{\nu_1}{2}\right)\Gamma\left(\frac{\nu_2}{2}\right)}
   \frac{x^{\frac{\nu_1}{2}-1}}{\left(\nu_1 x  + \nu_2 \right)^{\frac{\nu_1+\nu_2}{2}}},
\end{equation}
wobei $\Gamma$ die Gammafunktion bezeichne. Dann sagen wir, dass $X$ einer
$f$-Verteilung mit Freiheitsgradparametern $\nu_1$ und $\nu_2$ unterliegt und nennen $X$ eine
$f$-Zufallsvariable mit Freiheitsgradparametern $\nu_1$ und $\nu_2$. Wir kürzen dies
mit $X \sim f(\nu_1,\nu_2)$ ab. Die Wahrscheinlichkeitsdichtefunktion (WDF) einer
$f$-Zufallsvariable bezeichnen wir mit $f(x;\nu_1,\nu_2)$, die kumulative Verteilungsfunktion (KVF)
einer $f$-Zufallsvariable bezeichnen wir mit $F(x;\nu_1,\nu_2)$, und die inverse
kumulative Verteilungsfunktion einer $f$-Zufallsvariable bezeichnen wir mit $F^{-1}(x;\nu_1,\nu_2)$.
\end{definition}
Bemerkungen

* Im univariaten Fall ist die $F$-Statistik der Varianzanalyse bei Zutreffen der Nullhypothese $f$-verteilt
* Im multivariaten Fall ist z.B. die T$^2$-Statistik bei Zutreffen der Nullhypothese $f$-verteilt.

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

# x space
x_min   = 0
x_max   = 6
x_res   = 1e3
x       = seq(x_min, x_max, len = x_res)

# parameters of interest
nu_1    = c( 2,  2,  3,  3,  4,  4)
nu_2    = c(13, 26, 13, 26, 13, 26)

# plotting
matplot(x,
matrix(
c(
df(x,nu_1[1],nu_2[1]),
df(x,nu_1[2],nu_2[2]),
df(x,nu_1[3],nu_2[3]),
df(x,nu_1[4],nu_2[4]),
df(x,nu_1[5],nu_2[5]),
df(x,nu_1[6],nu_2[6])),
ncol = 6),
type         = "l",
lty          = c(1,2,1,2,1,2),
lwd          = c(2,1,2,1,2,1),
col          = c("gray20", "gray20", "gray50", "gray50", "gray80", "gray80"),
ylim         = c(0,1),
xlim         = c(x_min,x_max),
ylab         = " ",
xlab         = "x",
main         = TeX("$f(x;\\nu_1,\\nu_2)$"))
legend(
4.00,
1.00,
c(
TeX("$\\nu_1 = 2,  \\nu_2 = 13 "),
TeX("$\\nu_1 = 2,  \\nu_2 = 26"),
TeX("$\\nu_1 = 3,  \\nu_2 = 13 "),
TeX("$\\nu_1 = 3,  \\nu_2 = 26"),
TeX("$\\nu_1 = 4,  \\nu_2 = 13 "),
TeX("$\\nu_1 = 4,  \\nu_2 = 26")),
lty         = c(1,2,1,2,1,2),
lwd         = c(2,1,2,1,2,1),
col         = c("gray20", "gray20", "gray50", "gray50", "gray80",  "gray80"),
bty         = "n",
cex         = .8,
y.intersp   = 1.4)
dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_f_zentral_wdf.pdf"),
width  = 6,
height = 4)
```
\small
Wahrscheinlichkeitsdichtefunktionen von $f$-Verteilungen
\vspace{1mm}
```{r, echo = FALSE, out.width="80%"}
knitr::include_graphics("10_Abbildungen/mvda_10_f_zentral_wdf.pdf")
```

# Vorbemerkungen

\footnotesize
\begin{definition}[Nichtzentrale $f$-Zufallsvariable]
\justifying
$X$ sei eine Zufallsvariable mit Ergebnisraum $\mathbb{R}_{>0}$ und Wahrscheinlichkeitsdichtefunktion
\begin{multline}
p_X : \mathbb{R} \to \mathbb{R}_{>0}, x \mapsto \\
p_X(x)
:= \sum_{k=0}^\infty \frac{e^{-\delta/2}(\delta/2)^k}{\frac{\Gamma(\nu_2/2)\Gamma(\nu_1/2 + k)}{\Gamma(\nu_2/2 + \nu_1/2 + k)}k!}
    \left(\frac{\nu_1}{\nu_2}\right)^{\nu_1/2 + k}
    \left(\frac{\nu_2}{\nu_2+\nu_1x}\right)^{(\nu_1+\nu_2)/2 + k}
    x^{\nu_1/2 - 1 + k}
\end{multline}
wobei $\Gamma$ die Gammafunktion bezeichne. Dann sagen wir, dass $X$ einer
nichtzentralen $f$-Verteilung mit Nichtzentralitätsparameter $\delta$ und
Freiheitsgradparametern $\nu_1$ und $\nu_2$ unterliegt und nennen $X$ eine nichtzentrale
$f$-Zufallsvariable mit Nichtzentralitätsparameter $\delta$ und Freiheitsgradparametern
$\nu_1$ und $\nu_2$. Wir kürzen dies mit $X \sim f(\delta,\nu_1,\nu_2)$ ab.
Die Wahrscheinlichkeitsdichtefunktion (WDF) einer $f$-Zufallsvariable
bezeichnen wir mit $f(x;\delta,\nu_1,\nu_2)$, die kumulative Verteilungsfunktion (KVF)
einer nichtzentralen $f$-Zufallsvariable bezeichnen
wir mit $F(x;\delta,\nu_1,\nu_2)$, und die inverse kumulative Verteilungsfunktion
einer nichtzentralen $f$-Zufallsvariable
bezeichnen wir mit $F^{-1}(x;\delta,\nu_1,\nu_2)$.
\end{definition}
Bemerkungen

* Es gilt $f(0,\nu_1,\nu_2) = f(\nu_1,\nu_2)$.
* Im univariaten Fall ist die $F$-Statistik bei Nichtzutreffen der Nullhypothese nichtzentral $f$-verteilt
* Im multivariaten Fall ist z.B. die T$^2$-Statistik bei Nichtzutreffen der Nullhypothese nichtzentral $f$-verteilt.

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

# x space
x_min   = 0
x_max   = 6
x_res   = 1e3
x       = seq(x_min, x_max, len = x_res)

# parameters of interest
delta   = c(  0,  0,  4,  4,  8,  8)
nu_1    = c(  2,  2,  2,  2,  2,  2)
nu_2     = c(13, 26, 13, 26, 13, 26)

# plotting
matplot(x,
matrix(
c(
df(x,nu_1[1],nu_2[1],delta[1]),
df(x,nu_1[2],nu_2[2],delta[2]),
df(x,nu_1[3],nu_2[3],delta[3]),
df(x,nu_1[4],nu_2[4],delta[4]),
df(x,nu_1[5],nu_2[5],delta[5]),
df(x,nu_1[6],nu_2[6],delta[6])),
ncol = 6),
type         = "l",
lty          = c(1,2,1,2,1,2),
lwd          = c(2,1,2,1,2,1),
col          = c("gray20", "gray20", "gray50", "gray50", "gray80",  "gray80"),
ylim         = c(0,1),
xlim         = c(x_min,x_max),
ylab         = " ",
xlab         = "x",
main         = TeX("$f(x; \\delta, \\nu_1, \\nu_2)$"))
legend(
4.00,
1.00,
c(
TeX("$\\delta = 0, \\nu_1 = 2,  \\nu_2 = 13"),
TeX("$\\delta = 0, \\nu_1 = 2,  \\nu_2 = 26"),
TeX("$\\delta = 4, \\nu_1 = 2,  \\nu_2 = 13"),
TeX("$\\delta = 4, \\nu_1 = 2,  \\nu_2 = 26"),
TeX("$\\delta = 8, \\nu_1 = 2,  \\nu_2 = 13 "),
TeX("$\\delta = 8, \\nu_1 = 2,  \\nu_2 = 26")),
lty         = c(1,2,1,2,1,2),
lwd         = c(2,1,2,1,2,1),
col         = c("gray20", "gray20", "gray50", "gray50", "gray80",  "gray80"),
bty         = "n",
cex         = .8,
y.intersp   = 1.4)
dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_f_nichtzentral_wdf.pdf"),
width  = 6,
height = 4)
```

\small
WDFen von nichtzentralen $f$-Verteilungen
\vspace{1mm}
```{r, echo = FALSE, out.width="80%"}
knitr::include_graphics("10_Abbildungen/mvda_10_f_nichtzentral_wdf.pdf")
```

#
\vfill
\setstretch{2.4}
\large
Vorbemerkungen

**Einstichproben-T$^2$-Tests**

Zweistichproben-T$^2$-Tests

Univariates vs. multivariates Testen

Selbstkontrollfragen
\vfill

# Einstichproben-T$^2$-Tests

Anwendungszenario

\small
* **Eine Stichprobe experimenteller Einheiten**.

* Annahme unabhängiger und identisch nach $N(\mu,\Sigma)$ multivariat normalverteilter Daten.

* $\mu$ und $\Sigma$ unbekannt.

* Quantifizieren der Unsicherheit beim inferentiellen Vergleich von $\mu$ mit $\mu_0$ beabsichtigt.

\normalsize
Anwendungsbeispiele

\small
* Gruppenanalyse von BDI und Glukokortikoid Daten

$\quad\quad\quad$ $\circ$ $\mu \neq \mu_0$ als Evidenz für eine multivariate Abweichung von einem Normwert $\mu_0$.

* Gruppenanalyse von Kognitionstestdaten

$\quad\quad\quad$ $\circ$ $\mu \neq \mu_0$ als Evidenz für eine multivariate Abweichung von einem Normwert $\mu_0$.

# Einstichproben-T$^2$ Tests
Anwendungsbeispiel

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

# Datepräprozessierung
m           = 2                                       # Datendimension von Interesse
n           = 15                                      # Anzahl Datenpunkte pro Gruppe
Y           = cbind(1:n,                              # VP Index
                    D$X1[D$Gruppe == "ungenügend"],   # Y_{1j_1} Intelligenztestscore
                    D$X2[D$Gruppe == "ungenügend"])   # Y_{1j_2} Mathematiktestscore
colnames(Y) = c("VP",
                "IQ Test Score",
                "Math Test Score")                    # Variablennamen
knitr::kable(Y, "pipe", align = "c")                  # Tabelle
```
Abweichung des (IQ Test Score, Math Test Score) Erwartungswertparameters vom Normwert $\mu_0 = \begin{pmatrix} 60 \\ 60 \end{pmatrix}$?


# Einstichproben-T$^2$ Tests
Anwendungsbeispiel
\vspace{3mm}

\tiny
\setstretch{1}
```{r, message = F, warning = F}
# Deskriptivstatistik
library(foreign)                                              # R Paket
library(ellipse)                                              # R Paket
library(matlib)                                               # R Paket
mu_0   = matrix(c(60,60), nrow = 2)                           # Normwert
fname  = file.path(getwd(), "10_Daten", "studienerfolg.sav")  # Dateinamen
D      = read.spss(fname, to.data.frame = T)                  # Dateneinlesen
Y      = rbind(D$X1[D$Gruppe == "ungenügend"],                # Y_{i_1} IQ Test Score
               D$X2[D$Gruppe == "ungenügend"])                # Y_{i_2} Math Test Score
n      = ncol(Y)                                              # Anzahl Datenpunkte
j_n    = matrix(rep(1,n), nrow = n)                           # 1_n
I_n    = diag(n)                                              # I_n
J_n    = matrix(rep(1,n^2), nrow = n)                         # 1_{nn}
Y_bar  = (1/n)*(Y %*% j_n)                                    # Stichprobenmittel (cf.(3) Wahrscheinlichkeitstheorie)
C      = (1/(n-1))*(Y %*% (I_n-(1/n)*J_n) %*% t(Y))           # Stichprobenkovarianzmatrix (cf. ibid.)
D      = t(Y_bar - mu_0) %*% inv(C) %*% (Y_bar - mu_0)        # Mahalanobis Distanz

# Ausgabe
cat("Y_bar =", Y_bar,
    "\nD     =", D)
```

# Einstichproben-T$^2$ Tests
Anwendungsbeispiel
\vspace{1mm}


```{r, echo = F, eval = F}
library(ellipse)
library(matlib)
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,
cex.main    = 1)

# Gruppe 1
plot(
Y[1,],
Y[2,],
col  = "Gray70",
bg   = "Gray70",
pch  = 21,
xlab = TeX("IQ Test Score"),
ylab = TeX("Math Test Score"),
xlim = c(0,100),
ylim = c(0,100))
points(
Y_bar[1],
Y_bar[2],
col  = "White",
bg   = "gray70",
pch  = 24)
lines(
ellipse(C, level = 0.40, centre = Y_bar),
type = "l",
col = "Gray70")

# Normwert
points(
mu_0[1],
mu_0[2],
col  = "White",
bg   = "Gray30",
pch  = 24)

# Legende
legend(
"topleft",
c("Datenpunkte", "Stichprobenmittel", "Stichprobenkovarianz", "Normwert"),
lty = c(NA,NA,1,NA),
pch = c(19,17,NA,17),
col = c("gray70","gray70","gray70","gray30"),
bty = "n",
cex = .8)

# Print
dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_einstichproben_t2_deskription.pdf"),
width  = 4.5,
height = 4.5)
```

```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics("10_Abbildungen/mvda_10_einstichproben_t2_deskription.pdf")
```

# Einstichproben-T$^2$-Tests
Hypothesenszenarien
\small
\setstretch{1.1}

Einfache Nullhypothese, einfache Alternativhypothese $H_0:\mu = \mu_0, H_1:\mu = \mu_1$
\begin{itemize}
\item Theoretisch wichtiges Szenario (Neymann-Pearson Lemma)
\item Praktische Relevanz eher gering
\end{itemize}

Einfache Nullhypothese, zusammengesetzte Alternativhypothese $H_0:\mu = \mu_0, H_1:\mu \neq \mu_0$
\begin{itemize}
\item Zweiseitiger Einstichproben-T-Test mit ungerichteter Hypothese
\item Ungerichtete Fragestellung nach einem Unterschied
\end{itemize}

Zusammengesetzte Nullhypothese/Alternativhypothese $H_0:\mu \le \mu_0, H_1:\mu > \mu_0$
\begin{itemize}
\item Einseitiger Einstichproben-T-Test mit gerichteter Hypothese
\item Gerichtete Fragestellung nach einem positiven Unterschied
\end{itemize}

Zusammengesetzte Nullhypothese/Alternativhypothese $H_0:\mu\ge\mu_0,H_1:\mu<\mu_0$
\begin{itemize}
\item Gerichtete Fragestellung nach einem negativen Unterschied
\item Qualitativ äquivalente Theorie zum umgekehrten Fall
\end{itemize}

# Einstichproben-T$^2$-Tests
Im Folgenden näher betrachtetes Hypothesenszenario
\small
\setstretch{1.1}

\textcolor{lightgray}{Einfache Nullhypothese, einfache Alternativhypothese $H_0:\mu = \mu_0, H_1:\mu = \mu_1$}
\begin{itemize}
\item[\textcolor{lightgray}{$\bullet$}] \textcolor{lightgray}{Theoretisch wichtiges Szenario (Neymann-Pearson Lemma)}
\item[\textcolor{lightgray}{$\bullet$}] \textcolor{lightgray}{Praktische Relevanz eher gering}
\end{itemize}
Einfache Nullhypothese, zusammengesetzte Alternativhypothese $H_0:\mu = \mu_0, H_1:\mu \neq \mu_0$
\begin{itemize}
\item Zweiseitiger Einstichproben-T-Test mit ungerichteter Hypothese
\item Ungerichtete Fragestellung nach einem Unterschied
\end{itemize}
\textcolor{lightgray}{Zusammengesetzte Nullhypothese/Alternativhypothese $H_0:\mu \le \mu_0, H_1:\mu > \mu_0$}
\begin{itemize}
\item[\textcolor{lightgray}{$\bullet$}] \textcolor{lightgray}{Einseitiger Einstichproben-T-Test mit gerichteter Hypothese}
\item[\textcolor{lightgray}{$\bullet$}] \textcolor{lightgray}{Gerichtete Fragestellung nach einem positiven Unterschied}
\end{itemize}
\textcolor{lightgray}{Zusammengesetzte Nullhypothese/Alternativhypothese $H_0:\mu\ge\mu_0,H_1:\mu<\mu_0$}
\begin{itemize}
\item[\textcolor{lightgray}{$\bullet$}] \textcolor{lightgray}{Gerichtete Fragestellung nach einem negativen Unterschied}
\item[\textcolor{lightgray}{$\bullet$}] \textcolor{lightgray}{Qualitativ äquivalente Theorie zum umgekehrten Fall}
\end{itemize}

# Einstichproben-T$^2$-Tests
\setstretch{2}
Gliederung (vgl. [\textcolor{red}{(12) Hypothesentests}](https://bit.ly/32ZUAiy))

\small
(1) Statistisches Modell in klassischer Form
(2) Statistisches Modell in generativer Form
(3) Testhypothesen, Teststatistik, Test
(4) Analyse der Teststatistik
(5) Analyse der Testgütefunktion
(6) Testumfangkontrolle
(7) p-Werte
(8) Analyse der Powerfunktion

\normalsize
Zur Wiederholung des univariaten Falls, siehe [\textcolor{red}{(13) Einstichproben-T-Tests}](https://bit.ly/3n0f0zb).

# Einstichproben-T$^2$-Tests
\noindent (1) Statistisches Modell in klassischer Form

\small
$Y_1,...,Y_n \sim N(\mu,\Sigma)$ sei die Stichproben eines $m$-dimensionalen
Normalverteilungsmodells mit unbekanntem Erwartungswertparameter
$\mu \in \mathbb{R}^m$ und unbekanntem Kovarianzmatrixparameter
$\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d}$. Als Parameter von Interesse
betrachten wir $\theta = \mu$, so dass sich der Parameterraum von Interesse zu
$\Theta = \mathbb{R}^m$ ergibt.

\vspace{5mm}
\normalsize
\noindent (2) Statistisches Modell in generativer Form

\small
Es sei
\begin{equation}
Y_i = \mu + \varepsilon_i \mbox{ mit } \varepsilon_i \sim N(0_m,\Sigma) \mbox{ für } i = 1,...,n
\end{equation}
wobei

* $Y_i, i = 1,...,m$ beobachtbare Zufallsvektoren,
* $\mu \in \mathbb{R}^m$ den festen und identischen Erwartungswertparameter über Zufallsvektoren und
* $\varepsilon_i, i = 1,...,n$ unabhängige normalverteilte nicht beobachtbare Zufallsvektoren

bezeichnen.

# Einstichproben-T$^2$-Tests
\noindent (2) Statistisches Modell in generativer Form (fortgeführt)

\small
Die generative Form betont, dass im vorliegenden Modell beobachtete Daten durch einen
systematischen deterministischen Prozess (hier $\mu \in \mathbb{R}^m$) unter dem
additiven Einfluss einer Vielzahl unabhängiger und deshalb in ihrer Summe normalverteilter
Störprozesse (hier in der Summe $\varepsilon_i, i = 1,...,n$) erzeugt konzipiert
werden.

\tiny
\underline{Beweis}

Wir halten zunächst fest, dass wenn $\xi \sim N(\alpha,S)$ ein $m_1$-dimensionaler
normalverteilter Zufallsvektor mit Erwartungswertparameter $\alpha \in \mathbb{R}^{m_1}$
und Kovarianzmatrixparameter $S \in \mathbb{R}^{m_1\times m_1} \mbox{ p.d. }$ ist
und für $A \in \mathbb{R}^{m_2 \times m_1}$ und $b \in \mathbb{R}^{m_2}$ ein
$m_2$-dimensionaler Zufallsvektor definiert ist als
\begin{equation}
\zeta := A\xi + b,
\end{equation}
dann gilt, dass
\begin{equation}
\zeta \sim N\left(A\alpha + b, ASA^T\right)
\end{equation}
(vgl. @anderson_2003, Section 2.4).  Aus $\varepsilon \sim N(0_m,\Sigma)$ folgt
hier mit
\begin{equation}
Y_i = I_m\varepsilon + \mu
\end{equation}
dann aber sofort, dass
\begin{equation}
Y \sim N\left(I_m0_m + \mu,I_m\Sigma I_m^T \right) = N\left(\mu,\Sigma\right).
\end{equation}
$\hfill \Box$


# Einstichproben-T$^2$-Tests
\noindent (3) Testhypothesen, Teststatistik, Test

\footnotesize
Für ein $\mu_0 \in \mathbb{R}^m$ betrachten wir die einfache Nullhypothese und
die zusammengesetzte Alternativhypothese
\begin{equation}
H_0 := \mu = \mu_0 \Leftrightarrow \Theta_0 :=  \{\mu_0\}
\mbox{ und }
H_0 := \mu \neq \mu_0 \Leftrightarrow \Theta_1 :=  \mathbb{R}^m \setminus \{\mu_0\}
\end{equation}
Weiterhin betrachten wir die Einstichproben-T$^2$-Teststatistik
\begin{equation}
T^2 := n(\bar{Y} - \mu_0)^T C^{-1}(\bar{Y} - \mu_0)
\end{equation}
wobei $\bar{Y}$ und $C$ das Stichprobenmittel und die Stichprobenkovarianzmatrix
der $Y_1,...,Y_n$, respektive, bezeichnen (vgl. (3) Wahrscheinlichkeitstheorie).

* $T^2$ ist die mit dem Stichprobenumfang skalierte Mahalanobis Distanz von $\bar{Y}$ und $\mu_0$ hinsichtlich $C$.
* $T^2 \uparrow$ für  $||\bar{Y} - \mu_0|| \uparrow$, $C \downarrow$ und $n\uparrow$.

Schließlich
definieren wir den kritischen Wert-basierten Test
\begin{equation}
\phi(Y) := 1_{\{T^2 > k\}} := \begin{cases} 1 & T^2  > k \\ 0 & T^2 \le k \end{cases}.
\end{equation}
wobei wie üblich 1 den Vorgang des Ablehnens von $H_0$ und 0 den Vorgang des Nichtablehnens
von $H_0$ repräsentieren.


# Einstichproben-T$^2$-Tests
\noindent (4) Analyse der Teststatistik

\footnotesize
\begin{theorem}[Verteilung der Einstichproben-T$^2$-Teststatistik]
\justifying
\normalfont
Es seien $Y_1,...,Y_n \sim N(\mu,\Sigma)$ mit $\mu \in \mathbb{R}^m$ und $\Sigma \in \mathbb{R}^{m\times m} \mbox{p.d.}$,
\begin{equation}
\nu:= \frac{n-m}{(n-1)m}
\end{equation}
und für $\mu \in \mathbb{R}^m$ sei die Einstichproben-T$^2$-Teststatistik definiert als
\begin{equation}
T^2 := n(\bar{Y} - \mu_0)^T C^{-1} (\bar{Y} - \mu_0).
\end{equation}
Dann gilt
\begin{equation}\label{eq:nuT2}
\nu T^2 \sim f(\delta, m, n-m)
\end{equation}
wobei $f(\delta,m,n-m)$ die nichtzentrale $f$-Verteilung mit Nichtzentralitätsparameter
\begin{equation}
\delta := n(\mu - \mu_0)^T\Sigma^{-1}(\mu - \mu_0)
\end{equation}
sowie mit Freiheitsgradparametern $m$ und $n-m$ bezeichnet.
\end{theorem}
\footnotesize
Bemerkungen

* Für einen Beweis von \eqref{eq:nuT2} verweisen wir auf @anderson_2003 und @hotelling_1931.
* Für $\mu = \mu_0$ und damit $\delta = 0$ entspricht $f(m,n-m,\delta)$ der $f$-Verteilung  $f(m,n-m)$

# Einstichproben-T$^2$-Tests

\footnotesize
\begin{theorem}[WDF und KVF der Einstichproben-$T^2$-Teststatistik]
\normalfont
\justifying
Im  Einstichproben-$T^2$-Testszenario sei
\begin{equation}
\nu:= \frac{n-m}{(n-1)m}
\end{equation}
Dann ist eine WDF der Einstichproben-$T^2$-Teststatistik gegeben durch
\begin{equation}\label{eq:pT1}
p_{T^2} : \mathbb{R}_{\ge 0} \to \mathbb{R},
t^2 \mapsto p_{T^2}(t^2) := \nu f(\nu t^2;\delta, m,n-m)
\end{equation}
und eine WDF der Einstichproben-$T^2$-Teststatistik ist gegeben durch
\begin{equation}\label{eq:PT2}
P_{T^2} : \mathbb{R}_{\ge 0} \to [0,1],
t^2 \mapsto P_{T^2}(t^2) := F(\nu t^2;\delta, m,n-m)
\end{equation}
\end{theorem}

Bemerkungen

* $\nu T^2$ hat die WDF $f(\delta,m,n-m)$, $T^2$ dagegen hat die WDF $\nu f(\nu t^2;\delta, m,n-m)$.
* Für $m := 1$ ist $\nu = (n-1)/(n-1)\cdot 1 = 1$ und mit der Stichprobenvarianz $S^2$ gilt
\begin{equation}
T^2 = n\frac{(\bar{Y} - \mu_0)^2}{S^2} = \left(\sqrt{n}\frac{\bar{Y} - \mu_0}{S} \right)^2
\end{equation}
* Das Quadrat der univariate Einstichproben-T-Teststatistik $T := \sqrt{n}\frac{\bar{Y}-\mu_0}{S}$ ist also $f(\delta,1,n-1)$ verteilt.


# Einstichproben-T$^2$-Tests
\footnotesize
\underline{Beweis}

Wir halten zunächst fest, dass das Theorem zur univariate WDF Transformation bei
linear-affinen Abbildungen (vgl. [(7) Transformationen der Normalverteilung](https://bit.ly/3zmkSYJ))
besagt, dass für eine Zufallsvariable $X$ mit WDF $p_X$ und der Definition
$Y = f(X)$ mit  $f(X) := aX + b$  für  $a\neq 0$ eine WDF von $Y$ definiert ist
$p_Y(y) := (1/|a|)p_X((y-b)/a)$. Im vorliegenden Fall ist $X = \nu T^2$ mit WDF
$f(\delta,m,n-m)$ und $Y := T^2 = \frac{1}{\nu}\nu T^2$, also $a = 1/\nu$ und $b = 0$.
Mit $\nu > 0$ ergibt sich \eqref{eq:pT1} also aus
\begin{equation}
p_{T^2}(t^2) = \frac{1}{a}p_X\left(\frac{t^2}{a}\right) = \nu f(\nu t^2; m, n-m)
\end{equation}
\eqref{eq:PT2} folgt dann mit der Tatsache, dass WDFen bei kontinuierlichen
Zufallsvariablen die Ableitungen der entsprechenden KVF sind, sowie der Kettenregel
der Differentiation
\begin{align}
\begin{split}
\frac{d}{dt^2}P_{T^2}\left(t^2\right)
& = \frac{d}{dt^2}\left(F(\nu t^2;m,n-m,\delta)\right) \\
& = \frac{d}{dt^2}F(\nu t^2;m,n-m,\delta)\frac{d}{dt^2}\left(\nu t^2 \right) \\
& = \nu f(\nu t^2;m,n-m,\delta) \\
& = p_{T^2}(t^2)
\end{split}
\end{align}
$\hfill\Box$

# Einstichproben-T$^2$-Tests
\noindent (4) Analyse der Teststatistik
\vspace{3mm}


\tiny
\setstretch{1}
```{r}
# Modellparameter
m     = 2                                                     # Dimensionalität der Zufallsvektoren/Daten
n     = 15                                                    # Anzahl der Datenpunkte
mu_0  = matrix(c(1,1) , nrow = 2)                             # H0 Hypothesenparameter
mu    = matrix(c(2,2) , nrow = 2)                             # wahrer, aber unbekannter, Erwartungswertparameter
Sigma = matrix(c(0.5,0.3, 0.3,0.5), nrow = 2, byrow = TRUE)   # wahrer, aber unbekannter, Kovarianzmatrixparameter

# Simulation
library(MASS)                                                 # R Paket für multivariate Normalverteilungen
library(matlib)                                               # R Paket für Matrizenrechnung
nsim  = 1e4                                                   # Anzahl Simulationen/Datensatzrealisierungen
Yb    = matrix(rep(NaN,m*nsim), nrow = 2)                     # Stichprobenmittelarray
T2    = rep(NaN,nsim)                                         # T2 Statistik Array
j_n   = matrix(rep(1,n), nrow = n)                            # 1_n
I_n   = diag(n)                                               # I_n
J_n   = matrix(rep(1,n^2), nrow = n)                          # 1_{nn}
for(s in 1:nsim){                                             # Simulationsiterationen
    Y      = t(mvrnorm(n,mu,Sigma))                           # Y_i \sim N(\mu,\Sigma), i = 1,...,n
    Y_bar  = (1/n)*(Y %*% j_n)                                # Stichprobenmittel
    C      = (1/(n-1))*(Y %*% (I_n-(1/n)*J_n) %*% t(Y))       # Stichprobenkovarianzmatrix
    T2[s]  = n*t(Y_bar - mu_0) %*% inv(C) %*% (Y_bar - mu_0)  # T2 Statistik
    Yb[,s] = Y_bar                                            # Stichprobenmittel für Visualisierung
}
```


# Einstichproben-T$^2$-Tests
\vspace{1mm}
\noindent (4) Analyse der Teststatistik
\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)

# Stichprobenmittel per Simulation
plot(
Yb[1,],
Yb[2,],
xlim     = c(1,3),
ylim     = c(1,3),
xlab     = TeX("$\\bar{y}_1$"),
ylab     = TeX("$\\bar{y}_2$"),
pch      = 21,
col      = "white",
bg       = "gray60",
main     = TeX("$\\bar{Y}$ Simulationen"))


# Unnormalisierte T2-Teststatistik WDF
t2_min      = 0
t2_max      = 200
t2_res      = 1e3
t2          = seq(t2_min,t2_max,len = t2_res)
nu          = (n-m)/(m*(n-1))
delta       = n*t(mu - mu_0) %*% inv(Sigma) %*% (mu - mu_0)
p_T2        = nu*df(nu*t2,m,n-m,delta)
hist(
T2,
breaks      = 100,
col         = "gray90",
prob        = TRUE,
xlim        = c(t2_min, t2_max),
ylim        = c(0,0.03),
xlab        = TeX("$t^2$"),
ylab        = "",
main        = TeX("$p(t^2) := \\nu f(\\nu t^2;\\delta, m, n-m)$"))
fdir        =  file.path(getwd(), "10_Abbildungen")
lines(
t2,
p_T2,
lwd         = 2,
col         = "darkorange")

# Normalisierte T2-Teststatistik Sample
t2_min      = 0
t2_max      = 50
t2_res      = 1e3
t2          = seq(t2_min,t2_max,len = t2_res)
delta       = n*t(mu-mu_0) %*% inv(Sigma) %*%(mu-mu_0)
p_nm_T2     = df(t2,m,n-m,delta)
hist(
((n-m)/m)/(n-1)*T2,
breaks      = 100,
col         = "gray90",
prob        = TRUE,
xlim        = c(t2_min, t2_max),
ylim        = c(0,0.06),
xlab        = TeX("$t^2$"),
ylab        = "",
main        = TeX("$\\nu T^2 ~ f(m,n-m, \\delta)$") )
lines(
t2,
p_nm_T2,
lwd         = 2,
col         = "darkorange")

# Unnormalisierte T2-Teststatistik CDF
t2_min      = 0
t2_max      = 200
t2_res      = 1e3
t2          = seq(t2_min,t2_max,len = t2_res)
nu          = (n-m)/(m*(n-1))
P_T2        = pf(nu*t2,m,n-m,delta)
ECDF_T2     = ecdf(T2)
plot(
ECDF_T2,
verticals   = TRUE,
lwd         = 6,
xlim        = c(t2_min, t2_max),
col         = "gray90",
xlab        = TeX("$t^2$"),
ylab        = "",
main        = TeX("$P(t^2) := F(\\nu t^2; \\delta, m, n-m)$"))
lines(
t2,
P_T2,
lwd         = 2,
col         = "darkorange")
dev.copy2pdf(
file        = file.path(fdir, "mvda_10_einstichproben_t2_teststatistik.pdf"),
width       = 6,
height      = 6)
```


```{r, echo = FALSE, out.width = "60%"}
knitr::include_graphics("10_Abbildungen/mvda_10_einstichproben_t2_teststatistik.pdf")
```

# Einstichproben-T$^2$-Tests
\noindent (5) Analyse der Testgütefunktion

\footnotesize
\begin{theorem}[Testgütefunktion]
\justifying
\normalfont
$\phi$ sei der im obigen Testszenario definiert Test. Dann ist die Testgütefunktion
von $\phi$ gegeben durch
\begin{equation}
q_\phi : \mathbb{R}^m \to [0,1], \mu \mapsto q_\phi(\mu) := 1 - F(\nu k;\delta_\mu,m,n-m)
\end{equation}
wobei $F(\cdot;\delta_\mu, m,n-m)$ die KVF der nichtzentralen $f$-Verteilung mit
Freiheitsgradparametern$m$ und $n-m$ sowie mit Nichtzentralitätsparameter
\begin{equation}
\delta_\mu := n(\mu - \mu_0)^T\Sigma^{-1}(\mu - \mu_0)
\end{equation}
bezeichnet.
\end{theorem}

Bemerkungen

* $q_\phi$ kann zur Bestimmung kritischer Werte für einen erwünschten Testumfang genutzt werden.
* $q_\phi$ kann zur Bestimmung der Testpower genutzt werden.

# Einstichproben-T$^2$-Tests
\noindent \noindent (5) Analyse der Testgütefunktion

\footnotesize
\underline{Beweis}

Die Testgütefunktion des betrachteten Tests im vorliegenden Testszenario ist definiert als
\begin{equation}
q_\phi : \mathbb{R}^m \to [0,1], \mu \mapsto q_\phi(\mu) := \mathbb{P}_{\mu}(\phi = 1)
\end{equation}
Da die Wahrscheinlichkeiten für $\phi = 1$ und dafür, dass die zugehörige Teststatistik
im Ablehnungsbereich des Tests liegt, gleich sind, benötigen wir also zunächst die Verteilung
der Teststatistik. Wir haben oben aber bereits gesehen, dass
\begin{equation}
\frac{n-m}{m(n-1)}T^2 \sim f(m,n-m,\delta_\mu) \mbox{ mit } \delta_\mu := n(\mu -\mu_0)^T\Sigma^{-1}(\mu-\mu_0)
\end{equation}
gilt. Der Ablehnungsbereich des betrachteten Tests ist $A := ]k,\infty[$. Also ergibt sich
\begin{align}
\begin{split}
q_\phi(\mu)
&  = \mathbb{P}_\mu(\phi = 1) \\
&  = \mathbb{P}_\mu\left(T^2 \in \,\,]k,\infty[\right) \\
&  = \mathbb{P}_\mu\left(T^2 > k \right) \\
&  = 1 - \mathbb{P}_\mu\left(T^2 \le k \right) \\
&  = 1 - F(\nu k; \delta_\mu,m, m-n)
\end{split}
\end{align}

# Einstichproben-T$^2$-Tests
\noindent (5) Analyse der Testgütefunktion

\small
Beispiele
\footnotesize
\begin{equation*}
m := 2,
n := 15,
\Sigma := \begin{pmatrix} 1 & 0  \\ 0 & 1 \end{pmatrix},
\mu_0 := \begin{pmatrix} 1 \\ 1 \end{pmatrix}
\end{equation*}
\vspace{1mm}

\tiny
\setstretch{1}
```{r, eval = F}
# Modellparameter
library(matlib)                                                       # R Paket
m           = 2                                                       # m
n           = 15                                                      # n
nu          = (n-m)/((n-1)*m)                                         # \nu
Sigma       = diag(m)                                                 # \Sigma = I_2
iSigma      = inv(Sigma)                                              # \Sigma^{-1}

# Testparameter
mu_0        = matrix(c(1,1), nrow = 2)                                # \mu_0
k_all       = c(2,4,6)                                                # k <-> \phi
n_k         = length(k_all)                                           # Anzahl k Werte/Tests

# q_\phi(\mu) Evaluation
mu_min      = 0                                                       # \mu_i Minimum
mu_max      = 2                                                       # \mu_i Maximum
mu_res      = 1e3                                                     # \mu_i Auflösung
mu_i        = seq(mu_min,mu_max,len = mu_res)                         # mu_i
q_phi       = array(dim = c(mu_res, mu_res, length(k_all)))           # q_\phi Array
for(k in 1:n_k){
    for(i in 1:mu_res){
        for(j in 1:mu_res){
            mu           = matrix(c(mu_i[i],mu_i[j]), nrow = 2)       # \mu
            delta_mu     = n*t(mu - mu_0) %*% iSigma %*% (mu -mu_0)   # \delta_\mu
            q_phi[i,j,k] = 1 - pf(nu*k_all[k], m, n-m, delta_mu)}}}   # q_\phi(\mu)
```

# Einstichproben-T$^2$-Tests
\noindent (5) Analyse der Testgütefunktion

\small
Beispiele
\footnotesize
\begin{equation*}
m := 2,
n := 15,
\Sigma := \begin{pmatrix} 1.0 & 0.0  \\ 0.0 & 1.0 \end{pmatrix},
\mu_0  := \begin{pmatrix} 1.0 \\ 1.0 \end{pmatrix}
\end{equation*}
\vspace{2mm}
```{r, echo = F, eval = F}
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,3),
mgp         = c(2,1,0),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1.2,
cex.main    = 1)
cols        = hcl.colors(11, "YlOrRd", rev = T)
labs        = list(TeX("$q_{\\phi}(\\mu), k = 2$"),
                   TeX("$q_{\\phi}(\\mu), k = 4$"),
                   TeX("$q_{\\phi}(\\mu), k = 6$"))
for(k in 1:n_k){
    contour(
    mu_i,
    mu_i,
    q_phi[,,k],
    levels = seq(0,1,len = 11),
    col    = cols,
    xlab   = TeX("$\\mu_1$"),
    ylab   = TeX("$\\mu_2$"),
    main   = labs[k])
}

fdir        =  file.path(getwd(), "10_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_10_einstichproben_t2_gütefunktion_1.pdf"),
width       = 12,
height      = 4)

```

```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("10_Abbildungen/mvda_10_einstichproben_t2_gütefunktion_1.pdf")
```

# Einstichproben-T$^2$-Tests
\noindent (5) Analyse der Testgütefunktion

\small
Beispiele

\footnotesize
\begin{equation*}
m := 2,
n := 15,
\Sigma := \begin{pmatrix} 1.0 & 0.9  \\ 0.9 & 1.0 \end{pmatrix},
\mu_0  := \begin{pmatrix} 1.0 \\ 1.0 \end{pmatrix} 
\end{equation*}
\vspace{1mm}

```{r, echo = F, eval = F}
# Modellparameter
library(matlib)                                                       # R Paket
m           = 2                                                       # m
n           = 15                                                      # n
nu          = (n-m)/((n-1)*m)                                         # \nu
Sigma       = matrix(c(1,.9,.9,1), nrow = 2)                          # \Sigma
iSigma      = inv(Sigma)                                              # \Sigma^{-1}

# Testparameter
mu_0        = matrix(c(1,1), nrow = 2)                                # \mu_0
k_all       = c(2,4,6)                                                # k <-> \phi
n_k         = length(k_all)                                           # Anzahl k Werte/Tests

# q_\phi(\mu) Evaluation
mu_min      = 0                                                       # \mu_i Minimum
mu_max      = 2                                                       # \mu_i Maximum
mu_res      = 1e3                                                     # \mu_i Auflösung
mu_i        = seq(mu_min,mu_max,len = mu_res)                         # mu_i
q_phi       = array(dim = c(mu_res, mu_res, length(k_all)))           # q_\phi Array
for(k in 1:n_k){
    for(i in 1:mu_res){
        for(j in 1:mu_res){
            mu           = matrix(c(mu_i[i],mu_i[j]), nrow = 2)       # \mu
            delta_mu     = n*t(mu - mu_0) %*% iSigma %*% (mu -mu_0)   # \delta_\mu
            q_phi[i,j,k] = 1 - pf(nu*k_all[k], m, n-m, delta_mu)}}}   # q_\phi(\mu)

graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,3),
mgp         = c(2,1,0),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1.2,
cex.main    = 1)
cols        = hcl.colors(11, "YlOrRd", rev = T)
labs        = list(TeX("$q_{\\phi}(\\mu), k = 2$"),
                   TeX("$q_{\\phi}(\\mu), k = 4$"),
                   TeX("$q_{\\phi}(\\mu), k = 6$"))
for(k in 1:n_k){
    contour(
    mu_i,
    mu_i,
    q_phi[,,k],
    levels = seq(0,1,len = 11),
    col    = cols,
    xlab   = TeX("$\\mu_1$"),
    ylab   = TeX("$\\mu_2$"),
    main   = labs[k])
}
fdir        =  file.path(getwd(), "10_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_10_einstichproben_t2_gütefunktion_2.pdf"),
width       = 12,
height      = 4)

```

```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("10_Abbildungen/mvda_10_einstichproben_t2_gütefunktion_2.pdf")
```


# Einstichproben-T$^2$-Tests
\noindent (6) Testumfangkontrolle

\footnotesize
\begin{theorem}[Testumfangkontrolle]
\justifying
\normalfont
$\phi$ sei der im obigen Testszenario definierte Test. Dann ist $\phi$ ein
Level-$\alpha_0$-Test mit Testumfang $\alpha_0$, wenn der kritische Wert definiert
ist durch
\begin{equation}
k_{\alpha_0} := \nu^{-1}F^{-1}\left(1-\alpha_0; m, n-m \right)
\end{equation}
wobei $\nu := (n-m)/((n-1)m)$ und $F^{-1}(\cdot;m,n-m)$ die inverse KVF der
$f$-Verteilung mit Freiheitsgradparametern $m$ und $n-m$ ist.
\end{theorem}


# Einstichproben-T$^2$-Tests
\footnotesize
\underline{Beweis}

Damit der betrachtete Test ein Level-$\alpha_0$-Test ist, muss bekanntlich
$q_\phi(\mu)\le \alpha_0$ für alle $\mu \in \{\mu_0\}$, also hier $q_\phi(\mu_0)\le \alpha_0$
gelten. Weiterhin ist der Testumfang des betrachteten Tests durch
$\alpha = \max_{\mu \in \{\mu_0\}} q_\phi(\mu)$, also hier durch $\alpha = q_\phi(\mu_0)$
gegeben. Wir müssen also zeigen, dass die Wahl von $k_{\alpha_0}$ garantiert,
dass $\phi$ ein Level-$\alpha_0$-Test mit Testumfang $\alpha_0$ ist. Dazu merken
wird zunächst an, dass für $\mu = \mu_0$ gilt, dass
\begin{equation}
q_\phi(\mu_0) = 1 - F(\nu k;m,n-m,\delta_{\mu_0}) = 1 - F(\nu k;m,n-m, 0)  = 1 - F(\nu k;m,n-m)
\end{equation}
wobei$F(\nu k; \delta, m, n-m)$ und $F(\nu k;m,n-m)$ die KVF der nichtzentralen
$f$-Verteilung mit Nichtzentralitätsparameter $\delta$ und Freiheitsgradparametern
$m$ und $n-m$ sowie der $f$-Verteilung mit Freiheitsgradparametern
$m$ und $n-m$, respektive, bezeichnen. Sei nun also $k := k_{\alpha_0}$. Dann gilt
\begin{align}
\begin{split}
q_\phi(\mu_0)
& = 1 - F(\nu k_{\alpha_0};m,n-m) \\
& = 1 - F\left(\nu \nu^{-1}F^{-1}\left(1-\alpha_0; m, n-m \right);m,n-m\right) \\
& = 1 - F\left(F^{-1}\left(1-\alpha_0; m, n-m \right);m,n-m\right) \\
& = 1 - (1 - \alpha_0)  = \alpha_0
\end{split}
\end{align}
Es folgt also direkt, dass bei der Wahl von $k = k_{\alpha_0}, q_\phi(\mu_0) \le \alpha_0$
ist der betrachtete Test somit ein Level-$\alpha_0$-Test ist. Weiterhin folgt direkt,
dass der Testumfang des betrachteten Tests bei der Wahl von $k = k_{\alpha_0}$ gleich
$\alpha_0$ ist.


# Einstichproben-T$^2$-Tests
\noindent (6) Testumfangkontrolle

\vfill
\center
\small
Wahl von $k_{\alpha_0} := \nu^{-1}F^{-1}\left(1-\alpha_0; m, n-m \right)$ mit $m = 2, n = 15$ und $\alpha_0 := 0.05$
\vspace{5mm}

```{r,echo = F, eval = F}
dev.new()
graphics.off()
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
m           = 2                                                                  # Datendimensionalität
n           = 15                                                                 # Anzahl Datenpunkte
nu          = (n-m)/(m*(n-1))                                                    # Parameter
alpha_0     = 0.05                                                               # Signifikanzlevel
k_alpha_0   = (1/nu)*qf(1 - alpha_0, m,n-m)                                       # kritischer Wert
t2          = seq(0,16,length=1e4)                                               # T^2-Statistikwerte
Pt2         = pf(nu*t2,m,n-m)                                                    # T^2-Statistik KVF für H_0
pt2         = nu*df(nu*t2,m,n-m)                                                 # T^2-Statistik WDF für H_0

# KVF Perspektive
plot(                                                                            # original density function
t2,
Pt2,
type        = "l",
ylab        = " ",
xlab        = TeX("$t^2$"),
ylim        = c(0,1),
main        = TeX("$P(t^2) = F(\\nu t^2; m,n-m)$"))

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

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

arrows(
x0          = min(t2),
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 ,-0.08, TeX("$\\k_{\\alpha_0}$"), xpd = TRUE)
text(1.5       , 1.00, TeX("$1 - \\alpha_0$"), xpd = TRUE)


# WDF Perspektive
plot(
t2,
pt2,
type        = "l",
ylab        = " ",
xlab        = TeX("$t^2$"),
ylim        = c(0,.4),
main        = TeX("$p(t^2) = \\nu f(\\nu t^2; m,n-m)$"))

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


lines(
seq(k_alpha_0, max(t2), 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,-0.03, TeX("$\\k_{\\alpha_0}$") , xpd = TRUE)
text(10, .05, TeX("$P(T^2 >= k_{\\alpha_0}) = \\alpha_0$"), xpd = TRUE, cex = .8, col = "gray50")

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

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

# Einstichproben-T$^2$-Tests
\noindent (6) Testumfangkontrolle
\vspace{1mm}

\tiny
\setstretch{.9}
```{r}
# Modellparameter
m         = 2                                                    # Dimensionalität der Zufallsvektoren/Daten
n         = 15                                                   # Anzahl der Datenpunkte
nu        = (n-m)/(m*(n-1))                                      # Parameter
mu_0      = matrix(c(1,1) , nrow = 2)                            # H0 Hypothesenparameter
mu        = mu_0                                                 # w.a.u. Erwartungswertparameter bei Zutreffen von H0
Sigma     = matrix(c(0.5,0.3, 0.3,0.5), nrow = 2, byrow = TRUE)  # wahrer, aber unbekannter, Kovarianzmatrixparameter

# Testparameter
alpha_0   = 0.05                                                 # Signifikanzlevel
k_alpha_0 = (1/nu)*qf(1-alpha_0, m,n-m)                          # kritischer Wert

# Simulation der Testumfangkontrolle
library(MASS)                                                    # R Paket für multivariate Normalverteilungen
library(matlib)                                                  # R Paket für Matrizenrechnung
nsim  = 1e5                                                      # Testentscheidungsarray
phi   = rep(NaN,nsim)                                            # Testentscheidungsarray
j_n   = matrix(rep(1,n), nrow = n)                               # 1_n
I_n   = diag(n)                                                  # I_n
J_n   = matrix(rep(1,n^2), nrow = n)                             # 1_{nn}
for(s in 1:nsim){                                                # Simulationsiterationen
    Y      = t(mvrnorm(n,mu,Sigma))                              # Y_i \sim N(\mu,\Sigma), i = 1,...,n
    Y_bar  = (1/n)*(Y %*% j_n)                                   # Stichprobenmittel
    C      = (1/(n-1))*(Y %*% (I_n-(1/n)*J_n) %*% t(Y))          # Stichprobenkovarianzmatrix
    T2     = n*t(Y_bar - mu_0) %*% inv(C) %*% (Y_bar - mu_0)     # T^2 Statistik
    if(T2 > k_alpha_0){                                         # Test 1_{T^2 >= k_alpha_0}
        phi[s] = 1                                               # Ablehnen von H_0
    } else {
        phi[s] = 0                                               # Nicht Ablehnen von H_0
    }
}
cat("\nKritischer Wert              = ", k_alpha_0,
    "\nGeschätzter Testumfang alpha = ", mean(phi))              # Ausgabe
```
# Einstichproben-T$^2$-Tests
\noindent (6) Testumfangkontrolle

\small
Praktisches Vorgehen

* \justifying Man nimmt an, dass ein vorliegender Datensatz $y_1,...,y_n$ eine
Realisation von $Y_1,...,Y_n \sim N(\mu,\Sigma)$ mit unbekannten Parametern
$\mu \in \mathbb{R}^m$ und $\Sigma \in \mathbb{R}^{m \times m} \mbox{p.d.}$ ist.

* Man möchte entscheiden ob für ein $\mu_0 \in \mathbb{R}^m$ eher $H_0 : \mu = \mu_0$ oder $H_1: \mu \neq \mu_0$ zutrifft.

* Man wählt ein Signifikanzlevel $\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 = 2$ und $n=15$, also
Freiheitsgradparametern $2$ und $13$, dass $k_{0.05}=\nu^{-1}F^{-1}(1 - 0.05;2,13) \approx 8.2$.

* Anhand von $m,n, \mu_0, \bar{Y}$ und $C$ berechnet man die Realisierung der 
Einstichproben-T$^2$-Teststatistik
\begin{equation}
T^2 := n(\bar{Y} - \mu_0)^T C^{-1}(\bar{Y} - \mu_0)
\end{equation}

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

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



# Einstichproben-T$^2$-Tests
\noindent (7) p-Werte

\small
Bestimmung des p-Wertes

* \justifying Per Definition ist der p-Wert das kleinste Signifikanzlevel $\alpha_0$, bei welchem
man die Nullhypothese basierend auf einem vorliegenden Wert der Teststatistik ablehnen würde.

* Bei $T^2 = t^2$ würde $H_0$ für jedes $\alpha_0$ mit $t^2 \ge \nu^{-1}F^{-1}(1-\alpha_0;m,n-m)$
abgelehnt werden. Für diese $\alpha_0$ gilt, wie unten gezeigt
\begin{equation}
\alpha_0 \ge \mathbb{P}\left(T^2 \ge t^2\right)
\end{equation}

* Das kleinste $\alpha_0 \in [0,1]$ mit $\alpha_0 \ge \mathbb{P}\left(T^2 \ge t^2\right)$
ist dann $\alpha_0 = \mathbb{P}(T^2 \ge t^2)$, also folgt
\begin{equation}
\mbox{ p-Wert } = \mathbb{P}\left(T^2 \ge t^2\right) = 1 - F(\nu t^2;m,n-m)
\end{equation}

```{r, echo = F, eval = F}
# Beispiel p-Werte
m  = c(2,2,2,4)
n  = c(15,15,99,15)
t2 = c(7,9,7,7)
p  = rep(NaN, length(n))
for(i in 1:length(n)){
  nu   = (n[i]-m[i])/(m[i]*(n[i]-1))
  p[i] = 1 - pf(nu*t2[i], m[i], n[i]-m[i])
}
print(p)
```

* Zum Beispiel ergibt sich bei

\footnotesize
* $m = 2$ und $n=15$ der p-Wert für $t^2 = 7.00$ zu 0.071
* $m = 2$ und $n=15$ der p-Wert für $t^2 = 9.00$ zu 0.040
* $m = 2$ und $n=99$ der p-Wert für $t^2 = 7.00$ zu 0.035
* $m = 4$ und $n=15$ der p-Wert für $t^2 = 7.00$ zu 0.304


# Einstichproben-T$^2$-Tests
\noindent (7) p-Werte

\small
Bestimmung des p-Wertes

\footnotesize
* Es bleibt zu zeigen, dass gilt
\begin{align}
\begin{split}
t^2 \ & \ge \nu^{-1}F^{-1}(1-\alpha_0;m,n-m) \\
\Leftrightarrow
\nu t^2 & \ge F^{-1}(1-\alpha_0;m,n-m) \\
\Leftrightarrow
\alpha_0 & \ge \mathbb{P}\left(T^2 \ge t^2\right)
\end{split}
\end{align}

* Dies aber folgt aus
\begin{align}
\begin{split}
t^2
& \ge \nu^{-1}F^{-1}(1-\alpha_0;m,n-m) \\
\nu t^2
& \ge F^{-1}(1-\alpha_0;m,n-m) \\
F(\nu t^2; m,n-m)
& \ge F\left(F^{-1}(1-\alpha_0;m,n-m); m,n - m\right) \\
F(\nu t^2; m,n-m)
& \ge 1 -\alpha_0\\
\mathbb{P}\left(T^2 \le t^2\right)
& \ge 1-\alpha_0 \\
\alpha_0
& \ge 1-\mathbb{P}\left(T^2 \le t^2\right)  \\
\end{split}
\end{align}

# Einstichproben-T$^2$-Tests
\small
Anwendungsbeispiel

\tiny
\setstretch{1}
```{r, message = F}
# R Pakete
library(foreign)                                                # Dateneinlesen
library(matlib)                                                 # Matrixalgebra

# Datenpräprozessierung
fname     = file.path(getwd(), "10_Daten", "studienerfolg.sav") # Dateinamen
D         = read.spss(fname, to.data.frame = T)                 # Dateneinlesen
Y         = rbind(D$X1[D$Gruppe == "ungenügend"],               # Y_{i_1} IQ Test Score
                  D$X2[D$Gruppe == "ungenügend"])               # Y_{i_2} Math Test Score

# Testparameter
m         = nrow(Y)                                             # Dimensionalität der Zufallsvektoren/Daten
n         = ncol(Y)                                             # Anzahl der Datenpunkte
nu        = (n-m)/(m*(n-1))                                     # Parameter
mu_0      = matrix(c(60,60) , nrow = 2)                         # H0 Hypothesenparameter ("Normwert")
alpha_0   = 0.05                                                # Signifikanzlevel
k_alpha_0 = (1/nu)*qf(1-alpha_0,m,n-m)                          # kritischer Wert

# Testevaluation
j_n       = matrix(rep(1,n), nrow = n)                          # 1_n
I_n       = diag(n)                                             # I_n
J_n       = matrix(rep(1,n^2), nrow = n)                        # 1_{nn}
Y_bar     = (1/n)*(Y %*% j_n)                                   # Stichprobenmittel
C         = (1/(n-1))*(Y %*% (I_n-(1/n)*J_n) %*% t(Y))          # Stichprobenkovarianzmatrix
T2        = n*t(Y_bar - mu_0) %*% inv(C) %*% (Y_bar - mu_0)     # T^2 Statistik
if(T2 > k_alpha_0){                                            # Test 1_{T^2 >= k_alpha_0}
    phi = 1                                                     # Ablehnen von H_0
} else {
    phi = 0                                                     # Nicht Ablehnen von H_0
}
p         = 1 - pf(nu*T2,m,n-m)                                 # p-Wert
```


# Einstichproben-T$^2$-Tests
\small
Anwendungsbeispiel

\vspace{1mm}
\tiny
\setstretch{.9}
```{r, message = F}
# Ausgabe
cat("Y_bar   = ", Y_bar,
    "\nC       = ", C,
    "\nT^2     = ", T2,
    "\nalpha_0 = ", alpha_0,
    "\nk       = ", k_alpha_0,
    "\nphi     = ", phi,
    "\np       = ", p)
```

\small
\vspace{1mm}
Anwendungsbeispiel mit `MVTests::OneSampleHT2()`
\vspace{1mm}

\tiny
\setstretch{.9}
```{r, message = F, warning = F}
library(MVTests)                                 # R Pakete
Y     = D[1:15,2:3]                              # Dataframe von Interesse
phi   = OneSampleHT2(Y, mu_0, alpha_0)           # Einstichproben-T^2-Test

# Ausgabe
cat("Y_bar   = "  , phi$Descriptive[2,],
    "\nT^2     = ", phi$HT2,
    "\nalpha_0 = ", phi$alpha,
    "\nk       = " ,phi$F,
    "\np       = ", phi$p.value)

```

# Einstichproben-T$^2$-Tests
\noindent (8) Analyse der Powerfunktion
\vfill
\justifying

\small
Wir betrachten die Testgütefunktion
\begin{equation}
q_\phi : \mathbb{R}^m \to [0,1], \mu \mapsto q_\phi(\mu) := 1 - F(\nu k;\delta_\mu,m,n-m)
\end{equation}
bei kontrolliertem Testumfang, also für
\begin{equation}
k_{\alpha_0} := \nu^{-1}F^{-1}\left(1-\alpha_0; m, n-m \right)
\end{equation}
mit festem $\alpha_0$ als Funktion des Nichtzentralitätsparameters und des Stichprobenumfangs.
Namentlich hängt hier $k_{\alpha_0}$ auch von $n$ ab.

Es ergibt sich die bivariate reellwertige Funktion
\begin{equation}
\pi : \mathbb{R} \times \mathbb{N} \to [0,1],
(\delta_\mu,n) \mapsto
\pi(\delta_\mu,n) :=  1 - F(\nu k_{\alpha_0};\delta_\mu,m,n-m)
\end{equation}
Bei festgelegtem $\alpha_0$ hängt die Powerfunktion des Einstichproben-T$^2$-Tests
also vom unbekannten Wert $\delta_\mu$, von der Datendimensionalität $m$ und
von der Stichprobengröße $n$ ab. Wir evaluieren und visualisieren diese Abhängigkeiten
untenstehend.
\vfill

# Einstichproben-T$^2$-Tests
\noindent (8) Analyse der Powerfunktion

\vspace{1mm}
\tiny
\setstretch{1}
```{r}
# Szenariospezifikationen
a_0_all   = c(0.05,0.01)                                              # \alpha_0 raum
d_mu_min  = 0                                                         # \delta_\mu Minimum
d_mu_max  = 20                                                        # \delta_\mu Maximum
d_mu_res  = 30                                                        # \delta_\mu Auflösung
d_mu_all  = seq(d_mu_min, d_mu_max, len = d_mu_res)                   # \delta_\mu d Raum
n_min     = 5                                                         # n Minimum
n_max     = 20                                                        # n Maximum
n_res     = 30                                                        # n Auflösung
n_all     = seq(n_min,n_max, len = n_res)                             # n Raum
m_all     = c(2,4)                                                    # m Raum

# Evaluation der Powerfunktion
pi        = array(dim = c(d_mu_res, n_res, 2,2))                      # Powerfunktionsarray
for (a in 1:length(a_0_all)){
    for (l in 1:length(m_all)){                                       # m Iterationen
        for(i in 1:length(d_mu_all)){                                 # \delta_\mu Iterationen
            for(j in 1:length(n_all)){                                # n Iterationen
                m           = m_all[l]                                # Datendimensionalität
                n           = n_all[j]                                # Stichprobenumfang
                d_mu        = d_mu_all[i]                             # wahrer, aber unbekannter, Parameter
                nu          = (n-m)/(m*(n-1))                         # Parameter
                alpha_0     = a_0_all[a]                              # Signifikanzlevel
                k_alpha_0   = (1/nu)*qf(1-alpha_0,m,n-m)              # kritischer Wert
                pi[i,j,l,a] = 1 - pf(nu*k_alpha_0, m, n-m, d_mu)}}}}  # Powerfunktionswert
```

```{r, echo = F, eval = F}
graphics.off()
fdir        =  file.path(getwd(), "10_Abbildungen")
dev.new()
par(
family      = "sans",
mfcol       = c(2,2),
pty         = "m",
bty         = "l",
lwd         = 1,
las         = 1,
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 1)
for(i in 1:length(m_all)){
    for(j in 1:length(a_0_all)){
        persp(
        d_mu_all,
        n_all,
        pi[,,i,j],
        d           = 1,
        col         = "gray90",
        theta       = -45,
        phi         = 30,
        lwd         = .5,
        scale       = T,
        ticktype    = "detailed",
        xlab        = "delta_mu",
        ylab        = "n",
        zlab        = "",
        main        = sprintf("alpha_0 = %0.2f, m = %0.0f ", a_0_all[j], m_all[i]),
        zlim        = c(0,1),
        r           = 1.5)
    }
}
dev.copy2pdf(
file        = file.path(fdir, "mvda_10_einstichproben_t2_pi.pdf"),
width       = 9,
height      = 9)
```

# Einstichproben-T$^2$-Tests
\vspace{2mm}
\noindent (8) Analyse der Powerfunktion


```{r, echo = FALSE, out.width = "65%"}
knitr::include_graphics("10_Abbildungen/mvda_10_einstichproben_t2_pi.pdf")
```

# Einstichproben-T$^2$-Tests
\noindent (8) Analyse der Powerfunktion

\small
Praktisches Vorgehen

\footnotesize
Mit größerem $n$ steigt die Powerfunktion des Tests an

* Ein großer Stichprobenumfang ist besser als ein kleiner Stichprobenumfang.

* Kosten für die Erhöhung des Stichprobenumfangs werden aber nicht berücksichtigt.

$\Rightarrow$ Die Theorie statistischer Hypothesentests ist nicht besonders lebensnah.

\vspace{1mm}

Die Powerfunktion hängt vom wahren, aber unbekannten, Parameterwert $\delta_\mu = n(\mu-\mu_0)^T\Sigma^{-1}(\mu-\mu_0)$ ab.

$\Rightarrow$ Wenn man $\delta_\mu$ schon kennen würde, würde man den Test nicht durchführen.

\vspace{1mm}

Generell wird folgendes Vorgehen favorisiert

* Man legt das Signifikanzlevel $\alpha_0$ fest und evaluiert die Powerfunktion.

* Man wählt einen Mindestparameterwert $\delta_\mu^*$, den man mit $\pi(\delta_\mu,n) = \beta$ detektieren möchte.

* Ein konventioneller Wert ist $\beta = 0.8$.

* Man liest die für $\pi(\delta_\mu = \delta_\mu^*,n) = \beta$ nötige Stichprobengröße $n$ ab.


#

(7) Analyse der Powerfunktion

\small
Praktisches Vorgehen

\tiny
\setstretch{1}
```{r}
# Szenariospezifikation
n_min      = 5                                               # n Minimum
n_max      = 20                                              # n Maximum
n_res      = 1e2                                             # n Auflösung
n          = seq(n_min,n_max, len = n_res)                   # n Raum
alpha_0    = 0.05                                            # Signifikanzlevel

# Poweranalyse
m          = 2                                               # Datendimensionalität
d_mu_fix   = 12                                              # fester Nichtzentralitätsparameter
nu         = (n-m)/(m*(n-1))                                 # Parameter
k_alpha_0  = (1/nu)*qf(1-alpha_0,m,n-m)                      # kritischer Wert
pi_n       = 1 - pf(nu*k_alpha_0, m, n-m, d_mu_fix)          # Powerfunktionswert
beta       = 0.8                                             # gewünschter Powerfunktionswert
i          = 1                                               # Indexinitialisierung
n_min      = NaN                                             # minimales n Initialisierung
while(pi_n[i] < beta){                                       # Solange \pi(\delta_\mu*,n) < \beta
    n_min = n[i]                                             # Aufnahme des minimal nötigen ns
    i     = i + 1                                            # und Erhöhung des Indexes
}
cat("Minimal nötiges n =", ceiling(n_min))                   # Ausgabe
```
#

(7) Analyse der Powerfunktion

\small
Praktisches Vorgehen
\vspace{1mm}
```{r, echo = F, eval = F}
dev.new()
graphics.off()
fdir        =  file.path(getwd(), "10_Abbildungen")
par(
family      = "sans",
mfcol       = c(1,1),
pty         = "m",
bty         = "l",
lwd         = 1,
las         = 1,
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1.1,
cex.main    = 1.2)

plot(
n,
pi_n,
type        = "l",
lwd         = 1,
ylab        = " ",
ylim        = c(0,1),
xlab        = TeX("$n$"),
main        = TeX("$\\pi(d = 12,n)\\, für \\,\\alpha_0 = 0.05,\\, \\mu_0 = 0"))

lines(
5,
beta,
type        = "p",
pch         = 16,
xpd         = TRUE)

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

arrows(
x0          = min(n),
y0          = beta,
x1          = n_min,
y1          = beta,
col         = "darkorange",
angle       = 20,
length      = .1)

arrows(
x0          = n_min,
y0          = beta,
x1          = n_min,
y1          = 0,
col         = "darkorange",
angle       = 20,
length      = .1)
text(3 , 0.85 , TeX("$\\beta$") ,xpd = TRUE, cex = 1.2)
text(17 ,.05  , TeX("$n_{opt}") ,xpd = TRUE, cex = 1.2)

dev.copy2pdf(
file        = file.path(fdir, "mvda_10_einstichproben_t2_stichprobenumfang.pdf"),
width       = 6,
height      = 5)

```

```{r, echo = FALSE, out.width = "60%"}
knitr::include_graphics("10_Abbildungen/mvda_10_einstichproben_t2_stichprobenumfang.pdf")
```

#
\vfill
\setstretch{2.4}
\large
Vorbemerkungen

Einstichproben-T$^2$-Tests

**Zweistichproben-T$^2$-Tests**

Univariates vs. multivariates Testen

Selbstkontrollfragen
\vfill


# Zweistichproben-T$^2$-Tests

Zweistichproben-T$^2$-Tests **bei unabhängigen Stichproben**

Anwendungsszenario

\small

* **Zwei Stichproben** experimenteller Einheiten.

* Annahme unabhängiger und identisch nach $N(\mu_1,\Sigma_1)$ und $N(\mu_2,\Sigma_2)$ verteilter Daten.

* $\mu_1,\Sigma_1, \mu_2,\Sigma_2$ unbekannt.

* Quantifizieren der Unsicherheit beim inferentiellen Vergleich von $\mu_1$ und $\mu_2$ beabsichtigt.

\normalsize
Anwendungsbeispiele

\small
* Gruppenvergleich von BDI und Glukokortikoid Daten

$\quad\quad\quad$ $\circ$ $\mu_1 \neq \mu_2$ als Evidenz für multivariate Gruppenunterschiede


* Gruppenvergleich von Testdaten bei erfolgreichen vs. nicht-erfolgreichem Studienabschluss

$\quad\quad\quad$ $\circ$ $\mu_1 \neq \mu_2$ als Evidenz für multivariate Gruppenunterschiede


# Zweistichproben-T$^2$-Tests
\setstretch{1.6}

Mögliche Modellszenarien

* Annahme identischer Kovarianzmatrixparameter

* Annahme eines bekannten Varianzverhältnisses

* Keine Annahmen zu Varianzen

\vspace{3mm}

Mögliche Hypothesenszenarien

* $H_0: \mu_1 =   \mu_2$ und $H_1: \mu_1 \neq \mu_2$

* $H_0: \mu_1 \le \mu_2$ und $H_1: \mu_1 >    \mu_2$

* $H_0: \mu_1 \ge \mu_2$ und $H_1: \mu_1 <    \mu_2$

# Zweistichproben-T$^2$-Tests
\setstretch{1.6}

Hier betrachtetes Modellszenario

\begin{itemize}
\item Annahme identischer Kovarianzmatrixparameter 
\item \textcolor{lightgray}{Annahme eines bekannten Varianzverhältnisses}
\item \textcolor{lightgray}{Keine Annahmen zu Varianzen}
\end{itemize}
\vspace{3mm}

Hier betrachtetes Hypothesenszenario

\begin{itemize}
\item $H_0: \mu_1 =   \mu_2$ und $H_1: \mu_1 \neq \mu_2$
\item \textcolor{lightgray}{$H_0: \mu_1 \le \mu_2$ und $H_1: \mu_1 > \mu_2$}
\item \textcolor{lightgray}{$H_0: \mu_1 \ge \mu_2$ und $H_1: \mu_1 <    \mu_2$}
\end{itemize}

#  Zweistichproben-T$^2$-Tests
Anwendungsbeispiel

\tiny
\setstretch{.8}
```{r, echo = F, message = F}
# Einlesen des Datensatzes aus Rudolf & Buse (2019) Multivariate Verfahren
library("foreign")
D          = read.spss(file.path(getwd(), "10_Daten", "studienerfolg.sav"), to.data.frame = T)
D          = rbind(D[D$Gruppe == "ungenügend",2:3],
                   D[D$Gruppe == "gut",2:3])
D          = cbind(rep(NaN,30),rep(NaN,30), D)
D[,1]      = c(rep(1,15), rep(2,15))
D[,2]      = c(1:15,1:15)
colnames(D) = c("Gruppe",
                "VP",
                "IQ Test Score",
                "Math Test Score")
knitr::kable(D, "pipe", align = "c", row.names = F)
```

\footnotesize
Unterschiede zwischen den (IQ Test Score, Math Test Score) Erwartungswertparametern zwischen Gruppen?


#  Zweistichproben-T$^2$-Tests
Anwendungsbeispiel
\vspace{1mm}

\tiny
\setstretch{.9}
```{r, message = F}
# Deskriptivstatistik
library(foreign)                                                    # R Paket
library(ellipse)                                                    # R Paket
library(matlib)                                                     # R Paket
fname    = file.path(getwd(), "10_Daten", "studienerfolg.sav")      # Dateiname
D        = read.spss(fname, to.data.frame = T)                      # Dateneinlesen
m        = 2                                                        # Datendimensionalität
n        = 15                                                       # Datenpunkte pro Gruppe
n_1      = n                                                        # Datenpunkte Gruppe 1
n_2      = n                                                        # Datenpunkte Gruppe 2
j_n      = matrix(rep(1,n), nrow = n)                               # 1_n
I_n      = diag(n)                                                  # I_n
J_n      = matrix(rep(1,n^2), nrow = n)                             # 1_{nn}
Y_1      = t(as.matrix(D[D$Gruppe == "ungenügend", 2:3]))           # Daten Gruppe 1
Y_2      = t(as.matrix(D[D$Gruppe == "gut", 2:3]))                  # Daten Gruppe 2
Y_1_bar  = (1/n)*(Y_1 %*% j_n)                                      # Stichprobenmittel
Y_2_bar  = (1/n)*(Y_2 %*% j_n)                                      # Stichprobenmittel
C_1      = (1/(n_1-1))*(Y_1 %*% (I_n-(1/n)*J_n) %*% t(Y_1))         # Stichprobenkovarianzmatrix
C_2      = (1/(n_2-1))*(Y_2 %*% (I_n-(1/n)*J_n) %*% t(Y_2))         # Stichprobenkovarianzmatrix
C        = ((n_1-1)*C_1+(n_2-1)*C_2)/(n_1+n_2-2)                    # Gepoolte Stichprobenkovarianzmatrix
D        = t(Y_1_bar - Y_2_bar) %*% inv(C) %*% (Y_1_bar - Y_2_bar)  # Mahalanobis Distanz

# Ausgabe
cat("Y_1_bar  =", Y_1_bar,
    "\nY_2_bar  =", Y_2_bar,
    "\nD        =", D)
```

# Zweistichproben-T$^2$ Tests
Anwendungsbeispiel

```{r, echo = F, eval = F}
library(ellipse)
library(matlib)
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,
cex.main    = 1)

# Gruppe 1
plot(
Y_1[1,],
Y_1[2,],
col  = "Gray70",
bg   = "Gray70",
pch  = 21,
xlab = TeX("IQ Test Score"),
ylab = TeX("Math Test Score"),
xlim = c(0,100),
ylim = c(0,100))
points(
Y_1_bar[1],
Y_1_bar[2],
col  = "White",
bg   = "gray70",
pch  = 24)
lines(
ellipse(C, level = 0.40, centre = Y_1_bar),
type = "l",
col = "Gray70")

# Gruppe 2
points(
Y_2[1,],
Y_2[2,],
col  = "White",
bg   = "Gray30",
pch  = 21)
points(
Y_2_bar[1],
Y_2_bar[2],
col  = "White",
bg   = "Gray30",
pch  = 24)
lines(
ellipse(C, level = 0.40, centre = Y_2_bar),
type = "l",
col = "Gray30")

legend(
"topleft",
c("Datenpunkte Gruppe 1", "Stichprobenmittel Gruppe 1", "Gepoolte Stichprobenkovarianz",
  "Datenpunkte Gruppe 2", "Stichprobenmittel Gruppe 2", "Gepoolte Stichprobenkovarianz"),
lty = c(NA,NA,1,NA,NA,1),
pch = c(19,17,NA,19,17,NA),
col = c("gray70","gray70","gray70","gray30","gray30","gray30"),
bty = "n",
cex = .6)

dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_zweistichproben_t2_deskription.pdf"),
width  = 4.5,
height = 4.5)
```

```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics("10_Abbildungen/mvda_10_zweistichproben_t2_deskription.pdf")
```

# Zweistichproben-T$^2$-Tests
\setstretch{2}
Gliederung (vgl. [\textcolor{red}{(12) Hypothesentests}](https://bit.ly/32ZUAiy))

\small
(1) Statistisches Modell in klassischer Form
(2) Statistisches Modell in generativer Form
(3) Testhypothesen, Teststatistik, Test
(4) Analyse der Teststatistik
(5) Analyse der Testgütefunktion
(6) Testumfangkontrolle
(7) p-Werte

\normalsize
Zur Wiederholung des univariaten Falls, siehe [\textcolor{red}{(14) Zweistichproben-T-Tests}](https://bit.ly/3qXv2eB)

# Zweistichproben-T$^2$-Tests
\justifying
\noindent (1) Statistisches Modell in klassischer Form

\footnotesize
$Y_{11},...,Y_{1n_1} \sim N(\mu_1,\Sigma)$ sei eine Stichprobe eines multivariaten
Normalverteilungsmodells mit unbekanntem Erwartungswertparameter $\mu_1 \in \mathbb{R}^m$
und unbekanntem  Kovarianzmatrixparameter $\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d}$.
$Y_{21},..., Y_{2n_2} \sim N(\mu_2,\Sigma)$  sei eine weitere Stichprobe eines
multivariaten Normalverteilungsmodells mit unbekanntem Erwartungswertparameter
$\mu_2 \in \mathbb{R}^m$ und unbekanntem  Kovarianzmatrixparameter
$\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d}$. Die Kovarianzmatrixparameter
beider Stichproben werden also als identisch vorausgesetzt. Der Parameter von
Interesse ist $(\mu_1,\mu_2)$, der Parameterraum des Modells ist
$\Theta := \mathbb{R}^m \times \mathbb{R}^m$.

\vspace{1mm}
\normalsize
\noindent (2) Statistisches Modell in generativer Form

\footnotesize
Es sei
\begin{equation}
Y_{ij} = \mu_i + \varepsilon_{ij}
\mbox{ mit } \varepsilon_{ij} \sim N(0,\Sigma)
\mbox{ für } i = 1,2
\mbox{ und } j = 1,...,n_i
\end{equation}
wobei

* $i$ die Stichproben indiziert,
* $j$ die experimentellen Einheiten indiziert,
* $n_i$ die Stichprobengrößen sind,
* $Y_{ij}$ beobachtbare Zufallsvariablen sind,
* $\mu_i \in \mathbb{R}^m$ feste Erwartungswertparameter der Stichprobenvariablen sind
* $\Sigma \in \mathbb{R}^{m \times m} \mbox{ p.d. }$ der identische Kovarianzmatrixparameter über Stichproben, und
* $\varepsilon_{ij}$ unabhängige multivariat-normalverteilte nicht-beobachtbare Zufallsvariablen sind.

Der Zusammenhang zwischen klassischer und generativer Modellform ergibt sich analog zum Einstichprobenfall.

# Zweistichproben-T$^2$-Tests
\justifying
\noindent (3) Testhypothesen, Teststatistik, Test

\footnotesize
\setstretch{1.2}
Wir betrachten die einfache Nullhypothese
\begin{equation}
H_0 : \mu_1 = \mu_2    \Leftrightarrow \Theta_0 := \{(\mu_1,\mu_2) \in \mathbb{R}^m \times \mathbb{R}^m | \mu_1 = \mu_2\}
\end{equation}
und die zusammengesetzte Alternativhypothese
\begin{equation}
H_1 : \mu_1 \neq \mu_2 \Leftrightarrow \Theta_1 := \{(\mu_1,\mu_2) \in \mathbb{R}^m \times \mathbb{R}^m | \mu_1 \neq \mu_2\}
\end{equation}
Weiterhin betrachten die *Zweistichproben-T$^2$-Teststatistik*
\begin{equation}
T^2 := \frac{n_1n_2}{n_1 + n_2}\left(\bar{Y}_1 - \bar{Y}_2 \right)^TC^{-1}\left(\bar{Y}_1 - \bar{Y}_2 \right)
\end{equation}
wobei für $i=1,2$ und respektiven Stichprobenkovarianzmatrizen $C_1$ und $C_2$
\begin{equation}
\bar{Y}_i := \frac{1}{n_i}\sum_{j=1}^{n_i} Y_{ij}\mbox{ und } C := \frac{(n_1 - 1)C_1 + (n_2 - 1)C_2}{n_1 + n_2 - 2}
\end{equation}
die Stichprobenmittel und die *gepoolte Stichprobenkovarianzmatrix*, respektive, bezeichnen.

* $T^2$ ist die mit den Stichprobenumfängen skalierte Mahalanobis Distanz von $\bar{Y}_1$ und $\bar{Y}_2$ hinsichtlich $C$.
* Größere $T^2$ Werte ergeben sich für größere Abstände, geringere Kovarianzen und größere Stichprobenumfänge.

Schließlich definieren wir den kritischen Wert-basierten Test
\begin{equation}
\phi(X) := 1_{\{T^2 \ge k \}}.
\end{equation}


# Zweistichproben-T$^2$-Tests
\noindent (4) Analyse der Teststatistik

\setstretch{1.2}
\footnotesize
\begin{theorem}[Verteilung der Zweistichproben-T$^2$ Statistik]
\justifying
\normalfont
Für $i = 1,2$ seien $Y_{i1},...,Y_{in_i} \sim N(\mu_i,\Sigma)$
mit $\mu_i \in \mathbb{R}^{m}$ und $\Sigma \in \mathbb{R}^{m\times m} \mbox{ p.d.,}$,
\begin{equation}
\nu := \frac{n_1+n_2-m-1}{m(n_1+n_2-2)}
\end{equation}
und die Zweistichproben-$T^2$-Teststatistik sei definiert als
\begin{equation}
T^2 := \frac{n_1n_2}{n_1 + n_2}(\bar{Y}_1 - \bar{Y}_2)^T C^{-1}(\bar{Y}_1 - \bar{Y}_2)
\end{equation}
mit den Stichprobenmitteln $\bar{Y}_i$ und der gepoolten Stichprobenkovarianzmatrix $C$.
Dann gilt
\begin{equation}
\nu T^2 \sim f(\delta,m,n_1+n_2-m-1),
\end{equation}
wobei $f(\delta,m,n_1+n_2-m-1)$ die nichtzentrale $f$-Verteilung mit Nichtzentralitätsparameter
\begin{equation}
\delta := \frac{n_1n_2}{n_1+n_2} (\mu_1 - \mu_2)^T\Sigma^{-1}(\mu_1 - \mu_2)
\end{equation}
sowie mit Freiheitsgradparametern $m$ und $n_1+n_2-m-1$ bezeichnet.
\end{theorem}

\tiny
Bemerkungen

* Für einen Beweis verweisen wir auf @anderson_2003.
* Für $\mu_1 = \mu_2$ und damit $\delta = 0$ entspricht $f(\delta,m,n_1+n_2-m-1)$ der $f$-Verteilung $f(m,n_1+n_2-m-1)$.

# Zweistichproben-T$^2$-Tests
\noindent (4) Analyse der Teststatistik

\setstretch{1.2}
\footnotesize
\begin{theorem}[WDF und KVF der Zweistichproben-T$^2$ Statistik]
\justifying
\normalfont
Im Zweistichproben-$T^2$-Testszenario sei
\begin{equation}
\nu := \frac{n_1+n_2-m-1}{m(n_1+n_2-2)}
\end{equation}
Dann ist eine WDF der Zweistichproben-T$^2$ Statistik gegeben durch
\begin{equation}
p_{T^2} : \mathbb{R}_{\ge 0} \to \mathbb{R}, t^2 \mapsto p_{T^2}(t^2) := \nu f(\nu t^2; \delta, m,n_1+n_2-m-1 )
\end{equation}
und eine KVF der Zweistichproben-T$^2$ Statistik ist gegeben durch
\begin{equation}
P_{T^2} : \mathbb{R}_{\ge 0} \to \mathbb{R}, t^2 \mapsto P_{T^2}(t^2) := Ft(\nu t^2; \delta, m,n_1+n_2-m-1)
\end{equation}
\end{theorem}
Bemerkung

* Der Beweis erfolgt in Analogie zum Einstichproben-$T^2$-Testszenario.

# Zweistichproben-T$^2$-Tests
\vspace{1mm}
\noindent (4) Analyse der Teststatistik
\vspace{1mm}

\tiny
\setstretch{0.8}
```{r}
# Modellparameter
m     = 2                                                       # Dimensionalität der Zufallsvektoren/Daten
n     = 15                                                      # Anzahl Datenpunkte pro Stichprobe
n_1   = n                                                       # Anzahl Datenpunkte Stichprobe 1
n_2   = n                                                       # Anzahl Datenpunkte Stichprobe 1
mu_1  = matrix(c(1,1) , nrow = 2)                               # wahrer, aber unbekannter, Erwartungswertparameter
mu_2  = matrix(c(2,2) , nrow = 2)                               # wahrer, aber unbekannter, Erwartungswertparameter
Sigma = matrix(c(0.4,0.1, 0.1,0.4), nrow = 2, byrow = TRUE)     # wahrer, aber unbekannter, Kovarianzmatrixparameter

# Simulation
library(MASS)                                                   # R Paket für multivariate Normalverteilungen
library(matlib)                                                 # R Paket für Matrizenrechnung
nsim  = 1e4                                                     # Anzahl Simulationen/Datensatzrealisierungen
Y1b   = matrix(rep(NaN,m*nsim), nrow = 2)                       # Stichprobenmittelarray
Y2b   = matrix(rep(NaN,m*nsim), nrow = 2)                       # Stichprobenmittelarray
T2    = rep(NaN,nsim)                                           # T2 Statistik Array
j_n   = matrix(rep(1,n), nrow = n)                              # 1_n
I_n   = diag(n)                                                 # I_n
J_n   = matrix(rep(1,n^2), nrow = n)                            # 1_{nn}
for(s in 1:nsim){                                               # Simulationsiterationen
    Y_1      = t(mvrnorm(n,mu_1,Sigma))                         # Y_1j \sim N(\mu_1,\Sigma), j = 1,...,n_1
    Y_2      = t(mvrnorm(n,mu_2,Sigma))                         # Y_2j \sim N(\mu_2,\Sigma), j = 1,...,n_2
    Y_1_bar  = (1/n)*(Y_1 %*% j_n)                              # Stichprobenmittel
    Y_2_bar  = (1/n)*(Y_2 %*% j_n)                              # Stichprobenmittel
    C_1      = (1/(n_1-1))*(Y_1 %*% (I_n-(1/n)*J_n) %*% t(Y_1)) # Stichprobenkovarianzmatrix
    C_2      = (1/(n_2-1))*(Y_2 %*% (I_n-(1/n)*J_n) %*% t(Y_2)) # Stichprobenkovarianzmatrix

    # Gepoolte Stichprobenkovarianzmatrix und T2 Statistik
    C        = ((n_1-1)*C_1+(n_2-1)*C_2)/(n_1+n_2-2)
    T2[s]    = (((n_1*n_2)/(n_1+n_2))*t(Y_1_bar-Y_2_bar) %*% inv(C) %*% (Y_1_bar-Y_2_bar))

    # Stichprobenmittel für Visualisierung
    Y1b[,s]  = Y_1_bar
    Y2b[,s]  = Y_2_bar
}
```


# Zweistichproben-T$^2$-Tests
\vspace{1mm}
\noindent (4) Analyse der Teststatistik
\vspace{1mm}

```{r, echo = F, eval = F}
# figure parameters
dev.new()
graphics.off()
fdir        =  file.path(getwd(), "10_Abbildungen")
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)

# Stichprobenmittel per Simulation
plot(
Y1b[1,],
Y1b[2,],
xlim     = c(0,3),
ylim     = c(0,3),
xlab     = TeX("$\\bar{y}_1$"),
ylab     = TeX("$\\bar{y}_2$"),
pch      = 21,
col      = "white",
bg       = "gray70",
main     = TeX("$\\bar{Y}_1$ und $\\bar{Y}_2$  Simulationen"))
points(
Y2b[1,],
Y2b[2,],
pch      = 21,
col      = "white",
bg       = "gray30")

# Unnormalisierte T2-Teststatistik WDF
t2_min      = 0
t2_max      = 100
t2_res      = 1e3
nu          = (n_1+n_2-m-1)/(m*(n_1+n_2-2))
delta       = ((n_1*n_2)/(n_1+n_2))*t(mu_1 - mu_2) %*% inv(Sigma) %*% (mu_1-mu_2)
t2          = seq(t2_min,t2_max,len = t2_res)
p_T2        = nu*df(nu*t2,m,n_1+n_2-m-1,delta)
hist(
T2,
breaks      = 100,
col         = "gray90",
prob        = TRUE,
xlim        = c(t2_min, t2_max),
ylim        = c(0,0.04),
xlab        = TeX("$t^2$"),
ylab        = "",
main        = TeX("$p(t^2) := \\nu f(\\nu t^2; \\delta, m, n_1+n_2-m-1,)$"))
lines(
t2,
p_T2,
lwd         = 2,
col         = "darkorange")

# Normalisierte T2-Teststatistik Sample
t2_min      = 0
t2_max      = 50
t2_res      = 1e3
t2          = seq(t2_min,t2_max,len = t2_res)
p_nm_T2     = df(t2,m,n_1+n_2-m-1,delta)
hist(
nu*T2,
breaks      = 100,
col         = "gray90",
prob        = TRUE,
xlim        = c(t2_min, t2_max),
ylim        = c(0,0.1),
xlab        = TeX("$t^2$"),
ylab        = "",
main        = TeX("$\\nu T^2 ~ f(\\delta,m,n_1+n_2-m-1)$") )
lines(
t2,
p_nm_T2,
lwd         = 2,
col         = "darkorange")

# Unnormalisierte T2-Teststatistik CDF
t2_min      = 0
t2_max      = 200
t2_res      = 1e3
t2          = seq(t2_min,t2_max,len = t2_res)
P_T2        = pf(nu*t2,m,n_1+n_2-m-1,delta)
ECDF_T2     = ecdf(T2)
plot(
ECDF_T2,
verticals   = TRUE,
lwd         = 6,
xlim        = c(t2_min, t2_max),
col         = "gray90",
xlab        = TeX("$t^2$"),
ylab        = "",
main        = TeX("$P(t^2) := F(\\nu t^2;\\delta, m,n_1+n_2-m-1)$"))
lines(
t2,
P_T2,
lwd         = 2,
col         = "darkorange")


dev.copy2pdf(
file        = file.path(fdir, "mvda_10_zweistichproben_t2_teststatistik.pdf"),
width       = 6,
height      = 6)
```

```{r, echo = FALSE, out.width = "60%"}
knitr::include_graphics("10_Abbildungen/mvda_10_zweistichproben_t2_teststatistik.pdf")
```


# Zweistichproben-T$^2$-Tests
\noindent \noindent (5) Analyse der Testgütefunktion

\setstretch{1.2}
\footnotesize
\begin{theorem}[Testgütefunktion]
\justifying
\normalfont
$\phi$ sei der im obigen Testszenario definiert Test. Dann ist die Testgütefunktion von $\phi$
gegeben durch
\begin{equation}
q_\phi : \mathbb{R}^m \times \mathbb{R}^m \to [0,1], \mu \mapsto q_\phi(\mu_1,\mu_2) := 1 - F(\nu k;\delta_{\mu_{1,2}},m,n_1+n_2-m-1),
\end{equation}
wobei
\begin{equation}
\nu := \frac{n_1+n_2-m-1}{m(n_1+n_2-2)}, \quad
\delta_{\mu_{1,2}} := \frac{n_1n_2}{n_1+n_2} (\mu_1 - \mu_2)^T\Sigma^{-1}(\mu_1 - \mu_2)
\end{equation}
und $F(\cdot;\delta_{\mu_{1,2}}, m,n_1+n_2-m-1)$ die KVF der nichtzentralen
$f$-Verteilung mit Freiheitsgradparametern $m$ und $n_1+n_2-m-1$ sowie mit
Nichtzentralitätsparameter $\delta_{\mu_{1,2}}$ bezeichnet.
\end{theorem}

\footnotesize
Bemerkungen

* $q_\phi$ kann zur Bestimmung kritischer Werte für einen erwünschten Testumfang genutzt werden.
* $q_\phi$ kann zur Bestimmung der Testpower genutzt werden.

# Einstichproben-T$^2$-Tests
\noindent \noindent (5) Analyse der Testgütefunktion

\footnotesize
\setstretch{1.2}
\underline{Beweis}

Die Testgütefunktion des betrachteten Tests im vorliegenden Testszenario ist definiert als
\begin{equation}
q_\phi : \mathbb{R}^m \times \mathbb{R}^m \to [0,1], (\mu_1, \mu_2) \mapsto q_\phi(\mu_1, \mu_2) := \mathbb{P}_{\mu_1,\mu_2}(\phi = 1)
\end{equation}
Da die Wahrscheinlichkeiten für $\phi = 1$ und dafür, dass die zugehörige Teststatistik
im Ablehnungsbereich des Tests liegt, gleich sind, benötigen wir also zunächst die Verteilung
der Teststatistik. Wir haben oben aber bereits gesehen, dass für
\begin{equation}
\nu := \frac{n_1+n_2-m-1}{(m(n_1+n_2-2)}
\mbox{ und }
\delta_{\mu_{1,2}} := \frac{n_1n_2}{n_1+n_2} (\mu_1 - \mu_2)^T\Sigma^{-1}(\mu_1 - \mu_2)
\end{equation}
gilt, dass
\begin{equation}
\nu T^2 \sim f(\delta_{\mu_{1,2}}, m,n_1+n_2-m-1).
\end{equation}
Der Ablehnungsbereich des betrachteten Tests ergibt sich als $A := ]k,\infty[$
Also ergibt sich
\begin{align}
\begin{split}
q_\phi(\mu)
&  = \mathbb{P}_\mu(\phi = 1) \\
&  = \mathbb{P}_\mu\left(T^2 \in \,\,]k,\infty[\right) \\
&  = \mathbb{P}_\mu\left(T^2 > k \right) \\
&  = 1 - \mathbb{P}_\mu\left(T^2 \le k \right) \\
&  = 1 - F\left(\nu k;\delta_{\mu_{1,2}}, m,n_1+n_2-m-1 \right)
\end{split}
\end{align}

# Zweistichproben-T$^2$-Tests
\noindent (5) Analyse der Testgütefunktion

\small
Beispiele
\footnotesize
\begin{equation*}
m      := 2,
n_1    := 15,
n_2    := 15,
\mu_2  := \begin{pmatrix} 1 \\ 1 \end{pmatrix},
\Sigma := \begin{pmatrix} 1 & 0  \\ 0 & 1 \end{pmatrix},
\end{equation*}

\tiny
\setstretch{1}
```{r, eval = F}
# Modellparameter
library(matlib)                                                        # R Paket
m           = 2                                                        # m
n_1         = 15                                                       # n_1
n_2         = 15                                                       # n_2
nu          = (n_1+n_2-m-1)/(m*(n_1+n_2-2))                            # \nu
Sigma       = diag(m)                                                  # \Sigma = I_2
iSigma      = inv(Sigma)                                               # \Sigma^{-1}
mu_2        = matrix(c(1,1), nrow = 2)                                 # \mu_2

# Testparameter
k_all       = c(2,4,6)                                                 # k <-> \phi
n_k         = length(k_all)                                            # Anzahl k Werte/Tests

# q_\phi(\mu) Evaluation
mu_min      = 0                                                        # \mu_1 Minimum
mu_max      = 2                                                        # \mu_1 Maximum
mu_res      = 1e3                                                      # \mu_1 Auflösung
mu_i        = seq(mu_min,mu_max,len = mu_res)                          # mu_1i
q_phi       = array(dim = c(mu_res, mu_res, length(k_all)))            # q_\phi Array
for(k in 1:n_k){
    for(i in 1:mu_res){
        for(j in 1:mu_res){
            mu_1         = matrix(c(mu_i[i],mu_i[j]), nrow = 2)        # \mu_1
            delta        = (((n_1*n_2)/(n_1+n_2))*                     # \delta_{\mu_{1,2}}
                           t(mu_1 - mu_2) %*% iSigma %*% (mu_1-mu_2))
            q_phi[i,j,k] = 1-pf(nu*k_all[k],m,n_1+n_2-m-1,delta)}}}    # q_\phi(\mu)
```

# Zweistichproben-T$^2$-Tests
\noindent (5) Analyse der Testgütefunktion

\small
Beispiele
\footnotesize
\begin{equation*}
m      := 2,
n_1    := 15,
n_2    := 15,
\mu_2  := \begin{pmatrix} 1 \\ 1 \end{pmatrix},
\Sigma := \begin{pmatrix} 1 & 0  \\ 0 & 1 \end{pmatrix} 
\end{equation*}


```{r, echo = F, eval = F}
graphics.off()
dev.new()
par(
family        = "sans",
mfcol       = c(1,3),
mgp         = c(2,1,0),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1.2,
cex.main    = 1)
cols        = hcl.colors(11, "YlOrRd", rev = T)
labs        = list(TeX("$q_{\\phi}(\\mu_1, \\mu_2), k = 2$"),
                   TeX("$q_{\\phi}(\\mu_1, \\mu_2), k = 4$"),
                   TeX("$q_{\\phi}(\\mu_1, \\mu_2), k = 6$"))
for(k in 1:n_k){
    contour(
    mu_i,
    mu_i,
    q_phi[,,k],
    levels = seq(0,1,len = 11),
    col    = cols,
    xlab   = TeX("$\\mu_{1_1}$"),
    ylab   = TeX("$\\mu_{1_2}$"),
    main   = labs[k])
}
fdir        =  file.path(getwd(), "10_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "mvda_10_zweistichproben_t2_gütefunktion_1.pdf"),
width       = 12,
height      = 4)

```

```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("10_Abbildungen/mvda_10_zweistichproben_t2_gütefunktion_1.pdf")
```

# Zweistichproben-T$^2$-Tests
\vspace{2mm}
\noindent (6) Testumfangkontrolle
\footnotesize
\setstretch{1.1}
\begin{theorem}[Testumfangkontrolle]
\justifying
\normalfont
$\phi$ sei der im obigen Testszenario definierte Test. Dann ist $\phi$ ein
Level-$\alpha_0$-Test mit Testumfang $\alpha_0$, wenn der kritische Wert definiert
ist durch
\begin{equation}
k_{\alpha_0} := \frac{1}{\nu}F^{-1}\left(1-\alpha_0; m, n_1 + n_2 - m - 1 \right)
\end{equation}
wobei $F^{-1}(\cdot;m,n_1+n_2-m-1)$ die inverse KVF der $f$-Verteilung mit Freiheitsgradparametern
$m$ und $n_1 + n_2 - m - 1$ ist.
\end{theorem}

Bemerkung

* Der Beweis erfolgt analog zum Einstichprobenszenario.

# Zweistichproben-T$^2$-Tests
\noindent (6) Testumfangkontrolle

\vfill
\center
\footnotesize
Wahl von $k_{\alpha_0} := \frac{1}{\nu}F^{-1}\left(1-\alpha_0; m, n_1+n_2-m-1 \right)$ mit $m = 2, n_1 = 15, n_2 = 15,\alpha_0 := 0.05$
\vspace{5mm}

```{r,echo = F, eval = F}
dev.new()
graphics.off()
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
m           = 2                                                                  # Datendimensionalität
n_1         = 15                                                                 # n_1
n_2         = 15                                                                 # n_1
nu          = (n_1+n_2-m-1)/(m*(n_1+n_2-2))                                      # \nu
alpha_0     = 0.05                                                               # Signifikanzlevel
k_alpha_0   = (1/nu)*qf(1-alpha_0,m,n_1+n_2-m-1)                                 # kritischer Wert
t2          = seq(0,16,length=1e4)                                               # T^2-Statistikwerte
Pt2         = pf(nu*t2,m,n-m)                                                    # T^2-Statistik KVF für H_0
pt2         = nu*df(nu*t2,m,n-m)                                                 # T^2-Statistik WDF für H_0

# KVF Perspektive
plot(                                                                            # original density function
t2,
Pt2,
type        = "l",
ylab        = " ",
xlab        = TeX("$t^2$"),
ylim        = c(0,1),
main        = TeX("$P(t^2) = F(\\nu t^2; m, n_1 + n_2 -m - 1)$"))

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

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

arrows(
x0          = min(t2),
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, -.1 , TeX("$\\k_{\\alpha_0}$"), xpd = TRUE)
text(1.5       , 1.0 , TeX("$1 - \\alpha_0$"), xpd = TRUE)


# WDF Perspektive
plot(
t2,
pt2,
type        = "l",
ylab        = " ",
xlab        = TeX("$t^2$"),
ylim        = c(0,.4),
main        = TeX("$p(t^2) = \\nu f(\nu t^2; m, n_1 + n_2 -m - 1)$"))

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


lines(
seq(k_alpha_0, max(t2), 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,-.04, TeX("$\\k_{\\alpha_0}$") , xpd = TRUE)
text(10, .05, TeX("$P(T^2 >= k_{\\alpha_0}) = \\alpha_0$"), xpd = TRUE, cex = .8, col = "gray50")

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

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

# Zweistichproben-T$^2$-Tests
\noindent (6) Testumfangkontrolle
\vspace{1mm}

\tiny
\setstretch{.8}
```{r}
# Modellparameter bei Zutreffen von H0
m     = 2                                                          # Dimensionalität der Zufallsvektoren/Dat
n     = 15                                                         # Anzahl Datenpunkte pro Stichprobe
n_1   = n                                                          # Anzahl Datenpunkte Stichprobe 1
n_2   = n                                                          # Anzahl Datenpunkte Stichprobe 1
mu_1  = matrix(c(1,1) , nrow = 2)                                  # wahrer, aber unbekannter, Erwartungswertparameter
mu_2  = mu_1                                                       # wahrer, aber unbekannter, Erwartungswertparameter
Sigma = matrix(c(0.4,0.1, 0.1,0.4), nrow = 2, byrow = TRUE)        # wahrer, aber unbekannter, Kovarianzmatrixparameter

# Testparameter
nu          = (n_1+n_2-m-1)/(m*(n_1+n_2-2))                        # \nu
alpha_0     = 0.05                                                 # Signifikanzlevel
k_alpha_0   = (1/nu)*qf(1-alpha_0,m,n_1+n_2-m-1)                   # kritischer Wert

# Simulation
library(MASS)                                                      # R Paket für multivariate Normalverteilungen
library(matlib)                                                    # R Paket für Matrizenrechnung
nsim  = 1e4                                                        # Anzahl Simulationen/Datensatzrealisierungen
Y1b   = matrix(rep(NaN,m*nsim), nrow = 2)                          # Stichprobenmittelarray
Y2b   = matrix(rep(NaN,m*nsim), nrow = 2)                          # Stichprobenmittelarray
phi   = rep(NaN,nsim)                                              # phi Array
j_n   = matrix(rep(1,n), nrow = n)                                 # 1_n
I_n   = diag(n)                                                    # I_n
J_n   = matrix(rep(1,n^2), nrow = n)                               # 1_{nn}
for(s in 1:nsim){                                                  # Simulationsiterationen
    Y_1     = t(mvrnorm(n,mu_1,Sigma))                             # Y_1j \sim N(\mu_1,\Sigma), j = 1,...,n_1
    Y_2     = t(mvrnorm(n,mu_2,Sigma))                             # Y_2j \sim N(\mu_2,\Sigma), j = 1,...,n_2
    Y_1_bar = (1/n)*(Y_1 %*% j_n)                                  # Stichprobenmittel
    Y_2_bar = (1/n)*(Y_2 %*% j_n)                                  # Stichprobenmittel
    C_1     = (1/(n_1-1))*(Y_1 %*% (I_n-(1/n)*J_n) %*% t(Y_1))     # Stichprobenkovarianzmatrix
    C_2     = (1/(n_2-1))*(Y_2 %*% (I_n-(1/n)*J_n) %*% t(Y_2))     # Stichprobenkovarianzmatrix
    C       = ((n_1-1)*C_1+(n_2-1)*C_2)/(n_1+n_2-2)                # Gepoolte Stichprobenkovarianzmatrix
    T2      = (((n_1*n_2)/(n_1+n_2))*                              # T^2 Statistik
              t(Y_1_bar-Y_2_bar) %*% inv(C) %*% (Y_1_bar-Y_2_bar))
    if(T2 > k_alpha_0){phi[s] = 1} else {phi[s] = 0}              # Evaluation von phi
}
cat("\nKritischer Wert              =", k_alpha_0,
    "\nGeschätzter Testumfang alpha =", mean(phi))
```
# Einstichproben-T$^2$-Tests
\noindent (6) Testumfangkontrolle

\small
Praktisches Vorgehen

* \justifying Man nimmt an, dass zwei vorliegende Datensätze
$y_{11},...,y_{1n_1}$ und $y_{21},...,y_{2n_2}$ Realisationen von $m$-dimensionalen
Zufallsvektoren
\begin{equation}
Y_{11},...,Y_{1n_1} \sim N(\mu_1,\Sigma)
\mbox{ und }
Y_{21},...,Y_{2n_2} \sim N(\mu_2,\Sigma)
\end{equation}
mit unbekannten Parametern $\mu_1,\mu_2,\Sigma$ sind.

* Man möchte entscheiden, ob eher $H_0 : \mu_1 = \mu_2$ oder
$H_1: \mu_1 \neq \mu_2$ zutrifft.

* Man wählt ein Signifikanzlevel $\alpha_0$ und bestimmt den zugehörigen kritischen
Wert $k_{\alpha_0}$.

* Anhand von $m, n_1,n_2,\bar{Y}_1,\bar{Y}_2$ und der gepoolten
Stichprobenstandardabweichung $C$ berechnet man die Realisierung der
Zweistichproben-T$^2$-Teststatistik $t^2$.

* Wenn $t^2$ größer als $k_{\alpha_0}$ ist lehnt man die Nullhypothese ab,
andernfalls lehnt man sie nicht ab.

* Die oben entwickelte Theorie des Zweistichproben-T-$^2$Tests garantiert dann,
dass man in höchstens $\alpha_0 \cdot 100$ von $100$ Fällen die Nullhypothese
fälschlicherweise ablehnt.

# Zweistichproben-T$^2$-Tests
\noindent (7) p-Werte

\small
Bestimmung des p-Wertes

* \justifying Per Definition ist der p-Wert das kleinste Signifikanzlevel $\alpha_0$, bei welchem
man die Nullhypothese basierend auf einem vorliegenden Wert der Teststatistik ablehnen würde.

* Bei $T^2 = t^2$ würde $H_0$ für jedes $\alpha_0$ mit $t^2 \ge \frac{1}{\nu}F^{-1}(1-\alpha_0;m,n_1+n_2-m-1)$
abgelehnt werden. Für diese $\alpha_0$ gilt, wie analog im Eintsichprobenstzenario gezeigt
\begin{equation}
\alpha_0 \ge \mathbb{P}\left(T^2 \ge t^2\right)
\end{equation}

* Das kleinste $\alpha_0 \in [0,1]$ mit $\alpha_0 \ge \mathbb{P}\left(T^2 \ge t^2\right)$
ist dann $\alpha_0 = \mathbb{P}(T^2 \ge t^2)$, also folgt
\begin{equation}
\mbox{ p-Wert } = \mathbb{P}\left(T^2 \ge t^2\right) = 1 - F(\nu t^2;m,n_1+n_2-m-1)
\end{equation}

```{r, echo = F, eval = F}
# Beispiel p-Werte
m   = c(2,2,2,4)
n_1 = c(15,15,99,15)
n_2 = c(15,15,99,15)
t2  = c(5,7,5,7)
p   = rep(NaN, length(m))
for(i in 1:length(m)){
    nu   = (n_1[i]+n_2[i]-m[i]-1)/(m[i]*(n_1[i]+n_2[i]-2))
    p[i] = 1 - pf(nu*t2[i], m[i], n_1[i]+n_2[i]-m[i]-1)
}
print(p)
```

* Zum Beispiel ergibt sich bei
\begin{itemize}
\begin{small}
\item[$\circ$] $m = 2, n_1 = 15$ und $n_2=15$ der p-Wert für $t^2 = 5$ zu 0.11
\item[$\circ$] $m = 2, n_1 = 15$ und $n_2=15$ der p-Wert für $t^2 = 7$ zu 0.05
\item[$\circ$] $m = 2, n_1 = 99$ und $n_2=99$ der p-Wert für $t^2 = 5$ zu 0.09
\item[$\circ$] $m = 4, n_1 = 15$ und $n_2=15$ der p-Wert für $t^2 = 7$ zu 0.21
\end{small}
\end{itemize}

# Zweistichproben-T$^2$-Tests
\small
\vspace{1mm}
Anwendungsbeispiel

\tiny
\setstretch{1}
```{r, message = F}
# R Pakete
library(foreign)                                                 # Dateneinlesen
library(matlib)                                                  # Matrixalgebra

# Datenpräprozessierung
fname     = file.path(getwd(), "10_Daten", "studienerfolg.sav")  # Dateiname
D         = read.spss(fname, to.data.frame = T)                  # Dateneinlesen
Y_1       = t(as.matrix(D[D$Gruppe == "ungenügend", 2:3]))       # Daten Gruppe 1
Y_2       = t(as.matrix(D[D$Gruppe == "gut", 2:3]))              # Daten Gruppe 2

# Testparameter
m         = 2                                                    # Datendimensionalität
n         = 15                                                   # Datenpunkte pro Gruppe
n_1       = n                                                    # Datenpunkte Gruppe 1
n_2       = n                                                    # Datenpunkte Gruppe 2
nu        = (n_1+n_2-m-1)/(m*(n_1+n_2-2))                        # \nu
alpha_0   = 0.05                                                 # Signifikanzlevel
k_alpha_0 = (1/nu)*qf(1-alpha_0,m,n_1+n_2-m-1)                   # kritischer Wert

# Testevaluation
j_n      = matrix(rep(1,n), nrow = n)                            # 1_n
I_n      = diag(n)                                               # I_n
J_n      = matrix(rep(1,n^2), nrow = n)                          # 1_{nn}
Y_1_bar  = (1/n)*(Y_1 %*% j_n)                                   # Stichprobenmittel
Y_2_bar  = (1/n)*(Y_2 %*% j_n)                                   # Stichprobenmittel
C_1      = (1/(n_1-1))*(Y_1 %*% (I_n-(1/n)*J_n) %*% t(Y_1))      # Stichprobenkovarianzmatrix
C_2      = (1/(n_2-1))*(Y_2 %*% (I_n-(1/n)*J_n) %*% t(Y_2))      # Stichprobenkovarianzmatrix
C        = ((n_1-1)*C_1+(n_2-1)*C_2)/(n_1+n_2-2)                 # Gepoolte Stichprobenkovarianzmatrix
T2       = (((n_1*n_2)/(n_1+n_2))*                               # T^2 Statistik
           t(Y_1_bar-Y_2_bar) %*% inv(C) %*% (Y_1_bar-Y_2_bar))
if(T2 > k_alpha_0){phi = 1} else {phi = 0}                       # Test 1_{T^2 > k_alpha_0}
p         = 1 - pf(nu*T2,m,n_1+n_2-m-1)                          # p-Wert
```


# Zweistichproben-T$^2$-Tests
\small
Anwendungsbeispiel

\vspace{1mm}
\tiny
\setstretch{.8}
```{r, message = F}
# Ausgabe
cat("Y_1_bar = ", Y_1_bar,
    "\nY_2_bar = ", Y_2_bar,
    "\nC       = ", C,
    "\nT^2     = ", T2,
    "\nalpha_0 = ", alpha_0,
    "\nk       = ", k_alpha_0,
    "\nphi     = ", phi,
    "\np       = ", p)
```

\small
Anwendungsbeispiel mit `MVTests::TwoSamplesHT2()`
\vspace{1mm}

\tiny
\setstretch{.9}
```{r, message = F}
library(MVTests)                                 # R Pakete
Y     = D[c(1:15, 31:45), 2:3]                   # Dataframe von Interesse
G     = c(rep(1,15), rep(2,15))                  # Gruppenindikatoren
phi   = TwoSamplesHT2(Y,G, alpha_0)              # Zweistichproben-T^2-Test

# Ausgabe
cat("Y_1_bar = " , phi$Descriptive1[1,],
    "\nY_2_bar = ",  phi$Descriptive2[1,],
    "\nT^2     = ", phi$HT2,
    "\nalpha_0 = ", phi$alpha,
    "\nk       = " ,phi$F,
    "\np       = ", phi$p.value)

```

#
\vfill
\setstretch{2.4}
\large
Vorbemerkungen

Einstichproben-T$^2$-Tests

Zweistichproben-T$^2$-Tests

**Univariates vs. multivariates Testen**

Selbstkontrollfragen
\vfill


# Univariates vs. multivariates Testen
\setstretch{1.5}

\textcolor{darkblue}{Gegeben sei ein Anwendungsszenario mit $n$ Beobachtungen von $m$ Variablen}

Bei Durchführung von $m$ univariaten Tests entsteht ein multiples Testproblem

\small
* Induktion multipler Typ I und Typ II Fehlerraten (cf. @ostwald_2019)

* Familywise Error Rate (FWER) = $\mathbb{P}(\ge 1 {\mbox{ Typ I Fehler}})$

* Bei unabhängigen Variablen bietet zur FWER Kontrolle die *Bonferroni Korrektur* an.

* Bei Durchführung eines multivariaten Tests entsteht kein multiples Testproblem.

\normalsize
Multivariate Tests beziehen Variablenkovarianzen explizit mit ein.

\small
* Bei Durchführung von $m$ univariaten Tests werden Variablenkorrelationen aktiv ignoriert.

\normalsize
Sollten $m$ univariate Tests oder ein multivariater Test durchgeführt werden?

\small
* Je nachdem, ob die Daten als multi- oder univariate Realisierung konzipiert werden.

* Je nachdem, welche geometrische Form des Annahmebereiches gewünscht ist.

* Prinzipiell sollte im wissenschaftlichen Diskurs überhaupt nicht getestet werden.

\normalsize
\textcolor{darkblue}{Wir betrachten in der Folge Simulationsszenarien mit $m := 2$}.

# Univariates vs. multivariates Testen
\vspace{1mm}
\textcolor{darkblue}{Univariate vs. multivariate Einstichproben-T$^{(2)}$-Tests mit $\Sigma = \sigma^2I_2$}
\vspace{1mm}

\tiny
\setstretch{1}
```{r}
# R Pakete
library(MASS)                                               # R Paket für multivariate Normalverteilungen
library(matlib)                                             # R Paket für Matrizenrechnung

# Modellparameter
m     = 2                                                   # Dimensionalität der Zufallsvektoren/Daten
n     = 15                                                  # Anzahl der Datenpunkte
mu    = matrix(c(0,0) , nrow = 2)                           # Erwartungswertparameter
Sigma = matrix(c(1,0,0,1), nrow = 2, byrow = TRUE)          # Kovarianzmatrixparameter

# Testparameter
alpha = 0.05                                                # Signifikanzlevel
nu    = (n-m)/(m*(n-1))                                     # T^2-Test Parameter
k_T2  = (1/nu)*qf(1-alpha, m, n-m)                          # T^2-Test kritischer Wert
k_Tu  = qt(1-(1/2)*alpha  ,n-1)                             # T-Test kritischer Wert unkorrigiert
k_Tc  = qt(1-(1/2)*alpha/m,n-1)                             # T-Test kritischer Wert Bonferonni korrigiert

# Simulation der Testumfangkontrolle
nsim  = 1e4                                                 # Anzahl Simulation
phi   = matrix(rep(NaN,nsim*5), nrow = 5)                   # Testentscheidungsarray
j_n   = matrix(rep(1,n), nrow = n)                          # 1_n
I_n   = diag(n)                                             # I_n
J_n   = matrix(rep(1,n^2), nrow = n)                        # 1_{nn}
for(s in 1:nsim){                                           # Simulationsiterationen
    Y        = t(mvrnorm(n,mu,Sigma))                       # Y_i \sim N(\mu,\Sigma), i = 1,...,n
    Y_bar    = (1/n)*(Y %*% j_n)                            # Stichprobenmittel
    C        = (1/(n-1))*(Y %*% (I_n-(1/n)*J_n) %*% t(Y))   # Stichprobenkovarianzmatrix
    phi[1,s] =  n*t(Y_bar) %*% inv(C) %*% (Y_bar)  > k_T2   # T^2-Test mit \mu_0
    for(i in 1:m){                                          # T-Test Iterationen
      y_bar      = Y_bar[i]                                 # Stichprobenmittel
      sigma_hat  = sqrt(C[i,i])                             # Stichprobenstandardabweichung
      phi[i+1,s] = abs(sqrt(n)*y_bar/sigma_hat)  > k_Tu     # Unkorrigierter T-Test mit \mu_0 = 0
      phi[i+3,s] = abs(sqrt(n)*y_bar/sigma_hat)  > k_Tc}}   # Korrigierter   T-Test mit \mu_0 = 0
```

# Univariates vs. multivariates Testen
\vspace{1mm}
\textcolor{darkblue}{Univariate vs. multivariate Einstichproben-T$^{(2)}$-Tests mit $\Sigma = \sigma^2I_2$}
\vspace{2mm}

\small
```{r, echo = F}
cat("Kritischer Wert T^2-Test                                        =", k_T2         ,
    "\nKritischer Wert T-Test                                          =", k_Tu         ,
    "\nKritischer Wert T-Test Bonferroni                               =", k_Tc         ,
    "\n"                                                                                 ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T^2-Test              =", mean(phi[1,]),
    "\n"                                                   ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T-Test 1 unkorrigiert =", mean(phi[2,]),
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T-Test 2 unkorrigiert =", mean(phi[3,]),
    "\nGeschätzte FWER T-Tests unkorrigiert                            =", mean(phi[2,] == 1 | phi[3,] == 1),
    "\n"                                                                                 ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T-Test 1 Bonferroni   =", mean(phi[4,]),
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T-Test 2 Bonferroni   =", mean(phi[5,]),
    "\nGeschätzte FWER T-Tests Bonferroni                              =", mean(phi[4,] == 1 | phi[5,] == 1))
```


# Univariates vs. multivariates Testen
\vspace{1mm}
\textcolor{darkblue}{Univariate vs. multivariate Einstichproben-T$^{(2)}$-Tests mit $\Sigma \neq \sigma^2I_2$}
\vspace{2mm}

\small
```{r, echo = F}
# R Pakete
library(MASS)                                               # R Paket für multivariate Normalverteilungen
library(matlib)                                             # R Paket für Matrizenrechnung

# Modellparameter
m      = 2                                                  # Dimensionalität der Zufallsvektoren/Daten
n      = 15                                                 # Anzahl der Datenpunkte
mu    = matrix(c(0,0) , nrow = 2)                           # Erwartungswertparameter
Sigma = matrix(c(1,-.9,-.9,1), nrow = 2, byrow = TRUE)      # Kovarianzmatrixparameter

# Testparameter
alpha = 0.05                                                # Signifikanzlevel
nu    = (n-m)/(m*(n-1))                                     # T^2-Test Parameter
k_T2  = (1/nu)*qf(1-alpha, m, n-m)                          # T^2-Test kritischer Wert
k_Tu  = qt(1-(1/2)*alpha  ,n-1)                             # T-Test kritischer Wert unkorrigiert
k_Tc  = qt(1-(1/2)*alpha/m,n-1)                             # T-Test kritischer Wert Bonferonni korrigiert

# Simulation der Testumfangkontrolle
nsim  = 1e4                                                 # Anzahl Simulation
phi   = matrix(rep(NaN,nsim*5), nrow = 5)                   # Testentscheidungsarray
j_n   = matrix(rep(1,n), nrow = n)                          # 1_n
I_n   = diag(n)                                             # I_n
J_n   = matrix(rep(1,n^2), nrow = n)                        # 1_{nn}
for(s in 1:nsim){                                           # Simulationsiterationen
    Y        = t(mvrnorm(n,mu,Sigma))                       # Y_i \sim N(\mu,\Sigma), i = 1,...,n
    Y_bar    = (1/n)*(Y %*% j_n)                            # Stichprobenmittel
    C        = (1/(n-1))*(Y %*% (I_n-(1/n)*J_n) %*% t(Y))   # Stichprobenkovarianzmatrix
    phi[1,s] =  n*t(Y_bar) %*% inv(C) %*% (Y_bar)  > k_T2   # T^2-Test mit \mu_0
    for(i in 1:m){                                          # T-Test Iterationen
      y_bar      = Y_bar[i]                                 # Stichprobenmittel
      sigma_hat  = sqrt(C[i,i])                             # Stichprobenstandardabweichung
      phi[i+1,s] = abs(sqrt(n)*y_bar/sigma_hat)  > k_Tu     # Unkorrigierter T-Test mit \mu_0 = 0
      phi[i+3,s] = abs(sqrt(n)*y_bar/sigma_hat)  > k_Tc}}   # Korrigierter   T-Test mit \mu_0 = 0

cat("Kritischer Wert T^2-Test                                        =", k_T2         ,
    "\nKritischer Wert T-Test                                          =", k_Tu         ,
    "\nKritischer Wert T-Test Bonferroni                               =", k_Tc         ,
    "\n"                                                                                 ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T^2-Test              =", mean(phi[1,]),
    "\n"                                                   ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T-Test 1 unkorrigiert =", mean(phi[2,]),
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T-Test 2 unkorrigiert =", mean(phi[3,]),
    "\nGeschätzte FWER T-Tests unkorrigiert                            =", mean(phi[2,] == 1 | phi[3,] == 1),
    "\n"                                                                                 ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T-Test 1 Bonferroni   =", mean(phi[4,]),
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit T-Test 2 Bonferroni   =", mean(phi[5,]),
    "\nGeschätzte FWER T-Tests Bonferroni                              =", mean(phi[5,] == 1 | phi[5,] == 1))
```

\vspace{2mm}

* Kovariabilität von Variablen reduziert die FWER.
* Die Bonferroni FWER Korrektur wird *konservativ*, also $\mathbb{P}(\ge 1 {\mbox{ Typ I Fehler}}) < \alpha_0$.


# Univariates vs. multivariates Testen
\setstretch{1.5}
\textcolor{darkblue}{Univariate vs. multivariate Einstichproben-Z$^{(2)}$-Tests}

\small
Zur Visualisierung von Stichprobenmittel und Testentscheidung bieten sich (nur) $Z^2$-Test an.

Z$^2$-Test $\approx$ T$^2$-Test mit als bekannt vorausgesetzter Kovarianzmatrix bei $Y_i \sim N(\mu,\Sigma)$.

* Einstichproben-$Z^2$-Teststatistik:
\begin{equation}
Z^2 := n(\bar{Y} - \mu_0)^T\Sigma^{-1}(\bar{Y} - \mu_0)
\end{equation}

* Verteilung der Einstichproben-$Z^2$-Teststatistik bei $H_0 : \mu = \mu_0$:
\begin{equation}
Z^2 \sim \chi^2(m)
\end{equation}

* Kritischer Wert für Testumfangkontrolle: 
\begin{equation}
k_{\alpha_0} := \Xi^{2^{-1}}(1-\alpha_0;m)
\end{equation}
wobei $\Xi^{2^{-1}}$ die inversen KVF der $\chi^2$ Verteilung bezeichnet. 

# Univariates vs. multivariates Testen
\textcolor{darkblue}{Univariate vs. multivariate Einstichproben-Z$^{(2)}$-Tests mit $\Sigma = \sigma^2I_2$}

\setstretch{1}
\tiny
```{r}
# R Pakete
library(MASS)                                               # R Paket für multivariate Normalverteilungen
library(matlib)                                             # R Paket für Matrizenrechnung

# Modellparameter
m     = 2                                                   # Dimensionalität der Zufallsvektoren/Daten
n     = 15                                                  # Anzahl der Datenpunkte
mu    = matrix(c(0,0) , nrow = 2)                           # Erwartungswertparameter 
Sigma = matrix(c(1,0,0,1), nrow = 2, byrow = TRUE)          # Kovarianzmatrixparameter

# Testparameter
alpha = 0.05                                                # Signifikanzlevel 
k_Z2  = qchisq(1-alpha, m)                                  # Z^2-Test kritischer Wert
k_Zu  = qnorm(1-(1/2)*alpha)                                # Z-Test kritischer Wert unkorrigiert
k_Zc  = qnorm(1-(1/2)*alpha/m)                              # Z-Test kritischer Wert Bonferonni korrigiert

# Simulation der Testumfangkontrolle
nsim  = 2e3                                                 # Anzahl Simulation
YB    = matrix(rep(NaN,nsim*2), nrow  = 2)                  # Stichprobenmittelarray
phi   = matrix(rep(NaN,nsim*5), nrow  = 5)                  # Testentscheidungsarray
j_n   = matrix(rep(1,n), nrow = n)                          # 1_n
I_n   = diag(n)                                             # I_n
J_n   = matrix(rep(1,n^2), nrow = n)                        # 1_{nn}
for(s in 1:nsim){                                           # Simulationsiterationen
    Y        = t(mvrnorm(n,mu,Sigma))                       # Y_i \sim N(\mu,\Sigma), i = 1,...,n
    YB[,s]   = (1/n)*(Y %*% j_n)                            # Stichprobenmittel
    phi[1,s] = (n*t(YB[,s])%*%inv(Sigma)%*%YB[,s]) > k_Z2   # T^2-Test mit \mu_0 
    for(i in 1:m){                                          # T-Test Iterationen
      y_bar      = YB[i,s]                                  # Stichprobenmittel
      sigma      = sqrt(Sigma[i,i])                         # Stichprobenstandardabweichung
      phi[i+1,s] = abs(sqrt(n)*y_bar/sigma)  > k_Zu         # Unkorrigierter T-Test mit \mu_0 = 0
      phi[i+3,s] = abs(sqrt(n)*y_bar/sigma)  > k_Zc}}       # Korrigierter   T-Test mit \mu_0 = 0
```

```{r, echo = F, eval = F}
cat("Kritischer Wert Z^2-Test                                        =", k_Z2         ,
    "\nKritischer Wert Z-Test                                          =", k_Zu         ,
    "\nKritischer Wert Z-Test Bonferroni                               =", k_Zc         ,
    "\n"                                                                                 , 
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z^2-Test              =", mean(phi[1,]),
    "\n"                                                   ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z-Test 1 unkorrigiert =", mean(phi[2,]),
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z-Test 2 unkorrigiert =", mean(phi[3,]),
    "\nGeschätzte FWER Z-Tests unkorrigiert                            =", mean(phi[2,] == 1 | phi[3,] == 1),
    "\n"                                                                                 ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z-Test 1 Bonferroni   =", mean(phi[4,]),
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z-Test 2 Bonferroni   =", mean(phi[5,]),
    "\nGeschätzte FWER Z-Tests Bonferroni                              =", mean(phi[4,] == 1 | phi[5,] == 1))
```
# Univariates vs. multivariates Testen
\vspace{1mm}
\textcolor{darkblue}{Univariate vs. multivariate Einstichproben-Z$^{(2)}$-Tests mit $\Sigma = \sigma^2I_2$}
\vspace{1mm}

```{r, echo = F, eval = F}
library(ellipse)
library(matlib)
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,1),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 1)

# Z^2 Teststatistik 
res        = 1e2
Y_bar_i    = seq(-1,1,len = res)
Z2         = matrix(rep(NaN, res*res), nrow = res)
for(i in 1:res){
    for(j in 1:res){
        Y_bar   = c(Y_bar_i[i], Y_bar_i[j])
        Z2[i,j] = n*t(Y_bar) %*% inv(Sigma) %*% Y_bar    
    }
}

# Z^2 Test Entscheidungsgrenze
contour(
Y_bar_i,
Y_bar_i,
Z2,
level = k_Z2,
drawlabels = F,
xlab = TeX("$\\bar{Y}_1$"),
ylab = TeX("$\\bar{Y}_2$"),
main = "Kritische Werte im Stichprobenmittelraum")

# Z Test Entscheidungsgrenzen
abline(h =  k_Zu/sqrt(n), col = "Gray60", lty = 2)
abline(h = -k_Zu/sqrt(n), col = "Gray60", lty = 2)
abline(v =  k_Zu/sqrt(n), col = "Gray60", lty = 2)
abline(v = -k_Zu/sqrt(n), col = "Gray60", lty = 2)
abline(h =  k_Zc/sqrt(n), col = "Gray60", lty = 1)
abline(h = -k_Zc/sqrt(n), col = "Gray60", lty = 1)
abline(v =  k_Zc/sqrt(n), col = "Gray60", lty = 1)
abline(v = -k_Zc/sqrt(n), col = "Gray60", lty = 1)

legend(
-0.5,
-0.6,
c("Multivariater Test", "Univariater Test, unkorrigiert", "Univariater Test, korrigiert"),
lty = c(1,2,1),
col = c("Black", "Gray60","Gray60"),
bty = "n",
cex = 0.7)

dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_univariate_vs_multivariate_1.pdf"),
width  = 6,
height = 6)
```

```{r, echo = FALSE, out.width="65%"}
knitr::include_graphics("10_Abbildungen/mvda_10_univariate_vs_multivariate_1.pdf")
```

# Univariates vs. multivariates Testen
\vspace{1mm}
\textcolor{darkblue}{Univariate vs. multivariate Einstichproben-Z$^{(2)}$-Tests mit $\Sigma = \sigma^2I_2$}


```{r, echo = F, eval = F}
library(ellipse)
library(matlib)
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,1),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 1)

# Z^2 Teststatistik 
res        = 1e2
Y_bar_i    = seq(-1,1,len = res)
Z2         = matrix(rep(NaN, res*res), nrow = res)
for(i in 1:res){
    for(j in 1:res){
        Y_bar   = c(Y_bar_i[i], Y_bar_i[j])
        Z2[i,j] = n*t(Y_bar) %*% inv(Sigma) %*% Y_bar    
    }
}

# Z^2 Test Entscheidungsgrenze
contour(
Y_bar_i,
Y_bar_i,
Z2,
level = k_Z2,
drawlabels = F,
xlab = TeX("$\\bar{Y}_1$"),
ylab = TeX("$\\bar{Y}_2$"),
main = "(Korrigierte) H0 Ablehnungen im Stichprobenmittelraum")


# Z^2 Test Ablehnungen
points(
YB[1,phi[1,] == 1],
YB[2,phi[1,] == 1],
pch = 21,
col = "White",
bg  = "Gray80",
cex = 2)

# Z Test Ablehnungen
points(
YB[1,phi[4,] == 1],
YB[2,phi[4,] == 1],
pch = 21,
col = "White",
bg  = "Red",
cex = 1)

points(
YB[1,phi[5,] == 1],
YB[2,phi[5,] == 1],
pch = 21,
col = "White",
bg  = "Blue",
cex = 1)

points(
YB[1,phi[4,] == 1 & phi[5,] == 1],
YB[2,phi[4,] == 1 & phi[5,] == 1],
pch = 21,
col = "White",
bg  = "Black",
cex = 1)

# Z Test Entscheidungsgrenzen
abline(h =  k_Zc/sqrt(n))
abline(h = -k_Zc/sqrt(n))
abline(v =  k_Zc/sqrt(n))
abline(v = -k_Zc/sqrt(n))


dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_univariate_vs_multivariate_2.pdf"),
width  = 6,
height = 6)
```

```{r, echo = FALSE, out.width="65%"}
knitr::include_graphics("10_Abbildungen/mvda_10_univariate_vs_multivariate_2.pdf")
```
\vspace{-2mm}
\footnotesize
\center
\textcolor{lightgray}{$\bullet$} $\phi(Y) = 1$,
\textcolor{red}{$\bullet$} $\phi(Y_1) = 1$,
\textcolor{blue}{$\bullet$} $\phi(Y_2) = 1$,
\textcolor{black}{$\bullet$} $\phi(Y_1) = 1 \mbox{ und }\phi(Y_2) = 1$

# Univariates vs. multivariates Testen
\vspace{1mm}
\textcolor{darkblue}{Univariate vs. multivariate Einstichproben-Z$^{(2)}$-Tests mit $\Sigma \neq \sigma^2I_2$}
\vspace{1mm}

```{r, eval = F, echo = F}
# R Pakete
library(MASS)                                               # R Paket für multivariate Normalverteilungen
library(matlib)                                             # R Paket für Matrizenrechnung

# Modellparameter
m     = 2                                                   # Dimensionalität der Zufallsvektoren/Daten
n     = 15                                                  # Anzahl der Datenpunkte
mu    = matrix(c(0,0) , nrow = 2)                           # Erwartungswertparameter 
Sigma = matrix(c(1,.9,.9,1), nrow = 2, byrow = TRUE)          # Kovarianzmatrixparameter

# Testparameter
alpha = 0.05                                                # Signifikanzlevel 
k_Z2  = qchisq(1-alpha, m)                                  # Z^2-Test kritischer Wert
k_Zu  = qnorm(1-(1/2)*alpha)                                # Z-Test kritischer Wert unkorrigiert
k_Zc  = qnorm(1-(1/2)*alpha/m)                              # Z-Test kritischer Wert Bonferonni korrigiert

# Simulation der Testumfangkontrolle
nsim  = 2e3                                                 # Anzahl Simulation
YB    = matrix(rep(NaN,nsim*2), nrow  = 2)                  # Stichprobenmittelarray
phi   = matrix(rep(NaN,nsim*5), nrow  = 5)                  # Testentscheidungsarray
j_n   = matrix(rep(1,n), nrow = n)                          # 1_n
I_n   = diag(n)                                             # I_n
J_n   = matrix(rep(1,n^2), nrow = n)                        # 1_{nn}
for(s in 1:nsim){                                           # Simulationsiterationen
    Y        = t(mvrnorm(n,mu,Sigma))                       # Y_i \sim N(\mu,\Sigma), i = 1,...,n
    YB[,s]   = (1/n)*(Y %*% j_n)                            # Stichprobenmittel
    phi[1,s] = (n*t(YB[,s])%*%inv(Sigma)%*%YB[,s]) > k_Z2   # T^2-Test mit \mu_0 
    for(i in 1:m){                                          # T-Test Iterationen
      y_bar      = YB[i,s]                                  # Stichprobenmittel
      sigma      = sqrt(Sigma[i,i])                         # Stichprobenstandardabweichung
      phi[i+1,s] = abs(sqrt(n)*y_bar/sigma)  > k_Zu         # Unkorrigierter T-Test mit \mu_0 = 0
      phi[i+3,s] = abs(sqrt(n)*y_bar/sigma)  > k_Zc}}       # Korrigierter   T-Test mit \mu_0 = 0
```

```{r, echo = F, eval = F}
cat("Kritischer Wert Z^2-Test                                        =", k_Z2         ,
    "\nKritischer Wert Z-Test                                          =", k_Zu         ,
    "\nKritischer Wert Z-Test Bonferroni                               =", k_Zc         ,
    "\n"                                                                                 , 
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z^2-Test              =", mean(phi[1,]),
    "\n"                                                   ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z-Test 1 unkorrigiert =", mean(phi[2,]),
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z-Test 2 unkorrigiert =", mean(phi[3,]),
    "\nGeschätzte FWER Z-Tests unkorrigiert                            =", mean(phi[2,] == 1 | phi[3,] == 1),
    "\n"                                                                                 ,
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z-Test 1 Bonferroni   =", mean(phi[4,]),
    "\nGeschätzte Typ I Fehlerwahrscheinlichkeit Z-Test 2 Bonferroni   =", mean(phi[5,]),
    "\nGeschätzte FWER Z-Tests Bonferroni                              =", mean(phi[4,] == 1 | phi[5,] == 1))
```


```{r, echo = F, eval = F}
library(ellipse)
library(matlib)
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,1),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 1)

# Z^2 Teststatistik 
res        = 1e2
Y_bar_i    = seq(-1,1,len = res)
Z2         = matrix(rep(NaN, res*res), nrow = res)
for(i in 1:res){
    for(j in 1:res){
        Y_bar   = c(Y_bar_i[i], Y_bar_i[j])
        Z2[i,j] = n*t(Y_bar) %*% inv(Sigma) %*% Y_bar    
    }
}

# Z^2 Test Entscheidungsgrenze
contour(
Y_bar_i,
Y_bar_i,
Z2,
level = k_Z2,
drawlabels = F,
xlab = TeX("$\\bar{Y}_1$"),
ylab = TeX("$\\bar{Y}_2$"),
main = "Kritische Werte im Stichprobenmittelraum")

# Z Test Entscheidungsgrenzen
abline(h =  k_Zu/sqrt(n), col = "Gray60", lty = 2)
abline(h = -k_Zu/sqrt(n), col = "Gray60", lty = 2)
abline(v =  k_Zu/sqrt(n), col = "Gray60", lty = 2)
abline(v = -k_Zu/sqrt(n), col = "Gray60", lty = 2)
abline(h =  k_Zc/sqrt(n), col = "Gray60", lty = 1)
abline(h = -k_Zc/sqrt(n), col = "Gray60", lty = 1)
abline(v =  k_Zc/sqrt(n), col = "Gray60", lty = 1)
abline(v = -k_Zc/sqrt(n), col = "Gray60", lty = 1)

legend(
-0.5,
-0.6,
c("Multivariater Test", "Univariater Test, unkorrigiert", "Univariater Test, korrigiert"),
lty = c(1,2,1),
col = c("Black", "Gray60","Gray60"),
bty = "n",
cex = 0.7)

dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_univariate_vs_multivariate_3.pdf"),
width  = 6,
height = 6)
```

```{r, echo = FALSE, out.width="65%"}
knitr::include_graphics("10_Abbildungen/mvda_10_univariate_vs_multivariate_3.pdf")
```

# Univariates vs. multivariates Testen
\vspace{1mm}
\textcolor{darkblue}{Univariate vs. multivariate Einstichproben-Z$^{(2)}$-Tests mit $\Sigma = \sigma^2I_2$}


```{r, echo = F, eval = F}
library(ellipse)
library(matlib)
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,1),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 1)

# Z^2 Teststatistik 
res        = 1e2
Y_bar_i    = seq(-1,1,len = res)
Z2         = matrix(rep(NaN, res*res), nrow = res)
for(i in 1:res){
    for(j in 1:res){
        Y_bar   = c(Y_bar_i[i], Y_bar_i[j])
        Z2[i,j] = n*t(Y_bar) %*% inv(Sigma) %*% Y_bar    
    }
}

# Z^2 Test Entscheidungsgrenze
contour(
Y_bar_i,
Y_bar_i,
Z2,
level = k_Z2,
drawlabels = F,
xlab = TeX("$\\bar{Y}_1$"),
ylab = TeX("$\\bar{Y}_2$"),
main = "(Korrigierte) H0 Ablehnungen im Stichprobenmittelraum")


# Z^2 Test Ablehnungen
points(
YB[1,phi[1,] == 1],
YB[2,phi[1,] == 1],
pch = 21,
col = "White",
bg  = "Gray80",
cex = 2)

# Z Test Ablehnungen
points(
YB[1,phi[4,] == 1],
YB[2,phi[4,] == 1],
pch = 21,
col = "White",
bg  = "Red",
cex = 1)

points(
YB[1,phi[5,] == 1],
YB[2,phi[5,] == 1],
pch = 21,
col = "White",
bg  = "Blue",
cex = 1)

points(
YB[1,phi[4,] == 1 & phi[5,] == 1],
YB[2,phi[4,] == 1 & phi[5,] == 1],
pch = 21,
col = "White",
bg  = "Black",
cex = 1)

# Z Test Entscheidungsgrenzen
abline(h =  k_Zc/sqrt(n))
abline(h = -k_Zc/sqrt(n))
abline(v =  k_Zc/sqrt(n))
abline(v = -k_Zc/sqrt(n))


dev.copy2pdf(
file   = file.path(getwd(), "10_Abbildungen", "mvda_10_univariate_vs_multivariate_4.pdf"),
width  = 6,
height = 6)
```

```{r, echo = FALSE, out.width="65%"}
knitr::include_graphics("10_Abbildungen/mvda_10_univariate_vs_multivariate_4.pdf")
```
\vspace{-2mm}
\footnotesize
\center
\textcolor{lightgray}{$\bullet$} $\phi(Y) = 1$,
\textcolor{red}{$\bullet$} $\phi(Y_1) = 1$,
\textcolor{blue}{$\bullet$} $\phi(Y_2) = 1$,
\textcolor{black}{$\bullet$} $\phi(Y_1) = 1 \mbox{ und }\phi(Y_2) = 1$

#
\vfill
\setstretch{2.4}
\large
Vorbemerkungen

Einstichproben-T$^2$-Tests

Zweistichproben-T$^2$-Tests

Univariates vs. multivariates Testen

**Selbstkontrollfragen**
\vfill

# Selbstkontrollfragen
\footnotesize
\setstretch{1.2}
\begin{enumerate}
\justifying
\itemsep0mm
\item Erläutern Sie das Modell und die Standardprobleme der Frequentistischen Inferenz.
\item Erläutern Sie die Standardannahmen Frequentistischer Inferenz.
\item Definieren Sie den Begriff der Mahalanobis Distanz.
\item Erläutern Sie den Unterschied zwischen einer $f$-Verteilung und einer nichtzentralen $f$-Verteilung.
\item Beschreiben Sie das Anwendungsszenario für einen Einstichproben-T$^2$-Test.
\item Geben Sie das statistische Modell eines Einstichproben-T$^2$-Test in klassischer Form an.
\item Geben Sie das statistische Modell eines Einstichproben-T$^2$-Test in generativer Form an.
\item Erläutern Sie das statistische Modell eines Einstichproben-T$^2$-Test in generativer Form.
\item Definieren Sie die Testhypothesen eines Einstichproben-T$^2$-Test mit einfacher Nullhypothese und zusammengesetzter Alternativhypothese.
\item Definieren Sie die Einstichproben-T$^2$-Teststatistik.
\item Erläutern Sie, wann die Einstichproben-T$^2$-Teststatistik hohe Werte annimmt.
\item Geben Sie die WDF und KVF der Einstichproben-T$^2$-Teststatistik an.
\item Geben Sie den kritischen Wert für einen Level-$\alpha_0$-Einstichproben-T$^2$-Testmit einfacher Nullhypothese und zusammengesetzter Alternativhypothese und Testumfang $\alpha_0$ an.
\item Erläutern Sie das praktische Vorgehen bei der Durchführung eines Einstichproben-T$^2$-Tests.
\item Definieren Sie den Begriff des p-Wertes.
\item Geben Sie den p-Wert für einen Einstichproben-T$^2$-Test mit einfacher Nullhypothese und zusammengesetzter Alternativhypothese an und erläutern Sie die Komponenten des entsprechenden Ausdrucks.
\item Definieren Sie die Powerfunktion eines Einstichproben-T$^2$-Test mit einfacher Nullhypothese und zusammengesetzter Alternativhypothese und erläutern Sie ihre Komponenten.
\end{enumerate}



# References
\footnotesize