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


```{r, include = F}
source("7_R_common.R")
fdir        = file.path(getwd(), "7_Abbildungen")                                # Abbildungsverzeichnis
```


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

\vspace{2mm}

\huge
Allgemeines Lineares Modell
\vspace{6mm}

\large
BSc Psychologie SoSe 2022

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

#  {.plain}
\center
\huge
\vfill
\noindent (7) Modellevaluation
\vfill

# Überblick

\large

Naturwissenschaft \vspace{7mm}

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

# Überblick
\vspace{1mm}
\normalsize
Modellformulierung
\vspace{1mm}
\small
\begin{equation}
y = X\beta + \varepsilon, \varepsilon \sim N(0_n,\sigma^2I_n)
\end{equation}
\vspace{5mm}

\normalsize
Modellschätzung
\small
\begin{equation}
\hat{\beta} = (X^TX)^{-1} X^Ty,  \hat{\sigma}^2 = \frac{(y - X\hat{\beta})^T(y - X\hat{\beta})}{n-p}
\end{equation}
\vspace{4mm}

\normalsize
Modellevaluation
\small
\begin{equation}
T = \frac{c^T\hat{\beta} - c^T\beta_0}{\sqrt{\hat{\sigma}^2c^T(X^TX)^{-1}c}}, 
F = \frac{(\hat{\varepsilon}_1^T\hat{\varepsilon}_1 - \hat{\varepsilon}^T\hat{\varepsilon})/p_2}{\hat{\varepsilon}^T\hat{\varepsilon}/(n-p)}
\end{equation}

# Überblick
Standardprobleme Frequentistischer Inferenz
\small

\noindent (1) Parameterschätzung

Ziel der Parameterschätzung ist es, einen möglichst guten Tipp für die wahren,
aber unbekannten, Parameterwerte (oder eine Funktion derer) abzugeben,
typischerweise basierend auf der Beobachtung einer Datenrealisierung.
\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 Daten in einer möglichst sinnvollen Form zu entscheiden, ob ein
wahrer, aber unbekannter Parameterwert, sich in einer von zwei sich gegenseitig 
ausschließenden Untermengen des Parameterraumes, welche man als Hypothesen bezeichnet, 
liegt.


#  {.plain}
\center
```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("7_Abbildungen/alm_7_frequentistische_inferenz.pdf")
```
\center
\footnotesize
$\theta := (\beta,\sigma^2)$,
$\Theta := \mathbb{R}^p \times \mathbb{R}_{>0}$
$\mathbb{P}_\theta(y) := \mathbb{P}_{\beta,\sigma^2}(y)$ mit WDF $p_{\beta,\sigma^2}(\upsilon) := N(\upsilon;X\beta,\sigma^2I_n)$


# Überblick
\small
Standardannahmen Frequentistischer Inferenz

\footnotesize
\setstretch{1.2}
Gegeben sei ein statistisches Modell. Es wird angenommen, dass ein 
vorliegender Datensatz eine der möglichen Realisierungen der Daten des Modells ist.
Aus Frequentistischer Sicht kann man unendlich oft Datensätze basierend auf einem
Modell generieren und zu jedem Datensatz Schätzer oder Statistiken auswerten, z.B.
den Betaparameterschätzer
\vspace{1mm}

\begin{itemize}
\item[] Datensatz (1) : $y^{(1)} = \left(y_1^{(1)}, y_2^{(1)}, ...,y_n^{(1)}\right)^T$ 	mit $\hat{\beta}^{(1)} = (X^TX)^{-1}X^Ty^{(1)}$
\item[] Datensatz (2) : $y^{(2)} = \left(y_1^{(2)}, y_2^{(2)}, ...,y_n^{(2)}\right)^T$ 	mit $\hat{\beta}^{(2)} = (X^TX)^{-1}X^Ty^{(2)}$
\item[] Datensatz (3) : $y^{(3)} = \left(y_1^{(3)}, y_2^{(3)}, ...,y_n^{(3)}\right)^T$ 	mit $\hat{\beta}^{(3)} = (X^TX)^{-1}X^Ty^{(3)}$
\item[] Datensatz (4) : $y^{(4)} = \left(y_1^{(4)}, y_2^{(4)}, ...,y_n^{(4)}\right)^T$ 	mit $\hat{\beta}^{(4)} = (X^TX)^{-1}X^Ty^{(4)}$
\item[] Datensatz (5) : $y^{(5)} = ...$
\end{itemize}

\vspace{1mm}
Um die Qualität statistischer Methoden zu beurteilen betrachtet die Frequentistische
Statistik die Wahrscheinlichkeitsverteilungen von Schätzern und Statistiken
unter Annahme der Datenverteilung. Was zum Beispiel ist die Verteilung von
$\hat{\beta}^{(1)}$, $\hat{\beta}^{(2)}$, $\hat{\beta}^{(3)}$, $\hat{\beta}^{(4)}$, ... also die
Verteilung der Zufallsvariable $\hat{\beta} := (X^TX)^{-1}X^Ty$? Wenn eine statistische 
Methode im Sinne der Frequentistischen Standardannahmen "gut" ist, dann heißt das 
also, dass sie bei häufiger Anwendung "im Mittel gut" ist. Im Einzelfall, also 
im Normalfall nur eines vorliegenden Datensatzes, kann sie auch "schlecht" sein.




#  {.plain}
\vfill
\large
\setstretch{2.5}
Schätzerverteilungen

T-Statistiken

F-Statistiken

Selbstkontrollfragen
\vfill

#  {.plain}
\vfill
\large
\setstretch{2.5}
**Schätzerverteilungen**

T-Statistiken

F-Statistiken

Selbstkontrollfragen
\vfill

# Schätzerverteilungen
Überblick
\setstretch{2}
\footnotesize

* \justifying Unter den Standadannahmen der Frequentistischen Inferenz kann man fragen, wie $\hat{\beta}$ und $\hat{\sigma}^2$ verteilt sind.
* Das erste Theorem dieses Abschnitts besagt, dass $\hat{\beta}$ multivariat normalverteilt ist
* Das zweite Theorem dieses Abschnitts besagt, dass eine skalierte Version von $\hat{\sigma}^2$ $\chi^2$-verteilt ist.
* Diese Ergebnisse sind sehr ähnlich zu den Ergebnissen aus (7) Transformationen der Normalverteilung in WTFI.
* T- und F-Statistiken sind Funktionen von $\hat{\beta}$ und $\hat{\sigma}^2$.
* Die Verteilungen von $\hat{\beta}$ und $\hat{\sigma}^2$ bilden die Grundlage der Verteilungen der T- und F-Statistiken.
* In der konkreten Datenanalyse eines Datensatzes $y$ berechnet man genau ein $\hat{\beta}$ und ein $\hat{\sigma}^2$.
* Probabilistische Aussagen über die aus diesen $\hat{\beta}$ und $\hat{\sigma}^2$ berechneten T- und F-Statistiken
beziehen sich auf die Standardannahmen Frequentistischer Inferenz und weren vor diesem Hintergrund verstanden.


# Schätzerverteilungen
\footnotesize
\begin{theorem}[Frequentistische Verteilung des Betaparameterschätzers]
\justifying
\normalfont
Es sei
\begin{equation}
y = X\beta + \varepsilon \mbox{ mit } \varepsilon \sim N(0_n,\sigma^2I_n)
\end{equation}
das ALM in generativer Form. Weiterhin sei
\begin{equation}
\hat{\beta} := (X^TX)^{-1}X^Ty
\end{equation}
der Betaparameterschätzer. Dann gilt
\begin{equation}
\hat{\beta} \sim N(\beta,\sigma^2(X^T X)^{-1}).
\end{equation}
\end{theorem}

Bemerkungen

* \justifying Es gilt also wie bereits gesehen $\mathbb{E}(\hat{\beta}) = \beta$ und außerdem $\mathbb{C}(\hat{\beta}) = \sigma^2(X^TX)^{-1}$.
* Die Varianzen der Komponenten von $\hat{\beta}$ sind die Diagonalelemente von $\mathbb{C}(\hat{\beta})$, also
\begin{equation}
\mathbb{V}(\hat{\beta}_i) = (\sigma^2(X^TX)^{-1})_{ii} \mbox{ für } i = 1,...p.
\end{equation}
* Die Streuung von $\hat{\beta}$ hängt von $\sigma^2$ und der Designmatrix $X$ ab.
$\sigma^2$ ist ein experimentell nicht zu beinflussender wahrer, aber unbekannter, Parameter
$X$ dagegen kann so gewählt werden, um zum Beispiel die Diagonalelemente von
$\mathbb{C}(\hat{\beta})$ bei festem $\sigma^2$ zu minimieren.

# Schätzerverteilungen
\footnotesize
\underline{Beweis}

Das Theorem folgt direkt mit dem Theorem zur linearen Transformation von
multivariaten Normalverteilung aus Einheit (4) Normalverteilungen. Speziell
gilt hier:
\begin{equation}
\hat{\beta}
\sim
N\left(
(X^{T}X)^{-1}X^{T}X\beta,
(X^{T}X)^{-1}X^{T}(\sigma^{2}I_{n})((X^{T}X)^{-1}X^{T})^{T}
\right).
\end{equation}
Dabei können der Erwartungswertparameter wie folgt
\begin{equation}
(X^{T}X)^{-1}X^{T}X\beta = \beta,
\end{equation}
und der Kovarianzmatrixparamter wie folgt
\begin{align}
\begin{split}
(X^{T}X)^{-1}X^{T}(\sigma^{2}I_{n})((X^{T}X)^{-1}X^{T})^{T}
&  =(X^{T}X)^{-1}X^{T}(\sigma^{2}I_{n})X(X^{T}X)^{-1} 				\\
&  =\sigma^{2}(X^{T}X)^{-1}X^{T}X(X^{T}X)^{-1}						\\
&  =\sigma^{2}(X^{T}X)^{-1}.
\end{split}
\end{align}
Dabei hier die erste Gleichung aus der Tatsache, dass sowohl $X^{T}X$ als auch
ihre Inverse $(X^{T}X)^{-1}$ symmetrische Matrizen sind. Damit folgt dann aber sofort
\begin{equation}
\hat{\beta} \sim N(\beta,\sigma^2(X^T X )^{-1}).
\end{equation}

