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


# {.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 2026

\vspace{5mm}
Joram Soch


<!-- Sitzung 15: Kovarianzanalyse -->
# {.plain}

\vfill
\center
\huge
\textcolor{black}{(15) Kovarianzanalyse}
\vfill


<!-- Inhaltsverzeichnis -->
#

\large
\setstretch{3}
\vfill

Vorbemerkungen

Additive Kovarianzanalyse

Kovarianzanalyse mit Interaktion

Anwendung/Praxis

Selbstkontrollfragen

\vfill


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

\large
\setstretch{3}
\vfill

**Vorbemerkungen**

Additive Kovarianzanalyse

Kovarianzanalyse mit Interaktion

Anwendung/Praxis

Selbstkontrollfragen

\vfill


<!-- Überblick -->
# Vorbemerkungen

\setstretch{2}
\textcolor{darkblue}{Variabilität der Werte einer abhängigen Variable}

\small
* Die Werte einer AV unterscheiden sich im Allgemeinen zwischen experimentellen Einheiten.
* Die Schwankungen der Werte einer AV bezeichnet man als *Datenvariabilität*.
* Ziel jeder Datenanalyse ist die Erklärung von Datenvariabilität durch Zerlegung.

\normalsize
\textcolor{darkblue}{Datenvariabilität und Allgemeines Lineares Modell}

\small
* Eine spezielle Art der Quantifizierung von Datenvariabilität ist die Stichprobenvarianz.
* Einen additiven Zugang zur Dekomposition von Datenvariabilität bietet das ALM.
* Den folgenden Überlegungen liegt datenanalytisch die Kovarianzanalyse zugrunde.


<!-- Variabilitätsquellen (1) -->
# Vorbemerkungen

\large
\textcolor{darkblue}{Gesamtvariabilität = Primärvariabilität + Residualvariabilität}

\vspace{2mm}
\normalsize
\textcolor{darkblue}{Primärvariabilität}

Systematische Veränderung der AV, die allein auf Variation der UV zurückzuführen ist.

\vspace{2mm}
\large
\textcolor{darkblue}{Residualvariabilität = Sekundärvariabilität + Zufallsfehler}

\vspace{2mm}
\normalsize
\textcolor{darkblue}{Sekundärvariabilität}

Systematische Veränderung der AV, die auf die Wirkung von unkontrollierten Störvariablen, nicht aber auf die Variation der UV zurückzuführen ist.

\vspace{2mm}
\normalsize
\textcolor{darkblue}{Zufallsfehler}

Unsystematische Veränderung der AV, die weder auf die Variation der UV, noch auf den Einfluss von Störvariablen zurückzuführen ist.


<!-- Variabilitätsquellen (2) -->
# Vorbemerkungen

```{r, echo = FALSE, out.width = "100%"}
knitr::include_graphics("Abbildungen/variabilitätsquellen.pdf")
```


<!-- Parametrische vs. faktorielle Designs -->
# Vorbemerkungen

\vspace{2mm}
\textcolor{darkblue}{Faktorielle und parametrische ALM-Designs}

\small
Faktorielle ALM-Designs:
\vspace{-2mm}

* Designmatrizen mit $1$en und $0$en, manchmal $-1$en.
* Betaparameter repräsentieren Erwartungswerte für Gruppen.
* Betaparameterschätzer repräsentieren Stichprobenmittel für Gruppen.
* $\Rightarrow$ T-Tests, einfaktorielle Varianzanalyse, zweifaktorielle Varianzanalyse

Parametrische ALM-Designs:
\vspace{-2mm}

* Designmatrizen besitzen Spalten mit kontinuierlichen reellen Werten.
* Die Designmatrixsspalten werden *Regressoren*, *Prädiktoren*, oder *Kovariaten* genannt.
* Betaparameter repräsentieren Steigungsparameter.
* Betaparameterschätzer ergeben sich als normalisierte Kovarianzen von Regressoren und Daten.
* Es besteht ein enger Bezug zur Theorie der Korrelation.
* $\Rightarrow$ einfache lineare Regression, multiple lineare Regression

Faktoriell-parametrische ALM-Designs:
\vspace{-2mm}

* Designmatrizen mit mehreren faktoriellen und parametrischen Werten.
* Die parametrischen Regressoren werden oft als kontrollierte Kovariaten betrachtet.
* $\Rightarrow$ Kovarianzanalyse $\Leftrightarrow$ Allgemeines Lineares Modell

\center
\normalsize
\textcolor{darkblue}{Wir betrachten hier exemplarisch den Fall 1 Faktor + 1 parametrischer Regressor.}


<!-- Abschnitt 2: Additive Kovarianzanalyse -->
#

\large
\setstretch{3}
\vfill

Vorbemerkungen

**Additive Kovarianzanalyse**

Kovarianzanalyse mit Interaktion

Anwendung/Praxis

Selbstkontrollfragen

\vfill


<!-- Überblick -->
# Additive Kovarianzanalyse

\footnotesize
\textcolor{darkblue}{Modell der additiven einfaktoriellen Kovarianzanalyse}

Im Modell der additiven einfaktoriellen Kovarianzanalyse möchte man die Erwartungswertparameter $\mu_{ij}$ der insgesamt $n = \sum_{i=1}^p n_i$ Datenvariablen
\begin{equation}
y_{ij} \sim N(\mu_{ij}, \sigma^2 I_n)
\end{equation}

basierend auf den Leveln $i = 1,...,p$ des experimentellen Faktors und den entsprechenden Werten $x_{ij}$ für $j = 1,...,n_i$ der Kovariate modellieren.

Dazu wählt man für jeden Erwartungswertparameter $\mu_{ij}$ zunächst einen Parameter $\mu_0$, der einen Datenvariablen-unspezifischen Offset modelliert. Weiterhin wählt man einen Parameter $\alpha_i$, der den Beitrag des $i$ten Levels des experimentellen Faktors modelliert und schließlich einen Gewichtungsparameter $\beta_0$, der den Beitrag der Kovariatenwerten $x_{ij}$ zu $\mu_{ij}$ quantifiziert.

Das Modell für den Erwartungswertparameter von $y_{ij}$ für $i = 1,...,p$, $j = 1,...,n_i$ nimmt im Fall der additiven einfaktoriellen Kovarianzanalyse also folgende Form an:
\begin{equation}
\mu_{ij} := \mu_0 + \alpha_i + \beta_0x_{ij}
\quad \mbox{mit} \quad
\alpha_1 := 0,
\end{equation}

wobei $\alpha_1 := 0$ wie im Falle der einfaktoriellen Varianzanalyse eine Überparameterisierung des Modells verhindert.

Man beachte, dass die Form des $ij$ten Erwartungswertparameters damit eine einfache lineare Regression für jedes Level $i = 1,...,p$ des experimentellen Faktors definiert, wobei $\mu_0 + \alpha_i$ als Faktorlevel-spezifische Offsetparameter interpretiert werden können und der Steigungsparameter $\beta_0$ für alle Faktorlevel identisch ist.


<!-- Definition: Modell der additiven einfaktoriellen Kovarianzanalyse -->
# Additive Kovarianzanalyse

\footnotesize
\begin{definition}[Modell der additiven einfaktoriellen Kovarianzanalyse]

\justifying
$y_{ij}$ mit $i = 1,...,p$, $j = 1,...,n_i$ sei die Zufallsvariable, die den zum $i$ten Level des Faktors gehörenden $j$ten Datenpunkt modelliert und $x_{ij}$ sei der entsprechende Wert der Kovariate. Dann hat das \textit{Modell der additiven einfaktoriellen Kovarianzanalyse} die strukturelle Form
\begin{equation}
y_{ij} = \mu_{ij} + \varepsilon_{ij}
\quad \mbox{mit} \quad
\varepsilon_{ij} \sim N(0,\sigma^2)
\quad \mbox{u.i.v. für} \quad
i = 1,...,p, \; j = 1,...,n_i, \; \sigma^2 > 0
\end{equation}

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

wobei
\begin{equation}
\mu_{ij} := \mu_0 + \alpha_i + \beta_0x_{ij}
\quad \mbox{mit} \quad
\alpha_1 := 0 \; .
\end{equation}

\end{definition}

Bemerkungen

\begin{itemize}

\item Die Form des $ij$ten Erwartungswertparameters 
\begin{equation}
\mu_{ij} := \mu_0 + \alpha_i + \beta_0x_{ij}
\quad \mbox{mit} \quad
\alpha_1 := 0
\end{equation}
definiert eine einfache lineare Regression für jedes $i = 1,...,p$, wobei
\vspace{1mm}

