---
fontsize: 8pt
format:
  beamer:
    include-in-header: ../ALM_Header.tex
bibliography: ../ALM_Referenzen.bib
---

```{r, include = F}
source("../ALM_Common.R")
```

# {.plain}
<!-- Vorlesungstitel -->
\center
```{r, echo = F, out.width = "20%"}
knitr::include_graphics("../OvGU_Logo.png")
```

\vspace{2mm}
\huge
Allgemeines Lineares Modell

\vspace{6mm}
\large
BSc Psychologie, SoSe 2025

\vspace{5mm}
Joram Soch


<!-- Sitzung 11: Zweifaktorielle Varianzanalyse -->
# {.plain}

\vfill
\center
\huge
\textcolor{black}{(11) Zweifaktorielle Varianzanalyse}
\vfill


<!-- Inhaltsverzeichnis -->
#

\large
\setstretch{2.5}
\vfill

Anwendungsszenario

Modellformulierung

Modellschätzung

Modellevaluation

Anwendung/Praxis

Selbstkontrollfragen

\vfill


<!-- Abschnitt 1: Anwendungsszenario -->
#

\large
\setstretch{2.5}
\vfill

**Anwendungsszenario**

Modellformulierung

Modellschätzung

Modellevaluation

Anwendung/Praxis

Selbstkontrollfragen

\vfill


<!-- Grundbegriffe der ZVA -->
# Anwendungsszenario

\setstretch{1.5}
\textcolor{darkblue}{Randomisiertes zweifaktorielles Studiendesigns mit gekreuztem Design}

Zweifaktorielle Varianzanalyse (ZVA) = two-way analysis of variance (ANOVA):

\flushleft
\small
* Eine univariate abhängige Variable, bestimmt an randomisierten experimentellen Einheiten.
* Zwei diskrete unabhängige Variablen, die mindestens zweistufig sind.
* Die unabhängigen Variablen werden **Faktoren** genannt.
* Die Stufen der Faktoren werden auch **Level** der Faktoren genannt.
* Jedes Level eines Faktors wird mit allen Leveln des anderen Faktors kombiniert.
* Die Kombination zweier spezifischer Faktorlevel wird **Zelle** des Designs genannt.

\normalsize
Zweifaktorielle Studiendesigns werden üblicherweise anhand ihrer Faktorlevel bezeichnet:

\small
\center
\begin{tabular}{lll}
$2 \times 2$ ANOVA: & Faktor A mit Level 1,2     & Faktor B mit Level 1,2   \\
$2 \times 3$ ANOVA: & Faktor A mit Level 1,2     & Faktor B mit Level 1,2,3 \\
$4 \times 2$ ANOVA: & Faktor A mit Level 1,2,3,4 & Faktor B mit Level 1,2   \\
$3 \times 1$ ANOVA: & Faktor A mit Level 1,2,3   & Faktor B mit Level 1
\end{tabular}

\flushleft
* Generell sind **2 x 2 Designs** sehr populär, wir fokussieren auf diesen Fall.
* Die Zellen eines 2 x 2 Designs werden im Folgenden auch als **Gruppen** bezeichnet.


<!-- Konzept der ZVA -->
# Anwendungsszenario

\textcolor{darkblue}{Konzeptuelles Design}
\vspace{-2mm}

\center
```{r, echo = F, out.width = "70%"}
knitr::include_graphics("Abbildungen/zva_konzept.pdf")
```

\small
Die Zellen/Gruppen des Designs sind hier mit A1B1, A1B2, A2B1 und A2B2 bezeichnet.


<!-- Notation der ZVA -->
# Anwendungsszenario

\textcolor{darkblue}{Datennotation}
\vspace{-2mm}

\center
```{r, echo = F, out.width = "70%"}
knitr::include_graphics("Abbildungen/zva_datennotation.pdf")
```

\center
\small
$y_{ijk}$ bezeichnet die Datenvariable der $k$ten experimentellen Einheit ($k = 1,...,n_{ij}$)

im $i$ten Level von Faktor A und $j$ten Level von Faktor B ($i = 1,2, j = 1,2$).


<!-- Haupteffekte und Interaktionen -->
# Anwendungsszenario

\textcolor{darkblue}{Haupteffekte (\textit{main effects}) und Interaktionen (\textit{interactions})}

\small
Man unterscheidet intuitiv hinsichtlich der Gruppenmittelwerte *Haupteffekte* und *Interaktionen*:

\footnotesize
* \justifying Intuitiv spricht man vom Vorliegen eines *Haupteffekts von Faktor A*, wenn sich die Gruppenmittelwerte zwischen Level 1 und Level 2 von Faktor A, jeweils gemittelt über die zwei Level von Faktor B, unterscheiden.
* Intuitiv spricht man vom Vorliegen eines *Haupteffekts von Faktor B*, wenn sich die Gruppenmittelwerte zwischen Level 1 und Level 2 von Faktor B, jeweils gemittelt über die zwei Level von Faktor A, unterscheiden.
* Intuitiv spricht man vom Vorliegen einer *Interaktion der Faktoren A und B*, wenn der Unterschied der Gruppenmittelwerte von Faktor A zwischen Level 1 und 2 unterschiedlich für Level 1 und Level 2 von Faktor B ausgeprägt ist bzw. wenn der Unterschied der Gruppenmittelwerte von Faktor B zwischen Level 1 und 2 unterschiedlich für Level 1 und Level 2 von Faktor A ausgeprägt ist.

\small
Intuitiv beziehen sich Haupteffekte also auf (marginale) Unterschiede (Differenzen), während sich Interaktionen auf Unterschiede von Unterschieden (Differenzen von Differenzen) beziehen.

Das Vorhandensein einer Interaktion besagt lediglich, dass sich die Unterschiede der Gruppenmittelwerte zwischen den Leveln eines experimentellen Faktors in Abhängigkeit von den Leveln des anderen experimentellen Faktors ändern, es macht aber keine Aussage darüber, warum dies so ist.

\textcolor{darkblue}{$\Rightarrow$ Haupteffekte und Interaktionen sind Datenmuster, keine wissenschaftlichen Theorien.}


<!-- Anwendungsbeispiel -->
# Anwendungsszenario

\setstretch{1.5}
\textcolor{darkblue}{Anwendungsbeispiel}

BDI-Differenzwertanalyse für zwei CBT-Settings und CBT-Varianten

\vspace{2mm}
\center
```{r, echo = F, out.width = "70%"}
knitr::include_graphics("Abbildungen/beispielstudie.pdf")
```

```{r, echo = F, echo = F}
# Datensimulation
set.seed(1)                                                                         # Ergebnisreproduzierbarkeit
library(MASS)                                                                       # multivariate Normalverteilung
I        = 2                                                                        # Anzahl Level Faktor A (Therapiesetting)
J        = 2                                                                        # Anzahl Level Faktor B (Therapievariante)
n_ij     = 12                                                                       # balanciertes ANOVA Design
n        = I*J*n_ij                                                                 # Anzahl Datenpunkte
p        = 1 + (I-1)+(J-1)+(I*J-3)                                                  # Anzahl Parameter
D        = matrix(c(1,0,0,0,                                                        # Referenzgruppenregressor
                    1,1,0,0,                                                        # Haupteffekt A Regressor
                    1,0,1,0,                                                        # Haupteffekt B Regressor
                    1,1,1,1), nrow = p, byrow = TRUE)                               # Interaktionseffekt Regressor
C        = matrix(rep(1,n_ij),nrow = n_ij)                                          # prototypischer Zellenvektor für balancierte Designs
X        = kronecker(D,C)                                                           # Kroneckerprodukt für balancierte Designs
I_n      = diag(n)                                                                  # n x n Einheitsmatrix
beta     = matrix(c(10,4,0,0), nrow = p)                                            # \beta = (\mu_0,\alpha_2,\beta_2,\gamma_22)
sigsqr   = 8                                                                        # \sigma^2
y        = mvrnorm(1, X %*% beta, sigsqr*diag(n))                                   # Datenrealisierung
y        = round(y)                                                                 # diskrete BDI-Werte
D        = data.frame(Setting   = c(rep("F2F"  , n_ij*2),                           # face-to-face
                                    rep("ONL"  , n_ij*2)),                          # online
                      Variant   = c(rep("MND"  , n_ij),                             # CBT-M
                                    rep("EXC"  , n_ij),                             # CBT-E
                                    rep("MND"  , n_ij),                             # CBT-M
                                    rep("EXC"  , n_ij)),                            # CBT-E
                      dBDI      = y)                                                # BDI-Differenzwerte
fname    = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"                         # Dateiname
write.csv(D, file = fname, row.names = FALSE)                                       # Speicherung
```


<!-- Anwendungsbeispiel -->
# Anwendungsszenario

\setstretch{1.5}
\textcolor{darkblue}{Anwendungsbeispiel}

BDI-Differenzwertanalyse für zwei CBT-Settings und CBT-Varianten

\small
* Faktor (A) CBT-Setting  mit Level (1) face-to-face (F2F) und Level (2) online (ONL).
* Faktor (B) CBT-Variante mit Level (1) mindfulness (MND)  und Level (2) exercise (EXC).
* $n_{ij} := 12, \; 1 \le i,j \le 2$ Patient:innen in jeder Bedingung, $n = 48$ Patient:innen insgesamt.

\normalsize
2 x 2 randomisiertes zweifaktorielles Studiendesign $\Rightarrow$ 2 x 2 ANOVA

\small
* inferentielle Evidenz für einen Haupteffekt des CBT-Settings?
* inferentielle Evidenz für einen Haupteffekt der CBT-Variante?
* inferentielle Evidenz für Interaktion von CBT-Setting und CBT-Variante?


<!-- Abschnitt 2: Modellformulierung -->
#

\large
\setstretch{2.5}
\vfill

Anwendungsszenario

**Modellformulierung**

Modellschätzung

Modellevaluation

Anwendung/Praxis

Selbstkontrollfragen

\vfill


<!-- Einführung: Modell der additiven ZVA -->
# Modellformulierung

\vspace{2mm}
\textcolor{darkblue}{Modell der additiven ZVA}

\footnotesize
\setstretch{1.2}
In Analogie zur einfaktoriellen Varianzanalyse (EVA) möchte man in der additiven ZVA die Gruppenerwartungswerte $\mu_{ij}$ mit $i = 1,...,I$ für die Level von Faktor A und $j = 1,...,J$ für die Level von Faktor B als Summe eines gruppenunspezifischen Erwartungswertes und den Effekten der Level von Faktor A und der Level von Faktor B modellieren.

Wir bezeichnen den gruppenunspezifischen Erwartungswertparameter mit $\mu_0$, den Effekt von Level $i$ von Faktor A mit $\alpha_i$ und den Effekt von Level $j$ von Faktor B mit $\beta_j$ ($\beta_j$ bezeichnet hier also \underline{nicht} den $j$ten Eintrag des Betaparametervektors). Dann ergibt sich zum Beispiel für $I := J := 2$:
\begin{center}
\begin{tabular}{l|l}
$\mu_{11} := \mu_0 + \alpha_1 + \beta_1$ & $\mu_{12} := \mu_0 + \alpha_1 + \beta_2$ \\ \hline
$\mu_{21} := \mu_0 + \alpha_2 + \beta_1$ & $\mu_{22} := \mu_0 + \alpha_2 + \beta_2$ \\
\end{tabular}
\end{center}

Wie im Falle der EVA ist diese Darstellung der Gruppenerwartungswerte $\mu_{ij}$ allerdings überparameterisiert. Um eine eindeutige Darstellung der $\mu_{ij}$ zu gewährleisten, bietet sich auch hier die Restriktion an, den Effekt des ersten Levels jedes Faktors als Null zu definieren:
\begin{equation}
\alpha_1 := \beta_1 := 0 \; .
\end{equation}

Dies etabliert die Faktorlevelkombination A1B1 als Referenzgruppe. Es ergibt sich somit zum Beispiel für für $I := J := 2$:
\begin{center}
\begin{tabular}{l|l}
$\mu_{11} := \mu_0$            & $\mu_{12} := \mu_0 + \beta_2$            \\ \hline
$\mu_{21} := \mu_0 + \alpha_2$ & $\mu_{22} := \mu_0 + \alpha_2 + \beta_2$ \\
\end{tabular}
\end{center}

Auch bei dieser Effektdarstellung des Modells der additiven 2 x 2 ZVA mit Referenzgruppe ändern sich die Interpretation der Parameter $\mu_0,\alpha_2,\beta_2$ im Vergleich zum überparameterisierten Fall ohne Referenzgruppe: $\mu_0$ entspricht dem Erwartungswert der Faktorlevelkombination A1B1, $\alpha_2$ der Differenz beim Übergang von Level 1 zu Level 2 von Faktor A und  $\beta_2$ der Differenz beim Übergang von Level 1 zu Level 2 von Faktor B.


<!-- Definition: Modell der additiven ZVA mit Referenzgruppe -->
# Modellformulierung

\footnotesize
\begin{definition}[Modell der additiven ZVA mit Referenzgruppe]

\justifying
$y_{ijk}$ mit $i = 1,...,I$, $j = 1,...,J$, $k = 1,...,n_{ij}$ sei die Zufallsvariable, die den $k$ten Datenpunkt zum $i$ten Level von Faktor A und dem $j$ten Level von Faktor B in einem ZVA-Anwendungsszenario modelliert. Dann hat das \textit{Modell der additiven ZVA mit Referenzgruppe} die strukturelle Form
\begin{equation}
y_{ijk} = \mu_{ij} + \varepsilon_{ijk}
\quad \mbox{mit} \quad
\varepsilon_{ijk} \sim N(0,\sigma^2)
\quad \mbox{u.i.v. für} \quad
i = 1,...,I, \; j = 1,...,J, \; k = 1,...,n_{ij}
\end{equation}

und die Datenverteilungsform
\begin{equation}
y_{ijk} \sim N(\mu_{ij}, \sigma^2)
\quad \mbox{u.v. für} \quad
i = 1,...,I, \; j = 1,...,J, \; k = 1,...,n_{ij}
\end{equation}

mit
\begin{equation}
\mu_{ij} := \mu_0 + \alpha_i + \beta_j
\quad \mbox{für} \quad
i = 1,...,I, \; j = 1,...,J
\quad \mbox{mit} \quad
\alpha_1 := \beta_1 := 0
\end{equation}

und $\sigma^2 > 0$.

\end{definition}

Bemerkungen

* Das Modell der additiven ZVA modelliert ausschließlich Haupteffekte, keine Interaktionen.


<!-- Beispiele: Parameter -->
# Modellformulierung

\textcolor{darkblue}{Parameterbeispiele}

\small
\setstretch{1.5}
\noindent (1) Es sei $\mu_0 := 1, \alpha_2 := 1, \beta_2 := 0$. Dann gilt:

\begin{center}
\begin{tabular}{ll}
$\mu_{11} = \mu_0 + \alpha_1 + \beta_1 = 1 + 0 + 0 = 1$ &
$\mu_{12} = \mu_0 + \alpha_1 + \beta_2 = 1 + 0 + 0 = 1$ \\
$\mu_{21} = \mu_0 + \alpha_2 + \beta_1 = 1 + 1 + 0 = 2$ &
$\mu_{22} = \mu_0 + \alpha_2 + \beta_2 = 1 + 1 + 0 = 2$ \\
\end{tabular}
\end{center}

\vspace{1mm}
$\Rightarrow$ Haupteffekt von Faktor A, kein Haupteffekt von Faktor B

\vspace{1mm}
\noindent (2) Es sei $\mu_0 := 1, \alpha_2 := 0, \beta_2 := 1$. Dann gilt:

\begin{center}
\begin{tabular}{ll}
$\mu_{11} = \mu_0 + \alpha_1 + \beta_1 = 1 + 0 + 0 = 1$ &
$\mu_{12} = \mu_0 + \alpha_1 + \beta_2 = 1 + 0 + 1 = 2$ \\
$\mu_{21} = \mu_0 + \alpha_2 + \beta_1 = 1 + 0 + 0 = 1$ &
$\mu_{22} = \mu_0 + \alpha_2 + \beta_2 = 1 + 0 + 1 = 2$ \\
\end{tabular}
\end{center}

\vspace{1mm}
$\Rightarrow$ kein Haupteffekt von Faktor A, Haupteffekt von Faktor B

\vspace{1mm}
\noindent (3) Es sei $\mu_0 := 1, \alpha_2 := 1, \beta_2 := 1$. Dann gilt:

\begin{center}
\begin{tabular}{ll}
$\mu_{11} = \mu_0 + \alpha_1 + \beta_1 = 1 + 0 + 0 = 1$ &
$\mu_{12} = \mu_0 + \alpha_1 + \beta_2 = 1 + 0 + 1 = 2$ \\
$\mu_{21} = \mu_0 + \alpha_2 + \beta_1 = 1 + 1 + 0 = 2$ &
$\mu_{22} = \mu_0 + \alpha_2 + \beta_2 = 1 + 1 + 1 = 3$ \\
\end{tabular}
\end{center}

\vspace{1mm}
$\Rightarrow$ Haupteffekt von Faktor A, Haupteffekt von Faktor B


<!-- Beispiele: Visualisierung -->
# Modellformulierung

\vspace{2mm}
\textcolor{darkblue}{Parameterbeispiele}
\vspace{2mm}

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

# Parameterbeispiele
mu          = list(matrix(c(1,1,2,2) , ncol = 2),
                   matrix(c(1,2,1,2) , ncol = 2),
                   matrix(c(1,2,2,3) , ncol = 2))
labs        = list(TeX("$\\mu_0 = 1,\\, \\alpha_2 = 1,\\, \\beta_2 = 0$"),
                   TeX("$\\mu_0 = 1,\\, \\alpha_2 = 0,\\, \\beta_2 = 1$"),
                   TeX("$\\mu_0 = 1,\\, \\alpha_2 = 1,\\, \\beta_2 = 1$"))
for(i in 1:3){

    # Barplot
    B        = as.vector(mu[[i]])
    names(B) = c("A1B1", "A1B2", "A2B1", "A2B2")
    barplot(B,
        ylim        = c(0,3.5),
        main        = labs[[i]],
        col         = "gray90")

    # Lineplot
    matplot(1:2, mu[[i]],
        type        = "b",
        pch         = 21,
        ylim        = c(  0,   5),
        xlim        = c(0.5, 2.5),
        lty         = 1,
        col         = c("gray30", "gray70"),
        lwd         = c(3,1),
        main        = labs[[i]],
        xlab        = "Faktor B",
        ylab        = "",
        xaxt        = "n")
    legend(0.01, 6.00, c("Faktor A, Level 1", "Faktor A, Level 2"),
        lty         = 1,
        col         = c("gray30", "gray90"),
        lwd         = c(3,1),
        bty         = "n",
        seg.len     = 0.4,
        x.intersp   = 0.3,
        y.intersp   = 0.2,
        cex         = 1)
    text(1, -0.4, "Level 1", xpd = T)
    text(2, -0.4, "Level 2", xpd = T)
}

# PDF-Speicherung
dev.copy2pdf(
    file        = "Abbildungen/zva_beispiel_add.pdf",
    width       = 15,
    height      = 10)
```

\center
```{r, echo = F, out.width = "100%"}
knitr::include_graphics("Abbildungen/zva_beispiel_add.pdf")
```
\vfill


<!-- Theorem: Designmatrixform des Modells der additiven 2 x 2 ZVA mit Referenzgruppe -->
# Modellformulierung

\small
\vspace{2mm}
\setstretch{1.2}
\begin{theorem}[\small Designmatrixform des Modells der additiven 2 x 2 ZVA mit Referenzgruppe]

\justifying
\normalfont
Gegeben sei die strukturelle Form des Modells der additiven 2 x 2 ZVA mit Referenzgruppe. Dann hat dieses Modell die Designmatrixform
\begin{equation}
y = X\beta + \varepsilon
\quad \mbox{mit} \quad
\varepsilon \sim N(0_n,\sigma^2 I_n),
\end{equation}

wobei
\begin{equation}
y :=
\begin{pmatrix}
y_{111}      \\
\vdots       \\
y_{11n_{11}} \\
y_{121}      \\
\vdots       \\
y_{12n_{12}} \\
y_{211}      \\
\vdots       \\
y_{21n_{21}} \\
y_{221}      \\
\vdots       \\
y_{22n_{22}} \\
\end{pmatrix}, \quad
X =
\begin{pmatrix}
1       &  0       &  0       \\
\vdots  &  \vdots  &  \vdots  \\
1       &  0       &  0       \\
1       &  0       &  1       \\
\vdots  &  \vdots  &  \vdots  \\
1       &  0       &  1       \\
1       &  1       &  0       \\
\vdots  &  \vdots  &  \vdots  \\
1       &  1       &  0       \\
1       &  1       &  1       \\
\vdots  &  \vdots  &  \vdots  \\
1       &  1       &  1       \\
\end{pmatrix} \in \mathbb{R}^{n \times 3}, \quad
\beta :=
\begin{pmatrix}
\mu_0      \\
\alpha_{2} \\
\beta_{2}  \\
\end{pmatrix}
\in \mathbb{R}^{3}
\quad \mbox{und} \quad
\sigma^2 > 0.
\end{equation}

\end{theorem}

\vspace{-2mm}
Bemerkung
\vspace{-2mm}

* Das Theorem folgt direkt aus den Regel der Matrixmultiplikation.


<!-- Designmatrix des Modells -->
# Modellformulierung

Designmatrix des Modells der additiven ZVA mit Referenzgruppe ($n = 16$, $p = 3$)

```{r, echo = F, eval = F}
# Designmatrixerzeugung
I      = 2                                  # Anzahl Level Faktor A
J      = 2                                  # Anzahl Level Faktor B
n_ij   = 4                                  # Anzahl von Datenpunkten der (i,j)ten Gruppe
n      = I*J*n_ij                           # Anzahl Datenpunkte insgesamt
p      = 1 + (I-1)+(J-1)+(I*J-3)            # Anzahl Parameter insgesamt
D      = matrix(c(1,0,0,0,                  # prototypische Designmatrix für balancierte Designs
                  1,0,1,0,
                  1,1,0,0,
                  1,1,1,1),
                nrow  = p,
                byrow = TRUE)
C      = matrix(rep(1,n_ij), nrow = n_ij)   # prototypischer Zellenvektor für balancierte Designs
X      = kronecker(D,C)                     # Kroneckerprodukt-Designmatrix-Erzeugung für balancierte Designs
X      = X[,-4]                             # eliminiere vierte Spalte
Xp     = X