# Schätzerverteilungen
Beispiel (1) Unabhängige und identisch normalverteilte Zufallsvariablen
\vspace{1mm}

\footnotesize
Es sei
\begin{equation}
y \sim N(X\beta,\sigma^2 I_n)
\mbox{ mit }
X := 1_n \in \mathbb{R}^{n\times 1},
\beta := \mu \in \mathbb{R}
\mbox{ und } \sigma^2 > 0.
\end{equation}
das ALM Szenario unabhängiger und identisch normalverteilter Zufallsvariablen.
Wir haben bereits gesehen, dass $\hat{\beta} = \bar{y}$. Das Theorem zur Frequentistischen
Verteilung des Betaparameterschätzers impliziert damit
\begin{equation}
\bar{y} \sim N\left(\mu, \frac{\sigma^2}{n}\right).
\end{equation}
Das Stichprobenmittel von $n$ unabhängigen und identisch normalverteilten Zufallsvariablen
mit Erwartungswertparameter $\mu$ und Varianzparameter $\sigma^2$ ist also normalverteilt
mit Erwartungswertparameter $\mu$ und Varianzparameter  $\sigma^2/n$. Wir haben diese
Tatsache bereits in Einheit (7) Transformationen der Normalverteilungen in Wahrscheinlichkeitstheorie
und Frequentistische Inferenz unter dem Begriff der Mittelwertstransformation kennengelernt.

# Schätzerverteilungen
Beispiel (1) Unabhängige und identisch normalverteilte Zufallsvariablen
\vspace{2mm}

\setstretch{1}
\footnotesize
```{r}
# Multivariate Normalverteilung
library(MASS)

# Modellformulierung
n        = 12                                      # Anzahl von Datenpunkten
p        = 1                                       # Anzahl von Betparametern
X        = matrix(rep(1,n), nrow = n)              # Designmatrix
I_n      = diag(n)                                 # n x n Einheitsmatrix
beta     = 2                                       # wahrer, aber unbekannter, Betaparameter
sigsqr   = 1                                       # wahrer, aber unbekannter, Varianzparameter

# Frequentistische Simulation
nsim     = 1e4                                     # Anzahl Realisierungen des n-dimensionalen ZVs
beta_hat = rep(NaN,nsim)                           # \hat{\beta} Realisierungsarray
for(i in 1:nsim){
  y           = mvrnorm(1, X %*% beta, sigsqr*I_n) # eine Realisierung eines n-dimensionalen ZVs
  beta_hat[i] = solve(t(X) %*% X) %*% t(X) %*% y   # \hat{\beta} = (X^T)X^{-1}X^Ty
}
```

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

# figure setup
dev.new()
fig = par(
family      = "sans",
mfcol       = c(1,1),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(3,1,0),
xaxs        = "i",
yaxs        = "i",
font.main   = 1,
cex         = 1,
cex.main    = 1)

# density
b_hat_min  = 0.5
b_hat_max  = 3.5
b_hat_res  = 1e3
b_hat      = seq(b_hat_min, b_hat_max, len = b_hat_res)
p_beta_hat = dnorm(b_hat, beta, sqrt(sigsqr/n))

# histogram
hist(
beta_hat,
breaks = 50,
col   = "gray90",
prob  = TRUE,
xlab  = "",
ylab  = "",
xlim  = c(0.5,3.5),
ylim  = c(0,1.5),
main  = "")

# density
lines(
b_hat,
p_beta_hat,
lwd   = 2,
col   = "darkorange")

# print
dev.copy2pdf(
file   = file.path(getwd(), "7_Abbildungen", "alm_7_beta_hat_1.pdf"),
width  = 5,
height = 5)

```
# Schätzerverteilungen
Beispiel (1) Unabhängige und identisch normalverteilte Zufallsvariablen
\vspace{3mm}

\center
\small
$N\left(\hat{\mu}; \mu,\frac{\sigma^2}{n}\right)$
```{r, echo = FALSE, out.width = "50%"}
knitr::include_graphics("7_Abbildungen/alm_7_beta_hat_1.pdf")
```
\vfill

# Schätzerverteilungen
Beispiel (2) Einfache lineare Regression
\vspace{1mm}

\footnotesize
Es sei
\begin{equation}
y \sim N(X\beta, \sigma^2I_n) \mbox{ with}
\begin{pmatrix}
1 		& x_1		\\
\vdots 	& \vdots	\\
1 		& x_n
\end{pmatrix}
\in \mathbb{R}^{n\times 2},
\beta
\in \mathbb{R}^2,
\sigma^2 > 0.
\end{equation}
das ALM Szenario der einfachen linearen Regression. Wir haben bereits gesehen,
dass
\begin{equation}
\sigma^2(X^TX)^{-1}
=
\frac{\sigma^2}{s_{xx}}
\begin{pmatrix}
  \frac{s_{xx}}{n} + \bar{x}^2
& -\bar{x}
\\
  -\bar{x}
&  1
\end{pmatrix}
\mbox{ mit  }
s_{xx} := \sum_{i=1}^n (x_i - \bar{x})^2.
\end{equation}
Die Varianz des Offsetparameterschätzers hängt also sowohl von der Summe der
quadrierten Differenzen und dem Stichprobenmittel der unabhängigen Variablen $x_1,...,x_n$ ab,
wohingegen die Varianz des Steigungsparameterschätzers nur von der Summe der quadrierten
Differenzen der $x_1,...,x_n$ abhängt. Die Kovarianz von Offset- und Steigungsparameterschätzern
hängt vom Mittelwert der $x_1,...,x_n$ ab.

# Schätzerverteilungen
Beispiel (2) Einfache lineare Regression
\vspace{1mm}
\setstretch{1.2}
\footnotesize
```{r}
# Multivariate Normalverteilung und Matrizenrechnung
library(MASS)
library(matlib)

# Modellformulierung
n        = 10                                         # Anzahl von Datenpunkten
p        = 2                                          # Anzahl von Betparametern
x        = 1:n                                        # Prädiktorwerte
X        = matrix(c(rep(1,n),x), nrow = n)            # Designmatrix
I_n      = diag(n)                                    # n x n Einheitsmatrix
beta     = matrix(c(0,1), nrow = p)                   # wahrer,aber unbekannter,Betaparameter
sigsqr   = .5                                         # wahrer,aber unbekannter,Varianzparameter

# Frequentistische Simulation
nsim     = 10                                         # Anzahl Realisierungen n-dimensionaler ZV
y        = matrix(rep(NaN,n*nsim), nrow = n)          # y Realisierungsarray
beta_hat = matrix(rep(NaN,p*nsim), nrow = p)          # \hat{\beta} Realisierungsarray
for(i in 1:nsim){
  y[,i]        = mvrnorm(1, X %*% beta, sigsqr*I_n)   # eine Realisierung n-dimensionaler ZV
  beta_hat[,i] = inv(t(X) %*% X) %*% t(X) %*% y[,i]   # \hat{\beta} = (X^T)X^{-1}X^Ty
}
```

# Schätzerverteilungen
Beispiel (2) Einfache lineare Regression

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

# Graulevel
gfun        = colorRampPalette(c("grey20", "grey80"))
greys       = gfun(nsim)


# Datenrealisierungen
xlimits     = c(0,11)
ylimits     = c(-3,13)
plot(
x,
X %*% beta,
type        = "b",
lty         = 2,
pch         = 1,
col         = "white",
xlab        = "x",
ylab        = "y",
xlim        = xlimits,
ylim        = ylimits)
for(i in 1:nsim){
  lines(
  type = "b",
  x,
  y[,i],
  col = greys[i],
  pch = 16)
}

# Wahrscheinlichkeitsdichtefunktion
b_hat_min  = -2
b_hat_max  = 2
b_hat_res  = 1e3
b_hat_1    = seq(b_hat_min, b_hat_max, len = b_hat_res)
b_hat_2    = seq(b_hat_min, b_hat_max, len = b_hat_res)
b_hat      = expand.grid(b_hat_1,b_hat_2)
mu         = beta
Sigma      = sigsqr*inv(t(X) %*% X)
p_beta_hat = matrix(dmvnorm(as.matrix(b_hat), mu, Sigma), nrow = b_hat_res)

# Visualisierung
contour(
b_hat_1,
b_hat_2,
p_beta_hat,
xlim      = c(-1.5,1.5),
ylim      = c(0.5,1.5),
xlab      = TeX("$\\hat{\\beta}_1$"),
ylab      = TeX("$\\hat{\\beta}_2$"),
nlevels   = 5,
col       = "orange")
for(i in 1:nsim){
  points(
  type = "p",
  beta_hat[1,i],
  beta_hat[2,i],
  col = greys[i],
  pch = 16)
}


# Speichern
dev.copy2pdf(
file        = file.path(fdir, "alm_7_beta_hat_2.pdf"),
width       = 7,
height      = 3.5)
```

\vspace{4mm}
\center
$\quad\quad\quad\quad y \sim (X\beta,\sigma^2I_n)$ \hspace{3cm} $\hat{\beta} \sim N(\beta,\sigma^2(X^TX)^{-1})$

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

# Schätzerverteilungen
\footnotesize
\begin{theorem}[Frequentistische Verteilung des Varianzparameterschätzers]
\normalfont
\justifying
Es sei
\begin{equation}
y = X\beta + \varepsilon \mbox{ mit } \varepsilon \sim N(0_n,\sigma^2I_n)
\end{equation}
das ALM in generative Form. Weiterhin sei
\begin{equation}
\hat{\sigma}^2 = \frac{(y - X\hat{\beta})^T(y - X\hat{\beta})}{n-p}
\end{equation}
der Varianzparameterschätzer. Dann gilt
\begin{equation}
\frac{n-p}{\sigma^2}\hat{\sigma}^2 \sim \chi^2(n-p)
\end{equation}
\end{theorem}

Bemerkungen

* Wir verzichten auf einen Beweis.
* Da es sich bei $(y - X\hat{\beta})^T(y - X\hat{\beta})$ um eine Summe
normalverteilter Zufallsvariablen handelt, liegt die $\chi^2$-Verteilung im Lichte
der $\chi^2$ Transformation aus Einheit (7) Transformationen der Normalverteilung
in Wahrscheinlichkeitstheorie und Frequentistische Inferenz nahe.