\begin{itemize}
\begin{footnotesize}
\item $\mu_0 + \alpha_i$ ein Faktorlevel-spezifischer Offsetparameter ist und
\item der Steigungsparameter $\beta_0$ für alle Faktorlevel identisch ist.
\end{footnotesize}
\end{itemize}

\end{itemize}


<!-- Theorem: Designmatrixform der additiven einfaktoriellen Kovarianzanalyse -->
# Additive Kovarianzanalyse

\vspace{2mm}
\footnotesize
\begin{theorem}[\begin{small}Designmatrixform der additiven einfaktoriellen Kovarianzanalyse\end{small}]

\justifying
\normalfont
Gegeben sei die strukturelle Form eines additiven einfaktoriellen Kovarianzanalysemodells mit Anzahl der Faktorlevel $p$ und es sei $n := \sum_{i=1}^p n_i$ 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_{11}   \\
\vdots   \\
y_{1n_1} \\
y_{21}   \\
\vdots   \\
y_{2n_2} \\
         \\
\vdots   \\
         \\
y_{p1}   \\
\vdots   \\
y_{pn_p} \\
\end{pmatrix}, \quad
X := 
\begin{pmatrix}
1       & 0       &        & 0      & x_{11}   \\
\vdots  & \vdots  & \cdots & \vdots & \vdots   \\
1       & 0       &        & 0      & x_{1n_1} \\
1       & 1       &        & 0      & x_{21}   \\
\vdots  & \vdots  & \cdots & \vdots & \vdots   \\
1       & 1       &        & 0      & x_{2n_2} \\
        &         &        &        &          \\
\vdots  & \vdots  &        & \cdots & \vdots   \\
        &         &        &        &          \\
1       & 0       &        & 1      & x_{p1}   \\
\vdots  & \vdots  & \cdots & \vdots & \vdots   \\
1       & 0       &        & 1      & x_{pn_p} \\
\end{pmatrix}
\in \mathbb{R}^{n \times (p+1)}, \quad
\beta := 
\begin{pmatrix}
\mu_0    \\
\alpha_2 \\
\vdots   \\
\alpha_p \\
\beta_0  \\
\end{pmatrix} 
\in
\mathbb{R}^{p+1}
\quad \mbox{und} \quad
\sigma^2 > 0.
\end{equation}

\end{theorem}

\vspace{-1mm}
Bemerkungen
\vspace{-1mm}

* Das Theorem ergibt sich direkt mit den Regeln der Matrixmultiplikation.


<!-- Designmatrix des Modells -->
# Additive Kovarianzanalyse

Designmatrix der additiven einfaktoriellen Kovarianzanalyse ($n = 18$, $p = 3$, $p+1 = 4$)

```{r, echo = F, eval = F}
# Designmatrixerzeugung
set.seed(1)                                 # reproduzierbare Ergebnisse
e     = 0.5                                 # Wert für Null in Visualisierung
de    = 0.1                                 # vom Wert zu korrigierende Abweichung
n     = 18                                  # Anzahl von Datenpunkten
p     = 3                                   # Anzahl von Faktorleveln
m     = round(n/p)                          # Anzahl Datenpunkte pro Gruppe
Xt    = cbind(matrix(1, nrow = n, ncol = 1),# Varianzanalyse-Designmatrix
              kronecker(matrix(e, nrow = p, ncol = p) + (1-e)*diag(p),
                        matrix(1, nrow = m, ncol = 1)))
X     = Xt[,-2]                             # eliminiere zweite Spalte
x     = runif(n, -1, +1)                    # Werte der Kovariate
x[x>e-de & x<e] = x[x>e-de & x<e] - de
x[x>e & x<e+de] = x[x>e & x<e+de] + de
X     = cbind(X, matrix(x, nrow = n))       # n x p Designmatrix
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(-1, 1, length.out = 257),
    key       = NULL,
    axis.col  = NULL,
    axis.row  = NULL,
    asp       = 1)

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

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


<!-- Anwendungsbeispiel -->
# Additive Kovarianzanalyse

\normalsize
\setstretch{2}
\textcolor{darkblue}{Anwendungsbeispiel}

experimenteller Faktor/Studiengruppe ($p := 2, \; i = 1,2$)

* Face-to-face- vs. Online-Psychotherapie

parametrischer Regressor/Kovariate ($x_{ij}, \; j = 1,...,n_i$)  

* Dauer der Depressionssymptomatik zu Beginn der Intervention 

abhängige Variable/primäres Ergebnismaß ($y_{ij}, \; j = 1,...,n_i$)

* Prä-/Post-Differenzwerte von BDI-Scores

\vfill


<!-- Beispieldatensatz -->
# Additive Kovarianzanalyse

\vspace{2mm}
\textcolor{darkblue}{Beispieldatensatz} ($n = 32$; hier: $i = 1,2$, $j = 1,...,10$)

```{r, echo = F}
library(MASS)                                                                   # multivariate Normalverteilung
set.seed(1)                                                                     # reproduzierbare Ergebnisse
n_i     = c(16,16)                                                              # Anzahl Patient:innen pro Gruppe
n       = sum(n_i)                                                              # Gesamtzahl Datenpunkte
x       = round(rnorm(n,5,2), digits = 1)                                       # Depressionsdauerkovariate
p       = 3                                                                     # Anzahl Betaparameter
beta    = matrix(c(10,5,-1), nrow = p)                                          # wahre,  aber unbekannte  Betaparameter
sigsqr  = 4                                                                     # wahrer, aber unbekannter Varianzparameter
X       = matrix(c(rep(1,n_i[1]), rep(1,n_i[2]),                                # \mu_0    Regressor
                   rep(0,n_i[1]), rep(1,n_i[2]),                                # \alpha_2 Regressor
                   x), nrow = n)                                                # \beta_0  Regressor
y       = mvrnorm(1, X %*% beta, sigsqr*diag(n))                                # Realisierung der Daten
y       = round(y)                                                              # diskrete BDI-Werte
D       = data.frame(Therapie  = c(rep("F2F", n_i[1]),                          # Studiengruppe F2F
                                   rep("ONL", n_i[2])),                         # Studiengruppe ONL
                     Dauer     = x,                                             # Depressionsdauer vor Therapie
                     dBDI      = y)                                             # Delta-BDI-Werte
write.csv(D, file = "Daten/Kovarianzanalyse_A.csv")                             # Datenspeicherung
```

\footnotesize
\setstretch{1}
```{r, echo = F}
fname      = "Daten/Kovarianzanalyse_A.csv"
D          = read.table(fname, sep = ",", header = TRUE)
knitr::kable(D[c(1:10,(n_i[1]+1):(n_i[1]+10)),2:4], "pipe", align = "c")          # erste 10 Patient:innen pro Gruppe
```


<!-- Visualisierung des Beispieldatensatzes -->
# Additive Kovarianzanalyse

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

```{r, echo = F, eval = F}
# Laden des Datensatzes
fname      = "Daten/Kovarianzanalyse_A.csv"
D          = read.table(fname, sep = ",", header = TRUE)

# Deskriptivstatistikanalyse
thp        = c("F2F", "ONL")                                                    # Therapiebedingungen
y_thp      = list()                                                             # therapiespezifische dBDI-Werte
x_thp      = list()                                                             # therapiespezifische Dauer-Werte
for(i in 1:length(thp)){                                                        # Iteration über Therapiebedingungen
    y_thp[[i]]    = D$dBDI[D$Therapie == thp[i]]                                # therapiespezifische dBDI-Werte
    x_thp[[i]]    = D$Dauer[D$Therapie == thp[i]]                               # therapiespezifische Dauer-Werte
}

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

# Punktwolken
plot(x_thp[[1]], y_thp[[1]],
    pch        = 16,
    xlab       = "Dauer",
    ylab       = "dBDI",
    xlim       = c(-1,11),
    ylim       = c(-2,20))
points(x_thp[[2]], y_thp[[2]],
    pch        = 16,
    col        = "gray80")
legend("topright", c("F2F", "ONL"),
    lty         = 0,
    pch         = 16,
    col         = c("black", "gray80"),
    bty         = "n",
    cex         = 1,
    x.intersp   = 1)