# Abbildungsparameter
library(plot.matrix)
graphics.off()
dev.new()
par(
    family    = "sans",
    mar       = c(1,1,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  = 2)

# Designmatrix
plot(Xp,
    col       = gray(seq(0, 1, length.out = 256)),
    border    = "black",
    xlab      = '',
    ylab      = '',
    main      = '',
    breaks    = seq(-2, 1, length.out = 257),
    key       = NULL,
    axis.col  = NULL,
    axis.row  = NULL,
    asp       = 1)

# Speichern
dev.copy2pdf(
    file      = "Abbildungen/zva_designmatrix_add.pdf",
    width     = 10,
    height    = 10)
```

\vspace{2mm}
```{r, echo = F, out.width = "60%"}
knitr::include_graphics("Abbildungen/zva_designmatrix_add.pdf")
```


<!-- Einführung: Modell der ZVA mit Interaktion -->
# Modellformulierung

\vspace{2mm}
\textcolor{darkblue}{Modell der ZVA mit Interaktion}

\footnotesize
\setstretch{1.1}
In der ZVA mit Interaktion möchte man die Gruppenerwartungswerte $\mu_{ij}$ mit $i = 1,...,I$ für die Level von Faktor A und $j = 1,...,J$ für die Level von Faktor B als Summe eines gruppenunspezifischen Erwartungswertes, der Effekte der Level von Faktor A und Faktor B und der Interaktion der Level der Faktoren modellieren.

Wir bezeichnen den gruppenunspezifischen Erwartungswertparameter mit $\mu_0$, den Effekt von Level $i$ von Faktor A mit $\alpha_i$, den Effekt von Level $j$ von Faktor B mit $\beta_j$, und die Interaktion von Level $i$ von Faktor A mit Level $j$ von Faktor B mit $\gamma_{ij}$. Dann ergibt sich zum Beispiel für $I := J := 2$:
\begin{center}
\begin{tabular}{l|l}
$\mu_{11} := \mu_0 + \alpha_1 + \beta_1 + \gamma_{11}$ & $\mu_{12} := \mu_0 + \alpha_1 + \beta_2 + \gamma_{12}$ \\ \hline
$\mu_{21} := \mu_0 + \alpha_2 + \beta_1 + \gamma_{21}$ & $\mu_{22} := \mu_0 + \alpha_2 + \beta_2 + \gamma_{22}$ \\
\end{tabular}
\end{center}

Wie in der die additiven ZVA ist diese Darstellung der Gruppenerwartungswerte $\mu_{ij}$ multipel überparameterisiert. Um eine eindeutige Darstellung der $\mu_{ij}$ zu gewährleisten, bietet sich auch hier die Restriktion an, den Effekt des ersten Levels jedes Faktors und jeder Interaktion als Null zu definieren:
\begin{equation}
\alpha_1 := \beta_1 := \gamma_{i1} := \gamma_{1j} := 0
\quad \mbox{für} \quad
i = 1,...,I, \; j = 1,..,J
\end{equation}

Dies etabliert die Faktorlevelkombination A1B1 als Referenzgruppe. Es ergibt sich somit zum Beispiel für $I := J := 2$:
\begin{center}
\begin{tabular}{l|l}
$\mu_{11} := \mu_0$            & $\mu_{12} := \mu_0 + \beta_2$                          \\ \hline
$\mu_{21} := \mu_0 + \alpha_2$ & $\mu_{22} := \mu_0 + \alpha_2 + \beta_2 + \gamma_{22}$ \\
\end{tabular}
\end{center}

Auch bei dieser Effektdarstellung des Modells der 2 x 2 ZVA mit Interaktion und Referenzgruppe ändert sich die Interpretation der Parameter $\mu_0,\alpha_2,\beta_2,\gamma_{22}$ im Vergleich zum überparameterisierten Fall ohne Referenzgruppe: $\mu_0$ entspricht dem Erwartungswert der Faktorlevelkombination A1B1, $\alpha_2$ der Differenz beim Übergang von Level 1 zu Level 2 von Faktor A, $\beta_2$ der Differenz beim Übergang von Level 1 zu Level 2 von Faktor B und $\gamma_{22}$ der Differenz beim Übergang von Level 1 zu Level 2 von Faktor B im Unterschiede zum Übergang von Level 1 zu Level 2 von Faktor A.


<!-- Definition: Modell der ZVA mit Interaktion und Referenzgruppe -->
# Modellformulierung

\footnotesize
\begin{definition}[Modell der ZVA mit Interaktion und Referenzgruppe]

\justifying
$y_{ijk}$ mit $i = 1,...,I$, $j = 1,...,J$, $k = 1,...,n_{ij}$ sei die Zufallsvariable, die den $k$ten Datenpunkt zum $i$ten Level von Faktor A und dem $j$ten Level von Faktor B in einem ZVA-Anwendungsszenario modelliert. Dann hat das \textit{Modell der ZVA mit Interaktion und Referenzgruppe} die strukturelle Form
\begin{equation}
y_{ijk} = \mu_{ij} + \varepsilon_{ijk}
\quad \mbox{mit} \quad
\varepsilon_{ijk} \sim N(0,\sigma^2)
\quad \mbox{u.i.v. für} \quad
i = 1,...,I, \; j = 1,...,J, \; k = 1,...,n_{ij}
\end{equation}

und die Datenverteilungsform
\begin{equation}
y_{ijk} \sim N(\mu_{ij}, \sigma^2)
\quad \mbox{u.v. für} \quad
i = 1,...,I, \; j = 1,...,J, \; k = 1,...,n_{ij}
\end{equation}

mit
\begin{equation}
\mu_{ij} := \mu_0 + \alpha_i + \beta_j + \gamma_{ij}
\end{equation}

sowie
\begin{equation}
\alpha_1 := \beta_1 := \gamma_{i1} := \gamma_{1j} :=0
\quad \mbox{für} \quad
i = 1,...,I, \; j = 1,...,J
\end{equation}

und $\sigma^2 > 0$.

\end{definition}

Bemerkungen

* Das Modell der ZVA mit Interaktion beinhaltet durch die Terme $\gamma_{ij}$ für $i = 1,...,I$, $j = 1,...,J$ Interaktionseffekte, die Abweichungen vom Modell der additiven ZVA mit ausschließlich Haupteffekten repräsentieren.


<!-- Beispiele: Parameter -->
# Modellformulierung

\vspace{1mm}
\setstretch{1}
\textcolor{darkblue}{Parameterbeispiele}

\vspace{1mm}
\footnotesize
\noindent (1) Es sei $\mu_0 := 1, \alpha_2 := 0, \beta_2 := 0, \gamma_{22} = 2$. Dann gilt:
\begin{tiny}
\begin{tabular}{ll}
$\mu_{11} = \mu_0 + \alpha_1 + \beta_1 + \gamma_{11} = 1 + 0 + 0 + 0=  1$ &
$\mu_{12} = \mu_0 + \alpha_1 + \beta_2 + \gamma_{12} = 1 + 0 + 0 + 0 = 1$ \\\\
$\mu_{21} = \mu_0 + \alpha_2 + \beta_1 + \gamma_{21} = 1 + 0 + 0 + 0 = 1$ &
$\mu_{22} = \mu_0 + \alpha_2 + \beta_2 + \gamma_{22} = 1 + 0 + 0 + 2 = 3$ \\
\end{tabular}
\end{tiny}

\vspace{1mm}
$\Rightarrow$ kein Haupteffekt von Faktor A, kein Haupteffekt von Faktor B, Interaktion von A und B

\vspace{2mm}
\noindent (2) Es sei $\mu_0 := 1, \alpha_2 := 1, \beta_2 := 1, \gamma_{22} = -2$. Dann gilt:
\begin{tiny}
\begin{tabular}{ll}
$\mu_{11} = \mu_0 + \alpha_1 + \beta_1 + \gamma_{11} = 1 + 0 + 0 + 0=  1$ &
$\mu_{12} = \mu_0 + \alpha_1 + \beta_2 + \gamma_{12} = 1 + 0 + 1 + 0 = 2$ \\\\
$\mu_{21} = \mu_0 + \alpha_2 + \beta_1 + \gamma_{21} = 1 + 0 + 1 + 0 = 2$ &
$\mu_{22} = \mu_0 + \alpha_2 + \beta_2 + \gamma_{22} = 1 + 1 + 1 - 2 = 1$ \\
\end{tabular}
\end{tiny}

\vspace{1mm}
$\Rightarrow$ Haupteffekt von Faktor A, Haupteffekt von Faktor B, negative Interaktion von A und B

\vspace{2mm}
\noindent (3) Es sei $\mu_0 := 1, \alpha_2 := 1, \beta_2 := 0, \gamma_{22} = 1$. Dann gilt:
\begin{tiny}
\begin{tabular}{ll}
$\mu_{11} = \mu_0 + \alpha_1 + \beta_1 + \gamma_{11} = 1 + 0 + 0 + 0 = 1$ &
$\mu_{12} = \mu_0 + \alpha_1 + \beta_2 + \gamma_{12} = 1 + 0 + 0 + 0 = 1$ \\\\
$\mu_{21} = \mu_0 + \alpha_2 + \beta_1 + \gamma_{21} = 1 + 1 + 0 + 0 = 2$ &
$\mu_{22} = \mu_0 + \alpha_2 + \beta_2 + \gamma_{22} = 1 + 1 + 0 + 1 = 3$ \\
\end{tabular}
\end{tiny}

\vspace{1mm}
$\Rightarrow$ Haupteffekt von Faktor A, kein Haupteffekt von Faktor B, Interaktion von A und B

\vspace{2mm}
\noindent (4) Es sei $\mu_0 := 1, \alpha_2 := 0, \beta_2 := 1, \gamma_{22} = 1$. Dann gilt:
\begin{tiny}
\begin{tabular}{ll}
$\mu_{11} = \mu_0 + \alpha_1 + \beta_1 + \gamma_{11} = 1 + 0 + 0 + 0 = 1$ &
$\mu_{12} = \mu_0 + \alpha_1 + \beta_2 + \gamma_{12} = 1 + 0 + 1 + 0 = 2$ \\\\
$\mu_{21} = \mu_0 + \alpha_2 + \beta_1 + \gamma_{21} = 1 + 0 + 0 + 0 = 1$ &
$\mu_{22} = \mu_0 + \alpha_2 + \beta_2 + \gamma_{22} = 1 + 0 + 1 + 1 = 3$ \\
\end{tabular}
\end{tiny}

\vspace{1mm}
$\Rightarrow$ kein Haupteffekt von Faktor A, Haupteffekt von Faktor B, Interaktion von A und B


<!-- Beispiele: Visualisierung -->
# Modellformulierung

\vspace{2mm}
\textcolor{darkblue}{Parameterbeispiele}
\vspace{2mm}

```{r, echo = F, eval = F}
# Abbildungsparameter
graphics.off()
dev.new()
par(
    mfcol       = c(2,4),
    family      = "sans",
    pty         = "s",
    bty         = "l",
    lwd         = 1,
    las         = 1,
    mgp         = c(2,1,0),
    xaxs        = "i",
    yaxs        = "i",
    font.main   = 1,
    cex         = 1.1,
    cex.main    = 1)
mu          = list(matrix(c(1,1,1,3), ncol = 2),
                   matrix(c(1,2,2,1), ncol = 2),
                   matrix(c(1,1,2,3), ncol = 2),
                   matrix(c(1,2,1,3), ncol = 2))
labs        = list(TeX("$\\mu_0 = 1, \\, \\alpha_2 = 0, \\, \\beta_2 =  0,\\, \\gamma_{22} =  2$"),
                   TeX("$\\mu_0 = 1, \\, \\alpha_2 = 1, \\, \\beta_2 =  1,\\, \\gamma_{22} = -2$"),
                   TeX("$\\mu_0 = 1, \\, \\alpha_2 = 1, \\, \\beta_2 =  0,\\, \\gamma_{22} =  1$"),
                   TeX("$\\mu_0 = 1, \\, \\alpha_2 = 0, \\, \\beta_2 =  1,\\, \\gamma_{22} =  1$"))

# Parameterbeispiele
for(i in 1:4){

    # Barplot
    B        = as.vector(mu[[i]])
    names(B) = c("A1B1", "A1B2", "A2B1", "A2B2")
    barplot(B,
        ylim        = c(0,4),
        main        = labs[[i]],
        col         = "gray90")

    # Lineplot
    matplot(1:2, mu[[i]],
        type        = "b",
        pch         = 21,
        ylim        = c( 0,4.5),
        xlim        = c(.5,2.5),
        lty         = 1,
        col         = c("gray30", "gray70"),
        lwd         = c(3,1),
        main        = labs[[i]],
        xlab        = "Faktor B",
        ylab        = "",
        xaxt        = "n")
    legend(0.01, 6.00, c("Faktor A, Level 1", "Faktor A, Level 2"),
        lty         = 1,
        col         = c("gray30", "gray90"),
        lwd         = c(3,1),
        bty         = "n",
        seg.len     = 0.4,
        x.intersp   = 0.3,
        y.intersp   = 0.2,
        cex         = 1)
    text(1, -0.5, "Level 1", xpd = T)
    text(2, -0.5, "Level 2", xpd = T)
}

# PDF-Speicherung
dev.copy2pdf(
    file        = "Abbildungen/zva_beispiel_int.pdf",
    width       = 15,
    height      = 10)
```

\center
```{r, echo = F, out.width = "100%"}
knitr::include_graphics("Abbildungen/zva_beispiel_int.pdf")
```
\vfill


<!-- Theorem: Designmatrixform des Modells der 2 x 2 ZVA mit Interaktion und Referenzgruppe -->
# Modellformulierung

\small
\vspace{2mm}
\setstretch{1.2}
\begin{theorem}[\small Designmatrixform des Modells der 2 x 2 ZVA mit Interaktion und Referenzgruppe]

\justifying
\normalfont
Gegeben sei die strukturelle Form eines 2 x 2 ZVA-Modells mit Interaktion und Referenzgruppe und es sei $n := \sum_{i=1}^2 \sum_{j=1}^2 n_{ij}$ die Gesamtanzahl an Datenvariablen. Dann hat dieses Modell die Designmatrixform
\begin{equation}
y = X\beta + \varepsilon
\quad \mbox{mit} \quad
\varepsilon \sim N(0_n,\sigma^2 I_n),
\end{equation}

wobei
\begin{equation}
y :=
\begin{pmatrix}
y_{111}      \\
\vdots       \\
y_{11n_{11}} \\
y_{121}      \\
\vdots       \\
y_{12n_{12}} \\
y_{211}      \\
\vdots       \\
y_{21n_{21}} \\
y_{221}      \\
\vdots       \\
y_{22n_{22}} \\
\end{pmatrix}, \quad
X =
\begin{pmatrix}
1       &  0       &  0       &  0       \\
\vdots  &  \vdots  &  \vdots  &  \vdots  \\
1       &  0       &  0       &  0       \\
1       &  0       &  1       &  0       \\
\vdots  &  \vdots  &  \vdots  &  \vdots  \\
1       &  0       &  1       &  0       \\
1       &  1       &  0       &  0       \\
\vdots  &  \vdots  &  \vdots  &  \vdots  \\
1       &  1       &  0       &  0       \\
1       &  1       &  1       &  1       \\
\vdots  &  \vdots  &  \vdots  &  \vdots  \\
1       &  1       &  1       &  1       \\
\end{pmatrix} \in \mathbb{R}^{n \times 4}, \quad
\beta :=
\begin{pmatrix}
\mu_0       \\
\alpha_{2}  \\
\beta_{2}   \\
\gamma_{22}
\end{pmatrix}
\in \mathbb{R}^{4}
\quad \mbox{und} \quad
\sigma^2 > 0.
\end{equation}

\end{theorem}


<!-- Designmatrix des Modells -->
# Modellformulierung

Designmatrix des Modells der ZVA mit Interaktion ($n = 16$, $p = 4$)

```{r, echo = F, eval = F}
# Designmatrixerzeugung
I      = 2                                  # Anzahl Level Faktor A
J      = 2                                  # Anzahl Level Faktor B
n_ij   = 4                                  # Anzahl von Datenpunkten der (i,j)ten Gruppe
n      = I*J*n_ij                           # Anzahl Datenpunkte insgesamt
p      = 1 + (I-1)+(J-1)+(I*J-3)            # Anzahl Parameter insgesamt
D      = matrix(c(1,0,0,0,                  # prototypische Designmatrix für balancierte Designs
                  1,0,1,0,
                  1,1,0,0,
                  1,1,1,1),
                nrow  = p,
                byrow = TRUE)
C      = matrix(rep(1,n_ij), nrow = n_ij)   # prototypischer Zellenvektor für balancierte Designs
X      = kronecker(D,C)                     # Kroneckerprodukt-Designmatrix-Erzeugung für balancierte Designs
Xp     = X

# Abbildungsparameter
library(plot.matrix)
graphics.off()
dev.new()
par(
    family    = "sans",
    mar       = c(1,1,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  = 2)

# Designmatrix
plot(Xp,
    col       = gray(seq(0, 1, length.out = 256)),
    border    = "black",
    xlab      = '',
    ylab      = '',
    main      = '',
    breaks    = seq(-2, 1, length.out = 257),
    key       = NULL,
    axis.col  = NULL,
    axis.row  = NULL,
    asp       = 1)

# Speichern
dev.copy2pdf(
    file      = "Abbildungen/zva_designmatrix_int.pdf",
    width     = 10,
    height    = 10)
```

\vspace{2mm}
```{r, echo = F, out.width = "60%"}
knitr::include_graphics("Abbildungen/zva_designmatrix_int.pdf")
```


<!-- Beispiel: Designmatrix -->
# Modellformulierung

\textcolor{darkblue}{Beispiel}

\footnotesize
Es seien
\begin{equation}
I := 2, J := 2
\quad \mbox{und} \quad
n_{ij} := 4
\quad \mbox{für} \quad
i = 1,2, \; j = 1,2,
\quad \mbox{also} \quad
n = 16.
\end{equation}

Dann gilt
\begin{equation}
y = X\beta + \varepsilon
\quad \mbox{mit} \quad
\varepsilon \sim N(0_{16},\sigma^2 I_{16})
\end{equation}

mit
\footnotesize
\begin{equation}
y :=
\begin{pmatrix}
y_{111} \\
y_{112} \\
y_{113} \\
y_{114} \\
y_{121} \\
y_{122} \\
y_{123} \\
y_{124} \\
y_{211} \\
y_{212} \\
y_{213} \\
y_{214} \\
y_{221} \\
y_{222} \\
y_{223} \\
y_{224}
\end{pmatrix}, \quad
X :=
\begin{pmatrix}
1  &  0  &  0  &  0  \\
1  &  0  &  0  &  0  \\
1  &  0  &  0  &  0  \\
1  &  0  &  0  &  0  \\
1  &  0  &  1  &  0  \\
1  &  0  &  1  &  0  \\
1  &  0  &  1  &  0  \\
1  &  0  &  1  &  0  \\
1  &  1  &  0  &  0  \\
1  &  1  &  0  &  0  \\
1  &  1  &  0  &  0  \\
1  &  1  &  0  &  0  \\
1  &  1  &  1  &  1  \\
1  &  1  &  1  &  1  \\
1  &  1  &  1  &  1  \\
1  &  1  &  1  &  1  \\
\end{pmatrix}
\in \mathbb{R}^{16 \times 4}, \quad
\beta :=
\begin{pmatrix}
\mu_0       \\
\alpha_2    \\
\beta_2     \\
\gamma_{22} \\
\end{pmatrix}
\in \mathbb{R}^{4}
\quad \mbox{und} \quad
\sigma^2 > 0.
\end{equation}


<!-- Beispiel: Designmatrix -->
# Modellformulierung

\vspace{2mm}
\textcolor{darkblue}{Beispiel}
\vspace{2mm}

\tiny
\setstretch{0.85}
```{r, echo = T}
# Modellformulierung
library(MASS)                               # multivariate Normalverteilung
I      = 2                                  # Anzahl Level Faktor A
J      = 2                                  # Anzahl Level Faktor B
n_ij   = 4                                  # Anzahl von Datenpunkten der i,jten Gruppe
n      = I*J*n_ij                           # Anzahl Datenpunkte insgesamt
p      = 1 + (I-1)+(J-1)+(I*J-3)            # Anzahl Parameter insgesamt
D      = matrix(c(1,0,0,0,                  # prototypische Designmatrix für balancierte Designs
                  1,0,1,0,
                  1,1,0,0,
                  1,1,1,1),
                nrow  = p,
                byrow = TRUE)
C      = matrix(rep(1,n_ij), nrow = n_ij)   # prototypischer Zellenvektor für balancierte Designs
X      = kronecker(D,C)                     # Kroneckerprodukt-Designmatrix-Erzeugung für balancierte Designs
I_n    = diag(n)                            # n x n Einheitsmatrix
beta   = matrix(c(1,1,1,1), nrow = p)       # \beta = (\mu_0,\alpha_2,\alpha_3,\alpha_4)
sigsqr = 10                                 # \sigma^2

# Datenrealisierung
y      = mvrnorm(1, X %*% beta, sigsqr*I_n) # eine Realisierung des n-dimensionalen ZVs y
print(X)
```


<!-- Abschnitt 3: Modellschätzung -->
#

\large
\setstretch{2.5}
\vfill

Anwendungsszenario

Modellformulierung

**Modellschätzung**

Modellevaluation

Anwendung/Praxis

Selbstkontrollfragen

\vfill


<!-- Theorem: Betaparameterschätzung im additiven 2 x 2 ZVA-Modell mit Referenzgruppe -->
<!-- Beweis -->
# Modellschätzung

\footnotesize
\begin{theorem}[\small Betaparameterschätzung im additiven 2 x 2 ZVA-Modell mit Referenzgruppe]

\justifying
\normalfont
Gegeben sei die Designmatrixform eines balancierten additiven 2 x 2 ZVA-Modells mit Referenzgruppe. Dann ergibt sich für den Betaparameterschätzer
\begin{equation}
\hat{\beta} :=
\begin{pmatrix}
\hat{\mu}_0    \\
\hat{\alpha}_2 \\
\hat{\beta}_2  \\
\end{pmatrix} =
\begin{pmatrix}
\frac{3}{4}\bar{y}_{11} + \frac{1}{4}(\bar{y}_{12} + \bar{y}_{21}) - \frac{1}{4}\bar{y}_{22} \\
\frac{1}{2}(\bar{y}_{21} + \bar{y}_{22}) - \frac{1}{2}(\bar{y}_{11} + \bar{y}_{12})          \\
\frac{1}{2}(\bar{y}_{12} + \bar{y}_{22}) - \frac{1}{2}(\bar{y}_{11} + \bar{y}_{21})          \\
\end{pmatrix},
\end{equation}

wobei
\begin{equation}
\bar{y}_{ij} := \frac{1}{n_{ij}}\sum_{k = 1}^{n_{ij}} y_{ijk}
\quad \mbox{für} \quad
1 \le i,j \le 2
\end{equation}

das Stichprobenmittel der $i,j$ten Gruppe des 2 x 2 ZVA-Designs bezeichnet.

\end{theorem}

\small
\underline{Beweis}

Wir bestimmen zunächst $X^\mathrm{T}y, \; X^\mathrm{T}X$ und $(X^\mathrm{T}X)^{-1}$ bei konstantem $n_{ij}$ für $1 \le i,j \le 2$.


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
\begin{align}
\renewcommand{\arraystretch}{1.2}
\begin{split}
X^\mathrm{T}y & =
\setcounter{MaxMatrixCols}{20}
\begin{pmatrix}
1 & \cdots & 1 & 1 & \cdots & 1 & 1 & \cdots & 1 & 1 & \cdots & 1 \\
0 & \cdots & 0 & 0 & \cdots & 0 & 1 & \cdots & 1 & 1 & \cdots & 1 \\
0 & \cdots & 0 & 1 & \cdots & 1 & 0 & \cdots & 0 & 1 & \cdots & 1 \\
\end{pmatrix}
\begin{pmatrix*}[c]
y_{111} \\ \vdots \\ y_{11n_{11}} \\
y_{121} \\ \vdots \\ y_{12n_{12}} \\
y_{211} \\ \vdots \\ y_{21n_{21}} \\
y_{221} \\ \vdots \\ y_{22n_{22}} \\
\end{pmatrix*} \\
& =
\begin{pmatrix*}[l]
\sum_{i=1}^2 \sum_{j=1}^2 \sum_{k = 1}^{n_{ij}} y_{ijk} \\
\sum_{j=1}^2 \sum_{k = 1}^{n_{2j}} y_{2jk}              \\
\sum_{i=1}^2 \sum_{k = 1}^{n_{i2}} y_{i2k}              \\
\end{pmatrix*}
\end{split}
\end{align}


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
\begin{align}
\renewcommand{\arraystretch}{1}
\begin{split}
X^\mathrm{T}X & =
\setcounter{MaxMatrixCols}{20}
\begin{pmatrix}
1 & \cdots & 1 & 1 & \cdots & 1 & 1 & \cdots & 1 & 1 & \cdots & 1 \\
0 & \cdots & 0 & 0 & \cdots & 0 & 1 & \cdots & 1 & 1 & \cdots & 1 \\
0 & \cdots & 0 & 1 & \cdots & 1 & 0 & \cdots & 0 & 1 & \cdots & 1 \\
\end{pmatrix}
\begin{pmatrix}
1       &  0       &  0       \\
\vdots  &  \vdots  &  \vdots  \\
1       &  0       &  0       \\
1       &  0       &  1       \\
\vdots  &  \vdots  &  \vdots  \\
1       &  0       &  1       \\
1       &  1       &  0       \\
\vdots  &  \vdots  &  \vdots  \\
1       &  1       &  0       \\
1       &  1       &  1       \\
\vdots  &  \vdots  &  \vdots  \\
1       &  1       &  1       \\
\end{pmatrix} \\
&=
n_{ij}
\begin{pmatrix*}[r]
4 & 2 & 2 \\
2 & 2 & 1 \\
2 & 1 & 2 \\
\end{pmatrix*}
\end{split}
\end{align}


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
Ohne Beweis halten wir weiterhin fest, dass
\begin{align}
(X^\mathrm{T} X)^{-1} =
n_{ij}
\begin{pmatrix*}[r]
4 & 2 & 2 \\
2 & 2 & 1 \\
2 & 1 & 2 \\
\end{pmatrix*}^{-1}
=
\frac{1}{n_{ij}}
\begin{pmatrix*}[r]
 \frac{3}{4} & -\frac{1}{2} & -\frac{1}{2} \\
-\frac{1}{2} & 1            & 0            \\
-\frac{1}{2} & 0            & 1            \\
\end{pmatrix*}.
\end{align}

Es ergibt sich also
\begin{equation}
\hat{\beta} =
\begin{pmatrix}
\hat{\mu}_0    \\
\hat{\alpha}_2 \\
\hat{\beta}_2  \\
\end{pmatrix}
=
\frac{1}{n_{ij}}
\begin{pmatrix*}[r]
 \frac{3}{4} & -\frac{1}{2} & -\frac{1}{2} \\
-\frac{1}{2} & 1            & 0            \\
-\frac{1}{2} & 0            & 1            \\
\end{pmatrix*}
\begin{pmatrix*}[l]
\sum_{i=1}^2 \sum_{j=1}^2 \sum_{k = 1}^{n_{ij}} y_{ijk} \\
\sum_{j=1}^2 \sum_{k = 1}^{n_{2j}} y_{2jk}              \\
\sum_{i=1}^2 \sum_{k = 1}^{n_{i2}} y_{i2k}              \\
\end{pmatrix*}.
\end{equation}


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
Damit ergibt sich dann
\begin{align}
\begin{split}
\hat{\mu}_0 & =
\frac{1}{n_{ij}}
\left(
  \frac{3}{4}\sum_{i=1}^2 \sum_{j=1}^2 \sum_{k = 1}^{n_{ij}} y_{ijk}
- \frac{1}{2}\sum_{j=1}^2 \sum_{k = 1}^{n_{2j}} y_{2jk}
- \frac{1}{2}\sum_{i=1}^2 \sum_{k = 1}^{n_{i2}} y_{i2k}
\right) \\
& =
\frac{1}{n_{ij}}
\left(
  \frac{3}{4}\sum_{k = 1}^{n_{11}} y_{11k}
+ \frac{3}{4}\sum_{k = 1}^{n_{12}} y_{12k}
+ \frac{3}{4}\sum_{k = 1}^{n_{21}} y_{21k}
+ \frac{3}{4}\sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& +
  \frac{1}{n_{ij}}
\left(
- \frac{1}{2}\sum_{k = 1}^{n_{21}} y_{21k} - \frac{1}{2}\sum_{k = 1}^{n_{22}} y_{22k}
- \frac{1}{2}\sum_{k = 1}^{n_{12}} y_{12k} - \frac{1}{2}\sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& = \frac{1}{n_{ij}}
\left(
  \frac{3}{4}\sum_{k = 1}^{n_{11}} y_{11k}
+ \frac{1}{4}\sum_{k = 1}^{n_{12}} y_{12k}
+ \frac{1}{4}\sum_{k = 1}^{n_{21}} y_{21k}
- \frac{1}{4}\sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& =
  \frac{3}{4}\bar{y}_{11}
+ \frac{1}{4}(\bar{y}_{12} + \bar{y}_{21})
- \frac{1}{4}\bar{y}_{22}
\end{split}
\end{align}


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
sowie
\begin{align}
\begin{split}
\hat{\alpha}_2 & =
\frac{1}{n_{ij}}
\left(
- \frac{1}{2} \sum_{i=1}^2 \sum_{j=1}^2 \sum_{k = 1}^{n_{ij}} y_{ijk} + \sum_{j=1}^2 \sum_{k = 1}^{n_{2j}} y_{2jk}
\right) \\
& =
\frac{1}{n_{ij}}
\left(
- \frac{1}{2}\sum_{k = 1}^{n_{11}} y_{11k}
- \frac{1}{2}\sum_{k = 1}^{n_{12}} y_{12k}
- \frac{1}{2}\sum_{k = 1}^{n_{21}} y_{21k}
- \frac{1}{2}\sum_{k = 1}^{n_{22}} y_{22k}
+            \sum_{k = 1}^{n_{21}} y_{21k}
+            \sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& =
\frac{1}{n_{ij}}
\left(
- \frac{1}{2}\sum_{k = 1}^{n_{11}} y_{11k}
- \frac{1}{2}\sum_{k = 1}^{n_{12}} y_{12k}
+ \frac{1}{2}\sum_{k = 1}^{n_{21}} y_{21k}
+ \frac{1}{2}\sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& = \frac{1}{2}(\bar{y}_{21} + \bar{y}_{22}) - \frac{1}{2}(\bar{y}_{11} + \bar{y}_{12})
\end{split}
\end{align}

und analog für $\hat{\beta}_2$.
$\hfill\Box$


<!-- Beispiel: Betaparameterschätzer -->
# Modellschätzung

\vspace{2mm}
\textcolor{darkblue}{Beispiel}
\vspace{2mm}

\tiny
\setstretch{1.1}
```{r, echo = T}
# Datenreformatierung
fname      = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"
D          = read.table(fname, sep = ",", header = TRUE)     # Datensatz
A1B1       = D$dBDI[D$Setting == "F2F" & D$Variant == "MND"] # face-to-face, mindfulness
A1B2       = D$dBDI[D$Setting == "F2F" & D$Variant == "EXC"] # face-to-face, exercise
A2B1       = D$dBDI[D$Setting == "ONL" & D$Variant == "MND"] # online,       mindfulness
A2B2       = D$dBDI[D$Setting == "ONL" & D$Variant == "EXC"] # online,       exercise

# Datenmatrix für Gruppenmittelwerte
n_ij       = length(A1B1)                                    # Anzahl der Datenpunkte pro Gruppe
Y          = matrix(c(A1B1,A1B2,A2B1,A2B2), nrow = n_ij)     # n_ij x 4 Datenmatrix
y_bar      = colMeans(Y)                                     #    1 x 4 Zellenmittelwerte

# Modellschätzung
I          = 2                                               # Anzahl Level Faktor A (Therapiesetting)
J          = 2                                               # Anzahl Level Faktor B (Therapievariante)
n          = I*J*n_ij                                        # Anzahl Datenpunkte insgesamt
p          = 1 + (I-1)+(J-1)+(I*J-3)                         # Anzahl Parameter insgesamt
D          = matrix(c(1,0,0,                                 # prototypische Designmatrix für balancierte Designs
                      1,0,1,
                      1,1,0,
                      1,1,1),
                    nrow  = p,
                    byrow = TRUE)
C          = matrix(rep(1,n_ij),nrow = n_ij)                 # prototypischer Zellenvektor für balancierte Designs
X          = kronecker(D,C)                                  # Kroneckerprodukt-Designmatrix
y          = matrix(c(A1B1,A1B2,A2B1,A2B2), nrow = n)        # Datenvektor
beta_hat   = solve(t(X) %*% X) %*% t(X) %*% y                # Betaparameterschätzer
eps_hat    = y - X %*% beta_hat                              # Residuenvektor
sigsqr_hat = (t(eps_hat) %*% eps_hat) /(n-p)                 # Varianzparameterschätzer
```


<!-- Beispiel: Betaparameterschätzer -->
# Modellschätzung

Betaparameterschätzer und Stichprobenmittelwerte
\vspace{2mm}

\small
\setstretch{1.5}
```{r, echo = F}
# Ausgabe
cat(  "hat{beta}                                                : ", beta_hat,
    "\nhat{sigsqr}                                              : ", sigsqr_hat,
    "\nbar{y}_11, bar{y}_12, bar{y}_21, bar{y}_22               : ", y_bar,
    "\n3/4bar{y}_11 + 1/4(bar{y}_12 + bar{y}_21) - 1/4bar{y}_22 : ", (3/4)*y_bar[1]+(1/4)*(y_bar[2]+y_bar[3])-(1/4)*y_bar[4],
    "\n1/2(bar{y}_21 + bar{y}_22) - 1/2(bar{y}_11 + bar{y}_12)  : ", (1/2)*(y_bar[3]+y_bar[4])-(1/2)*(y_bar[1] + y_bar[2]),
    "\n1/2(bar{y}_12 + bar{y}_22) - 1/2(bar{y}_11 + bar{y}_21)  : ", (1/2)*(y_bar[2]+y_bar[4])-(1/2)*(y_bar[1] + y_bar[3]))
```
\vfill


<!-- Theorem: Betaparameterschätzung im 2 x 2 ZVA-Modell mit Interaktion und Referenzgruppe -->
# Modellschätzung

\footnotesize
\begin{theorem}[\small Betaparameterschätzung im 2 x 2 ZVA-Modell mit Interaktion und Referenzgruppe]

\justifying
\normalfont
Gegeben sei die Designmatrixform eines balancierten 2 x 2 ZVA-Modells mit Interaktion und Referenzgruppe. Dann ergibt sich für den Betaparameterschätzer
\begin{equation}
\hat{\beta} :=
\begin{pmatrix}
\hat{\mu}_0       \\
\hat{\alpha}_2    \\
\hat{\beta}_2     \\
\hat{\gamma}_{22} \\
\end{pmatrix} =
\begin{pmatrix}
\bar{y}_{11}                                              \\
\bar{y}_{21} - \bar{y}_{11}                               \\
\bar{y}_{12} - \bar{y}_{11}                               \\
\bar{y}_{11} + \bar{y}_{22} - \bar{y}_{12} - \bar{y}_{21} \\
\end{pmatrix},
\end{equation}

wobei
\begin{equation}
\bar{y}_{ij} := \frac{1}{n_{ij}}\sum_{k = 1}^{n_{ij}} y_{ijk}
\quad \mbox{für} \quad
1 \le i,j \le 2
\end{equation}

das Stichprobenmittel der $i,j$ten Gruppe des 2 x 2 ZVA-Designs bezeichnet.

\end{theorem}

\small
\underline{Beweis}

Wir bestimmen zunächst $X^\mathrm{T}y, \; X^\mathrm{T}X$ und $(X^\mathrm{T}X)^{-1}$ bei konstantem $n_{ij}$ für $1 \le i,j \le 2$.


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
\begin{align}
\renewcommand{\arraystretch}{1.2}
\begin{split}
X^\mathrm{T}y & =
\setcounter{MaxMatrixCols}{20}
\begin{pmatrix}
1 & \cdots & 1 & 1 & \cdots & 1 & 1 & \cdots & 1 & 1 & \cdots & 1 \\
0 & \cdots & 0 & 0 & \cdots & 0 & 1 & \cdots & 1 & 1 & \cdots & 1 \\
0 & \cdots & 0 & 1 & \cdots & 1 & 0 & \cdots & 0 & 1 & \cdots & 1 \\
0 & \cdots & 0 & 0 & \cdots & 0 & 0 & \cdots & 0 & 1 & \cdots & 1 \\
\end{pmatrix}
\begin{pmatrix*}[c]
y_{111} \\ \vdots \\ y_{11n_{11}} \\
y_{121} \\ \vdots \\ y_{12n_{12}} \\
y_{211} \\ \vdots \\ y_{21n_{21}} \\
y_{221} \\ \vdots \\ y_{22n_{22}} \\
\end{pmatrix*} \\
& =
\begin{pmatrix*}[l]
\sum_{i=1}^2 \sum_{j=1}^2 \sum_{k = 1}^{n_{ij}} y_{ijk} \\
\sum_{j=1}^2 \sum_{k = 1}^{n_{2j}} y_{2jk}              \\
\sum_{i=1}^2 \sum_{k = 1}^{n_{i2}} y_{i2k}              \\
\sum_{k = 1}^{n_{22}} y_{22k}                           \\
\end{pmatrix*}
\end{split}
\end{align}


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
\begin{align}
\renewcommand{\arraystretch}{1}
\begin{split}
X^\mathrm{T}X & =
\setcounter{MaxMatrixCols}{20}
\begin{pmatrix}
1 & \cdots & 1 & 1 & \cdots & 1 & 1 & \cdots & 1 & 1 & \cdots & 1  \\
0 & \cdots & 0 & 0 & \cdots & 0 & 1 & \cdots & 1 & 1 & \cdots & 1  \\
0 & \cdots & 0 & 1 & \cdots & 1 & 0 & \cdots & 0 & 1 & \cdots & 1  \\
0 & \cdots & 0 & 0 & \cdots & 0 & 0 & \cdots & 0 & 1 & \cdots & 1  \\
\end{pmatrix}
\begin{pmatrix}
1       &  0       &  0       &  0       \\
\vdots  &  \vdots  &  \vdots  &  \vdots  \\
1       &  0       &  0       &  0       \\
1       &  0       &  1       &  0       \\
\vdots  &  \vdots  &  \vdots  &  \vdots  \\
1       &  0       &  1       &  0       \\
1       &  1       &  0       &  0       \\
\vdots  &  \vdots  &  \vdots  &  \vdots  \\
1       &  1       &  0       &  0       \\
1       &  1       &  1       &  1       \\
\vdots  &  \vdots  &  \vdots  &  \vdots  \\
1       &  1       &  1       &  1       \\
\end{pmatrix} \\
& =
n_{ij}
\begin{pmatrix*}[r]
4 & 2 & 2 & 1 \\
2 & 2 & 1 & 1 \\
2 & 1 & 2 & 1 \\
1 & 1 & 1 & 1
\end{pmatrix*}
\end{split}
\end{align}


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
Ohne Beweis halten wir weiterhin fest, dass
\begin{align}
(X^\mathrm{T} X)^{-1} =
n_{ij}
\begin{pmatrix*}[r]
4 & 2 & 2 & 1 \\
2 & 2 & 1 & 1 \\
2 & 1 & 2 & 1 \\
1 & 1 & 1 & 1
\end{pmatrix*}^{-1}
=
\frac{1}{n_{ij}}
\begin{pmatrix*}[r]
 1 & -1 & -1 &  1 \\
-1 &  2 &  1 & -2 \\
-1 &  1 &  2 & -2 \\
 1 & -2 & -2 &  4
\end{pmatrix*}.
\end{align}

Es ergibt sich also
\begin{equation}
\hat{\beta} =
\begin{pmatrix}
\hat{\mu}_0       \\
\hat{\alpha}_2    \\
\hat{\beta}_2     \\
\hat{\gamma}_{22} \\
\end{pmatrix}
=
\frac{1}{n_{ij}}
\begin{pmatrix*}[r]
 1 & -1 & -1 &  1 \\
-1 &  2 &  1 & -2 \\
-1 &  1 &  2 & -2 \\
 1 & -2 & -2 &  4 \\
\end{pmatrix*}
\begin{pmatrix*}[l]
\sum_{i=1}^2 \sum_{j=1}^2 \sum_{k = 1}^{n_{ij}} y_{ijk} \\
\sum_{j=1}^2 \sum_{k = 1}^{n_{2j}} y_{2jk}              \\
\sum_{i=1}^2 \sum_{k = 1}^{n_{i2}} y_{i2k}              \\
\sum_{k = 1}^{n_{22}} y_{22k}                           \\
\end{pmatrix*}.
\end{equation}


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
Damit ergibt sich dann
\begin{align}
\begin{split}
\hat{\mu}_0 & =
\frac{1}{n_{ij}}
\left(
  \sum_{i=1}^2 \sum_{j=1}^2 \sum_{k = 1}^{n_{ij}} y_{ijk}
- \sum_{j=1}^2 \sum_{k = 1}^{n_{2j}} y_{2jk}
- \sum_{i=1}^2 \sum_{k = 1}^{n_{i2}} y_{i2k}
+ \sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& =
\frac{1}{n_{ij}}
\left(
  \sum_{k = 1}^{n_{11}} y_{11k}
+ \sum_{k = 1}^{n_{12}} y_{12k}
+ \sum_{k = 1}^{n_{21}} y_{21k}
+ \sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& +
\frac{1}{n_{ij}}
\left(
- \sum_{k = 1}^{n_{21}} y_{21k} - \sum_{k = 1}^{n_{22}} y_{22k}
- \sum_{k = 1}^{n_{12}} y_{12k} - \sum_{k = 1}^{n_{22}} y_{22k}
+ \sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& = \frac{1}{n_{11}} \sum_{k = 1}^{n_{11}} y_{11k} \\
& = \bar{y}_{11}
\end{split}
\end{align}


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
sowie
\begin{align}
\begin{split}
\hat{\alpha}_2 & =
\frac{1}{n_{ij}}
\left(
-  \sum_{i=1}^2 \sum_{j=1}^2 \sum_{k = 1}^{n_{ij}} y_{ijk}
+ 2\sum_{j=1}^2 \sum_{k = 1}^{n_{2j}} y_{2jk}
+ 1\sum_{i=1}^2 \sum_{k = 1}^{n_{i2}} y_{i2k}
- 2 \sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& =
\frac{1}{n_{ij}}
\left(
- \sum_{k = 1}^{n_{11}} y_{11k}
- \sum_{k = 1}^{n_{12}} y_{12k}
- \sum_{k = 1}^{n_{21}} y_{21k}
- \sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& +
\frac{1}{n_{ij}}
\left(
  2\sum_{k = 1}^{n_{21}} y_{21k} +2\sum_{k = 1}^{n_{22}} y_{22k}
+  \sum_{k = 1}^{n_{12}} y_{12k} + \sum_{k = 1}^{n_{22}} y_{22k}
- 2\sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& = \frac{1}{n_{ij}}\left(\sum_{k = 1}^{n_{21}} y_{21k} - \sum_{k = 1}^{n_{11}} y_{11k}\right) \\
& = \bar{y}_{21} - \bar{y}_{11}
\end{split}
\end{align}

und analog für $\hat{\beta}_2$.


<!-- Beweis -->
# Modellschätzung

\small
\underline{Beweis (fortgeführt)}

\footnotesize
Schließlich ergibt sich
\begin{align}
\begin{split}
\hat{\gamma}_{22} & =
\frac{1}{n_{ij}}
\left(
   \sum_{i=1}^2 \sum_{j=1}^2 \sum_{k = 1}^{n_{ij}} y_{ijk}
- 2\sum_{j=1}^2 \sum_{k = 1}^{n_{2j}} y_{2jk}
- 2\sum_{i=1}^2 \sum_{k = 1}^{n_{i2}} y_{i2k}
+ 4 \sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& =
\frac{1}{n_{ij}}
\left(
  \sum_{k = 1}^{n_{11}} y_{11k}
+ \sum_{k = 1}^{n_{12}} y_{12k}
+ \sum_{k = 1}^{n_{21}} y_{21k}
+ \sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& +
\frac{1}{n_{ij}}
\left(
- 2\sum_{k = 1}^{n_{21}} y_{21k} - 2\sum_{k = 1}^{n_{22}} y_{22k}
- 2\sum_{k = 1}^{n_{12}} y_{12k} - 2\sum_{k = 1}^{n_{22}} y_{22k}
+ 4\sum_{k = 1}^{n_{22}} y_{22k}
\right) \\
& =
\frac{1}{n_{ij}}
\left(
  \sum_{k = 1}^{n_{11}} y_{11k}
+ \sum_{k = 1}^{n_{22}} y_{22k}
- \sum_{k = 1}^{n_{12}} y_{12k}
- \sum_{k = 1}^{n_{21}} y_{21k}
\right) \\
& = \bar{y}_{11} + \bar{y}_{22} - \bar{y}_{12} - \bar{y}_{21}.
\end{split}
\end{align}
$\hfill\Box$


<!-- Beispiel: Betaparameterschätzer -->
# Modellschätzung

\vspace{2mm}
\textcolor{darkblue}{Beispiel}
\vspace{2mm}

\tiny
\setstretch{1.1}
```{r, echo = T}
# Datenreformatierung
fname      = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"
D          = read.table(fname, sep = ",", header = TRUE)     # Datensatz
A1B1       = D$dBDI[D$Setting == "F2F" & D$Variant == "MND"] # face-to-face, mindfulness
A1B2       = D$dBDI[D$Setting == "F2F" & D$Variant == "EXC"] # face-to-face, exercise
A2B1       = D$dBDI[D$Setting == "ONL" & D$Variant == "MND"] # online,       mindfulness
A2B2       = D$dBDI[D$Setting == "ONL" & D$Variant == "EXC"] # online,       exercise

# Datenmatrix für Gruppenmittelwerte
n_ij       = length(A1B1)                                    # Anzahl der Datenpunkte pro Gruppe
Y          = matrix(c(A1B1,A1B2,A2B1,A2B2), nrow = n_ij)     # n_ij x 4 Datenmatrix
y_bar      = colMeans(Y)                                     #    1 x 4 Zellenmittelwerte

# Modellschätzung
I          = 2                                               # Anzahl Level Faktor A (Therapiesetting)
J          = 2                                               # Anzahl Level Faktor B (Therapievariante)
n          = I*J*n_ij                                        # Anzahl Datenpunkte insgesamt
p          = 1 + (I-1)+(J-1)+(I*J-3)                         # Anzahl Parameter insgesamt
D          = matrix(c(1,0,0,0,                               # prototypische Designmatrix für balancierte Designs
                      1,0,1,0,
                      1,1,0,0,
                      1,1,1,1),
                    nrow  = p,
                    byrow = TRUE)
C          = matrix(rep(1,n_ij), nrow = n_ij)                # prototypischer Zellenvektor für balancierte Designs
X          = kronecker(D,C)                                  # Kroneckerprodukt-Designmatrix
y          = matrix(c(A1B1,A1B2,A2B1,A2B2), nrow = n)        # Datenvektor
beta_hat   = solve(t(X) %*% X) %*% t(X) %*% y                # Betaparameterschätzer
eps_hat    = y - X %*% beta_hat                              # Residuenvektor
sigsqr_hat = (t(eps_hat) %*% eps_hat) /(n-p)                 # Varianzparameterschätzer
```


<!-- Beispiel: Betaparameterschätzer -->
# Modellschätzung

Betaparameterschätzer und Stichprobenmittelwerte
\vspace{2mm}

\small
\setstretch{1.5}
```{r, echo = F}
# Ausgabe
cat(  "hat{beta}                                     : ", beta_hat,
    "\nhat{sigsqr}                                   : ", sigsqr_hat,
    "\nbar{y}_11, bar{y}_12, bar{y}_21, bar{y}_22    : ", y_bar,
    "\nbar{y}_11                                     : ", y_bar[1],
    "\nbar{y}_21 - bar{y}_11                         : ", y_bar[3]-y_bar[1],
    "\nbar{y}_12 - bar{y}_11                         : ", y_bar[2]-y_bar[1],
    "\nbar{y}_11 + bar{y}_22 - bar{y}_12 + bar{y}_21 : ", y_bar[1]+y_bar[4]-y_bar[3]-y_bar[2])
```
\vfill


<!-- Abschnitt 4: Modellevaluation -->
#

\large
\setstretch{2.5}
\vfill

Anwendungsszenario

Modellformulierung

Modellschätzung

**Modellevaluation**

Anwendung/Praxis

Selbstkontrollfragen

\vfill


<!-- Überblick: Modellevaluation -->
# Modellevaluation

\textcolor{darkblue}{Überblick}

\small
Wie bei der EVA kann auch bei der ZVA eine Modellevaluationstheorie mithilfe einer Quadratsummenzerlegung entwickelt werden (siehe Einheit (10) in *Allgemeines Lineares Modell*). Mit zunehmender Designkomplexität wird eine solche Theorie allerdings zunehmend unübersichtlich. Prinzipiell können alle Quadratsummenzerlegung-basierten F-Statistiken auf Likelihood-basierte Modellvergleiche zurückgeführt werden. Allerdings fehlt dazu in diesem Kurs bisher eine allgemeine Kontrasttheorie. Wir beschränken und im Folgenden deshalb auf die

(1) Evaluation der Haupteffekte im additiven Modell der 2 x 2 ZVA mit Referenzgruppe und die
(2) Evaluation der Interaktion im Modell der 2 x 2 ZVA mit Interaktion und Referenzgruppe.
\vfill


<!-- Inhaltsverzeichnis: Modellevaluation -->
# Modellevaluation

\setstretch{2}

\textcolor{darkblue}{Evaluation der Haupteffekte im additiven Modell der 2 x 2 ZVA mit Referenzgruppe}

(1) Statistische Modelle und Teststatistiken

(2) Testhypothesen und Tests

(3) Testumfangkontrollen und p-Werte

\textcolor{darkblue}{Evaluation der Interaktion im Modell der 2 x 2 ZVA mit Interaktion und Referenzgruppe}

(1) Statistisches Modell und Teststatistik

(2) Testhypothese und Test

(3) Testumfangkontrolle und p-Wert


<!-- Theorem: Teststatistiken für Haupteffekte -->
# Modellevaluation: Haupteffekte (1) Teststatistiken

\vspace{2mm}
\footnotesize
\begin{theorem}[Teststatistiken für Haupteffekte]

\justifying
\normalfont
Gegeben sei die Designmatrixform des additiven Modells der 2 x 2 ZVA mit Referenzgruppe
\begin{equation}
y = X\beta + \varepsilon
\quad \mbox{mit} \quad
\varepsilon \sim N(0_n,\sigma^2I_n) \; ,
\end{equation}

wobei die Spalten von $X$ bezeichnet seien durch
\begin{equation}
X := \begin{pmatrix} X_{\mu_0} & X_{\alpha_2} & X_{\beta_2} \end{pmatrix} \in \mathbb{R}^{n \times 3} \; .
\end{equation}

Dann gelten:

\begin{itemize}
\item[(A)] Eine \textit{F-Teststatistik für den Haupteffekt von Faktor A}, $F_A$, ist die F-Statistik unter der Partitionierung
\begin{equation}
X       := \begin{pmatrix} X_{\mu_0} & X_{\beta_2} & X_{\alpha_2} \end{pmatrix}, \;
\beta   := \begin{pmatrix} \mu_0 \\ \beta_2 \\ \alpha_2 \end{pmatrix} \quad \mbox{und} \quad
X_0     := \begin{pmatrix} X_{\mu_0} & X_{\beta_2} \end{pmatrix}, \;
\beta_0 := \begin{pmatrix} \mu_0 \\ \beta_2 \end{pmatrix} \; .
\end{equation}

\item[(B)] Eine \textit{F-Teststatistik für den Haupteffekt von Faktor B}, $F_B$, ist die F-Statistik unter der Partitionierung
\begin{equation}
X       := \begin{pmatrix} X_{\mu_0} & X_{\alpha_2} & X_{\beta_2} \end{pmatrix}, \;
\beta   := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \beta_2  \end{pmatrix} \quad \mbox{und} \quad
X_0     := \begin{pmatrix} X_{\mu_0} & X_{\alpha_2} \end{pmatrix}, \;
\beta_0 := \begin{pmatrix} \mu_0 \\ \alpha_2 \end{pmatrix} \; .
\end{equation}
\end{itemize}

\end{theorem}

Bemerkung
\vspace{-2mm}

* Der Beweis ergibt sich aus dem Theorem zur Verteilung der F-Statistik bei Partitionierung eines ALM (siehe Einheit (8) in *Allgemeines Lineares Modell*).


<!-- Definition: Testhypothesen und Tests -->
# Modellevaluation: Haupteffekte (2) Testhypothesen

\footnotesize
\begin{definition}[Testhypothesen und Tests]

\justifying
Gegeben sei das Modell der additiven 2 x 2 ZVA mit Referenzgruppe. Die F-Teststatistiken für die Haupteffekte von Faktor A und B seien mit $F_A$ und $F_B$ bezeichnet und wie oben definiert. Dann gilt:

\begin{itemize}
\item[(A)] Der kritische-Wert-basierte Test
\begin{equation}
\phi_A(y) := 1_{\{F_A \ge k\}}
\quad \mbox{mit Nullhypothese} \quad
H_0^A : \alpha_2 = 0
\end{equation}
definiert den \textit{F-Test des Haupteffekts von Faktor A}.

\item[(B)] Der kritische-Wert-basierte Test
\begin{equation}
\phi_B(y) := 1_{\{F_B \ge k\}}
\quad \mbox{mit Nullhypothese} \quad
H_0^B : \beta_2 = 0
\end{equation}
definiert den \textit{F-Test des Haupteffekts von Faktor B}.
\end{itemize}

\end{definition}

Bemerkung

* Die F-Tests beziehen sich auf die im vorangegangenen Theorem definierten F-Statistiken.


<!-- Theorem: Testumfangkontrolle und p-Werte -->
# Modellevaluation: Haupteffekte (3) Testumfangkontrolle

\footnotesize
\begin{theorem}[Testumfangkontrolle und p-Werte]

\justifying
\normalfont
Mit obigen Definitionen und der KVF $\varphi(\cdot; n_1, n_2)$ der $f$-Verteilung gelten:

\begin{itemize}
\item[(A)] $\phi_A$ ist ein Level-$\alpha_0$-Test mit Testumfang $\alpha_0$, wenn der kritische Wert definiert ist durch
\begin{equation}
k_{\alpha_0}^A := \varphi^{-1}(1-\alpha_0; 1,n-3).
\end{equation}

Der zu einem beobachteten Wert $f_A$ von $F_A$ assoziierte p-Wert ist gegeben durch
\begin{equation}
\mbox{p-Wert} := 1 - \varphi(f_A; 1,n-3).
\end{equation}

\item[(B)] $\phi_B$ ist ein Level-$\alpha_0$-Test mit Testumfang $\alpha_0$, wenn der kritische Wert definiert ist durch
\begin{equation}
k_{\alpha_0}^B := \varphi^{-1}(1-\alpha_0; 1,n-3).
\end{equation}

Der zu einem beobachteten Wert $f_B$ von $F_B$ assoziierte p-Wert ist gegeben durch
\begin{equation}
\mbox{p-Wert} := 1 - \varphi(f_B; 1,n-3).
\end{equation}
\end{itemize}

\end{theorem}

Bemerkung

* Wir verzichten auf einen Beweis.


<!-- Inhaltsverzeichnis: Modellevaluation -->
# Modellevaluation

\setstretch{2}

\textcolor{darkblue}{Evaluation der Haupteffekte im additiven Modell der 2 x 2 ZVA mit Referenzgruppe}

(1) Statistische Modelle und Teststatistiken

(2) Testhypothesen und Tests

(3) Testumfangkontrollen und p-Werte

\textcolor{darkblue}{Evaluation der Interaktion im Modell der 2 x 2 ZVA mit Interaktion und Referenzgruppe}

(1) Statistisches Modell und Teststatistik

(2) Testhypothese und Test

(3) Testumfangkontrolle und p-Wert


<!-- Theorem: Teststatistik für die Interaktion -->
# Modellevaluation: Interaktion (1) Teststatistik

\vspace{2mm}
\footnotesize
\begin{theorem}[Teststatistik für die Interaktion]

\justifying
\normalfont
Gegeben sei die Designmatrixform des additiven Modells der 2 x 2 ZVA mit Interaktion Referenzgruppe
\begin{equation}
y = X\beta + \varepsilon
\quad \mbox{mit} \quad
\varepsilon \sim N(0_n,\sigma^2I_n) \; ,
\end{equation}

wobei die Spalten von $X$ bezeichnet seien durch
\begin{equation}
X := \begin{pmatrix} X_{\mu_0} & X_{\alpha_2} & X_{\beta_2} & X_{\gamma_{22}} \end{pmatrix} \in \mathbb{R}^{n \times 4} \; .
\end{equation}

Dann gilt: Eine \textit{F-Teststatistik für die Interaktion von Faktor A und B} wird als $F_{A \times B}$ bezeichnet und ist die F-Statistik unter der Partitionierung
\begin{align}
\begin{split}
X       &:= \begin{pmatrix} X_{\mu_0} & X_{\alpha_2} & X_{\beta_2} & X_{\gamma_{22}} \end{pmatrix}, \;
\beta    := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \beta_2 \\ \gamma_{22} \end{pmatrix}, \;
p        := 4 \quad \mbox{und} \quad \\
X_0     &:= \begin{pmatrix} X_{\mu_0} & X_{\alpha_2} & X_{\beta_2} \end{pmatrix}, \;
\beta_0  := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \beta_2 \end{pmatrix}, \;
p_0      := 3 \; .
\end{split}
\end{align}

\end{theorem}

Bemerkung
\vspace{-2mm}

* Der Beweis ergibt sich aus dem Theorem zur Verteilung der F-Statistik bei Partitionierung eines ALM (siehe Einheit (8) in *Allgemeines Lineares Modell*).


<!-- Definition: Testhypothese und Test -->
<!-- Theorem: Testumfangkontrolle und p-Wert -->
# Modellevaluation: Interaktion (2) Testhypothese (3) Testumfangskontrolle

\footnotesize
\begin{definition}[Testhypothese und Test]

\justifying
Die F-Teststatistik für die Interaktion von Faktor A und Faktor B sei mit $F_{A \times B}$ bezeichnet und wie oben definiert. Dann definiert der kritische-Wert-basierte Test
\begin{equation}
\phi_{A\times B}(y) := 1_{\{F_{A \times B} \ge k\}}
\quad \mbox{mit Nullhypothese} \quad
H_0^{A\times B} : \gamma_{22} = 0
\end{equation}
den \textit{F-Test der Interaktion von Faktor A und Faktor B}.

\end{definition}

\vspace{4mm}
\begin{theorem}[Testumfangkontrolle und p-Wert]

\justifying
\normalfont
Mit der obigen Definition und der KVF $\varphi(\cdot; n_1, n_2)$ der $f$-Verteilung gilt, dass $\phi_{A\times B}$ ein Level-$\alpha_0$-Test mit Testumfang $\alpha_0$ ist, wenn der kritische Wert definiert ist durch
\begin{equation}
k_{\alpha_0}^{A \times B} := \varphi^{-1}(1-\alpha_0; 1,n-4).
\end{equation}

Der zu einem beobachteten Wert $f_{A \times B}$ von $F_{A \times B}$ assoziierte 
p-Wert ist gegeben durch
\begin{equation}
\mbox{p-Wert} := 1 - \varphi(f_{A \times B}; 1,n-4).
\end{equation}

\end{theorem}


<!-- Abschnitt 5: Anwendung/Praxis -->
#

\large
\setstretch{2.5}
\vfill

Anwendungsszenario

Modellformulierung

Modellschätzung

Modellevaluation

**Anwendung/Praxis**

Selbstkontrollfragen

\vfill


<!-- Beispieldatensatz -->
# Anwendung/Praxis

\small
\vspace{2mm}
\textcolor{darkblue}{Daten einlesen:} $j = 1,...,4$ für jede Gruppe

\footnotesize
\vspace{2mm}
```{r, echo = T}
fname       = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"
D           = read.table(fname, sep = ",", header = TRUE)
```

\center
```{r, echo = F}
library(knitr)
n_ij  = 12
kable(D[c(1:4,                                                                    
         (1:4)+sum(n_ij*1),                                                       
         (1:4)+sum(n_ij*2),                                                       
         (1:4)+sum(n_ij*3)),],
      align  = "c",
      "pipe")
```


<!-- Visualisierung: Box-Plot -->
# Anwendung/Praxis

\vspace{2mm}
\textcolor{darkblue}{Visualisierung der Daten: Box-Plot}
\vspace{2mm}

```{r, echo = F, eval = F}
# Datenselektion
fname    = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"                 # Dateiname
D        = read.table(fname, sep = ",", header = TRUE)                      # Daten einlesen
Y        = data.frame(                                                      # Dataframe anlegen
               F2FMND = D$dBDI[D$Setting == "F2F" & D$Variant == "MND"],    # A1B1
               F2FEXC = D$dBDI[D$Setting == "F2F" & D$Variant == "EXC"],    # A1B2
               ONLMND = D$dBDI[D$Setting == "ONL" & D$Variant == "MND"],    # A2B1
               ONLEXC = D$dBDI[D$Setting == "ONL" & D$Variant == "EXC"])    # A2B2

# Abbildungsparameter
library(latex2exp)
dev.new()
par(
    family      = "sans",
    pty         = "m",
    bty         = "l",
    lwd         = 1,
    las         = 1,
    font.main   = 1,
    cex         = 1.1,
    cex.main    = 1.2)

# Boxplot und Speichern
boxplot(Y,
    ylim        = c(0,20))
dev.copy2pdf(
    file        = "Abbildungen/zva_boxplot.pdf",
    width       = 6,
    height      = 4)
```

\center
```{r, echo = F, out.width = "100%"}
knitr::include_graphics("Abbildungen/zva_boxplot.pdf")
```


<!-- Visualisierung: Balkendiagramm -->
# Anwendung/Praxis

\vspace{2mm}
\textcolor{darkblue}{Visualisierung der Daten: Balkendiagramm}
\vspace{2mm}

```{r, echo = F, eval = F}
# Abbildungsparameter
library(latex2exp)                              # LaTeX-Annotations
dev.new()                                       # Figure-Initialisierung
par(                                            # für Details siehe ?par
    family      = "sans",                       # Serif-freier Fonttyp
    pty         = "m",                          # maximale Abbildungsregion
    bty         = "l",                          # L-förmige Box
    lwd         = 1,                            # Liniendicke
    las         = 1,                            # horizontale Achsenbeschriftung
    font.main   = 1,                            # Titel nicht fett
    cex         = 1.1,                          # Textvergrößerungsfaktor
    cex.main    = 1.2)                          # Titeltextvergrößerungsfaktor

# Balkendiagramm
groupmeans  = colMeans(Y)                       # Gruppenmittelwerte
groupstds   = apply(Y,2,sd)                     # Gruppenstandardabweichungen
x = barplot(groupmeans,
    ylim        = c(0,20),                      # Y-Achsen-Limits
    col         = "gray90")                     # Balkenfarbe
arrows(
    x0          = x,                            # arrow start x-ordinate
    y0          = groupmeans - groupstds,       # arrow start y-ordinate
    x1          = x,                            # arrow end   x-ordinate
    y1          = groupmeans + groupstds,       # arrow end   y-ordinate
    code        = 3,                            # Pfeilspitzen beiderseits
    angle       = 90,                           # Pfeilspitzenwinkel -> Linie
    length      = 0.05)                         # Linienlänge

# PDF-Speicherung
dev.copy2pdf(
    file        = "Abbildungen/zva_barplot.pdf",
    width       = 6,
    height      = 4)
```

\center
```{r, echo = F, out.width = "100%"}
knitr::include_graphics("Abbildungen/zva_barplot.pdf")
```


<!-- Visualisierung: Liniendiagramm -->
# Anwendung/Praxis

\vspace{2mm}
\textcolor{darkblue}{Visualisierung der Daten: Liniendiagramm}
\vspace{-2mm}

```{r, echo = F, eval = F}
# Vorbereitung
groupmeans  = matrix(colMeans(Y), nrow = 2)
groupstds   = matrix(apply(Y, 2, sd), nrow = 2)
cols        = c("gray30", "gray70")
lwds        = c(2,2)
print(groupmeans)
print(groupstds)

# Abbildungsparameter
library(latex2exp)                                  # LaTeX-Annotations
dev.new()                                           # Figure-Initialisierung
par(
    family      = "sans",
    pty         = "s",
    bty         = "l",
    lwd         = 1,
    las         = 1,
    mgp         = c(2,1,0),
    xaxs        = "i",
    yaxs        = "i",
    font.main   = 1,
    cex         = 1.5,
    cex.main    = 1)

# Lineplots
x               = 1:2
matplot(x, groupmeans,
    type        = "b",
    pch         = 21,
    xlim        = c(0.5, 2.5),
    ylim        = c(  5, 20),
    lty         = 1,
    col         = cols,
    lwd         = lwds,
    xlab        = "CBT Variant",
    ylab        = TeX("$\\Delta BDI$"),
    xaxt        = "n")
for(i in 1:2){
    arrows(
        x0      = x,                                # arrow start x-ordinate
        y0      = groupmeans[,i] - groupstds[,i],   # arrow start y-ordinate
        x1      = x,                                # arrow end   x-ordinate
        y1      = groupmeans[,i] + groupstds[,i],   # arrow end   y-ordinate
        col     = cols[[i]],                        # Farbe
        lwd     = lwds[[i]],                        # Linienstärke
        code    = 3,                                # Pfeilspitzen beiderseits
        angle   = 90,                               # Pfeilspitzenwinkel -> Linie
        length  = 0.05)                             # Linienlänge
}
legend(0.6, 20, c("CBT Setting F2F", "CBT Setting ONL"),
    lty         = 1,
    col         = c("gray30", "gray90"),
    lwd         = lwds,
    bty         = "n",
    cex         = 1,
    x.intersp   = 0.3,
    y.intersp   = 2,
    seg.len     = 1.5)
text(1, 4, "MND", xpd = T, cex = 1)
text(2, 4, "EXC", xpd = T, cex = 1)

# PDF-Speicherung
dev.copy2pdf(
    file        = "Abbildungen/zva_lineplot.pdf",
    width       = 7.5,
    height      = 7.5)
```

\center
```{r, echo = F, out.width = "80%"}
knitr::include_graphics("Abbildungen/zva_lineplot.pdf")
```


<!-- Anwendungsbeispiel: vollständig in R -->
# Anwendung/Praxis

\vspace{4mm}

```{r, echo = F}
# Datenreformatierung
fname      = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"
D          = read.table(fname, sep = ",", header = TRUE)     # Datensatz
A1B1       = D$dBDI[D$Setting == "F2F" & D$Variant == "MND"] # face-to-face, mindfulness
A1B2       = D$dBDI[D$Setting == "F2F" & D$Variant == "EXC"] # face-to-face, exercise
A2B1       = D$dBDI[D$Setting == "ONL" & D$Variant == "MND"] # online,       mindfulness
A2B2       = D$dBDI[D$Setting == "ONL" & D$Variant == "EXC"] # online,       exercise
```

\tiny
\setstretch{0.5}
```{r, echo = T}
# Modellevaluation
I          = 2                                           # Anzahl Level Faktor A (Therapiesetting)
J          = 2                                           # Anzahl Level Faktor B (Therapievariante)
n_ij       = length(A1B1)                                # balanciertes ZVA-Design
n          = I*J*n_ij                                    # Anzahl Datenpunkte
p          = 3                                           # Anzahl Parameter vollständiges Modell
y          = matrix(c(A1B1,A1B2,A2B1,A2B2), nrow = n)    # Datenvektor
D          = matrix(c(1,0,0,  1,0,1,  1,1,0,  1,1,1),
                    nrow = I*J, byrow = TRUE)            # prototypische Designmatrix
C          = matrix(rep(1,n_ij), nrow = n_ij)            # prototypischer Zellenvektor für balancierte Designs
X          = kronecker(D,C)                              # Kroneckerprodukt-Designmatrix
XH         = list(X[,c(1,3,2)], X)                       # Modellvarianten
alpha_0    = 0.05                                        # Signifikanzlevel
Eff        = rep(NaN,2)                                  # F-Teststatistik Arrayinitialisierung
k_alpha_0  = rep(NaN,2)                                  # kritischer Wert Arrayinitialisierung
phi        = rep(NaN,2)                                  # Testwert Arrayinitialisierung
p_vals     = rep(NaN,2)                                  # p-Wert Arrayinitialisierung
for(i in 1:2){                                           # Iteration über Modellvarianten
    X            = XH[[i]]                               # Designmatrix vollständiges Modell
    X_0          = X[,-3]                                # Designmatrix reduziertes Modell
    p            = ncol(X)                               # Anzahl Parameter vollständiges Modell
    p_0          = ncol(X_0)                             # Anzahl Parameter reduziertes Modell
    p_1          = p - p_0                               # Anzahl zusätzlicher Parameter im vollständigen Modell
    beta_hat_0   = solve(t(X_0)%*%X_0) %*% t(X_0)%*%y    # Betaparameterschätzer reduziertes Modell
    beta_hat     = solve(t(X)  %*%X )  %*% t(X)  %*%y    # Betaparameterschätzer vollständiges Modell
    eps_hat_0    = y - X_0 %*% beta_hat_0                # Residuenvektor reduziertes Modell
    eps_hat      = y - X   %*% beta_hat                  # Residuenvektor vollständiges Modell
    eh0_eh0      = t(eps_hat_0) %*% eps_hat_0            # residuelle QS reduziertes Modell
    eh_eh        = t(eps_hat)   %*% eps_hat              # esiduelle QS vollständiges Modell
    sigsqr_hat   = eh_eh/(n-p)                           # Varianzparameterschätzer vollständiges Modell
    Eff[i]       = ((eh0_eh0-eh_eh)/p_1)/sigsqr_hat      # F-Statistik
    k_alpha_0[i] = qf(1-alpha_0, p_1, n-p)               # kritischer Wert
    if(Eff[i] >= k_alpha_0[i]){ phi[i] = 1 }             # H_0 ablehnen
    else {                      phi[i] = 0 }             # H_0 nicht ablehnen
    p_vals[i]    = 1 - pf(Eff[i], p_1,n-p)               # p-Wert
}
data.frame("f"=Eff, "k"= k_alpha_0, "phi"=phi, "p-Wert"= p_vals, row.names = c("Setting", "Variant"))
```


<!-- Anwendungsbeispiel: mit R-Funktionen -->
# Anwendung/Praxis

\setstretch{0.9}

\small
\textcolor{darkblue}{Haupteffekte der additiven 2 x 2 ZVA mit R's aov() Funktion}

\vspace{1mm}
\tiny
```{r, echo = T}
# Daten einlesen
fname      = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"
D          = read.table(fname, sep = ",", header = TRUE)

# R's aov-Funktion
res.aov    = aov(dBDI ~ Setting + Variant, data = D)          # Modellformulierung und Modellschätzung
summary(res.aov)                                              # Modellevaluation
```

\vspace{2mm}

\small
\textcolor{darkblue}{Haupteffekte der additiven 2 x 2 ZVA mit R's lm() und anova() Funktionen}

\vspace{1mm}
\tiny
```{r, echo = T}
# Dateneinlesen
fname      = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"
D          = read.table(fname, sep = ",", header = TRUE)

# R's lm() und anova() Funktionen
glm        = lm(dBDI ~ Setting + Variant, data = D)           # Modellformulierung und Modellschätzung
anova(glm)                                                    # Modellevaluation
```


<!-- Anwendungsbeispiel: vollständig in R -->
# Anwendung/Praxis

\vspace{4mm}

```{r, echo = F}
# Datenreformatierung
fname      = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"
D          = read.table(fname, sep = ",", header = TRUE)     # Datensatz
A1B1       = D$dBDI[D$Setting == "F2F" & D$Variant == "MND"] # face-to-face, mindfulness
A1B2       = D$dBDI[D$Setting == "F2F" & D$Variant == "EXC"] # face-to-face, exercise
A2B1       = D$dBDI[D$Setting == "ONL" & D$Variant == "MND"] # online,       mindfulness
A2B2       = D$dBDI[D$Setting == "ONL" & D$Variant == "EXC"] # online,       exercise
```

\tiny
\setstretch{1.1}
```{r, echo = T}
# Modellevaluation
I          = 2                                           # Anzahl Level Faktor A (Therapiesetting)
J          = 2                                           # Anzahl Level Faktor B (Therapievariante)
n_ij       = length(A1B1)                                # balanciertes ZVA-Design
n          = I*J*n_ij                                    # Anzahl Datenpunkte
p          = 4                                           # Anzahl Parameter vollständiges Modell
y          = matrix(c(A1B1,A1B2,A2B1,A2B2), nrow = n)    # Datenvektor
D          = matrix(c(1,0,0,0,  1,0,1,0,  1,1,0,0,  1,1,1,1),
                    nrow = I*J, byrow = TRUE)            # prototypische Designmatrix
C          = matrix(rep(1,n_ij), nrow = n_ij)            # prototypischer Zellenvektor für balancierte Designs
X          = kronecker(D,C)                              # Kroneckerprodukt-Designmatrix
alpha_0    = 0.05                                        # Signifikanzlevel
# X        = X                                           # Designmatrix vollständiges Modell
X_0        = X[,-4]                                      # Designmatrix reduziertes Modell
p          = ncol(X)                                     # Anzahl Parameter vollständiges Modell
p_0        = ncol(X_0)                                   # Anzahl Parameter reduziertes Modell
p_1        = p - p_0                                     # Anzahl zusätzlicher Parameter im vollständigen Modell
beta_hat_0 = solve(t(X_0)%*%X_0) %*% t(X_0)%*%y          # Betaparameterschätzer reduziertes Modell
beta_hat   = solve(t(X)  %*%X )  %*% t(X)  %*%y          # Betaparameterschätzer vollständiges Modell
eps_hat_0  = y - X_0 %*% beta_hat_0                      # Residuenvektor reduziertes Modell
eps_hat    = y - X   %*% beta_hat                        # Residuenvektor vollständiges Modell
eh0_eh0    = t(eps_hat_0) %*% eps_hat_0                  # residuelle QS reduziertes Modell
eh_eh      = t(eps_hat)   %*% eps_hat                    # residuelle QS vollständiges Modell
sigsqr_hat = eh_eh/(n-p)                                 # Varianzparameterschätzer vollständiges Modell
f          = ((eh0_eh0-eh_eh)/p_1)/sigsqr_hat            # F-Statistik
k_alpha_0  = qf(1-alpha_0, p_1, n-p)                     # kritischer Wert
if(f >= k_alpha_0){phi = 1} else {phi = 0}               # Test
p_val      = 1 - pf(f, p_1,n-p)                          # p-Wert
data.frame("f"=f, "k"=k_alpha_0, "phi"=phi, "p-Wert"=p_val, row.names=c("Setting x Variant"))
```


<!-- Anwendungsbeispiel: mit R-Funktionen -->
# Anwendung/Praxis

\setstretch{0.8}

\small
\textcolor{darkblue}{Interaktion im Modell der 2 x 2 ZVA mit Interaktion mit R's aov() Funktion}

\vspace{1mm}
\tiny
```{r, echo = T}
# Dateneinlesen
fname      = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"
D          = read.table(fname, sep = ",", header = TRUE)

# R's aov Funktion
res.aov    = aov(dBDI ~ Setting + Variant + Setting:Variant, data = D)    # Modellformulierung und Modellschätzung
summary(res.aov)                                                          # Modellevaluation
```

\vspace{2mm}

\small
\textcolor{darkblue}{Interaktion im Modell der 2 x 2 ZVA mit Interaktion mit R's lm() und anova() Funktionen}

\vspace{1mm}
\tiny
```{r, echo = T}
# Dateneinlesen
fname      = "Daten/Zweifaktorielle_Varianzanalyse_Daten.csv"
D          = read.table(fname, sep = ",", header = TRUE)

# R's lm() und anova() Funktionen
glm        = lm(dBDI ~ Setting + Variant + Setting:Variant, data = D)     # Modellformulierung und Modellschätzung
anova(glm)                                                                # Modellevaluation
```


<!-- Datensatz: Küchenrollen -->
# Anwendung/Praxis

\small
Beispiel: Induzieren Bewegungen im Uhrzeigersinn Offenheit für neue Erfahrungen (= "psychische Zustände zeitlichen Fortschreitens und eine Orientierung auf die Zukunft und Neuartigkeit")? (@wagenmakers_turning_2015)

\vspace{2mm}
\begin{minipage}{0.58\textwidth}
  \includegraphics[width=\linewidth]{Abbildungen/kitchen_rolls_abstract.png}
\end{minipage}
\hfill
\begin{minipage}{0.40\textwidth}
  \includegraphics[width=\linewidth]{Abbildungen/kitchen_rolls_experiment.png}
\end{minipage}


<!-- Datensatz: Küchenrollen -->
# Anwendung/Praxis

\small
Datensatz: mittlere NEO-PI-R-Scores, getrennt nach Rotationsrichtung ($n = 102$)

```{r, echo = F, eval = T}
# Daten einlesen
filename = "Daten/kitchen_rolls.csv"
D        = read.csv(filename)

# Daten aufräumen
levs       = c("clock", "counter")
D$Rotation = factor(D$Rotation, levels = levs)
```

```{r, echo = F, eval = F}
# Daten visualisieren
library(vioplot)
par(
    family      = "sans",
    mar         = c(3,3,1,1),
    pty         = "m",
    bty         = "l",
    lwd         = 1,
    las         = 1,
    mgp         = c(2,1,0),
    xaxs        = "i",
    yaxs        = "i",
    font.main   = 1.5,
    cex         = 1.5)

# Violinplot
vioplot(D$mean_NEO ~ D$Rotation, D,
    col         = "gray80",
    rectCol     = "black",
    lineCol     = "white",
    colMed      = "black",
    border      = "black",
    pchMed      = 16,
    plotCentre  = "lines",
    xlab        = "rotation",
    ylab        = "mean NEO",
    ylim        = c(-1,2.5),
    drawRect    = FALSE)

# Datenpunkte
stripchart(D$mean_NEO ~ D$Rotation, D,
    method      = "jitter",
    xaxt        = "n",
    vertical    = TRUE,
    pch         = 19,
    col         = "black",
    add         = TRUE,
    cex         = 1)

# Mittelwerte
dx = 0.3
for (i in 1:length(levs)) {
    segments(
        x0      = i - dx,
        x1      = i + dx,
        y0      = median(D$mean_NEO[D$Rotation==levs[i]]),
        y1      = median(D$mean_NEO[D$Rotation==levs[i]]),
        col     = "black",
        lwd     = 2)
}

# Speichern
dev.copy2pdf(
    file        = "Abbildungen/kitchen_rolls_violin_zstt.pdf",
    width       = 9,
    height      = 9)
```

\vspace{2mm}
```{r, echo = F, out.width = "60%"}
knitr::include_graphics("Abbildungen/kitchen_rolls_violin_zstt.pdf")
```


<!-- Zweistichproben-T-Test -->
# Anwendung/Praxis

\small
Zweistichproben-T-Test: mittlere NEO-PI-R-Scores, im UZS vs. gegen den UZS ($n_1 = 48$, $n_2 = 54$)

\vspace{2mm}
```{r, echo = T, eval = T}
# Zweistichproben-T-Test
y1     = D$mean_NEO[D$Rotation=="clock"]                    # Gruppe 1
y2     = D$mean_NEO[D$Rotation=="counter"]                  # Gruppe 2
n1     = length(y1)                                         # Anzahl Gruppe 1
n2     = length(y2)                                         # Anzahl Gruppe 2
n      = n1 + n2                                            # Anzahl Datenpunkte
p      = 2                                                  # Anzahl Regressoren
y      = matrix(c(y1,y2), nrow = n)                         # Datenvektor
X      = matrix(c(rep(1,n1), rep(0,n2),                     # Designmatrix
                  rep(0,n1), rep(1,n2)), ncol = p)
```

```{r, echo = F, eval = T}
# Zweistichproben-T-Test
b_hat  = solve(t(X) %*% X) %*% t(X) %*% y                   # Betaparameterschätzer
e_hat  = y - X %*% b_hat                                    # Residuenvektor
s2_hat = (t(e_hat) %*% e_hat)/(n-p)                         # Varianzparameterschätzer
c      = matrix(c(1,-1), nrow = p)                          # Kontrastvektor
cTb0   = 0                                                  # Nullhypothese
t      = (t(c)%*%b_hat - cTb0)/sqrt(s2_hat * t(c)%*%solve(t(X) %*% X)%*%c)  # T-Statistik
d      = t/sqrt(n)                                          # Cohen's d 
pval   = 2*(1 - pt(abs(t), n-1))                            # p-Wert
cat(  "Datenpunkte pro Gruppe          : ", round(c(n1,n2), digits = 3),
    "\nBetaparameterschätzer           : ", round(b_hat,    digits = 3),
    "\nVarianzparameterschätzer        : ", round(s2_hat,   digits = 3),
    "\nZweistichproben-T-Teststatistik : ", round(t,        digits = 3),
    "\nCohen's d                       : ", round(d,        digits = 3),
    "\np-Wert                          : ", round(pval,     digits = 3), "\n\n")
```


<!-- Zweistichproben-T-Test -->
# Anwendung/Praxis

\small
Zweistichproben-T-Test: mittlere NEO-PI-R-Scores, Frauen vs. Männer ($n_1 = 77$, $n_2 = 25$)

\vspace{2mm}
```{r, echo = T, eval = T}
# Zweistichproben-T-Test
y1     = D$mean_NEO[D$Sex=="F"]                             # Gruppe 1
y2     = D$mean_NEO[D$Sex=="M"]                             # Gruppe 2
n1     = length(y1)                                         # Anzahl Gruppe 1
n2     = length(y2)                                         # Anzahl Gruppe 2
n      = n1 + n2                                            # Anzahl Datenpunkte
p      = 2                                                  # Anzahl Regressoren
y      = matrix(c(y1,y2), nrow = n)                         # Datenvektor
X      = matrix(c(rep(1,n1), rep(0,n2),                     # Designmatrix
                  rep(0,n1), rep(1,n2)), ncol = p)
```

```{r, echo = F, eval = T}
# Zweistichproben-T-Test
b_hat  = solve(t(X) %*% X) %*% t(X) %*% y                   # Betaparameterschätzer
e_hat  = y - X %*% b_hat                                    # Residuenvektor
s2_hat = (t(e_hat) %*% e_hat)/(n-p)                         # Varianzparameterschätzer
c      = matrix(c(1,-1), nrow = p)                          # Kontrastvektor
cTb0   = 0                                                  # Nullhypothese
t      = (t(c)%*%b_hat - cTb0)/sqrt(s2_hat * t(c)%*%solve(t(X) %*% X)%*%c)  # T-Statistik
d      = t/sqrt(n)                                          # Cohen's d 
pval   = 2*(1 - pt(abs(t), n-1))                            # p-Wert
cat(  "Datenpunkte pro Gruppe          : ", round(c(n1,n2), digits = 3),
    "\nBetaparameterschätzer           : ", round(b_hat,    digits = 3),
    "\nVarianzparameterschätzer        : ", round(s2_hat,   digits = 3),
    "\nZweistichproben-T-Teststatistik : ", round(t,        digits = 3),
    "\nCohen's d                       : ", round(d,        digits = 3),
    "\np-Wert                          : ", round(pval,     digits = 3), "\n\n")
```


<!-- Datensatz: Küchenrollen -->
# Anwendung/Praxis

\small
Datensatz: mittlere NEO-PI-R-Scores, getrennt nach Rotation und Geschlecht ($n = 102$)

```{r, echo = F, eval = F}
# Daten umsortieren
D$RotSex    = rep(NA, nrow(D))
D$RotSex[D$Rotation=="clock"   & D$Sex=="F"] = "clock-female"
D$RotSex[D$Rotation=="clock"   & D$Sex=="M"] = "clock-male"
D$RotSex[D$Rotation=="counter" & D$Sex=="F"] = "counter-female"
D$RotSex[D$Rotation=="counter" & D$Sex=="M"] = "counter-male"
levs_rs     = c("clock-female", "clock-male", "counter-female", "counter-male")
D$RotSex    = factor(D$RotSex, levels = levs_rs)

# Daten visualisieren
library(vioplot)
par(
    family      = "sans",
    mar         = c(3,3,1,1),
    pty         = "m",
    bty         = "l",
    lwd         = 1,
    las         = 1,
    mgp         = c(2,1,0),
    xaxs        = "i",
    yaxs        = "i",
    font.main   = 1.5,
    cex         = 1.5)

# Violinplot
vioplot(D$mean_NEO ~ D$RotSex, D,
    col         = "gray80",
    rectCol     = "black",
    lineCol     = "white",
    colMed      = "black",
    border      = "black",
    pchMed      = 16,
    plotCentre  = "lines",
    xlab        = "rotation and gender",
    ylab        = "mean NEO",
    ylim        = c(-1,2.5),
    drawRect    = FALSE)

# Datenpunkte
stripchart(D$mean_NEO ~ D$RotSex, D,
    method      = "jitter",
    xaxt        = "n",
    vertical    = TRUE,
    pch         = 19,
    col         = "black",
    add         = TRUE,
    cex         = 1)

# Mittelwerte
dx = 0.3
for (i in 1:length(levs_rs)) {
    segments(
        x0      = i - dx,
        x1      = i + dx,
        y0      = median(D$mean_NEO[D$RotSex==levs_rs[i]]),
        y1      = median(D$mean_NEO[D$RotSex==levs_rs[i]]),
        col     = "black",
        lwd     = 2)
}

# Speichern
dev.copy2pdf(
    file        = "Abbildungen/kitchen_rolls_violin_zva.pdf",
    width       = 16,
    height      = 9)
```

\vspace{2mm}
```{r, echo = F, out.width = "100%"}
knitr::include_graphics("Abbildungen/kitchen_rolls_violin_zva.pdf")
```


<!-- Zweifaktorielle Varianzanalyse -->
# Anwendung/Praxis

\small
$2 \times 2$ ZVA mit Interaktion: mittlere NEO-Scores, Rotation und Geschlecht ($n = 102$)

\vspace{2mm}
```{r, echo = T, eval = T}
# Zweifaktorielle Varianzanalyse
y11    = D$mean_NEO[D$Rotation=="clock"   & D$Sex=="F"]  # Gruppe 11
y12    = D$mean_NEO[D$Rotation=="clock"   & D$Sex=="M"]  # Gruppe 12
y21    = D$mean_NEO[D$Rotation=="counter" & D$Sex=="F"]  # Gruppe 21
y22    = D$mean_NEO[D$Rotation=="counter" & D$Sex=="M"]  # Gruppe 22
n      = length(y11) + length(y12) +
         length(y21) + length(y22)                       # Anzahl Datenpunkte
p      = 4                                               # Anzahl Regressoren
y      = matrix(c(y11, y12, y21, y22), nrow = n)         # Datenvektor
```

```{r, echo = F, eval = T}
# Zweifaktorielle Varianzanalyse
n11    = length(y11)                                     # Anzahl Datenpunkte Gr. 11
n12    = length(y12)                                     # Anzahl Datenpunkte Gr. 12
n21    = length(y21)                                     # Anzahl Datenpunkte Gr. 21
n22    = length(y22)                                     # Anzahl Datenpunkte Gr. 22
X      = matrix(c(rep(1,n),                              # Designmatrix vollständiges Modell
                  c(rep(0,n11+n12), rep(1,n21+n22)),
                  c(rep(0,n11), rep(1,n12), rep(0,n21), rep(1,n22)),
                  c(rep(0,n11+n12+n21), rep(1,n22))), nrow = n)
X_0    = X[,-4]                                          # Designmatrix reduziertes Modell

# Zweifaktorielle Varianzanalyse
p          = ncol(X)                                     # Anzahl Parameter vollständiges Modell
p_0        = ncol(X_0)                                   # Anzahl Parameter reduziertes Modell
p_1        = p - p_0                                     # Anzahl zusätzlicher Parameter im vollständigen Modell
beta_hat_0 = solve(t(X_0)%*%X_0) %*% t(X_0)%*%y          # Betaparameterschätzer reduziertes Modell
beta_hat   = solve(t(X)  %*%X )  %*% t(X)  %*%y          # Betaparameterschätzer vollständiges Modell
eps_hat_0  = y - X_0 %*% beta_hat_0                      # Residuenvektor reduziertes Modell
eps_hat    = y - X   %*% beta_hat                        # Residuenvektor vollständiges Modell
eh0_eh0    = t(eps_hat_0) %*% eps_hat_0                  # residuelle QS reduziertes Modell
eh_eh      = t(eps_hat)   %*% eps_hat                    # residuelle QS vollständiges Modell
sigsqr_hat = eh_eh/(n-p)                                 # Varianzparameterschätzer vollständiges Modell
f          = ((eh0_eh0-eh_eh)/p_1)/sigsqr_hat            # F-Statistik
p_val      = 1 - pf(f, p_1,n-p)                          # p-Wert
cat(  "Datenpunkte pro Gruppe            : ", c(n11,n12,n21,n22),
    "\nBetaparameterschätzer             : ", round(beta_hat,   digits = 3),
    "\nBetaparameterschätzer red. Mod.   : ", round(beta_hat_0, digits = 3),
    "\nresiduelle Quadratsumme           : ", round(eh_eh,      digits = 3),
    "\nresiduelle Quadratsumme red. Mod. : ", round(eh0_eh0,    digits = 3),
    "\nF-Teststatistik                   : ", round(f,          digits = 3),
    "\np-Wert                            : ", round(p_val,      digits = 3), "\n\n")
```


<!-- p-Werte (1) -->
# Anwendung/Praxis

```{r, echo = F, out.width = "80%"}
knitr::include_graphics("Abbildungen/pvalue1.png")
```

\vspace{2mm}
\center
\footnotesize
Webcomic zum p-Wert (Quelle: [*xkcd*: "Significant" (#882)](https://xkcd.com/882/); Lizenz: CC-BY-NC 2.5).


<!-- p-Werte (2) -->
# Anwendung/Praxis

\vspace{-2mm}
```{r, echo = F, out.width = "50%"}
knitr::include_graphics("Abbildungen/pvalue2.png")
```

\vspace{0.5mm}
\center
\footnotesize
Webcomic zum p-Wert (Quelle: [*xkcd*: "Significant" (#882)](https://xkcd.com/882/); Lizenz: CC-BY-NC 2.5).


<!-- p-Werte (3) -->
# Anwendung/Praxis

```{r, echo = F, out.width = "80%"}
knitr::include_graphics("Abbildungen/pvalue3.png")
```

\vspace{2mm}
\center
\footnotesize
Webcomic zum p-Wert (Quelle: [*xkcd*: "Significant" (#882)](https://xkcd.com/882/); Lizenz: CC-BY-NC 2.5).


<!-- Anhang: Selbstkontrollfragen -->
#

\large
\setstretch{2.5}
\vfill

Anwendungsszenario

Modellformulierung

Modellschätzung

Modellevaluation

**Selbstkontrollfragen**

\vfill


<!-- Selbstkontrollfragen -->
# Selbstkontrollfragen

\footnotesize
\begin{enumerate}

\justifying
\item Erläutern Sie das Anwendungsszenario der zweifaktoriellen Varianzanalyse (ZVA).
\item Aus wie vielen Datenpunkten besteht ein Datensatz eines 3 x 4 ZVA-Designs mit 10 Datenpunkten pro Zelle?
\item Erläutern Sie die intuitive Bedeutung eines Haupteffektes in einem ZVA-Design.
\item Erläutern Sie die intuitive Bedeutung einer Interaktion ein einem ZVA-Design.
\item Geben Sie die Definition des additiven Modells der ZVA mit Referenzgruppe (RG) wieder.
\item Erläutern Sie die Bedeutung der Parameter $\mu_0,\alpha_2$ und $\beta_2$ im additiven Modell der ZVA mit RG.
\item Bestimmen Sie $\mu_{ij}$ für $\mu_0 := 2, \alpha_2 := -1, \beta_2 := 3$ im additiven Modell der ZVA mit RG.
\item Geben Sie die Designmatrixform des Modells einer additiven 2 x 2 ZVA mit RG für $n_{ij} := 1$ an.
\item Geben Sie die Designmatrixform des Modells einer additiven 2 x 2 ZVA mit RG für $n_{ij} := 3$ an.
\item Geben Sie die Definition des Modells der ZVA mit Interaktion und Referenzgruppe wieder.
\item Erläutern Sie die Bedeutung der Parameter $\mu_0, \alpha_2, \beta_2$ und $\gamma_{22}$ im Modell der ZVA mit Interaktion und RG.
\item Bestimmen Sie $\mu_{ij}$ für $\mu_0 := 2, \alpha_2 := 1, \beta_2 := -1, \gamma_{22} := 3$ im Modell der ZVA mit Interaktion und RG.
\item Geben Sie die Designmatrixform des Modells einer 2 x 2 ZVA mit Interaktion und RG für $n_{ij} := 1$ an.
\item Geben Sie die Designmatrixform des Modells einer 2 x 2 ZVA mit Interaktion und RG für $n_{ij} := 3$ an.
\item Geben Sie das Theorem zur Betaparameterschätzung im 2 x 2 ZVA-Modell mit Interaktion und RG wieder.
\item Wie muss das Model der 2 x 2 ZVA mit Interaktion und Referenzgruppe partitioniert werden, um eine F-Teststatistik für die Interaktion von Faktor A und B zu bestimmen?

\end{enumerate}
\vfill


<!-- Literatur -->
# Referenzen