# Schätzerverteilungen
Beispiel (1) Unabhängige und identisch normalverteilte Zufallsvariablen
\vspace{1mm}

\footnotesize
Es sei
\begin{equation}
y \sim N(X\beta,\sigma^2 I_n)
\mbox{ mit }
X := 1_n \in \mathbb{R}^{n\times 1},
\beta := \mu \in \mathbb{R}
\mbox{ und } \sigma^2 > 0.
\end{equation}
das ALM Szenario unabhängiger und identisch normalverteilter Zufallsvariablen.
Wir haben bereits gesehen, dass in diesem Fall $\hat{\beta}$ mit dem
Stichprobenmittel $\bar{y}$ identisch ist und dass $\hat{\sigma}^2$ mit der
Stichprobenvarianz $s^2_y$ übereinstimmt.

In Einheit (11) Konfidenzintervalle von  Wahrscheinlichkeitstheorie und
Frequentistische Inferenz hatten wir für den Fall von $n$ unabhängig und identisch
normalverteilten Zufallsvariablen die Statistik
\begin{equation}
U := \frac{n-1}{\sigma^2}s^2
\end{equation}
definiert und festgehalten, dass
\begin{equation}
U \sim \chi^2(n-1).
\end{equation}
Offenbar ist $U$ für $p = 1$ mit der im obigen Theorem betrachten Zufallsvariable
$\frac{n-p}{\sigma^2}\hat{\sigma}^2$ identisch.

# Schätzerverteilungen
Beispiel (2) Einfache lineare Regression
\vspace{2mm}
\setstretch{1.2}
\tiny
```{r}
# Libraries
library(MASS)                                        # multivariate Normalverteilung
library(matlib)                                      # Matrizenrechnung

# Modellformulierung
n          = 10                                      # Anzahl von Datenpunkten
p          = 2                                       # Anzahl von Betparametern
x          = 1:n                                     # Prädiktorwerte
X          = matrix(c(rep(1,n),x), nrow = n)         # Designmatrix
I_n        = diag(n)                                 # n x n Einheitsmatrix
beta       = matrix(c(0,1), nrow = p)                # wahrer,aber unbekannter,Betaparameter
sigsqr     = .5                                      # wahrer,aber unbekannter,Varianzparameter

# Frequentistische Simulation
nsim       = 1e3                                      # Anzahl Realisierungen n-dimensionaler ZV
y          = matrix(rep(NaN,n*nsim), nrow = n)       # y Realisierungsarray
beta_hat   = matrix(rep(NaN,p*nsim), nrow = p)       # \hat{\beta} Realisierungsarray
sigsqr_hat = rep(NaN, nsim)                          # \hat{\sigma}^2 Realisierungsarray
for(i in 1:nsim){
  y[,i]         = mvrnorm(1, X %*% beta, sigsqr*I_n) # eine Realisierung n-dimensionaler ZV
  beta_hat[,i]  = inv(t(X) %*% X) %*% t(X) %*% y[,i] # \hat{\beta}    = (X^T)X^{-1}X^Ty
  eps_hat       = y[,i] - X %*% beta_hat[,i]         # \hat{\eps}     = y-X\hat{\beta}
  sigsqr_hat[i] = (t(eps_hat) %*% eps_hat)/(n-p)     # \hat{\sigma}^2 =\hat{\eps}^T\hat{\eps}/(n-p)
}
U = ((n-p)/sigsqr)*sigsqr_hat                        # \chi^2 verteilte Zufallsvariable
```

# Schätzerverteilungen
Beispiel (2) Einfache lineare Regression

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


# Datenrealisierungen
xlimits     = c(0,11)
ylimits     = c(-5,15)
plot(
x,
X %*% beta,
type        = "b",
lty         = 2,
pch         = 1,
col         = "white",
xlab        = "x",
ylab        = "y",
xlim        = xlimits,
ylim        = ylimits)
for(i in 1:nsim){
  lines(
  type = "l",
  x,
  y[,i],
  col = "grey")
}

# density
xlimits     = c(0,30)
ylimits     = c(0,.12)
u_min  = xlimits[1]
u_max  = xlimits[2]
u_res  = 1e3
u      = seq(u_min, u_max, len = u_res)
p_u    = dchisq(u,n-p)

# histogram
hist(
U,
col   = "gray90",
prob  = TRUE,
xlab  = TeX("$((n-p)/\\sigma^2)\\hat{\\sigma}^2$"),
ylab  = "",
xlim  = xlimits,
ylim  = ylimits,
main   = "")

# density
lines(
u,
p_u,
lwd   = 2,
col   = "darkorange")

# Speichern
dev.copy2pdf(
file        = file.path(fdir, "alm_7_sigsqr_hat_1.pdf"),
width       = 7,
height      = 3.5)
```
\vspace{4mm}
\center
$\quad\quad  y \sim (X\beta,\sigma^2I_n)$ \hspace{3cm} $\frac{n-p}{\sigma^2}\hat{\sigma}^2 \sim \chi^2(n-p)$

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


#  {.plain}
\vfill
\large
\setstretch{2.5}
Schätzerverteilungen

**T-Statistiken**

F-Statistiken

Selbstkontrollfragen
\vfill

# T-Statistiken
Überblick
\setstretch{1.6}

\footnotesize
* \justifying T-Statistiken quantifizieren die geschätzten Effekte in $\hat{\beta}$
in bezug auf die Residualvariabilität $\hat{\sigma}^2$. T-Statistiken sind also als
Signal-zu-Rauschen Verhältnisse (signal-to-noise ratios) zu verstehen.
\vspace{2mm}

* \justifying T-Statistiken erlauben die Frequentistische Evaluation von linearen Kombinationen
der $\hat{\beta}$ Komponenten. Eine einfache Linearkombination von $\hat{\beta}$
ist das Skalarprodukt meinem Einheitsvektor. Solch eine Linearkombination "extrahiert"
eine Komponente von $\hat{\beta}$ und erlaubt damit ihre Frequentistische Evaluation.
\vspace{2mm}

* \justifying T-Teststatistiken erlauben weiterhin die Evaluation von Linearkombinationen
von $\hat{\beta}$ im Sinne Frequentistischer Hypothesentests. Wir betrachten hier 
zunächst die nur die funktionale Form von T-Teststatistiken und ihre Frequentistische
Verteilung. Der Einsatz von T-Teststatistiken im Rahmen von Einstich- und Zweistichproben
T-Tests ist Inhalt von (9) T-Tests.


# T-Statistiken
\footnotesize
\begin{definition}[$t$-Zufallsvariable]
\justifying
$T$ sei eine Zufallsvariable mit Ergebnisraum $\mathbb{R}$ und WDF
\begin{equation}
p : \mathbb{R} \to \mathbb{R}_{>0}, t \mapsto p(t)
:= \frac{\Gamma\left(\frac{n+1}{2}\right)}{\sqrt{n\pi}\Gamma\left(\frac{n}{2}\right)}
\left(1 + \frac{t^2}{n} \right)^{-\frac{n+1}{2}},
\end{equation}
wobei $\Gamma$ die Gammafunktion bezeichne. Dann sagen wir, dass $T$ einer
$t$-Verteilung mit $n$ Freiheitsgraden unterliegt und nennen $T$ eine $t$-Zufallsvariable
mit $n$ Freiheitsgraden. Wir kürzen dies mit $T \sim t(n)$ ab. Die WDF einer
$t$-Zufallsvariable bezeichnen wir mit
\begin{equation}
T(t;n) := \frac{\Gamma\left(\frac{n+1}{2}\right)}{\sqrt{n\pi}\Gamma\left(\frac{n}{2}\right)}
\left(1 + \frac{t^2}{n} \right)^{-\frac{n+1}{2}}.
\end{equation}
\end{definition}

Bemerkungen

* Die Verteilung ist um 0 symmetrisch
* Steigendes $n$ verschiebt Wahrscheinlichkeitsmasse aus den Ausläufen zum Zentrum.
* Ab $n = 30$ gilt $T(t;n) \approx N(0,1)$.


# T-Statistiken

```{r, echo = F, eval = F}
dev.new()                                                                        # new figure
fig     = par(                                                                   # figure parameters
            family     = "sans",                                                 # font family
            pty        = "m",                                                    # square plots
            bty        = "l",                                                    # plot box, o, l, 7, c, or ]
            lwd        = 1,                                                      # line width
            las        = 1,                                                      # 0: axis parallel, 1: horizontal, 2: axis perpendicular, 3: vertical
            mgp        = c(2,1,0),                                               # margin line in mex unit
            xaxs       = "i",                                                    # "internal" (tight) x-axis style
            yaxs       = "i",                                                    # "internal" (tight) y-axis style
            font.main  = 1,                                                      # title font type
            cex        = 1.1,
            cex.main   = 1.5                                                     # title  magnification factor
)

# t space
t_min   = -5                                                                     # minimum t-value
t_max   = 5                                                                      # maximum t-value
t_res   = 1e3                                                                    # t-space resolution
t       = seq(t_min,t_max, len = t_res)                                          # t-space

# parameters of interest
n       = c(2,3,5,10,30)                                                        # degrees of freedom


# plotting
matplot(t,
        matrix(c(dt(t,n[1]),
                 dt(t,n[2]),
                 dt(t,n[3]),
                 dt(t,n[4]),
                 dt(t,n[5])),
               ncol = 5),
        type         = "l",                                                      # line plot
        lty          = 1,                                                        # solid line
        lwd          = 2,                                                        # line width
        col          = c("gray10", "gray20", "gray50", "gray70", "gray90"),       # line colors
        ylim         = c(0,.4),                                                  # y-axis limits
        xlim         = c(t_min,t_max),                                           # x-axis limits
        ylab         = " ",                                                      # y-axis label
        xlab         = "t",                                                      # x-axis label
        main         = TeX("$\\T(t;n)$"))                                        # main title

legend( 2,                                                                       # x-ordinate
        .4,                                                                      # y-ordinate
        c("n = 2", "n = 3", "n = 5", "n = 10", "n = 30"),                        # legend text
        lty         = 1,                                                         # line type
        lwd         = 2,                                                         # line width
        col         =   c("gray10", "gray20", "gray50", "gray70", "gray90"),     # line colors
        bty         = "n",                                                       # no legend box
        cex         = 1.1,                                                       # character expansion (fontsize)
        y.intersp   = 1.6)                                                       # y-direction character spacing