# Balkendiagramm                    
groupmeans = c(mean(y_thp[[1]]), mean(y_thp[[2]]))
groupstds  = c(sd(y_thp[[1]])  , sd(y_thp[[2]]))
names(groupmeans) = c("F2F", "ONL")
x = barplot(groupmeans,                                      
    ylim        = c(-2,20),
    xlim        = c(-.4,3),
    col         = c("black", "gray90"),
    ylab        = "dBDI",
    xlab        = "Therapie")
arrows(
    x0          = x,                                    
    y0          = groupmeans - groupstds,               
    x1          = x,                                    
    y1          = groupmeans + groupstds,               
    code        = 3,                                   
    angle       = 90,                                  
    length      = 0.05)   

# Speicherung
dev.copy2pdf(
    file        = "Abbildungen/kovarianzanalyse_A1.pdf",
    width       = 9,
    height      = 4.5)

# Ausgabe
cat(  "Mittelwert und Standardabweichung dBDI  F2F : ", groupmeans[1],    groupstds[1], 
    "\nMittelwert und Standardabweichung dBDI  ONL : ", groupmeans[2],    groupstds[2],
    "\nMittelwert und Standardabweichung Dauer F2F : ", mean(x_thp[[1]]), sd(x_thp[[1]]),
    "\nMittelwert und Standardabweichung Dauer ONL : ", mean(x_thp[[2]]), sd(x_thp[[2]]))
```

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


<!-- Modellformulierung -->
# Additive Kovarianzanalyse

\textcolor{darkblue}{Modellformulierung}

\normalsize
Modell 1 | keine Berücksichtigung der Kovariate
\vspace{-2mm}

\footnotesize
\begin{align}
&y_{ij} \sim N(\mu_{ij}, \sigma^2)
\quad \mbox{mit} \quad
\mu_{ij} := \mu_0 + \alpha_i
\quad \mbox{und} \quad
\alpha_1 := 0
\quad \mbox{und} \quad
\sigma^2 > 0 \\ 
\Leftrightarrow \quad
&y \sim N(X\beta,\sigma^2 I_n)
\quad \mbox{mit} \quad
X := 
\begin{pmatrix}
1      & 0      \\
\vdots & \vdots \\
1      & 0      \\
1      & 1      \\
\vdots & \vdots \\
1      & 1      \\
\end{pmatrix}, \quad
\beta :=
\begin{pmatrix}
\mu_0    \\
\alpha_2 \\
\end{pmatrix}
\quad \mbox{und} \quad
\sigma^2 > 0
\end{align}

\normalsize
Modell 2 | additive Berücksichtigung der Kovariate
\vspace{-2mm}

\footnotesize
\begin{align}
&y_{ij} \sim N(\mu_{ij}, \sigma^2)
\quad \mbox{mit} \quad
\mu_{ij} := \mu_0 + \alpha_i + \beta_0 x_{ij}
\quad \mbox{und} \quad
\alpha_1 := 0
\quad \mbox{und} \quad
\sigma^2 > 0 \\ 
\Leftrightarrow \quad
&y \sim N(X\beta,\sigma^2 I_n)
\quad \mbox{mit} \quad
X := 
\begin{pmatrix}
1      & 0      & x_{11}   \\
\vdots & \vdots & \vdots   \\
1      & 0      & x_{1n_1} \\
1      & 1      & x_{21}   \\
\vdots & \vdots & \vdots   \\
1      & 1      & x_{2n_2} \\
\end{pmatrix}, \quad
\beta := 
\begin{pmatrix}
\mu_0    \\
\alpha_2 \\
\beta_0  \\
\end{pmatrix}
\quad \mbox{und} \quad
\sigma^2 > 0
\end{align}


<!-- Implementation -->
# Additive Kovarianzanalyse

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

\tiny
\setstretch{1}
```{r, echo = T}
fname     = "Daten/Kovarianzanalyse_A.csv"                              # Dateiname
D         = read.table(fname, sep = ",", header = TRUE)                 # Datensatz
n_i       = c(sum(D$Therapie == "F2F"), sum(D$Therapie == "ONL"))       # Anzahl Datenpunkte pro Gruppe
y         = D$dBDI                                                      # Daten
n         = length(y)                                                   # Gesamtanzahldatenpunkte
XS        = list()                                                      # Modell 1 und 2 Liste
CS        = list()                                                      # Kontrastgewichtsvektorenliste
XS[[1]]   = matrix(c(rep(1,n_i[1]),rep(1,n_i[2]),                       # \mu_0    Regressor                            
                     rep(0,n_i[1]),rep(1,n_i[2])), nrow = n)            # \alpha_2 Regressor
XS[[2]]   = matrix(c(rep(1,n_i[1]),rep(1,n_i[2]),                       # \mu_0    Regressor                       
                     rep(0,n_i[1]),rep(1,n_i[2]),                       # \alpha_2 Regressor
                     D$Dauer), nrow = n)                                # x        Regressor
CS[[1]]   = matrix(c(0,1)  , nrow = 2)                                  # Kontrastgewichtsvektor \alpha_2
CS[[2]]   = matrix(c(0,1,0), nrow = 3)                                  # Kontrastgewichtsvektor \alpha_2
B         = list()                                                      # Betaparameterschätzerliste
S         = rep(NaN,2)                                                  # Varianzparameterschätzervektor
T         = rep(NaN,2)                                                  # T-Statistik-Vektor
for(i in 1:2){                                                          # Iteration über Modelle
    X           = XS[[i]]                                               # Designmatrix
    p           = ncol(X)                                               # Anzahl Betaparameter
    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
    c           = CS[[i]]                                               # Kontrastgewichtsvektor
    t_num       = t(c) %*% beta_hat                                     # Zähler der Zweistichproben-T-Teststatistik
    t_den       = sqrt(sigsqr_hat %*% t(c) %*% solve(t(X) %*% X) %*% c) # Nenner der Zweistichproben-T-Teststatistik
    t           = t_num/t_den                                           # Wert   der Zweistichproben-T-Teststatistik
    B[[i]]      = beta_hat                                              # Betaparameterschätzer
    S[i]        = sigsqr_hat                                            # Varianzparameterschätzer
    T[i]        = t                                                     # T-Statistik
}
```


<!-- Ergebnisse -->
# Additive Kovarianzanalyse

\textcolor{darkblue}{Ergebnisse}

\small
```{r, echo = F}
cat(  "Modell 1 | Betaparameterschätzer    : ", round(B[[1]],2),
    "\nModell 2 | Betaparameterschätzer    : ", round(B[[2]],2),
    "\nModell 1 | Varianzparameterschätzer : ", round(S[[1]],2),
    "\nModell 2 | Varianzparameterschätzer : ", round(S[[2]],2),
    "\nModell 1 | T-Statistik ONL-Therapie : ", round(T[[1]],2),
    "\nModell 2 | T-Statistik ONL-Therapie : ", round(T[[2]],2))