dev.copy2pdf(                                                                    # export to PDF
    file   = file.path(getwd(), "7_Abbildungen", "alm_7_t_wdf.pdf"),            # filename
    width  = 6,                                                                  # PDF width
    height = 5                                                                   # PDF height
)
```

\vfill
\vspace{5mm}
```{r, echo = FALSE, out.width="60%"}
knitr::include_graphics("7_Abbildungen/alm_7_t_wdf.pdf")
```
\vspace{-3mm}
\footnotesize


# $T$-Transformation
\footnotesize
\begin{theorem}[$T$-Transformation]
\justifying
\normalfont
$Z \sim N(0,1)$ sei eine  $Z$-Zufallsvariable, $U \sim \chi^2(n)$ sei eine
$\chi^2$-Zufallsvariable mit  $n$ Freiheitsgraden, und $Z$ und $U$ seien
unabhängige Zufallsvariablen. Dann ist die Zufallsvariable
\begin{equation}
T := \frac{Z}{\sqrt{U/n}}
\end{equation}
eine $t$-verteilte Zufallsvariable mit $n$ Freiheitsgraden, es gilt also $T \sim t(n)$.
\end{theorem}


Bemerkungen

* Das Theorem geht auf @student_1908 zurück.
* Das Theorem ist das zentrale Resultat der Frequentistischen Statistik.
* @zabell_2008 gibt hierzu einen historischen Überblick.
* Die T-Statistik und T-Tests sind wichtige Anwendungsfälle.


# T-Statistiken

\footnotesize
\begin{definition}[Nichtzentrale $t$-Zufallsvariable]
\justifying
$T$ sei eine Zufallsvariable mit Ergebnisraum $\mathbb{R}$ und WDF
\begin{multline}
p : \mathbb{R} \to \mathbb{R}_{>0}, t \mapsto p(t) :=
\frac{1}{2^{\frac{n-1}{2}}\Gamma\left(\frac{n}{2} \right)(n \pi)^{\frac{1}{2}}} \\
\times \int_{0}^\infty \tau^{\frac{n-1}{2}} \exp\left(-\frac{\tau}{2}\right)
\exp\left(-\frac{1}{2}\left(t \left(\frac{\tau}{n}\right)^{\frac{1}{2}} - \delta \right)^2 \right)\,d\tau.
\end{multline}
Dann sagen wir, dass $T$ einer
nichtzentralen $t$-Verteilung mit Nichtzentralitätsparameter $\delta$ und
Freiheitsgradparameter $n$ unterliegt und nennen $T$ eine \textit{nichtzentrale
$t$-Zufallsvariable mit Nichtzentralitätsparameter $\delta$ und
Freiheitsgradparameter $n$}. Wir kürzen dies mit $t(\delta, n)$ ab. Die WDF einer
nichtzentralen $t$-Zufallsvariable bezeichnen wir mit
$t(T;\delta,n)$. Die KVF und inverse KVF einer nichtzentralen $t$-Zufallsvariable
bezeichnen wir mit $\psi(\cdot; \delta, n)$ und $\psi^{-1}(\cdot; \delta, n)$, respektive.
\end{definition}

Bemerkungen

* Eine nichtzentrale $t$-Zufallsvariable mit $\delta = 0$ ist eine $t$-Zufallsvariable.
* Es gilt also $t(T;0,n) = t(T;n)$.
* Die funktionale Form der WDF findet sich zum Beispiel in @lehmann_1986, Seite 254, Gl. (80).


# T-Statistiken

Wahrscheinlichkeitsdichtefunktionen nichtzentraler $t$-Verteilungen
\vspace{4mm}

```{r, echo = F, eval = F}
options(warn=-1)                                                                 # warning off
t_min     = -5                                                                   # Minimum T-Wert
t_max     = 30                                                                   # Maximum T-Wert
t_res     = 1e3                                                                  # T-Wert Auflösung
t         = seq(t_min, t_max, len = t_res)                                       # T-Raum
delta     = c(0,5,15)                                                            # Nichtzentralitätsparameter
n         = c(5, 30)                                                             # Freiheitsgrade
p         = cbind(
            matrix(dt(t, n[1], delta[1]),nrow=length(t)),
            matrix(dt(t, n[2], delta[1]),nrow=length(t)),
            matrix(dt(t, n[1], delta[2]),nrow=length(t)),
            matrix(dt(t, n[2], delta[2]),nrow=length(t)),
            matrix(dt(t, n[1], delta[3]),nrow=length(t)),
            matrix(dt(t, n[2], delta[3]),nrow=length(t)))

# Visualisierung
dev.new()
graphics.off()
par(
family      = "sans",
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)
matplot(
t,
p,
type        = "l",
lty         = c(1,2,1,2,1,2),
col         = c("gray10", "gray10", "gray50", "gray50", "gray70", "gray70"),
lwd         = 2,
xlab        = "T",
ylab        = "",
ylim        = c(0,.4),
main        = TeX("$t(T;\\delta,n)$"))
legend(
18,
.4,
c(TeX("$\\delta = 0 , n = 5$"),
  TeX("$\\delta = 0 , n = 30$"),
  TeX("$\\delta = 5 , n = 5$"),
  TeX("$\\delta = 5 , n = 30$"),
  TeX("$\\delta = 15, n = 5$"),
  TeX("$\\delta = 15, n = 30$")),
lty         = c(1,2,1,2,1,2),
col         = c("gray10", "gray10", "gray50", "gray50", "gray70", "gray70"),
lwd         = 2,
bty         = "n",
seg.len     = 1.75)
fdir        =  file.path(getwd(), "7_Abbildungen")
dev.copy2pdf(
file        = file.path(fdir, "alm_7_nichtzentrale_t_verteilung.pdf"),
width       = 7,
height      = 4.5)
```

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

# T-Statistiken
\footnotesize
\begin{theorem}[Nichtzentrale T-Transformation]
\normalfont
\justifying
$X \sim N(\mu,1)$ sei eine normalverteilte Zufallsvariable, $U \sim \chi^2(n)$
sei eine $\chi^2$ Zufallsvariable mit Freiheitsgradparameter $n$, und $X$ und
$U$ seien unabhängige Zufallsvariablen. Dann ist die Zufallsvariable
\begin{equation}
T := \frac{X}{\sqrt{U/n}}
\end{equation}
eine nichtzentrale $t$-Zufallsvariable mit Nichtzentralitätsparameter $\mu$ und
Freiheitsgradparameter $n$, es gilt also $T \sim t(\mu,n)$
\end{theorem}

Bemerkung

* Wir verzichten auf einen Beweis.

# T-Statistiken
\footnotesize
\begin{definition}[T-Statistik]
Es sei
\begin{equation}
y = X\beta + \varepsilon \mbox{ mit } \varepsilon \sim N(0_n,\sigma^2I_n)
\end{equation}
das ALM in generativer Form. Weiterhin seien
\begin{equation}
\hat{\beta} := (X^TX)^{-1}X^Ty \mbox{ und } \hat{\sigma}^2 := \frac{(y - X\hat{\beta})^T(y - X\hat{\beta})}{n-p}
\end{equation}
die Betaparameter- und Varianzparameterschätzer, respektive. Dann ist für einen
\textit{Kontrastgewichtsvektor} $c \in \mathbb{R}^p$ die \textit{T-Statistik} definiert als
\begin{equation}
T := \frac{c^T\hat{\beta}}{\sqrt{\hat{\sigma}^2c^T(X^TX)^{-1}c}}.
\end{equation}
\end{definition}

Bemerkungen

* Die T-Statistik hängt via $\hat{\beta}$ und $\hat{\sigma}^2$ von den Daten $y$ ab.
* Der Kontrastgewichtsvektor projiziert $\hat{\beta}$ auf einen Skalar $c^T\hat{\beta}$.
* Die folgende Intuition zur T-Statistik ist im Allgemeinen hilfreich
\begin{equation}
T = \frac{\mbox{Geschätzte Effektstärke}}{\mbox{Geschätzte stichprobenumfangskalierte Datenvariabilität}}
\end{equation}
* Intuitiv repräsentiert eine T-Statistik also ein Signal-zu-Rauschen Verhältnis.

# T-Statistiken
Beispiel (1) Unabhängige und identisch normalverteilte Zufallsvariablen
\vspace{2mm}
\footnotesize

Es sei
\begin{equation}
y \sim N(X\beta,\sigma^2 I_n)
\mbox{ mit }
X := 1_n \in \mathbb{R}^{n\times 1},
\beta := \mu \in \mathbb{R}
\mbox{ und } \sigma^2 > 0.
\end{equation}
das ALM Szenario unabhängiger und identisch normalverteilter Zufallsvariablen und
es sei $c := 1$. Dann gilt für die T-Statistik
\begin{equation}
T
= \frac{c^T\hat{\beta}}{\sqrt{\hat{\sigma}^2c^T(X^TX)^{-1}c}}
= \frac{1^T\bar{y}}{\sqrt{s^2_y 1^T (1_n^T 1_n)^{-1}1}}
= \sqrt{n}\frac{\bar{y}}{s_y}
\end{equation}
was der Einstichproben-T-Teststatistik für den Fall $\mu_0 = 0$ entspricht (vgl.
Einheit (13) Einstichproben-T-Tests in Wahrscheinlichkeitstheorie und Frequentistische
Inferenz und Einheit (9) T-Tests in Allgemeines Lineares Modell). Die hier betrachtete
T-Statistik nimmt hohe Werte für hohe Werte von $\bar{y}$ (Effekt), kleine Werte
von $s_y^2$ (Datenvariabilität) und hohe Werte von $n$ (Stichprobenumfang) an.

Eine beliebte Definition in diesem Zusammenhang ist \textit{Cohen's $d$} als \textit{Effektstärkenmaß}.
Es gilt
\begin{equation}
d := \frac{\bar{y}}{s_y},
\end{equation}
so dass
\begin{equation}
T = \sqrt{n}d \mbox{ bzw. } d = \frac{1}{\sqrt{n}} T.
\end{equation}
Cohen's $d$ ist also ein stichprobenumfangunabhängiges Signal-zu-Rauschen Verhältnis.

# T-Statistiken
\vspace{1mm}
\footnotesize
\begin{theorem}[T-Teststatistik]
\normalfont
\justifying
Es sei
\begin{equation}
y = X\beta + \varepsilon \mbox{ mit } \varepsilon \sim N(0_n,\sigma^2I_n)
\end{equation}
das ALM in generativer Form. Weiterhin seien
\begin{equation}
\hat{\beta} := (X^TX)^{-1}X^Ty \mbox{ und } \hat{\sigma}^2 := \frac{(y - X\hat{\beta})^T(y - X\hat{\beta})}{n-p}
\end{equation}
die Betaparameter- und Varianzparameterschätzer, respektive. Schließlich sei für
einen \textit{Kontrastgewichtsvektor} $c \in \mathbb{R}^p$ und
einen \textit{Nullhypothesenbetaparameter} $\beta_0 \in \mathbb{R}^p$
die \textit{T-Teststatistik} definiert als
\begin{equation}
T := \frac{c^T\hat{\beta} - c^T\beta_0}{\sqrt{\hat{\sigma}^2 c^T(X^TX)^{-1}c}}.
\end{equation}
Dann gilt
\begin{equation}
T \sim t(\delta, n-p) \mbox{ mit } \delta := \frac{c^T\beta - c^T\beta_0}{\sqrt{\sigma^2 c^T(X^TX)^{-1}c}}
\end{equation}
\end{theorem}
\vspace{-2mm}

Bemerkungen
\vspace{-2mm}

* Wir verzichten auf einen Beweis.
* $T$ ist eine Funktion der Parameterschätzer, $\delta$ ist eine Funktion der wahren, aber unbekannten, Parameter
* Für $c^T\beta = c^T\beta_0$, also bei Zutreffen der Nullhypothese, gilt $\delta = 0$ und damit $T \sim t(n-p)$.
* Für $c^T\beta \neq c^T\beta_0$ kann die Verteilung von $T$ zur Herleitung von Powerfunktionen benutzt werden..

# T-Statistiken
Simulation (1) Unabhängig und identische normalverteilte Zufallsvariablen

\small
Wahre, aber unbekannte, Hypothesenszenarien $c^T\beta = c^T\beta_0$ und $c^T\beta \neq c^T\beta_0$

\vspace{2mm}
\tiny
\setstretch{1}
```{r}
# Libraries
library(MASS)                                                  # multivariate Normalverteilung