```

In Modell 2 ist der Varianzparameterschätzer kleiner als in Modell 1.

$\Rightarrow$ Der parametrische Regressor erklärt Datenvarianz.

Die Betaparameterschätzer für den ONL-Therapieeffekt sind in Modell 1 und Modell 2 ähnlich.

$\Rightarrow$ Die T-Statistik (Signal-zu-Rauschen-Verhältnis) für den ONL-Therapieffekt ist Modell 2 höher.

Prinzipiell kann man die beobachteten Daten für den geschätzten Dauer Effekt durch

\begin{equation}
\hat{y}_{ij}^{\tiny\mbox{adj}} = y_{ij} - \hat{\beta}_0x_{ij}
\end{equation}

"korrigieren" und erhält so einen "adjustierten Datensatz".


<!-- Visualisierung -->
# Additive Kovarianzanalyse

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

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

# Punktwolken
plot(D$Dauer[D$Therapie == "F2F"], D$dBDI[D$Therapie == "F2F"],
    pch        = 16,
    xlab       = "Dauer",
    ylab       = "dBDI",
    xlim       = c(-1,11),
    ylim       = c(-2,20))
points(D$Dauer[D$Therapie == "ONL"], D$dBDI[D$Therapie == "ONL"],
    pch        = 16,
    col        = "gray80")

# prädizierte Daten
X              = XS[[2]]
beta_hat       = B[[2]]
y_hat          = X %*% beta_hat
lines(D$Dauer[D$Therapie == "F2F"], y_hat[1:n_i[1]],
    col        = "black")
lines(D$Dauer[D$Therapie == "ONL"], y_hat[(n_i[1]+1):n],
    col        = "gray80")
legend("topright", c("F2F", "ONL"),
    lty        = 0,
    pch        = 16,
    col        = c("black", "gray80"),
    bty        = "n",
    cex        = 1,
    x.intersp  = 1)

# adjustierte Daten
y_adj       = y - X[,3]*beta_hat[3]
groupmeans  = c(mean(y_adj[1:n_i[1]]), mean(y_adj[(n_i[1]+1):n]))
groupstds   = c(sd(y_adj[1:n_i[1]]), sd(y_adj[(n_i[1]+1):n]))
names(groupmeans) = c("F2F", "ONL")
x = barplot(groupmeans,                                      
    ylim        = c(-2,20),
    xlim        = c(-.4,3),
    col         = c("black", "gray90"),
    ylab        = "dBDI",
    xlab        = "Therapie",
    main        = "Adjustierte Daten")
arrows(
    x0          = x,                                    
    y0          = groupmeans - groupstds,               
    x1          = x,                                    
    y1          = groupmeans + groupstds,               
    code        = 3,                                   
    angle       = 90,                                  
    length      = 0.05)   

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

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


<!-- Abschnitt 3: Kovarianzanalyse mit Interaktion -->
#

\large
\setstretch{3}
\vfill

Vorbemerkungen

Additive Kovarianzanalyse

**Kovarianzanalyse mit Interaktion**

Anwendung/Praxis

Selbstkontrollfragen

\vfill


<!-- Überblick -->
# Kovarianzanalyse mit Interaktion

\footnotesize
\vspace{1mm}
\textcolor{darkblue}{Modell der einfaktoriellen Kovarianzanalyse mit Interaktion}

Wie im Modell der additiven einfaktoriellen Kovarianzanalyse möchte man im Modell der einfaktoriellen Kovarianzanalyse mit Interaktion die Erwartungswertparameter $\mu_{ij}$ der insgesamt $n = \sum_{i=1}^p n_i$ Datenvariablen
\begin{equation}
y_{ij} \sim N(\mu_{ij}, \sigma^2 I_n) 
\end{equation}

\vspace{-1mm}
basierend auf den Leveln $i = 1,...,p$ des experimentellen Faktors und den entsprechenden Werten $x_{ij}$ für $j = 1,...,n_i$ der Kovariate modellieren. Über das Szenario der additiven einfaktoriellen Kovarianzanalyse hinaus möchte man dabei explizit einen Faktorlevel-spezifischen Einfluss der Kovariate auf die Erwartungswertparameter zulassen. 

Dazu wählt man für jeden Erwartungswertparameter $\mu_{ij}$ zunächst wieder einen Parameter $\mu_0$, der einen Datenvariablen-unspezifischen Offset modelliert. Weiterhin wählt man einen Parameter $\alpha_i$, der den Beitrag des $i$ten Levels des experimentellen Faktors modelliert. Den Wichtungsparameter des Beitrags der Kovariaten $x_{ij}$ zu $\mu_{ij}$ modelliert man nun hier mithilfe eines Faktorlevel-unspezifischen Parameter $\beta_0$ \underline{und} eines Faktorlevel-spezifischen Parameters $\gamma_i$ für $i = 1,...,p$. Das Modell für den Erwartungswertparameter von $y_{ij}$ für $i = 1,...,p$, $j = 1,...,n_i$ nimmt im Fall der einfaktoriellen Kovarianzanalyse mit Interaktion damit folgende Form an:
\begin{equation}
\mu_{ij} := \mu_0 + \alpha_i + (\beta_0 + \gamma_i)x_{ij}
\quad \mbox{mit} \quad
\alpha_1 := \gamma_1 := 0,
\end{equation}

\vspace{-1mm}
wobei $\alpha_1 := \gamma_1 := 0$ wie im Fall der zweiaktoriellen Varianzanalyse mit Interaktion eine Überparmeterisierung des Modells verhindert.

Man beachte, dass die Form des $ij$ten Erwartungswertparameters damit eine einfache lineare Regression für jedes Level $i = 1,...,p$ des experimentellen Faktors  definiert, bei $\mu_0 + \alpha_i$ als Faktorlevel-spezifische Offsetparameter interpretiert werden können und sich die Steigungsparameter $\beta_0 + \gamma_i$ über die Faktorlevel unterscheiden können. Es werden also Unterschiede (zwischen den Faktorleveln) von Unterschieden (zwischen den Erwartungswertparametern eines Faktorlevels in Abhängigkeit vom Wert der Kovariate), d.h. Interaktionen, mit modelliert.


<!-- Definition: Modell der einfaktoriellen Kovarianzanalyse mit Interaktion -->
# Kovarianzanalyse mit Interaktion

\footnotesize
\begin{definition}[Modell der einfaktoriellen Kovarianzanalyse mit Interaktion]

\justifying
$y_{ij}$ mit $i = 1,...,p$, $j = 1,...,n_i$ sei die Zufallsvariable, die den zum $i$ten Level des Faktors gehörenden $j$ten Datenpunkt modelliert und $x_{ij}$ sei der entsprechende Wert der Kovariate. Dann hat das \textit{Modell einfaktoriellen Kovarianzanalyse mit Interaktion} die strukturelle Form
\begin{equation}
y_{ij} = \mu_{ij} + \varepsilon_{ij}
\quad \mbox{mit} \quad
\varepsilon_{ij} \sim N(0,\sigma^2)
\quad \mbox{u.i.v. für} \quad
i = 1,...,p, \; j = 1,...,n_i, \; \sigma^2 > 0 
\end{equation}

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

wobei
\begin{equation}
\mu_{ij} := \mu_0 + \alpha_i + (\beta_0 + \gamma_i)x_{ij}
\quad \mbox{mit} \quad
\alpha_1 := \gamma_1 := 0.
\end{equation}

\end{definition}

Bemerkungen

\begin{itemize}

\item Die Form des $ij$ten Erwartungswertparameters 
\begin{equation}
\mu_{ij} := \mu_0 + \alpha_i + (\beta_0 + \gamma_i)x_{ij}
\quad \mbox{mit} \quad
\alpha_1 := \gamma_1 := 0
\end{equation}
definiert eine einfache lineare Regression für jedes $i = 1,...,p$, wobei
\vspace{1mm}

\begin{itemize}
\begin{footnotesize}
\item $\mu_0   + \alpha_i$ ein Faktorlevel-spezifischer Offsetparameter ist und
\item $\beta_0 + \gamma_i$ ein Faktorlevel-spezifischer Steigungsparameter ist.
\end{footnotesize}
\end{itemize}

\end{itemize}


<!-- Theorem: Designmatrixform der einfaktoriellen Kovarianzanalyse mit Interaktion -->
# Kovarianzanalyse mit Interaktion

\vspace{2mm}
\footnotesize
\begin{theorem}[\begin{small}Designmatrixform der einfaktoriellen Kovarianzanalyse mit Interaktion\end{small}]

\justifying
\normalfont
Gegeben sei die strukturelle Form eines einfaktoriellen Kovarianzanalysemodells mit Interaktion mit Anzahl der Faktorlevel $p$ und es sei $n := \sum_{i=1}^N n_i$ 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
\tiny
\begin{equation}
y := 
\begin{pmatrix}
y_{11}   \\
\vdots   \\
y_{1n_1} \\
y_{21}   \\
\vdots   \\
y_{2n_2} \\
         \\
\vdots   \\
         \\
y_{p1}   \\
\vdots   \\
y_{pn_p} \\
\end{pmatrix}, \quad
X := 
\begin{pmatrix}
1       & 0       &        & 0      & x_{11}   & 0        &        & 0        \\
\vdots  & \vdots  & \cdots & \vdots & \vdots   & \vdots   & \cdots & \vdots   \\
1       & 0       &        & 0      & x_{1n_1} & 0        &        & 0        \\
1       & 1       &        & 0      & x_{21}   & x_{21}   &        & 0        \\
\vdots  & \vdots  & \cdots & \vdots & \vdots   & \vdots   & \cdots & \vdots   \\
1       & 1       &        & 0      & x_{2n_2} & x_{2n_2} &        & 0        \\
        &         &        &        &          &          &        &          \\
\vdots  & \vdots  &        & \cdots & \vdots   & \vdots   &        & \cdots   \\
        &         &        &        &          &          &        &          \\
1       & 0       &        & 1      & x_{p1}   & 0        &        & x_{p1}   \\
\vdots  & \vdots  & \cdots & \vdots & \vdots   & \vdots   & \cdots & \vdots   \\
1       & 0       &        & 1      & x_{pn_p} & 0        &        & x_{pn_p} \\
\end{pmatrix}
\in \mathbb{R}^{n \times (2p)}, \quad
\beta := 
\begin{pmatrix}
\mu_0    \\
\alpha_2 \\
\vdots   \\
\alpha_p \\
\beta_0  \\
\gamma_2 \\
\vdots   \\
\gamma_p \\
\end{pmatrix} 
\in
\mathbb{R}^{2p}
\quad \mbox{und} \quad
\sigma^2 > 0.
\end{equation}

\end{theorem}

\footnotesize
\vspace{-1mm}
Bemerkungen
\vspace{-1mm}

* Das Theorem ergibt sich direkt mit den Regeln der Matrixmultiplikation.


<!-- Designmatrix des Modells -->
# Kovarianzanalyse mit Interaktion

Designmatrix der einfaktoriellen Kovarianzanalyse mit Interaktion ($n = 18$, $p = 3$, $2p = 6$)

```{r, echo = F, eval = F}
# Designmatrixerzeugung
set.seed(3)                                 # reproduzierbare Ergebnisse
e     = 0.5                                 # Wert für Null in Visualisierung
de    = 0.1                                 # vom Wert zu korrigierende Abweichung
n     = 18                                  # Anzahl von Datenpunkten
p     = 3                                   # Anzahl von Faktorleveln
m     = round(n/p)                          # Anzahl Datenpunkte pro Gruppe
Xt    = cbind(matrix(1, nrow = n, ncol = 1),# Varianzanalyse-Designmatrix
              kronecker(matrix(e, nrow = p, ncol = p) + (1-e)*diag(p),
                        matrix(1, nrow = m, ncol = 1)))