# Modellformulierung
n          = 12                                                # Anzahl von Datenpunkten
p          = 1                                                 # Anzahl von Betaparametern
X          = matrix(c(rep(1,n)), nrow = n)                     # Designmatrix
I_n        = diag(n)                                           # Einheitsmatrix
beta       = c(0,1)                                            # wahre , aber unbekannte , Betaparameter
nscn       = length(beta)                                      # Anzahl wahrer, aber unbekannter, Hypothesenszenarien
sigsqr     = 1                                                 # wahrer, aber unbekannter, Varianzparameter
c          = 1                                                 # Kontrastvektor von Interessse
beta_0     = 0                                                 # Nullhypothesenbetaparameter

# Frequentistische Simulation
nsim       = 1e4                                               # Anzahl Simulationen
delta      = rep(NaN, nscn)                                    # Anzahl Nichtzentralitätsparameter
Tee        = matrix(rep(NaN, nscn*nsim), ncol = nscn)          # T-Teststatistik Realisierungsarray
for(s in 1:nscn){                                              # Hypothesenszenarien
  delta[s]    = ((t(c) %*% beta[s] - t(c) %*% beta_0)/         # Nichtzentralitätsparameter
                sqrt(sigsqr*t(c)%*%solve(t(X)%*%X)%*%c))
  for(i in 1:nsim){                                            # Simulationsiterationen
    y          = mvrnorm(1, X %*% beta[s], sigsqr*I_n)         # y
    beta_hat   = solve(t(X) %*% X) %*% t(X) %*% y              # \hat{\beta}
    eps_hat    = y - X %*% beta_hat                            # \hat{\eps}
    sigsqr_hat = (t(eps_hat) %*% eps_hat)/(n-p)                # \hat{\sigma}^2
    Tee[i,s]   = ((t(c) %*% beta_hat - t(c) %*% beta_0)/       # T
                  sqrt(sigsqr_hat*t(c)%*%solve(t(X)%*%X)%*%c))
  }
}
```

# T-Statistiken
Simulation (1) Unabhängig und identische normalverteilte Zufallsvariablen

\footnotesize
Wahre, aber unbekannte, Hypothesenszenarien $c^T\beta = c^T\beta_0$ und $c^T\beta \neq c^T\beta_0$

\vspace{8mm}

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

# Visualisierung
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2.5,1,0),
xaxs        = "i",
yaxs        = "i",
xpd         = TRUE,
font.main   = 1,
cex         = 1,
cex.main    = 1.2)

# T-Teststatistik Ergebnisraum
xlims  = c(-5,12)
t_min  = xlims[1]
t_max  = xlims[2]
t_res  = 1e3
t      = seq(t_min, t_max, len = t_res)
lab    = c(TeX("$\\c^T beta = \\c^T beta_0$"), TeX("$\\c^T beta \\neq \\c^T beta_0$") )

# T-Teststatistiken
for(s in 1:nscn){
  p_t    = dt(t,n-p, delta[s])
  hist(
  Tee[,s],
  breaks = 50,
  col    = "gray90",
  prob   = TRUE,
  xlab   = TeX("$T$"),
  ylab   = "",
  xlim   = xlims,
  ylim   = c(0,.4),
  main   = lab[s])
  lines(
  t,
  p_t,
  type  = "l",
  lwd   = 2,
  col   = "darkorange")
}

# Speichern
dev.copy2pdf(
file        = file.path(fdir, "alm_7_T_Teststatistik_1.pdf"),
width       = 8,
height      = 4)
```


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

# T-Statistiken
Simulation (1) Einfache lineare Regression

\small
Wahre, aber unbekannte, Hypothesenszenarien $c^T\beta = c^T\beta_0$ und $c^T\beta \neq c^T\beta_0$

\vspace{2mm}
\tiny
\setstretch{1}
```{r}
# Libraries
library(MASS)                                                  # multivariate Normalverteilung
library(matlib)                                                # Matrizenrechnung

# Modellformulierung
n          = 10                                                # Anzahl von Datenpunkten
p          = 2                                                 # Anzahl von Betaparametern
x          = 1:n                                               # Prädiktorwerte
X          = matrix(c(rep(1,n),x), ncol = p)                   # Designmatrix
I_n        = diag(n)                                           # Einheitsmatrix
beta       = matrix(c(1,0,1,1), nrow = 2)                      # wahre , aber unbekannte , Betaparameter
nscn       = ncol(beta)                                        # Anzahl wahrer, aber unbekannter, Hypothesenszenarien
sigsqr     = 1                                                 # wahrer, aber unbekannter, Varianzparameter
c          = matrix(c(0,1), nrow = 2)                          # Kontrastvektor von Interessse
beta_0     = matrix(c(0,0), nrow = 2)                          # Nullhypothesenbetaparameter

# Frequentistische Simulation
nsim       = 1e4                                               # Anzahl Simulationen
delta      = rep(NaN, nscn)                                    # Anzahl Nichtzentralitätsparameter
Tee        = matrix(rep(NaN, nscn*nsim), ncol = nscn)          # T-Teststatistik Realisierungsarray
for(s in 1:nscn){                                              # Hypothesenszenarien
  delta[s]    = ((t(c) %*% beta[,s] - t(c) %*% beta_0)/        # Nichtzentralitätsparameter
                sqrt(sigsqr*t(c)%*%inv(t(X)%*%X)%*%c))
  for(i in 1:nsim){                                            # Simulationsiterationen
    y          = mvrnorm(1, X %*% beta[,s], sigsqr*I_n)        # y
    beta_hat   = inv(t(X) %*% X) %*% t(X) %*% y                # \hat{\beta}
    eps_hat    = y - X %*% beta_hat                            # \hat{\eps}
    sigsqr_hat = (t(eps_hat) %*% eps_hat)/(n-p)                # \hat{\sigma}^2
    Tee[i,s]   = ((t(c) %*% beta_hat - t(c) %*% beta_0)/       # T
                  sqrt(sigsqr_hat*t(c)%*%inv(t(X)%*%X)%*%c))
  }
}
```

# T-Statistiken
Simulation (1) Einfache lineare Regression

\small
Wahre, aber unbekannte, Hypothesenszenarien $c^T\beta = c^T\beta_0$ und $c^T\beta \neq c^T\beta_0$

\vspace{8mm}

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

# Visualisierung
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2.5,1,0),
xaxs        = "i",
yaxs        = "i",
xpd         = TRUE,
font.main   = 1,
cex         = 1,
cex.main    = 1.2)

# T-Teststatistik Ergebnisraum
xlims  = c(-5,20)
t_min  = xlims[1]
t_max  = xlims[2]
t_res  = 1e3
t      = seq(t_min, t_max, len = t_res)
lab    = c(TeX("$\\c^T beta = \\c^T beta_0$"), TeX("$\\c^T beta \\neq \\c^T beta_0$") )

# T-Teststatistiken
for(s in 1:nscn){
  p_t    = dt(t,n-p, delta[s])
  hist(
  Tee[,s],
  breaks = 50,
  col    = "gray90",
  prob   = TRUE,
  xlab   = TeX("$T$"),
  ylab   = "",
  xlim   = xlims,
  ylim   = c(0,.4),
  main   = lab[s])
  lines(
  t,
  p_t,
  type  = "l",
  lwd   = 2,
  col   = "darkorange")
}

# Speichern
dev.copy2pdf(
file        = file.path(fdir, "alm_7_T_Teststatistik_2.pdf"),
width       = 8,
height      = 4)
```

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


#  {.plain}
\vfill
\large
\setstretch{2.5}
Schätzerverteilungen

T-Statistiken

**F-Statistiken**

Selbstkontrollfragen
\vfill

# F-Statistiken
Überblick
\setstretch{1.6}

\footnotesize
* \justifying F-Statistiken können als Modellvergleichskriterien zwischen ALMs verstanden werden.
Im Gegensatz zu T-Statistiken kann das Ziel der Berechnung von F-Statistiken sein,
nicht nur Linearkombinationen von Betaparameterschätzwerten probabilistisch zu evaluieren,
sondern die Modellanpassung an einen Datensatz insgesamt zu evaluieren.

* Wir motivieren die funktionale Form der F-Statistik als Residualquadratsummendifferenz
hier vor dem Hintergrund eine Likelihood-Quotienten zum Modellvergleich. Die Modellvergleichskapazität 
von F-Statistiken ist allerdings etwas beschränkt, da sich die 
F-Statistik auf ALMs und insbesondere geschachtelte ALMs bezieht, in denen ein Modell
Bestandteil eines anderen Modells ist.

* F-Statistiken bilden die Grundlage für Hypothesentests im Rahmen der Varianzanalyse 
(Einheiten (10) Einfaktorielle Varianzanalyse, (11) Zweifaktorielle Varianzanalyse, (13) Kovarianzanalyse).
Der Einsatz von F-Statistiken ist aber \textit{perse} nicht auf Varianzanalysen beschränkt.

# F-Statistiken
\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 $\varphi(x;\nu_1,\nu_2)$, und die inverse
kumulative Verteilungsfunktion einer $f$-Zufallsvariable bezeichnen wir mit $\varphi^{-1}(x;\nu_1,\nu_2)$.
\end{definition}

Bemerkungen

* $f$-Zufallsvariablen sind nach Ronald A. Fisher benannt.
* George W. Snedecor hat die KVF der $F$-Verteilung wohl 1934 basierend auf Arbeiten von Fisher tabuliert.

# F-Statistiken
\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(), "7_Abbildungen", "alm_7_f_wdf.pdf"),
width  = 6,
height = 4)
```
\small
Wahrscheinlichkeitsdichtefunktionen von $f$-Verteilungen
\vspace{1mm}
```{r, echo = FALSE, out.width="80%"}
knitr::include_graphics("7_Abbildungen/alm_7_f_wdf.pdf")
```

# F-Statistiken
\small
\begin{theorem}[$F$-Transformation]
\justifying
\normalfont
$V \sim \chi^2(\nu_1)$ und $W \sim \chi^2(\nu_2)$ seien zwei unabhängige
$\chi^2$-Zufallfsvariablen mit $\nu_1$ und $\nu_2$ Freiheitsgraden, respektive.
Dann ist die Zufallsvariable
\begin{equation}
F := \frac{V/\nu_1}{W/\nu_2}
\end{equation}
eine $f$-verteilte Zufallsvariable mit $\nu_1,\nu_2$ Freiheitsgraden, es gilt also $F \sim f(\nu_1,\nu_2)$.
\end{theorem}

\footnotesize
Bemerkungen

* Wir verzichten auf einen Beweis.
* Das Theorem kann bewiesen werden, in dem man zunächst ein Transformationstheorem
für Quotienten von Zufallsvariablen mithilfe des multivariaten Transformationstheorems
und Marginalisierung herleitet und dieses Theorem dann auf die WDF von $\chi^2$-verteilten
ZVen anwendet. Dabei ist die Regel zur Integration durch Substitution von zentraler Bedeutung.

# F-Statistiken
\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 $\varphi(x;\delta,\nu_1,\nu_2)$, und die inverse kumulative Verteilungsfunktion
einer nichtzentralen $f$-Zufallsvariable
bezeichnen wir mit $\varphi^{-1}(x;\delta,\nu_1,\nu_2)$.
\end{definition}
Bemerkungen

* Es gilt $f(0,\nu_1,\nu_2) = f(\nu_1,\nu_2)$.