X1    = Xt[,-2]                             # eliminiere zweite Spalte
X2    = X1                                  # dupliziere Designmatrix
x     = runif(n, -1, +1)                    # Werte der Kovariate
x[x>e-de & x<e] = x[x>e-de & x<e] - de
x[x>e & x<e+de] = x[x>e & x<e+de] + de
for (i in 1:n) { X2[i,X2[i,]==1] = x[i] }   # Interaktions-Designmatrix
X     = cbind(X1, X2)                       # n x p Designmatrix
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(-1, 1, length.out = 257),
    key       = NULL,
    axis.col  = NULL,
    axis.row  = NULL,
    asp       = 1)

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

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


<!-- Anwendungsbeispiel -->
# Kovarianzanalyse mit Interaktion

\normalsize
\setstretch{2}
\textcolor{darkblue}{Anwendungsbeispiel}

experimenteller Faktor/Studiengruppe ($p := 2, \; i = 1,2$)

* Face-to-face- vs. Online-Psychotherapie

parametrischer Regressor/Kovariate ($x_{ij}, \; j = 1,...,n_i$)  

* Digitalaffinität der Studienteilnehmer

abhängige Variable/primäres Ergebnismaß ($y_{ij}, \; j = 1,...,n_i$)

* Prä-/Post-Differenzwerte von BDI-Scores

\vfill


<!-- Beispieldatensatz -->
# Kovarianzanalyse mit Interaktion

\vspace{2mm}
\textcolor{darkblue}{Beispieldatensatz} ($n = 32$; hier: $i = 1,2$, $j = 1,...,10$)

\tiny
\setstretch{1}
```{r, echo = F}
library(MASS)                                                                   # multivariate Normalverteilung
set.seed(1)                                                                     # reproduzierbare Ergebnisse
n_i     = c(16,16)                                                              # Anzahl Patient:innen pro Gruppe
n       = sum(n_i)                                                              # Gesamtzahl Datenpunkte
x       = round(rnorm(n,5,2), digits = 1)                                       # Digitalaffinitätskovariate
p       = 4                                                                     # Anzahl Betaparameter
beta    = matrix(c(5,-3,0,1), nrow = p)                                         # wahre,  aber unbekannte  Betaparameter
sigsqr  = 4                                                                     # wahrer, aber unbekannter Varianzparameter 
X       = matrix(c(rep(1,n_i[1]), rep(1,n_i[2]),                                # \mu_0    Regressor
                   rep(0,n_i[1]), rep(1,n_i[2]),                                # \alpha_2 Regressor
                   x,                                                           # \beta_2  Regressor 
                   x*(c(rep(0,n_i[1]), rep(1,n_i[2])))), nrow = n)              # \gamma_2 Regressor                                
y       = mvrnorm(1, X %*% beta, sigsqr*diag(n))                                # Realisierung der Daten
y       = round(y)                                                              # diskrete BDI-Werte
D       = data.frame(Therapie  = c(rep("F2F", n_i[1]),                          # Studiengruppe F2F
                                   rep("ONL", n_i[2])),                         # Studiengruppe ONL
                     Dig.Affin = x,                                             # Digitalaffinität bei Therapie
                     dBDI      = y)                                             # Delta-BDI-Werte
write.csv(D, file = "Daten/Kovarianzanalyse_B.csv")                             # Datenspeicherung
``` 

\footnotesize
\setstretch{1}
```{r, echo = F}
fname      = "Daten/Kovarianzanalyse_B.csv"
D          = read.table(fname, sep = ",", header = TRUE)
knitr::kable(D[c(1:10,(n_i[1]+1):(n_i[1]+10)),2:4], "pipe", align = "c")        # erste 10 Patient:innen pro Gruppe
```


<!-- Visualisierung des Beispieldatensatzes -->
# Kovarianzanalyse mit Interaktion

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

```{r, echo = F, eval = F}
# Laden des Datensatzes
fname      = "Daten/Kovarianzanalyse_B.csv"
D          = read.table(fname, sep = ",", header = TRUE)

# Deskriptivstatistikanalyse
thp        = c("F2F", "ONL")                                                    # Therapiebedingungen
y_thp      = list()                                                             # therapiespezifische dBDI-Werte
x_thp      = list()                                                             # therapiespezifische Affinitäts-Werte
for(i in 1:length(thp)){                                                        # Iteration über Therapiebedingungen
    y_thp[[i]]    = D$dBDI[D$Therapie == thp[i]]                                # therapiespezifische dBDI-Werte
    x_thp[[i]]    = D$Dig.Affin[D$Therapie == thp[i]]                           # therapiespezifische Affinitäts-Werte
}

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

# Punktwolken
plot(x_thp[[1]], y_thp[[1]],
    pch         = 16,
    xlab        = "Digitalaffinität",
    ylab        = "dBDI",
    xlim        = c(-1,11),
    ylim        = c(0,12))
points(x_thp[[2]], y_thp[[2]],
    pch         = 16,
    col         = "gray80")
legend("topleft", c("F2F", "ONL"),
    lty         = 0,
    pch         = 16,
    col         = c("black", "gray80"),
    bty         = "n",
    cex         = 1,
    x.intersp   = 1)

# Balkendiagramm                        
groupmeans = c(mean(y_thp[[1]]), mean(y_thp[[2]]))
groupstds  = c(sd(y_thp[[1]]), sd(y_thp[[2]]))
names(groupmeans) = c("F2F", "ONL")
x = barplot(groupmeans,                                      
    ylim        = c(0,12),
    xlim        = c(-.4,3),
    col         = c("black", "gray90"),
    ylab        = "dBDI",
    xlab        = "Therapie")
arrows(
    x0          = x,                                    
    y0          = groupmeans - groupstds,               
    x1          = x,                                    
    y1          = groupmeans + groupstds,               
    code        = 3,                                   
    angle       = 90,                                  
    length      = 0.05)   

# Speicherung
dev.copy2pdf(
    file        = "Abbildungen/kovarianzanalyse_B1.pdf",
    width       = 9,
    height      = 4.5)
```

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


<!-- Modellformulierung -->
# Kovarianzanalyse mit Interaktion

\textcolor{darkblue}{Modellformulierung}

\normalsize
Modell 1 | additive Berücksichtigung der Kovariate
\vspace{-2mm}

\footnotesize
\begin{align}
&y_{ij} \sim N(\mu_{ij}, \sigma^2)
\quad \mbox{mit} \quad
\mu_{ij} := \mu_0 + \alpha_i + \beta_0 x_{ij}
\quad \mbox{und} \quad
\alpha_1 := 0
\quad \mbox{und} \quad
\sigma^2 > 0 \\ 
\Leftrightarrow \quad
&y \sim N(X\beta,\sigma^2 I_n)
\quad \mbox{mit} \quad
X := 
\begin{pmatrix}
1      & 0      & x_{11}   \\
\vdots & \vdots & \vdots   \\
1      & 0      & x_{1n_1} \\
1      & 1      & x_{21}   \\
\vdots & \vdots & \vdots   \\
1      & 1      & x_{2n_2} \\
\end{pmatrix}, \quad
\beta := 
\begin{pmatrix}
\mu_0    \\
\alpha_2 \\
\beta_0  \\
\end{pmatrix}
\quad \mbox{und} \quad
\sigma^2 > 0
\end{align}