# F-Statistiken
\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(), "7_Abbildungen", "alm_7_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("7_Abbildungen/alm_7_f_nichtzentral_wdf.pdf")
```

# F-Statistiken
\small
\begin{theorem}[Nichtzentrale $F$-Transformation]
\justifying
\normalfont
$V \sim \chi^2(\delta,\nu_1)$ und $W \sim \chi^2(\nu_2)$ seien eine nichtzentrale $\chi^2$
Zufallsvariable und eine  $\chi^2$-Zufallfsvariablen mit $\nu_1$ und $\nu_2$ Freiheitsgraden,
respektive, und $V$ und $W$ seien unabhängig. Dann ist die Zufallsvariable
\begin{equation}
F := \frac{V/\nu_1}{W/\nu_2}
\end{equation}
eine nichtzentral $f$-verteilte Zufallsvariable mit Nichtzentralitätsparameter
$\delta$  und $\nu_1,\nu_2$ Freiheitsgraden, es gilt also $F \sim f(\delta, \nu_1,\nu_2)$.
\end{theorem}

\footnotesize
Bemerkungen

* Wir verzichten auf einen Beweis.


# F-Statistiken
\small
Motivation der F-Statistik $\vert$ Modellpartitionierung und Likelihood Ratios

\footnotesize
\begin{definition}[Vollständiges und reduziertes Modell]
Für $p > 1$ mit $p = p_1 + p_2$ seien
\begin{equation}
X := \begin{pmatrix} X_1 & X_2 \end{pmatrix}  \in \mathbb{R}^{n \times p}
\mbox{ mit }
X_1 \in \mathbb{R}^{n \times p_1}
\mbox{ und }
X_2 \in \mathbb{R}^{n \times p_2},
\end{equation}
sowie
\begin{equation}
\beta := \begin{pmatrix} \beta_1 \\ \beta_2 \end{pmatrix} \in \mathbb{R}^p
\mbox{ mit }
\beta_1 \in \mathbb{R}^{p_1}
\mbox{ und }
\beta_2 \in \mathbb{R}^{p_2}
\end{equation}
Partitionierungen einer $n \times p$ Designmatrix und eines $p$-dimensionalen
Betaparametervektors. Dann nennen wir
\begin{equation}
y = X\beta + \varepsilon \mbox{ mit } \varepsilon \sim N(0_n,\sigma^2I_n)
\end{equation}
das \textit{vollständige Modell} und
\begin{equation}
y = X_1\beta_1 + \varepsilon_1 \mbox{ mit } \varepsilon_1 \sim N(0_n,\sigma^2I_n)
\end{equation}
das \textit{reduzierte Modell} und sprechen von einer \textit{Partionierung eines
(vollständigen) Modells}.
\end{definition}
Bemerkungen

* Man sagt auch, dass das reduzierte Modell im vollständigen Modell \textit{geschachtelt (nested)} ist.

# F-Statistiken
\small
Motivation der F-Statistik $\vert$ Modellpartitionierung und Likelihood Ratios

\footnotesize
\begin{definition}[Likelihood-Quotient von vollständigem und reduzierten Modell]
Für $p = p_1 + p_2, p > 1$ sei eine Partitionierung eines vollständigen ALMs
gegeben und $\sigma^2 > 0$ als bekannt voraussgesetzt. Weiterhin seien
\begin{equation}
L : \mathbb{R}^{p} \to \mathbb{R}_{>0}, \beta \mapsto L(\beta) := p_{\beta}(\upsilon) = N(\upsilon;X\beta,\sigma^2 I_n)
\end{equation}
und
\begin{equation}
L_1 : \mathbb{R}^{p_1} \to \mathbb{R}_{>0}, \beta_1 \mapsto L(\beta_1) :=  p_{\beta_1}(\upsilon) = N(\upsilon;X_1\beta_1,\sigma^2 I_n)
\end{equation}
die Likelihood-Funktionen von vollständigem und reduzierten Modell, respektive.
Dann ist für einen Datenvektor $\upsilon \in \mathbb{R}^n$ der \textit{Likelihood-Quotient
von vollständigem und reduziertem Modell} gegeben als
\begin{equation}
\qoppa
:= \frac{\max_{\beta \in \mathbb{R}^p} p_{\beta}(\upsilon)}{\max_{\beta_1 \in \mathbb{R}^{p_1}} p_{\beta_1}(\upsilon)}
 = \frac{\max_{\beta \in \mathbb{R}^p} L(\beta)}{\max_{\beta_1 \in \mathbb{R}^{p_1}} L_1(\beta_1)}
 = \frac{L(\hat{\beta})}{L_1(\hat{\beta}_1)}
\end{equation}
\end{definition}

# F-Statistiken
\small
Motivation der F-Statistik $\vert$ Modellpartitionierung und Likelihood Ratios

\footnotesize
Bemerkungen

* \justifying Der Likelihood-Quotient setzt die Wahrscheinlichkeitsdichten eines beobachteten
Datenssatzes $\upsilon \in \mathbb{R}^n$ unter dem vollständigem und dem reduzierten
Modell nach Optimierung der jeweiligen Modellbetaparameter ins Verhältnis.

* \justifying Ein hoher Wert des Likelihood-Quotienten entspricht einer höhereren Wahrscheinlichkeit(sdichte)
des beobachteten Datensatzes $\upsilon \in \mathbb{R}^n$ unter dem optimierten vollständigem Modell
und vice versa.

* \justifying Die Wahrscheinlichkeiten beobachtbarerer Daten nach Modellschätzung unter verschiedenen Modellen
zu betrachten ist ein allgemeines Vorgehen zum Modellvergleich. Letztlich erlaubt dieses Vorgehen,
verschiedene wissenschaftliche Theorien (= Modelle) über die Genese beobachtbarer Daten quantitativ
zu vergleichen und die damit verbundene Unsicherheit zu quantifizieren.

* \justifying Modellvergleiche sind ein zentrales Thema in der Bayesianischen Statistik, die die Logik von
Likelihood-Quotienten zum Beispiel unter den Begriffen der Bayes Factors oder der des Bayesian
Information Criterions auf arbiträre probabilistische Modelle, d.h. nicht nur geschachtelte
ALMs generalisiert.

* \justifying Allerdings sind, wie hier gesehen, Modellvergleiche auch im Rahmen der Frequentistischen Statistik
möglich und sinnvoll, Modellvergleiche sind also kein Alleinstellungsmerkmal der Bayesianischen gegenüber
der Frequentistischen Statistik.


# F-Statistiken
\small
Motivation der F-Statistik $\vert$ Modellpartitionierung und Likelihood Ratios
\footnotesize
\begin{theorem}[Likelihood-Quotient und Residualquadratsummendifferenz]
\normalfont
\justifying
Für $p = p_1 + p_2, p > 1$ sei eine Partitionierung eines vollständigen ALMs
gegeben und $\sigma^2 > 0$ als bekannt voraussgesetzt. Weiterhin seien
\begin{equation}
\hat{\varepsilon}   := y - X\hat{\beta} \mbox{ und }
\hat{\varepsilon}_1 := y - X_1\hat{\beta}_1
\end{equation}
die Residuenvektoren des vollständigen und es reduzierten Modells, respektive.
Dann gilt
\begin{equation}
\ln \qoppa =
\frac{\hat{\varepsilon}_1^T\hat{\varepsilon}_1 - \hat{\varepsilon}^T\hat{\varepsilon}}{2\sigma^2}
\end{equation}
\end{theorem}

Bemerkungen

* $\ln \qoppa$ und die F-Statistik sind nicht identisch, aber eng verwandt.
* $\hat{\varepsilon}_1^T\hat{\varepsilon}_1$ ist die Residualquadratsumme (RQS), die bei Schätzung des reduzierten Modells von $y$ resultiert.
* $\hat{\varepsilon}^T\hat{\varepsilon}$ ist die RQS, die bei Schätzung des vollständigen Modells anhand desselben $y$ resultiert.
* Für $\beta_2 := 0_{p_2}$ sind vollständiges und reduziertes Modell identisch, es muss also $\hat{\varepsilon}_1^T\hat{\varepsilon}_1 \ge \hat{\varepsilon}^T\hat{\varepsilon}$ gelten.
* $\hat{\varepsilon}_1^T\hat{\varepsilon}_1 - \hat{\varepsilon}^T\hat{\varepsilon}$ misst, inwieweit für $\beta_2 \neq 0_{p_2}$ die Spalten von $X_2$ die RQS verringern.
* Ist $\hat{\varepsilon}_1^T\hat{\varepsilon}_1$ groß, aber $\hat{\varepsilon}^T\hat{\varepsilon}$ klein, so ist $\ln \qoppa$ groß und der Beitrag der $X_2$-Spalten zur Datenerklärung groß.
* Sind $\hat{\varepsilon}_1^T\hat{\varepsilon}_1$ und $\hat{\varepsilon}^T\hat{\varepsilon}$ gleich groß, also $\ln \qoppa$ klein, ist und der Beitrag der $X_2$-Spalten zur Datenerklärung klein.
* "Groß" und "klein" von $\hat{\varepsilon}_1^T\hat{\varepsilon}_1$ und $\hat{\varepsilon}^T\hat{\varepsilon}$ werden in Relation zur doppelten Fehlervektorvarianz $\sigma^2$ gemessen.

# F-Statistiken
\footnotesize
\begin{definition}[F-Statistik]
Für $X \in \mathbb{R}^{n \times p}, \beta \in \mathbb{R}^p$ und $\sigma^2 > 0$
sei ein ALM der Form
\begin{equation}
y = X\beta + \varepsilon \mbox{ mit } \varepsilon \sim N(0_n,\sigma^2I_n)
\end{equation}
mit der Partitionierung
\begin{equation}
X      = \begin{pmatrix} X_1     & X_2      \end{pmatrix},
X_1      \in \mathbb{R}^{n\times p_1},
X_2      \in \mathbb{R}^{n\times p_2},
\mbox{ und }
\beta := \begin{pmatrix} \beta_1 \\ \beta_2 \end{pmatrix},
\beta_1 \in \mathbb{R}^{p_1},
\beta_2 \in \mathbb{R}^{p_2},
\end{equation}
mit $p = p_1 + p_2$ gegeben. Weiterhin seien mit
\begin{equation}
\hat{\beta}_1 := (X_1^TX_1)^{-1} X_1^Ty \mbox{ und } \hat{\beta} := (X^TX)^{-1}X^Ty
\end{equation}
die Residuenvektoren
\begin{equation}
\hat{\varepsilon}_1 := y - X_1\hat{\beta}_1 \mbox{ und } \hat{\varepsilon} := y - X\hat{\beta}
\end{equation}
definiert. Dann ist die F-Statistik definiert als
\begin{equation}
F := \frac{\frac{\hat{\varepsilon}_1^T\hat{\varepsilon}_1 - \hat{\varepsilon}^T\hat{\varepsilon}}{p_2}}{\frac{\hat{\varepsilon}^T\hat{\varepsilon}}{n-p}}
\end{equation}
\end{definition}

# F-Statistiken
\footnotesize
Bemerkungen

* \justifying Der Zähler der F-Statistik
\begin{equation}
\frac{\hat{\varepsilon}_1^T\hat{\varepsilon}_1 - \hat{\varepsilon}^T\hat{\varepsilon}}{p_2}
\end{equation}
misst, inwieweit die $p_2$ Regressoren in $X_2$ die Residualquadratsumme reduzieren und zwar
im Verhältnis zur Anzahl dieser Regressoren. Das heißt, dass bei gleicher Größe der
Residualquadratsummenreduktion (und gleichem Nenner) ein größerer $F$ Wert resultiert,
wenn diese durch weniger zusätzliche Regressoren resultiert, also $p_2$ klein ist (und vice versa).
Im Sinne der Anzahl der Spalten von $X$ und der entsprechenden Komponenten von $\beta$
favorisiert die $F$-Statistik also weniger "komplexe" Modelle.
* Für den Nenner der F-Statistik gilt
\begin{equation}
\frac{\hat{\varepsilon}^T\hat{\varepsilon}}{n-p} = \hat{\sigma}^2,
\end{equation}
wobei $\hat{\sigma}^2$ hier der aufgrund des vollständigen Modells geschätzte
Schätzer von $\sigma^2$ ist. Werden die Daten tatsächlich unter dem reduzierten Modell
generiert, so kann das vollständige Modell dies durch $\hat{\beta}_{2} \approx 0_{p_2}$
abbilden und erreicht eine ähnliche $\sigma^2$ Schätzung wie das reduzierte Modell.
Werden die Daten de-facto unter dem vollständigem Modell generiert, so ist
$\hat{\varepsilon}^T\hat{\varepsilon}/(n-p)$ ein besserer Schätzer von $\sigma^2$
als $\hat{\varepsilon}^T_1\hat{\varepsilon_1}/(n-p)$, da sich für diesen Datenvariabilität,
die nicht durch die $p_1$ Regressoren in $X_1$ erklärt wird, in der Schätzung von $\sigma^2$
widerspiegeln würde. Der Nenner der F-Statistik ist also in beiden Fällen der sinnvollere
Schätzer von $\sigma^2$.
* Zusammengenommen misst die F-Statistik also die Residualquadratsummenreduktion
durch die $p_2$ Regressoren in $X_2$ gegenüber den  $p_1$ Regressoren in $X_1$
pro Datenvariabilitäts ($\sigma^2$)- und Regressor ($p_2$)-Einheit.


# F-Statistiken
\vspace{1mm}
\setstretch{.9}
\small
Beispiel (1)  Einfache lineare Regression
\footnotesize
\begin{equation}
X = \begin{pmatrix} X_1 & X_2 \end{pmatrix}, X_1 := 1_{n}, X_2 := (x_1,...,x_n)^T,
\beta := \begin{pmatrix} \beta_1 \\ \beta_2 \end{pmatrix},
\beta_1 = \mbox{Interzept,  }
\beta_2 = \mbox{Steigung }
\end{equation}
\tiny
```{r}
# Libraries
library(MASS)                                             # Multivariate Normalverteilung

# Modellformulierung
nmod   = 2                                                # Anzahl Modelle
n      = 10                                               # Anzahl Datenpunkte
p      = 2                                                # Anzahl Betaparameter
p_1    = 1                                                # Anzahl Betaparameter reduziertes Modell
p_2    = 1                                                # Anzahl zusätzlicher Betaparameter vollständiges Modell
p      = p_1 + p_2                                        # Anzahl Betaparameter im vollständigem Modell
x      = 1:n                                              # Prädiktorwerte
X      = matrix(c(rep(1,n),x), nrow = n)                  # Designmatrix des vollständigen Modells
X_1    = X[,1]                                            # Designmatrix des reduzierten Modells
I_n    = diag(n)                                          # n x n Einheitsmatrix
beta   = matrix(c(1,0,1,.5), nrow = 2)                    # wahre , aber unbekannte , Betaparameter
nscn   = ncol(beta)                                       # Anzahl wahrer, aber unbekannter, Hypothesenszenarien
sigsqr = 1                                                # wahrer, aber unbekannter, Varianzparameter

# Modellsimulation und Evaluierung
Eff    = matrix(rep(NaN, nscn), nrow = nscn)              # F-Statistiik Realisierungsarray
for(s in 1:nscn){                                         # Szenarieniterationen
  y               = mvrnorm(1, X %*%beta[,s], sigsqr*I_n) # Datenrealisierung
  beta_hat_1      = solve(t(X_1)%*%X_1)%*%t(X_1)%*%y      # Betaparameterschätzer reduziertes Modell
  beta_hat        = solve(t(X)  %*%X  )%*%t(X)  %*%y      # Betaparameterschätzer vollständiges Modell
  eps_1_hat       = y-X_1%*%beta_hat_1                    # Residuenvektor reduziertes Modell
  eps_hat         = y-X%*%beta_hat                        # Residuenvektor vollständiges Modell
  eps_1_eps_1_hat = t(eps_1_hat) %*% eps_1_hat            # RQS reduziertes Modell
  eps_eps_hat     = t(eps_hat)   %*% eps_hat              # RQS vollständiges Modell
  Eff[s]          = (((eps_1_eps_1_hat-eps_eps_hat)/p_2)/ # F-Statistik
                       (eps_eps_hat/(n-p)))}
```

```{r, echo = F}
cat("F-Statistik für beta_2  = 0_{p_2}:", Eff[1],
    "\nF-Statistik für beta_2 != 0_{p_2}:", Eff[2])