\normalsize
Modell 2 | additive und interaktive Berücksichtigung der Kovariate
\footnotesize

\begin{align}
&y_{ij} \sim N(\mu_{ij}, \sigma^2)
\quad \mbox{mit} \quad
\mu_{ij} := \mu_0 + \alpha_i + (\beta_0 + \gamma_i) x_{ij}
\quad \mbox{und} \quad
\alpha_1 := \gamma_1 := 0
\quad \mbox{und} \quad
\sigma^2 > 0 \\ 
\Leftrightarrow \quad
&y \sim N(X\beta,\sigma^2 I_n)
\quad \mbox{mit} \quad
X := 
\begin{pmatrix}
1      & 0      & x_{11}   & 0        \\
\vdots & \vdots & \vdots   & \vdots   \\
1      & 0      & x_{1n_1} & 0        \\
1      & 1      & x_{21}   & x_{21}   \\
\vdots & \vdots & \vdots   & \vdots   \\
1      & 1      & x_{2n_2} & x_{2n_2} \\
\end{pmatrix}, \quad
\beta := 
\begin{pmatrix}
\mu_0    \\
\alpha_2 \\
\beta_0  \\
\gamma_2 \\
\end{pmatrix} 
\quad \mbox{und} \quad
\sigma^2 > 0
\end{align}


<!-- Implementation -->
# Kovarianzanalyse mit Interaktion

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

\tiny
\setstretch{1}
```{r, echo = F}
fname     = "Daten/Kovarianzanalyse_B.csv"                            # Dateiname
D         = read.table(fname, sep = ",", header = TRUE)               # Datensatz
n_i       = c(sum(D$Therapie == "F2F"), sum(D$Therapie == "ONL"))     # Anzahl Datenpunkte pro Gruppes
y         = D$dBDI                                                    # Daten
n         = length(y)                                                 # Gesamtanzahl Datenpunkte
XS        = list()                                                    # Model 1 und 2 Liste
CS        = list()                                                    # Kontrastgewichtsvektorenliste
XS[[1]]   = matrix(c(rep(1,n_i[1]), rep(1,n_i[2]),                    # \mu_0    Regressor
                     rep(0,n_i[1]), rep(1,n_i[2]),                    # \alpha_2 Regressor
                     D$Dig.Affin), nrow = n)                          # \beta_2  Regressor
XS[[2]]   = matrix(c(rep(1,n_i[1]), rep(1,n_i[2]),                    # \mu_0    Regressor
                     rep(0,n_i[1]), rep(1,n_i[2]),                    # \alpha_2 Regressor
                     D$Dig.Affin,                                     # \beta_2  Regressor
                     D$Dig.Affin*(c(rep(0,n_i[1]),rep(1,n_i[2])))),   # \gamma_2 Regressor
                   nrow = n)                                          # Anzahl Zeilen der Designmatrix
CS[[1]]   = matrix(c(0,1,0),   nrow = 3)                              # Kontrastgewichtsvektor \alpha_2
CS[[2]]   = matrix(c(0,1,0,0), nrow = 4)                              # Kontrastgewichtsvektor \alpha_2
B         = list()                                                    # Betaparameterschätzerliste
S         = rep(NaN,2)                                                # Varianzparameterschätzervektor
T         = rep(NaN,2)                                                # T-Statistik-Vektor
for(i in 1:2){                                                        # Iteration über Modelle
    X           = XS[[i]]                                             # Designmatrix
    p           = ncol(X)                                             # Anzahl Betaparameter
    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
    c           = CS[[i]]                                             # Kontrastgewichtsvektor
    t_num       = t(c) %*% beta_hat                                   # Zähler der Zweistichproben-T-Teststatistik
    t_den       = sqrt(sigsqr_hat %*% t(c) %*% solve(t(X) %*% X) %*%c)# Nenner der Zweistichproben-T-Teststatistik
    t           = t_num/t_den                                         # Wert der   Zweistichproben-T-Teststatistik
    B[[i]]      = beta_hat                                            # Betaparameterschätzer
    S[i]        = sigsqr_hat                                          # Varianzparameterschätzer
    T[i]        = t                                                   # T-Statistik
}
```

```{r, echo = F, out.width = "90%", fig.align = "center"}
knitr::include_graphics("Abbildungen/kovarianzanalyse_interaktion.png")
```


<!-- Ergebnisse -->
# Kovarianzanalyse mit Interaktion

\textcolor{darkblue}{Ergebnisse}
\small

```{r, echo = F}
cat(  "Modell 1 | Betaparameterschätzer    : ", round(B[[1]],2),
    "\nModell 2 | Betaparameterschätzer    : ", round(B[[2]],2),
    "\nModell 1 | Varianzparameterschätzer : ", round(S[[1]],2),
    "\nModell 2 | Varianzparameterschätzer : ", round(S[[2]],2),
    "\nModell 1 | T-Statistik ONL-Therapie : ", round(T[[1]],2),
    "\nModell 2 | T-Statistik ONL-Therapie : ", round(T[[2]],2))
```

\small
In Modell 1 wird der ONL-Therapieeffekt positiv, in Modell 2 negativ geschätzt.

Entsprechend ist die T-Statistik für ONL-Therapie in Modell 1 positiv, in Modell 2 negativ.

Die wahren, aber unbekannten Betaparameterwerte sind $\beta := (5,-3,0,1)^\mathrm{T}$.

$\Rightarrow$ Die Abhängigkeit von Digitalaffinität kommt bei ONL zum Tragen, bei F2F nicht.

$\Rightarrow$ Der Haupteffekt der Therapieart ist schwierig zu beurteilen.

$\Rightarrow$ ONL-Therapie ist bei niedriger Digitalaffinität nicht so wirksam wie F2F-Therapie.

$\Rightarrow$ ONL-Therapie ist bei hoher Digitalaffinität wirksamer als F2F-Therapie.


<!-- Visualisierung -->
# Kovarianzanalyse mit Interaktion

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

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

# Punktwolken
plot(D$Dig.Affin[D$Therapie == "F2F"], D$dBDI[D$Therapie == "F2F"],
    pch        = 16,
    xlab       = "Digitalaffinität",
    ylab       = "dBDI",
    xlim       = c(-1,11),
    ylim       = c(0,12))
points(D$Dig.Affin[D$Therapie == "ONL"], D$dBDI[D$Therapie == "ONL"],
    pch        = 16,
    col        = "gray80")

# prädizierte Daten
X              = XS[[2]]
beta_hat       = B[[2]]
y_hat          = X %*% beta_hat
lines(D$Dig.Affin[D$Therapie == "F2F"], y_hat[1:n_i[1]],
    col        = "black")
lines(D$Dig.Affin[D$Therapie == "ONL"], y_hat[(n_i[1]+1):n],
    col        = "gray80")
legend("topleft", c("F2F", "ONL"),
    lty        = 0,
    pch        = 16,
    col        = c("black", "gray80"),
    bty        = "n",
    cex        = 1,
    x.intersp  = 1)

# adjustierte Daten
y_adj       = y - X %*% matrix(c(0, 0, beta_hat[3], beta_hat[4]), nrow = 4)
groupmeans  = c(mean(y_adj[1:n_i[1]]), mean(y_adj[(n_i[1]+1):n]))
groupstds   = c(sd(y_adj[1:n_i[1]]), sd(y_adj[(n_i[1]+1):n]))
names(groupmeans) = c("F2F", "ONL")
x = barplot(groupmeans,                                      
    ylim        = c(0,12),
    xlim        = c(-.4,3),
    col         = c("black", "gray90"),
    ylab        = "dBDI",
    xlab        = "Therapie",
    main        = "Adjustierte Daten")
arrows(                                             
    x0          = x,                                    
    y0          = groupmeans - groupstds,               
    x1          = x,                                    
    y1          = groupmeans + groupstds,               
    code        = 3,                                   
    angle       = 90,                                  
    length      = 0.05)   

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

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


#
<!-- Abschnitt 4: Anwendung/Praxis -->
\setcounter{framenumber}{32}

\large
\setstretch{3}
\vfill

Vorbemerkungen

Additive Kovarianzanalyse

Kovarianzanalyse mit Interaktion

**Anwendung/Praxis**

Selbstkontrollfragen

\vfill


<!-- Datensatz: Geburtsgewichte -->
# Anwendung/Praxis

\small
Datensatz: Geburtsgewicht und Muttergewicht für Mutter-Kind-Paare ($n = 189$; hier: $i = 1,...,25$)

\setstretch{1.0}

\footnotesize
\vspace{2mm}
```{r, echo = T}
filename = "Daten/birth_weights.csv"        # Dateiname
D        = read.csv(filename)               # Datenframe
```

\tiny
\vspace{-2mm}
```{r, echo = F}
knitr::kable(D[c(1:25),],
             align = "rrrcccc",
             "pipe")
```


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

\small
Datensatz: Geburtsgewicht und Muttergewicht, getrennt nach Raucherstatus ($n = 189$)

```{r, echo = F, eval = F}
# Daten extrahieren
filename = "Daten/birth_weights.csv"                # Dateiname
D        = read.csv(filename)                       # Datenframe
y        = D$Birth_weight                           # abhängige Variable
f        = D$Smoker                                 # kategorische UV
x        = D$Weight                                 # parametrische UV

# Modellschätzung
y_1      = y[f=="no"]                               # Datenpunkte Gr. 1
y_2      = y[f=="yes"]                              # Datenpunkte Gr. 2
x_1      = x[f=="no"]                               # Kovariatenwerte Gr. 1
x_2      = x[f=="yes"]                              # Kovariatenwerte Gr. 2
n_1      = length(y_1)                              # Anzahl Datenpunkte Gr. 1
n_2      = length(y_2)                              # Anzahl Datenpunkte Gr. 2
n        = n_1 + n_2                                # Anzahl Datenpunkte
p        = 4                                        # Anzahl Regressoren
y        = matrix(c(y_1, y_2), nrow = n)            # Datenvektor
x        = matrix(c(x_1, x_2), nrow = n)            # Kovariatenvektor
X        = matrix(c(rep(1,n_1), rep(1,n_2),         # \mu_0    Regressor
                    rep(0,n_1), rep(1,n_2),         # \alpha_2 Regressor
                    x,                              # \beta_0  Regressor
                    x*(c(rep(0,n_1), rep(1,n_2)))), # \gamma_2 Regressor
                  nrow = n)                         # Anzahl Zeilen der Designmatrix
beta_hat = solve(t(X) %*% X) %*% t(X) %*% y         # Betaparameterschätzer vollständiges Modell

# Visualisierung
library(latex2exp)
par(
    family     = "sans",
    mfcol      = c(1,2),
    pty        = "s",
    bty        = "l",
    lwd        = 1,
    las        = 1,
    mgp        = c(2,1,0),
    xaxs       = "i",
    yaxs       = "i",
    font.main  = 1,
    cex        = 1,
    cex.main   = 1)

# Punktwolken
plot(x_1, y_1,
    pch        = 16,
    xlab       = "Muttergewicht",
    ylab       = "Geburtsgewicht",
    xlim       = c(35,115),
    ylim       = c(1,5))
points(x_2, y_2,
    pch        = 16,
    col        = "gray80")

# prädizierte Daten
y_hat          = X %*% beta_hat
lines(x_1, y_hat[1:n_1],
    col        = "black")
lines(x_2, y_hat[(n_1+1):n],
    col        = "gray80")
legend("topleft", c("Nichtraucherin", "Raucherin"),
    lty        = 0,
    pch        = 16,
    col        = c("black", "gray80"),
    bty        = "n",
    cex        = 1,
    x.intersp  = 1)

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

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


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

\small
Zweistichproben-T-Test: Geburtsgewicht, Nichtraucherinnen vs. Raucherinnen ($n_1 = 115$, $n_2 = 74$)

\vspace{2mm}
```{r, echo = T, eval = T}
# Zweistichproben-T-Test
y1     = D$Birth_weight[D$Smoker=="no"]                     # Gruppe 1
y2     = D$Birth_weight[D$Smoker=="yes"]                    # 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")
```


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

\small
Kovarianzanalyse: T-Test auf Haupteffekt von Raucherstatus

\vspace{2mm}
\footnotesize
```{r, echo = T, eval = T}
# Kovarianzanalyse
D$Weight    = D$Weight - mean(D$Weight)                    # Kovariate mittelwertzentrieren
y_1         = D$Birth_weight[D$Smoker=="no"]               # Datenpunkte Gr. 1
y_2         = D$Birth_weight[D$Smoker=="yes"]              # Datenpunkte Gr. 2
x_1         = D$Weight[D$Smoker=="no"]                     # Kovariatenwerte Gr. 1
x_2         = D$Weight[D$Smoker=="yes"]                    # Kovariatenwerte Gr. 2
n_1         = length(y_1)                                  # Anzahl Datenpunkte Gr. 1
n_2         = length(y_2)                                  # Anzahl Datenpunkte Gr. 2
n           = n_1 + n_2                                    # Anzahl Datenpunkte
p           = 4                                            # Anzahl Regressoren
y           = matrix(c(y_1, y_2), nrow = n)                # Datenvektor
x           = matrix(c(x_1, x_2), nrow = n)                # Kovariatenvektor
X           = matrix(c(rep(1,n_1), rep(1,n_2),             # \mu_0    Regressor
                       rep(0,n_1), rep(1,n_2),             # \alpha_2 Regressor
                       x,                                  # \beta_0  Regressor
                       x*(c(rep(0,n_1), rep(1,n_2)))),     # \gamma_2 Regressor
                     nrow = n)                             # Anzahl Zeilen der Designmatrix
c           = matrix(c(0,1,0,0), nrow = p)                 # Kontrastgewichtsvektor
```

```{r, echo = F, eval = T}
# Evaluation der T-Teststatistik
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
c           = matrix(c(0,1,0,0), nrow = p)                     # Kontrastgewichtsvektor
t_num       = t(c) %*% beta_hat                                # Zähler der T-Statistik
ste         = sqrt(sigsqr_hat*t(c) %*% solve(t(X) %*%X ) %*% c)# Nenner der T-Statistik
tee         = t_num/ste                                        # T-Statistik
pval        = 2*(1 - pt(abs(tee),n-p))                         # p-Wert

# Ausgabe der Ergebnisse
R           = data.frame(t_num, ste, tee, pval)
rownames(R) = c("Smoking")
colnames(R) = c("Estimate", "Std. Error", "t-value", "Pr(>|t|)")
print(R)
```


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

\small
Kovarianzanalyse: F-Test auf Effekt der Kovariate Muttergewicht

\vspace{2mm}
\footnotesize
```{r, echo = T, eval = T}
# Kovarianzanalyse
D$Weight    = D$Weight - mean(D$Weight)                    # Kovariate mittelwertzentrieren
y_1         = D$Birth_weight[D$Smoker=="no"]               # Datenpunkte Gr. 1
y_2         = D$Birth_weight[D$Smoker=="yes"]              # Datenpunkte Gr. 2
x_1         = D$Weight[D$Smoker=="no"]                     # Kovariatenwerte Gr. 1
x_2         = D$Weight[D$Smoker=="yes"]                    # Kovariatenwerte Gr. 2
n_1         = length(y_1)                                  # Anzahl Datenpunkte Gr. 1
n_2         = length(y_2)                                  # Anzahl Datenpunkte Gr. 2
n           = n_1 + n_2                                    # Anzahl Datenpunkte
p           = 4                                            # Anzahl Regressoren
y           = matrix(c(y_1, y_2), nrow = n)                # Datenvektor
x           = matrix(c(x_1, x_2), nrow = n)                # Kovariatenvektor
X           = matrix(c(rep(1,n_1), rep(1,n_2),             # \mu_0    Regressor
                       rep(0,n_1), rep(1,n_2),             # \alpha_2 Regressor
                       x,                                  # \beta_0  Regressor
                       x*(c(rep(0,n_1), rep(1,n_2)))),     # \gamma_2 Regressor
                     nrow = n)                             # Anzahl Zeilen der Designmatrix
p_0         = 2                                            # Anzahl Parameter reduziertes Modell
```

```{r, echo = F, eval = T}
# Evaluation der F-Teststatistik
p_0         = 2                                                # Anzahl Parameter reduziertes Modell
p_1         = p - p_0                                          # Anzahl zusätzlicher Parameter im vollst. Modell
X_0         = X[,1:p_0]                                        # Designmatrix reduzierters 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 Quadratsumme reduziertes Modell
eh_eh       = t(eps_hat)   %*% eps_hat                         # residuelle Quadratsumme 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
pval        = 1 - pf(f,p_1,n-p)                                # p-Wert