```

# F-Statistiken
\footnotesize
\begin{theorem}[F-Statistik]
\justifying
\normalfont
Für $X \in \mathbb{R}^{n \times p}, \beta \in \mathbb{R}^p$ und $\sigma^2 > 0$
sei ein ALM der Form
\begin{equation}
y = X\beta + \varepsilon \mbox{ mit } \varepsilon \sim N(0_n,\sigma^2I_n)
\end{equation}
mit der Partitionierung
\begin{equation}
X      = \begin{pmatrix} X_1     & X_2      \end{pmatrix},
X_1      \in \mathbb{R}^{n\times p_1},
X_2      \in \mathbb{R}^{n\times p_2},
\mbox{ und }
\beta := \begin{pmatrix} \beta_1 \\ \beta_2 \end{pmatrix},
\beta_1 \in \mathbb{R}^{p_1},
\beta_2 \in \mathbb{R}^{p_2},
\end{equation}
mit $p = p_1 + p_2$ gegeben. Schließlich sei
\begin{equation}
K := \begin{pmatrix} 0_{p_1} \\ 1_{p_2} \end{pmatrix} \in \mathbb{R}^p
\end{equation}
ein Kontrastgewichtsvektor. Dann gilt
\begin{equation}
F \sim f(\delta, p_2, n-p) \mbox{ mit } \delta := \frac{K^T\beta \left(K^T(X^TX)^{-1}K\right)^{-1}K^T \beta}{\sigma^2}
\end{equation}
\end{theorem}

Bemerkungen

* Wir verzichten auf einen Beweis.
* $F$ ist eine Funktion der Parameterschätzer, $\delta$ ist eine Funktion der wahren, aber unbekannten, Parameter
* Diese Verteilung von $F$ kann zum Nullhypothesentesten und zur Powerfunktionsevaluation genutzt werden.


# F-Statistiken
\tiny
\setstretch{.8}

```{r}
# Libraries
library(MASS)                                               # Multivariate Normalverteilung

# Modellformulierung
nmod   = 2                                                  # Anzahl Modelle
n      = 10                                                 # Anzahl Datenpunkte
p      = 2                                                  # Anzahl Betaparameter
p_1    = 1                                                  # Anzahl Betaparameter im reduzierten Modell
p_2    = 1                                                  # Anzahl zusätzlicher Betaparameter im vollständigen Modell
p      = p_1 + p_2                                          # Anzahl Betaparameter im vollständigem Modell
x      = 1:n                                                # Prädiktorwerte
X      = matrix(c(rep(1,n),x), nrow = n)                    # Designmatrix des vollständigen Modells
X_1    = X[,1]                                              # Designmatrix des reduzierten Modells
I_n    = diag(n)                                            # n x n Einheitsmatrix
beta   = matrix(c(1,0,1,.5), nrow = 2)                      # wahre , aber unbekannte , Betaparameter
nscn   = ncol(beta)                                         # Anzahl wahrer, aber unbekannter, Hypothesenszenarien
sigsqr = 1                                                  # wahrer, aber unbekannter, Varianzparameter
K      = matrix((c(0,1)), nrow = 2)

# Frequentistische Simulation
nsim   = 1e4                                                # Anzahl Realisierungen des n-dimensionalen ZVs
delta  = rep(NaN,nscn)                                      # Nichtzentralitätsparameterarray
Eff    = matrix(rep(NaN, nscn*nsim), nrow = nscn)           # F-Statistiik Realisierungsarray
for(s in 1:nscn){                                           # Szenarieniterationen
  delta[s] = (t(t(K)%*%beta[,s])%*%                         # Nichtzentralitätsparameter
              solve(t(K)%*%solve(t(X)%*%X)%*%K) %*%
               (t(K)%*%beta[,s])/sigsqr)
  for(i in 1:nsim){                                         # Simulationsiterationen
    y               = mvrnorm(1, X %*%beta[,s], sigsqr*I_n) # Datenrealisierung
    beta_hat_1      = solve(t(X_1)%*%X_1)%*%t(X_1)%*%y      # Betaparameterschätzer reduziertes Modell
    beta_hat        = solve(t(X)  %*%X  )%*%t(X)  %*%y      # Betaparameterschätzer vollständiges Modell
    eps_1_hat       = y-X_1%*%beta_hat_1                    # Residuenvektor reduziertes Modell
    eps_hat         = y-X%*%beta_hat                        # Residuenvektor vollständiges Modell
    eps_1_eps_1_hat = t(eps_1_hat) %*% eps_1_hat            # RQS reduziertes Modell
    eps_eps_hat     = t(eps_hat)   %*% eps_hat              # RQS vollständiges Modell
    Eff[s,i]        = (((eps_1_eps_1_hat-eps_eps_hat)/p_2)/ # F-Statistik
                         (eps_eps_hat/(n-p)))
  }
}
```

# F-Statistiken
\vspace{1mm}
Beispiel (1)  Einfache lineare Regression
\vspace{2mm}

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

# Visualisierung
graphics.off()
dev.new()
par(
family      = "sans",
mfcol       = c(1,2),
pty         = "s",
bty         = "l",
lwd         = 1,
las         = 1,
mgp         = c(2.5,1,0),
xaxs        = "i",
yaxs        = "i",
xpd         = TRUE,
font.main   = 1,
cex         = 1,
cex.main    = 1.2)

# F-Statistik Ergebnisraum
lab         = c(TeX("$\\beta = (1,0)^T$"), TeX("$\\beta = (1,0.5)^T$") )
xlimits     = list(c(.05,6), c(0,80))
ylimits     = list(c(0,2), c(0,0.05))
breaks      = c(600,200)

# F-Teststatistiken
for(s in 1:nscn){
  xlims  = xlimits[[s]]
  f_min  = xlims[1]
  f_max  = xlims[2]
  f_res  = 1e3
  f      = seq(f_min, f_max, len = f_res)
  p_f    = df(f,p_2,n-p, delta[s])
  hist(
  Eff[s,],
  breaks = breaks[s],
  col    = "gray90",
  prob   = TRUE,
  xlab   = TeX("$F$"),
  ylab   = "",
  xlim   = xlims,
  ylim   = ylimits[[s]],
  main   = lab[s])
  lines(
  f,
  p_f,
  type  = "l",
  lwd   = 2,
  col   = "darkorange")
}

# Speichern
dev.copy2pdf(
file        = file.path(fdir, "alm_7_F_Statistik.pdf"),
width       = 8,
height      = 4)
```
\vfill
```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("7_Abbildungen/alm_7_F_Statistik.pdf")
```
\vfill

# F-Statistiken
Schlussbemerkungen

\footnotesize
* \justifying Die Theorie von T- und F-Statistiken kann unter dem Gesichtspunkt von linearen
Transformationen von Betaparametern verallgemeinert und integriert werden. Dabei
betrachtet man Hypothesen der Form
\begin{equation}
K\beta = \beta_0,
\end{equation}
wobei $K \in \mathbb{R}^{s \times p}$ eine beliebige Matrix ist. Zum Beispiel ergibt
sich für $K \in \mathbb{R}^{1 \times p}$ der hier betrachtete Fall von T-Statistiken
und für $K = I_p \in \mathbb{R}^{p\times p}$ der hier betrachtete Fall von F-Statistiken.
\vspace{2mm}

* \justifying Man kann weiterhin zeigen, dass im Wesentlichen gilt
\begin{equation}
F = T^2,
\end{equation}
dass also auch das Quadrat einer T-Statistik $f$-verteilt ist und T-Statistiken damit 
spezielle F-Statistiken sind. 
\vspace{2mm}

* \justifying Dennoch wird in der Anwendung sehr stark zwischen T- und F-Statistiken unterschieden
und es ist sinnvoll, sich der unterschiedlichen Anwenungsfälle von T- und F-Statistiken
bewusst zu sein. Eine sehr gute Einführung in die allgemeine Kontrasttheorie gibt
@searle_1971, Chapter 3.6.

#  {.plain}
\vfill
\large
\setstretch{2.5}
Schätzerverteilungen

T-Statistiken

F-Statistiken

**Selbstkontrollfragen**
\vfill

# Selbstkontrollfragen
\footnotesize
\setstretch{2.2}
1. Geben Sie das Theorem zur Frequentistischen Verteilung des Betaparameterschätzers wieder.
2. Geben Sie die Verteilung des Betaparameterschätzers im Szenarion von $n$ u.i.v. Zufallsvariablen wieder.
3. Geben Sie das Theorem zur Frequentistischen Verteilung des Varianzparameterschätzers wieder.
4. Geben Sie die Verteilung des skalierten Varianzparameterschätzers bei $n$ u.i.v. Zufallsvariablen wieder.
5. Skizzieren Sie die WDFen von $t$-Zufallsvariablen mit 2, 10 und 30 Freiheitsgraden.
6. Skizzieren Sie die WDFen von nichtzentralen $t$-Zufallsvariablen mit Nichtzentralitätsparametern 0,5 und 15.
7. Geben Sie die Definition der T-Statistik wieder.
8. Erläutern Sie die Definition der T-Statistik.
9. Warum kann die T-Statistik als Signal-zu-Rauschen Verhältnis interpretiert werden.
10. Geben Sie die Form der T-Statistik im Szenario von $n$ u.i.v. Zufallsvariablen wieder.
11. Erläutern Sie den Zusammenhang der T-Statistik und Cohen's $d$.
12. Geben Sie die Definition der T-Teststatistik wieder.

# Selbstkontrollfragen
\footnotesize
\setstretch{2.2}
13. Erläutern Sie die Definition der T-Teststatistik.
14. Geben Sie das Theorem zur Verteilung der T-Teststatistik wieder.
15. Geben Sie die Definition eines vollständigem und eines reduziertem ALMs wieder.
16. Geben Sie die Definition des Likelihood-Quotienten eines vollständigen und eines reduzierten ALMs wieder.
17. Geben Sie das Theorem zu Likelihood-Quotient und Residualquadratsummendifferenz wieder.
18. Definieren Sie die F-Statistik.
19. Erläutern Sie den Zusammenhang zwischen F-Statistik und Likelihood-Quotient.
20. Definieren Sie und erläutern Sie den Zähler der F-Statistik
21. Definieren Sie und erläutern Sie den Nenner der F-Statistik.
22. Erläutern Sie die F-Statistik.
23. Geben Sie das Theorem zur Verteilung der F-Statistik wieder.


# References
\footnotesize
\setstretch{2.2}