# Ausgabe der Ergebnisse
R           = data.frame(p_1, n-p, f, pval)
rownames(R) = c("Weight")
colnames(R) = c("numerator DFs", "denominator DFs", "F-value", "Pr(>F)")
print(R)
```


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

\small
Kovarianzanalyse: T-Test auf Interaktion von Raucherstatus und Muttergewicht

\vspace{2mm}
\footnotesize
```{r, echo = T, eval = T}
# Kovarianzanalyse
D$Weight    = D$Weight - mean(D$Weight)                    # Kovariate mittelwertzentrieren
y_1         = D$Birth_weight[D$Smoker=="no"]               # Datenpunkte Gr. 1
y_2         = D$Birth_weight[D$Smoker=="yes"]              # Datenpunkte Gr. 2
x_1         = D$Weight[D$Smoker=="no"]                     # Kovariatenwerte Gr. 1
x_2         = D$Weight[D$Smoker=="yes"]                    # Kovariatenwerte Gr. 2
n_1         = length(y_1)                                  # Anzahl Datenpunkte Gr. 1
n_2         = length(y_2)                                  # Anzahl Datenpunkte Gr. 2
n           = n_1 + n_2                                    # Anzahl Datenpunkte
p           = 4                                            # Anzahl Regressoren
y           = matrix(c(y_1, y_2), nrow = n)                # Datenvektor
x           = matrix(c(x_1, x_2), nrow = n)                # Kovariatenvektor
X           = matrix(c(rep(1,n_1), rep(1,n_2),             # \mu_0    Regressor
                       rep(0,n_1), rep(1,n_2),             # \alpha_2 Regressor
                       x,                                  # \beta_0  Regressor
                       x*(c(rep(0,n_1), rep(1,n_2)))),     # \gamma_2 Regressor
                     nrow = n)                             # Anzahl Zeilen der Designmatrix
c           = matrix(c(0,0,0,1), nrow = p)                 # Kontrastgewichtsvektor
```

```{r, echo = F, eval = T}
# Evaluation der T-Teststatistik
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
c           = matrix(c(0,0,0,1), nrow = p)                     # Kontrastgewichtsvektor
t_num       = t(c) %*% beta_hat                                # Zähler der T-Statistik
ste         = sqrt(sigsqr_hat*t(c) %*% solve(t(X) %*%X ) %*% c)# Nenner der T-Statistik
tee         = t_num/ste                                        # T-Statistik
pval        = 2*(1 - pt(abs(tee),n-p))                         # p-Wert

# Ausgabe der Ergebnisse
R           = data.frame(t_num, ste, tee, pval)
rownames(R) = c("Smoking x Weight")
colnames(R) = c("Estimate", "Std. Error", "t-value", "Pr(>|t|)")
print(R)
```


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

\small
Datensatz: Geburtsgewicht und Muttergewicht, getrennt nach Raucherstatus ($n = 189$)

```{r, echo = F, eval = F}
# Daten extrahieren
filename = "Daten/birth_weights.csv"                # Dateiname
D        = read.csv(filename)                       # Datenframe
y        = D$Birth_weight                           # abhängige Variable
f        = D$Smoker                                 # kategorische UV
x        = D$Weight                                 # parametrische UV

# Modellschätzung
y_1      = y[f=="no"]                               # Datenpunkte Gr. 1
y_2      = y[f=="yes"]                              # Datenpunkte Gr. 2
x_1      = x[f=="no"]                               # Kovariatenwerte Gr. 1
x_2      = x[f=="yes"]                              # Kovariatenwerte Gr. 2
n_1      = length(y_1)                              # Anzahl Datenpunkte Gr. 1
n_2      = length(y_2)                              # Anzahl Datenpunkte Gr. 2
n        = n_1 + n_2                                # Anzahl Datenpunkte
p        = 4                                        # Anzahl Regressoren
y        = matrix(c(y_1, y_2), nrow = n)            # Datenvektor
x        = matrix(c(x_1, x_2), nrow = n)            # Kovariatenvektor
X        = matrix(c(rep(1,n_1), rep(1,n_2),         # \mu_0    Regressor
                    rep(0,n_1), rep(1,n_2),         # \alpha_2 Regressor
                    x,                              # \beta_0  Regressor
                    x*(c(rep(0,n_1), rep(1,n_2)))), # \gamma_2 Regressor
                  nrow = n)                         # Anzahl Zeilen der Designmatrix
beta_hat = solve(t(X) %*% X) %*% t(X) %*% y         # Betaparameterschätzer vollständiges Modell

# Visualisierung
library(latex2exp)
par(
    family     = "sans",
    mfcol      = c(1,2),
    pty        = "s",
    bty        = "l",
    lwd        = 1,
    las        = 1,
    mgp        = c(2,1,0),
    xaxs       = "i",
    yaxs       = "i",
    font.main  = 1,
    cex        = 1,
    cex.main   = 1)

# Punktwolken
plot(x_1, y_1,
    pch        = 16,
    xlab       = "Muttergewicht",
    ylab       = "Geburtsgewicht",
    xlim       = c(35,115),
    ylim       = c(1,5))
points(x_2, y_2,
    pch        = 16,
    col        = "gray80")

# prädizierte Daten
y_hat          = X %*% beta_hat
lines(x_1, y_hat[1:n_1],
    col        = "black")
lines(x_2, y_hat[(n_1+1):n],
    col        = "gray80")
legend("topleft", c("Nichtraucherin", "Raucherin"),
    lty        = 0,
    pch        = 16,
    col        = c("black", "gray80"),
    bty        = "n",
    cex        = 1,
    x.intersp  = 1)

# adjustierte Daten
y_adj       = y - X %*% matrix(c(0, 0, beta_hat[3], beta_hat[4]), nrow = 4)
groupmeans  = c(mean(y_adj[1:n_1]), mean(y_adj[(n_1+1):n]))
groupstds   = c(sd(y_adj[1:n_1]), sd(y_adj[(n_1+1):n]))
names(groupmeans) = c("Nichtraucherin", "Raucherin")
bp = barplot(groupmeans,
    ylim        = c(0,5),
    xlim        = c(-0.5,3),
    col         = c("black", "gray90"),
    ylab        = "Geburtsgewicht",
    xlab        = "Studiengruppe",
    main        = "adjustierte Daten")
arrows(
    x0          = bp,
    y0          = groupmeans - groupstds,
    x1          = bp,
    y1          = groupmeans + groupstds,
    code        = 3,
    angle       = 90,
    length      = 0.05)

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

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


<!-- Anhang: Selbstkontrollfragen -->
#

\large
\setstretch{3}
\vfill

Vorbemerkungen

Additive Kovarianzanalyse

Kovarianzanalyse mit Interaktion

Anwendung/Praxis

**Selbstkontrollfragen**

\vfill


<!-- Selbstkontrollfragen -->
# Selbstkontrollfragen

\footnotesize
\setstretch{2}

1. Erläutern Sie die Begriffe der Gesamt-/Primär-/Sekundär- und Residualvaribilität.
1. Geben Sie jeweils zwei Beispiele für parametrische bzw. faktorielle ALM-Designs an.
1. Geben Sie die Definition des Modells der additiven einfaktoriellen Kovarianzanalyse (EKA) wieder.
1. Erläutern Sie die Bedeutung der Parameter $\mu_0, \alpha_i$ und $\beta_0$ im Modell der additiven EKA.
1. Geben Sie die Designmatrixform des Modells der additiven EKA für $p := 2$ und $n_1 = n_2 = 3$ an.
1. Erläutern Sie den Unterschied zwischen dem Modell der einfaktoriellen Varianzanalyse und dem Modell der additiven einfaktoriellen Kovarianzanalyse.
1. Geben Sie die Definition des Modells der einfaktoriellen Kovarianzanalyse (EKA) mit Interaktion wieder.
1. Erläutern Sie die Bedeutung der Parameter $\mu_0, \alpha_i, \beta_0$ und $\gamma_i$ im Modell der EKA mit Interaktion.
1. Geben Sie die Designmatrixform des Modells der EKA mit Interaktion für $p := 2$ und $n_1 = n_2 = 3$ an.
1. Erläutern Sie den Unterschied zwischen dem Modell der additiven einfaktoriellen Kovarianzanalyse und dem Modell der einfaktoriellen Kovarianzanalyse mit Interaktion.
