--- fontsize: 8pt format: beamer: include-in-header: ../ALM_Header.tex bibliography: ../ALM_Referenzen.bib --- ```{r, include = F} source("../ALM_Common.R") ``` # {.plain} \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 # {.plain} \vfill \center \huge \textcolor{black}{(10) Einfaktorielle Varianzanalyse} \vfill # \large \setstretch{2.5} \vfill Anwendungsszenario Modellformulierung Modellschätzung Modellevaluation Anwendung/Praxis Selbstkontrollfragen \vfill # \large \setstretch{2.5} \vfill **Anwendungsszenario** Modellformulierung Modellschätzung Modellevaluation Anwendung/Praxis Selbstkontrollfragen \vfill # Anwendungsszenario \textbf{\textcolor{darkblue}{Zwei oder mehr Gruppen}} (Stichproben) randomisierter experimenteller Einheiten. Annahme der unabhängigen und identischen Normalverteilung $N(\mu_i,\sigma^2)$ der Daten. $\mu_i$ und $\sigma^2$ unbekannt. Inferentieller Test der Nullhypothese identischer Gruppenerwartungswerte beabsichtigt. \center \vspace{1mm} $\Downarrow$ \vspace{1mm} Generalisierung des Zweistichproben-T-Tests bei unabhängigen Stichproben \linebreak mit einfacher Nullhypothese für mehr als zwei Stichproben \vfill # Anwendungsszenario \textcolor{darkblue}{Anwendungsbeispiel} ```{r, echo = F, out.width = "70%"} knitr::include_graphics("Abbildungen/beispielstudie.pdf") ``` \small \vspace{2mm} Pre-Post-BDI-Differenzwertanalyse für drei Gruppen von Patient:innen (F2F, ONL, WLC): \vspace{-2mm} * inferentielle Evidenz für Gruppenerwartungswertunterschiede? * empirische Evidenz für unterschiedliche Therapiewirksamkeit? ```{r, echo = F, eval = F} # Initialisierung set.seed(1) # Zufallszahlengenerator initialisieren # Simulationsparameter n_c = 3 # Anzahl Gruppen n_i = 40 # Anzahl Proband:innen pro Gruppe n = n_c*n_i # Gesamtanzahl Datenpunte # Simulationsparameter mu_pre = c(31, 32, 29) # \mu pre F2F, ONL, WLC mu_pos = c(27, 26, 30) # \mu post F2F, ONL, WLC sigsqr = 10 # \sigma^2 # Datensimulation D = data.frame("ID" = 1:n) # Dataframe-Initialisierung und ID-Variable D$Setting = c(rep("F2F",n_i), rep("ONL", n_i), rep("WLC", n_i)) # Bedingung D$PreBDI = c(round(rnorm(n_i, mu_pre[1], sqrt(sigsqr))), # PreBDI round(rnorm(n_i, mu_pre[2], sqrt(sigsqr))), round(rnorm(n_i, mu_pre[3], sqrt(sigsqr)))) D$PostBDI = c(round(rnorm(n_i, mu_pos[1], sqrt(sigsqr))), # PostBDI round(rnorm(n_i, mu_pos[2], sqrt(sigsqr))), round(rnorm(n_i, mu_pos[3], sqrt(sigsqr)))) D$dBDI = -(D$PostBDI - D$PreBDI) # -(PostBDI - PreBDI) = PreBDI - PostBDI # Datenspeicherung write.csv(D, file = "Daten/Einfaktorielle_Varianzanalyse_Daten.csv") ``` # \large \setstretch{2.5} \vfill Anwendungsszenario **Modellformulierung** Modellschätzung Modellevaluation Anwendung/Praxis Selbstkontrollfragen \vfill # Modellformulierung \footnotesize \begin{definition}[EVA-Modell in Erwartungswertparameterdarstellung] \justifying $y_{ij}$ mit $i= 1,...,p$, welche die Gruppen indizieren, und $j = 1,...,n_i$, welche die experimentellen Einheiten innerhalb der Gruppen indizieren, seien Zufallsvariablen, die die Datenpunkte eines Einfaktoriellen-Varianzanalyse-Szenarios (EVA) modellieren. Dann hat das \textit{EVA-Modell in Erwartungswertparameterdarstellung} die strukturelle Form \begin{equation} y_{ij} = \mu_i + \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 \quad \mbox{mit} \quad \mu_i \in \mathbb{R}, \; \sigma^2 > 0 \end{equation} und die Datenverteilungsform \begin{equation} \label{eq:eva_klassisch} y_{ij} \sim N(\mu_i,\sigma^2) \quad \mbox{u.v. für} \quad i = 1,...,p, \; j = 1,...,n_i \quad \mbox{mit} \quad \mu_i \in \mathbb{R}, \; \sigma^2 > 0. \end{equation} Wenn $n_i := m$ für alle $i = 1,...,p$, heißt das EVA-Szenario \textit{balanciert}. \end{definition} Bemerkungen * Die Äquivalenz der beiden Modellformulierungen folgt aus den Ergebnissen zu Transformationen der Normalverteilung (vgl. Einheit (7) in *Wahrscheinlichkeitstheorie und Frequentistische Inferenz* und Einheit (4) in *Allgemeines Lineares Modell*). * Es handelt sich um die Generalisierung des Zweistichproben-T-Test-Modells für unabhängige Stichproben unter Annahme identischer Varianzparameter von $p = 2$ auf ein beliebiges $p \in \mathbb{N}$. * Bei balancierten Varianzanalyseszenarien besteht jede Datengruppe aus der gleichen Anzahl von Datenpunkten. # Modellformulierung \textcolor{darkblue}{Motivation der Effektdarstellung} \footnotesize Die Erwartungswertparameterdarstellung des EVA Modells ist ein valides ALM, das sich in dieser Form auch in der Literatur findet (z.B. @georgii2009, Kapitel 12.4). Im Sinne der Konsistenz mit mehrfaktoriellen Varianzanalysemodellen bietet sich jedoch eine Reparameterisierung der Erwartungswertparameter an. Kern dieser Reparameterisierung ist es, den Erwartungswertparameter der $i$ten Gruppe als Summe eines *gruppenübergreifenden Erwartungswertparameters* $\mu_0 \in \mathbb{R}$ und eines *gruppenspezifischen Effektparameters* $\alpha_i \in \mathbb{R}$ zu modellieren, \begin{equation} \label{eq:uber} \mu_i := \mu_0 + \alpha_i \quad \mbox{für} \quad i = 1,...,p. \end{equation} \vspace{-2mm} Dabei modelliert $\alpha_i$ den Unterschied (die Differenz) zwischen dem $i$ten Erwartungswertparameter $\mu_i$ und dem gruppenübergreifenden Erwartungswertparameter $\mu_0$, \begin{equation} \alpha_i = \mu_i - \mu_0 \quad \mbox{für} \quad i = 1,...,p. \end{equation} \vspace{-2mm} Allerdings hat die in dieser Form vorgenommene Reparameterisierung einen entscheidenen Nachteil: es werden $p$ Erwartungswertparameter $\mu_i, i = 1,...,p$ durch die $p + 1$ Parameter $\mu_0$ und $\alpha_i, i = 1,...,p$ dargestellt. Diese Darstellung ist im Allgemeinen nicht eindeutig. Zum Beispiel können die Erwartungswertparameter \begin{equation} \mu_1 = 3, \; \mu_2 = 5, \; \mu_3 = 6 \end{equation} sowohl durch \begin{equation} \mu_0 = 0 \quad \mbox{und} \quad \alpha_1 = 3, \; \alpha_2 = 5, \; \alpha_3 = 6 \end{equation} als auch durch \begin{equation} \mu_0 = 1 \quad \mbox{und} \quad \alpha_1 = 2, \; \alpha_2 = 4, \; \alpha_3 = 5 \end{equation} dargestellt werden. Man sagt in diesem Kontext auch, dass das EVA-Modell mit \eqref{eq:uber} *überparameterisiert* ist. # Modellformulierung \textcolor{darkblue}{Motivation der Effektdarstellung} \footnotesize Datenanalytisch hat die Überparameterisierung eines Varianzanalysemodells den Nachteil, dass aus $p$ geschätzten Erwartungswertparametern $p + 1$ Betaparameterschätzer bestimmt werden müssten, was wie gesehen nicht eindeutig erfolgen kann. Um diese Probleme in der Effektparameterdarstellung des EVA-Modells zu umgehen und sie konsistent auf mehrfaktorielle Varianzanalysemodelle zu übertragen, bietet sich die Einführung einer Nebenbedingung an: \begin{equation} \alpha_1 := 0 \; . \end{equation} Es wird also ein Effektparameter von vornherein als "gleich Null" angenommen. Für die gruppenspezifischen Erwartungswertparameter ergibt sich damit \begin{align} \begin{split} \mu_1 & := \mu_0 \\ \mu_i & := \mu_0 + \alpha_i \quad \mbox{für} \quad i = 2,...,p. \end{split} \end{align} Hierbei wird die erste Gruppe nun als *Referenzgruppe* bezeichnet und die $\alpha_i$ modellieren die Differenz zwischen dem Erwartungswertparameter der $i$ten Gruppe und dem Erwartungswertparameter der ersten Gruppe: \begin{equation} \alpha_i = \mu_i - \mu_0 = \mu_i - \mu_1 \quad \mbox{für} \quad i = 1,...,p. \end{equation} $\mu_0$ ist also kein gruppenübergreifender Erwartungswertparameter mehr, sondern identisch mit dem Erwartungswertparameter der ersten Gruppe. Welche tatsächliche experimentelle Gruppe dabei als "erste Gruppe" definiert wird, ist unerheblich. Entscheidend ist, dass die entsprechenden Erwartungswertparameterschätzer $\hat{\mu}_0$,$\hat{\alpha}_2, ..., \hat{\alpha}_p$ korrekt als (1) Erwartungswertparameterschätzer der Referenzgruppe ($\hat{\mu}_0$) und (2) geschätzte Differenz zwischen dem Erwartungswertparameter der Referenzgruppe und dem Erwartungswertparameter der $i$ten Gruppe ($\hat{\alpha}_2, ..., \hat{\alpha}_p$) verstanden werden. Wir formalisieren das oben Gesagte in folgendem Theorem. # Modellformulierung \footnotesize \begin{theorem}[EVA-Modell in Effektdarstellung mit Referenzgruppe I] \justifying \normalfont Gegeben sei das EVA-Modell in Erwartungswertparameterdarstellung. Dann können die Zufallsvariablen, die die Datenpunkte des EVA-Szenarios modellieren, äquivalent in der strukturellen Form \begin{align} \label{eq:eva_effekt_1} \begin{split} y_{1j} & = \mu_0 + \varepsilon_{1j} \hphantom{+ \alpha_i} \quad \mbox{mit} \quad \varepsilon_{1j} \sim N(0,\sigma^2) \quad \mbox{u.i.v. für} \quad j = 1,...,n_1 \\ y_{ij} & = \mu_0 + \alpha_i + \varepsilon_{ij} \quad \mbox{mit} \quad \varepsilon_{ij} \sim N(0,\sigma^2) \quad \mbox{u.i.v. für} \quad i = 2,...,p, \; j = 1,...,n_i \end{split} \end{align} mit $\alpha_i := \mu_i - \mu_1$ für $i = 2,...,p$ und in der entsprechenden Datenverteilungsform \begin{align}\label{eq:eva_effekt_2} \begin{split} y_{1j} & \sim N(\mu_0,\sigma^2) \hphantom{+ \alpha_i} \quad \mbox{u.i.v. für} \quad j = 1,...,n_i \quad \mbox{mit} \quad \mu_0 \in \mathbb{R}, \; \sigma^2 > 0 \\ y_{ij} & \sim N(\mu_0 + \alpha_i,\sigma^2) \quad \mbox{u.i.v. für} \quad i = 2,..., p, \; j = 1,...,n_i \quad \mbox{mit} \quad \alpha_i \in \mathbb{R}, \; \sigma^2 > 0 \\ \end{split} \end{align} geschrieben werden. Wir nennen \eqref{eq:eva_effekt_1} und \eqref{eq:eva_effekt_2} die strukturelle und die Datenverteilungsform des \textit{EVA-Modells in Effektdarstellung mit Referenzgruppe}. \end{theorem} \footnotesize \underline{Beweis} Es gilt \begin{equation} \mu_i = \mu_0 + \mu_i - \mu_0 \; . \end{equation} Die Parameterisierungen mit $\mu_i$ und mit $\mu_0 + \mu_i - \mu_0$ sind also gleich und damit äquivalent. Dann folgt aber auch \begin{equation} \mu_i = \mu_0 + (\mu_i - \mu_0) =: \mu_0 + \alpha_i \quad \mbox{für} \quad i = 1,...,p \; . \end{equation} Mit $\alpha_1 := 0$ gilt dann $\mu_1 = \mu_0$ und $\mu_i = \mu_0 + \alpha_i$ für $i = 2,...,p$, wie im Theorem behauptet. $\hfill\Box$ # Modellformulierung \vspace{2mm} \footnotesize \begin{theorem}[EVA-Modell in Effektdarstellung mit Referenzgruppe II] \justifying \normalfont Gegeben sei die strukturelle Form des EVA-Modells in Effektdarstellung 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 $n := \sum_{i=1}^p n_i$ und \begin{equation*} \renewcommand{\arraystretch}{0.8} y := \begin{pmatrix*}[c] y_{11} \\ \vdots \\ y_{1n_1} \\ y_{21} \\ \vdots \\ y_{2n_2} \\ \\ \vdots \\ \\ y_{p1} \\ \vdots \\ y_{pn_p} \end{pmatrix*} \in \mathbb{R}^n, \quad X := \begin{pmatrix} 1 & 0 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 0 \\ 1 & 1 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 1 & & 0 \\ & & & \\ \vdots & \vdots & \cdots & \vdots \\ & & & \\ 1 & 0 & & 1 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 1 \\ \end{pmatrix} \in \mathbb{R}^{n \times p}, \quad \beta := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \vdots \\ \alpha_p \end{pmatrix} \in \mathbb{R}^{p} \quad \mbox{und} \quad \sigma^2 > 0. \end{equation*} \end{theorem} \vspace{-2mm} \underline{Beweis} Das Theorem ergibt sich direkt mit den Regeln der Matrixmultiplikation. $\hfill\Box$ # Modellformulierung Designmatrix des EVA-Modells in Effektdarstellung mit Referenzgruppe ($n = 12$, $p = 4$) ```{r, echo = F, eval = F} # Designmatrixerzeugung m = 3 # Anzahl Datenpunkte i-te Gruppe p = 4 # Anzahl Gruppen n = p*m # Gesamtanzahl Datenpunkte Xt = cbind(matrix(1, nrow = n, ncol = 1), # n x p Designmatrix kronecker(diag(p), matrix(1, nrow = m, ncol = 1))) X = Xt[,-2] # eliminiere zweite 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/eva_designmatrix.pdf", width = 10, height = 10) ``` \vspace{2mm} ```{r, echo = F, out.width = "60%"} knitr::include_graphics("Abbildungen/eva_designmatrix.pdf") ``` # Modellformulierung \textcolor{darkblue}{Beispiel} \small Es seien \begin{equation} p = 4 \quad \mbox{und} \quad n_i := 3 \quad \mbox{für} \quad i = 1,...,p, \quad \mbox{also} \quad n = 12. \end{equation} Dann gilt \begin{equation} y = X\beta + \varepsilon \quad \mbox{mit} \quad \varepsilon \sim N(0_{12},\sigma^2 I_{12}) \end{equation} mit \footnotesize \begin{equation} y := \begin{pmatrix} y_{11} \\ y_{12} \\ y_{13} \\ y_{21} \\ y_{22} \\ y_{23} \\ y_{31} \\ y_{32} \\ y_{33} \\ y_{41} \\ y_{42} \\ y_{43} \end{pmatrix}, \quad X := \begin{pmatrix} 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ \end{pmatrix} \in \mathbb{R}^{12 \times 4}, \quad \beta := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \alpha_3 \\ \alpha_4 \\ \end{pmatrix} \in \mathbb{R}^{4} \quad \mbox{und} \quad \sigma^2 > 0. \end{equation} # Modellformulierung \vspace{2mm} \textcolor{darkblue}{Beispiel} \vspace{1mm} \tiny \setstretch{1} ```{r, echo = T} # Modellformulierung library(MASS) # multivariate Normalverteilung m = 3 # Anzahl Datenpunkte i-te Gruppe p = 4 # Anzahl Gruppen n = p*m # Gesamtanzahl Datenpunkte Xt = cbind(matrix(1, nrow = n, ncol = 1), # n x p Designmatrix kronecker(diag(p), matrix(1, nrow = m, ncol = 1))) X = Xt[,-2] # eliminiere zweite Spalte I_n = diag(n) # n x n Einheitsmatrix beta = matrix(c(2,0,1,-1), nrow = p) # \beta = (\mu_0,\alpha_2,\alpha_3,\alpha_4)^T sigsqr = 5 # \sigma^2 # Datenrealisierung y = mvrnorm(1, X %*% beta, sigsqr*I_n) # eine Realisierung des n-dimensionalen ZVs y ``` ```{r, echo = F} print(X) ``` # \large \setstretch{2.5} \vfill Anwendungsszenario Modellformulierung **Modellschätzung** Modellevaluation Anwendung/Praxis Selbstkontrollfragen \vfill # Modellschätzung \footnotesize \begin{theorem}[Betaparameterschätzer im EVA-Modell] \justifying \normalfont Gegeben sei die Designmatrixform der EVA in Effektdarstellung mit Referenzgruppe. Dann ergibt sich für den Betaparameterschätzer \begin{equation} \hat{\beta} := \begin{pmatrix} \hat{\mu}_0 \\ \hat{\alpha}_2 \\ \vdots \\ \hat{\alpha}_p \end{pmatrix} = \begin{pmatrix} \bar{y}_1 \\ \bar{y}_2 - \bar{y}_1 \\ \vdots \\ \bar{y}_p - \bar{y}_1 \end{pmatrix} \; , \end{equation} wobei \begin{equation} \bar{y}_i := \frac{1}{n_i}\sum_{j=1}^{n_i}y_{ij} \end{equation} das Stichprobenmittel der $i$ten Gruppe bezeichent. \end{theorem} Bemerkungen * \justifying Der Erwartungswertparameter der ersten Gruppe wird durch das Stichprobenmittel der ersten Gruppe geschätzt; der Effektparameter der $i$ten Gruppe wird durch die Differenz des Stichprobenmittels der $i$ten und der ersten Gruppe geschätzt. # Modellschätzung \vspace{2mm} \small \underline{Beweis} Wir halten zunächst fest, dass \tiny \begin{align*} \begin{split} X^\mathrm{T} X & = \setcounter{MaxMatrixCols}{20} \begin{pmatrix} 1 & \cdots & 1 & 1 & \cdots & 1 & & \cdots & & 1 & \cdots & 1 \\ 0 & \cdots & 0 & 1 & \cdots & 1 & & \cdots & & 0 & \cdots & 0 \\ & \vdots & & & \vdots & & & \vdots & & & \vdots & \\ 0 & \cdots & 0 & 0 & \cdots & 0 & & \cdots & & 1 & \cdots & 1 \\ \end{pmatrix} \begin{pmatrix} 1 & 0 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 0 \\ 1 & 1 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 1 & & 0 \\ & & & \\ \vdots & \vdots & \cdots & \vdots \\ & & & \\ 1 & 0 & & 1 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 1 \\ \end{pmatrix} \\ & = \begin{pmatrix} n & n_2 & n_3 & \cdots & n_p \\ n_2 & n_2 & 0 & \cdots & 0 \\ n_3 & 0 & n_3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ n_p & 0 & 0 & \cdots & n_p \\ \end{pmatrix}. \end{split} \end{align*} # Modellschätzung \small \underline{Beweis (fortgeführt)} Die Inverse von $X^\mathrm{T}X$ ist \begin{equation} \renewcommand{\arraystretch}{1.5} (X^\mathrm{T}X)^{-1} = \begin{pmatrix*}[c] \frac{1}{n_1} & -\frac{1}{n_1} & \cdots & -\frac{1}{n_1} \\ -\frac{1}{n_1} & \frac{n_1+n_2}{n_1 n_2} & \cdots & \frac{1}{n_1} \\ \vdots & \vdots & \ddots & \vdots \\ -\frac{1}{n_1} & \frac{1}{n_1} & \cdots & \frac{n_1+n_p}{n_1 n_p} \end{pmatrix*}. \end{equation} So gilt zum Beispiel für $p = 3$, dass \begin{equation} X^\mathrm{T} X = \begin{pmatrix} n & n_2 & n_3 \\ n_2 & n_2 & 0 \\ n_3 & 0 & n_3 \\ \end{pmatrix} \quad \mbox{und} \quad \renewcommand{\arraystretch}{1.5} (X^\mathrm{T}X)^{-1} = \begin{pmatrix*}[r] \frac{1}{n_1} & -\frac{1}{n_1} & -\frac{1}{n_1} \\ -\frac{1}{n_1} & \frac{n_1+n_2}{n_1 n_2} & \frac{1}{n_1} \\ -\frac{1}{n_1} & \frac{1}{n_1} & \frac{n_1+n_3}{n_1 n_3} \end{pmatrix*}. \end{equation} # Modellschätzung \small \underline{Beweis (fortgeführt)} Wir halten weiterhin fest, dass \tiny \begin{align} \begin{split} \renewcommand{\arraystretch}{1} X^\mathrm{T}y & = \setcounter{MaxMatrixCols}{20} \begin{pmatrix} 1 & \cdots & 1 & 1 & \cdots & 1 & & \cdots & & 1 & \cdots & 1 \\ 0 & \cdots & 0 & 1 & \cdots & 1 & & \cdots & & 0 & \cdots & 0 \\ & \vdots & & & \vdots & & & \vdots & & & \vdots & \\ 0 & \cdots & 0 & 0 & \cdots & 0 & & \cdots & & 1 & \cdots & 1 \\ \end{pmatrix} \begin{pmatrix*}[c] y_{11} \\ \vdots \\ y_{1n_1} \\ y_{21} \\ \vdots \\ y_{2n_2} \\ \\ \vdots \\ \\ y_{p1} \\ \vdots \\ y_{pn_p} \\ \end{pmatrix*} \\ & = \begin{pmatrix} \sum_{i=1}^p \sum_{j=1}^{n_i} y_{ij} \\ \sum_{j=1}^{n_2} y_{2j} \\ \vdots \\ \sum_{j=1}^{n_p} y_{pj} \\ \end{pmatrix}. \end{split} \end{align} # Modellschätzung \vspace{2mm} \small \underline{Beweis (fortgeführt)} Es ergibt sich also \tiny \begin{equation} \renewcommand{\arraystretch}{2} \hat{\beta} = (X^\mathrm{T} X)^{-1}X^\mathrm{T}y = \begin{pmatrix*}[c] \frac{1}{n_1} & -\frac{1}{n_1} & \cdots & -\frac{1}{n_1} \\ -\frac{1}{n_1} & \frac{n_1+n_2}{n_1 n_2} & \cdots & \frac{1}{n_1} \\ \vdots & \vdots & \ddots & \vdots \\ -\frac{1}{n_1} & \frac{1}{n_1} & \cdots & \frac{n_1+n_p}{n_1 n_p} \end{pmatrix*} \begin{pmatrix} \sum_{i=1}^p \sum_{j=1}^{n_i} y_{ij} \\ \sum_{j=1}^{n_2} y_{2j} \\ \vdots \\ \sum_{j=1}^{n_p} y_{pj} \\ \end{pmatrix}. \end{equation} \footnotesize \small Für die erste Komponente von $\hat{\beta}$ ergibt sich damit \tiny \begin{align} \begin{split} \hat{\beta}_1 & = \frac{1}{n_1}\sum_{i=1}^p \sum_{j=1}^{n_i} y_{ij} - \frac{1}{n_1}\sum_{j=1}^{n_2} y_{2j} - \cdots - \frac{1}{n_1}\sum_{j=1}^{n_p} y_{pj} \\ & = \frac{1}{n_1} \left( \left( \sum_{j=1}^{n_1} y_{1j} + \sum_{j=1}^{n_2} y_{2j} + \cdots + \sum_{j=1}^{n_p} y_{pj} \right) - \sum_{j=1}^{n_2} y_{2j} - \cdots - \sum_{j=1}^{n_p} y_{pj} \right) \\ & = \frac{1}{n_1}\sum_{j=1}^{n_1} y_{1j} \\ & = \bar{y}_1. \end{split} \end{align} # Modellschätzung \vspace{2mm} \small \underline{Beweis (fortgeführt)} Für die zweite Komponente von $\hat{\beta}$ und analog für alle weiteren ergibt sich \tiny \begin{align} \begin{split} \hat{\beta}_2 & = - \frac{1}{n_1} \sum_{i=1}^p \sum_{j=1}^{n_i} y_{ij} + \frac{n_1 + n_2}{n_1n_2} \sum_{j=1}^{n_2} y_{2j} + \frac{1}{n_1} \sum_{j=1}^{n_3} y_{3j} + \cdots + \frac{1}{n_1} \sum_{j=1}^{n_p} y_{pj} \\ & = \frac{n_1 + n_2}{n_1n_2} \sum_{j=1}^{n_2} y_{2j} -\frac{1}{n_1} \left( \left( \sum_{j=1}^{n_1} y_{1j} + \sum_{j=1}^{n_2} y_{2j} + \cdots + \sum_{j=1}^{n_p} y_{pj} \right) - \sum_{j=1}^{n_3} y_{3j} - \cdots - \sum_{j=1}^{n_p} y_{pj} \right) \\ & = \frac{n_1 + n_2}{n_1n_2} \sum_{j=1}^{n_2} y_{2j} -\frac{1}{n_1} \sum_{j=1}^{n_1} y_{1j} -\frac{1}{n_1} \sum_{j=1}^{n_2} y_{2j} \\ & = \frac{n_1+n_2}{n_1n_2} \sum_{j=1}^{n_2} y_{2j} -\frac{n_2}{n_1n_2} \sum_{j=1}^{n_2} y_{2j} -\frac{1}{n_1} \sum_{j=1}^{n_1} y_{1j} \\ & = \frac{n_1}{n_1n_2} \sum_{j=1}^{n_2} y_{2j} -\frac{1}{n_1} \sum_{j=1}^{n_1} y_{1j} \\ & = \frac{1}{n_2} \sum_{j=1}^{n_2} y_{2j} -\frac{1}{n_1} \sum_{j=1}^{n_1} y_{1j} = \bar{y}_2 - \bar{y}_1. \end{split} \end{align} $\hfill\Box$ # Modellschätzung \footnotesize \begin{theorem}[Varianzparameterschätzer im EVA-Modell] \normalfont \justifying Gegeben sei die Designmatrixform der EVA in Effektdarstellung mit Referenzgruppe. Dann ergibt sich für den Varianzparameterschätzer \begin{equation} \hat{\sigma}^2 = \frac{\sum_{i=1}^p \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_i)^2}{\sum_{i=1}^p n_i - p} =: s_{1 \ldots p}^2 \end{equation} wobei \begin{equation} \bar{y}_i := \frac{1}{n_i}\sum_{j=1}^{n_i}y_{ij} \end{equation} das Stichprobenmittel der $i$ten Gruppe bezeichent. \end{theorem} Bemerkungen * \justifying Der Varianzparameterschätzer ergibt sich wiederum als gebündelte Stichprobenvarianz, generalisiert auf eine beliebige Anzahl von Gruppen $p$. Wir erinnern daran, dass die gebündelte Stichprobenvarianz $s_{1 \ldots p}^2$ im Allgemeinen nicht der Stichprobenvarianz des Datenvektors $s_y^2$ entspricht. # Modellschätzung \footnotesize \underline{Beweis} Der Varianzparameterschätzer ist gegeben durch \begin{align} \begin{split} \hat{\sigma}^{2} & = \frac{(y-X\hat{\beta})^\mathrm{T}(y - X\hat{\beta})}{n - p} \\ & = \frac{1}{n - p} \left( \begin{pmatrix*}[c] y_{11} \\ \vdots \\ y_{1n_1} \\ y_{21} \\ \vdots \\ y_{2n_2} \\ \\ \vdots \\ \\ y_{p1} \\ \vdots \\ y_{pn_p} \\ \end{pmatrix*} - \begin{pmatrix} 1 & 0 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 0 \\ 1 & 1 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 1 & & 0 \\ & & & \\ \vdots & \vdots & \cdots & \vdots \\ & & & \\ 1 & 0 & & 1 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 1 \\ \end{pmatrix} \begin{pmatrix} \bar{y}_1 \\ \bar{y}_2 - \bar{y}_1 \\ \vdots \\ \bar{y}_p - \bar{y}_1 \end{pmatrix} \right)^\mathrm{T} \left( \begin{pmatrix*}[c] y_{11} \\ \vdots \\ y_{1n_1} \\ y_{21} \\ \vdots \\ y_{2n_2} \\ \\ \vdots \\ \\ y_{p1} \\ \vdots \\ y_{pn_p} \\ \end{pmatrix*} - \begin{pmatrix} 1 & 0 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 0 \\ 1 & 1 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 1 & & 0 \\ & & & \\ \vdots & \vdots & \cdots & \vdots \\ & & & \\ 1 & 0 & & 1 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 1 \\ \end{pmatrix} \begin{pmatrix} \bar{y}_1 \\ \bar{y}_2 - \bar{y}_1 \\ \vdots \\ \bar{y}_p - \bar{y}_1 \end{pmatrix} \right) \end{split} \end{align} # Modellschätzung \footnotesize \underline{Beweis (fortgeführt)} Es ergibt sich also \begin{align} \begin{split} \hat{\sigma}^{2} & = \frac{1}{n_1+\ldots+n_p - p} \begin{pmatrix*}[c] y_{11} - \bar{y}_1 \\ \vdots \\ y_{1n_1} - \bar{y}_1 \\ y_{21} - \bar{y}_1 - (\bar{y}_2 - \bar{y}_1) \\ \vdots \\ y_{2n_2} - \bar{y}_1 - (\bar{y}_2 - \bar{y}_1) \\ \\ \vdots \\ \\ y_{p1} - \bar{y}_1 - (\bar{y}_p - \bar{y}_1) \\ \vdots \\ y_{pn_p} - \bar{y}_1 - (\bar{y}_p - \bar{y}_1) \\ \end{pmatrix*}^\mathrm{T} \begin{pmatrix*}[c] y_{11} - \bar{y}_1 \\ \vdots \\ y_{1n_1} - \bar{y}_1 \\ y_{21} - \bar{y}_1 - (\bar{y}_2 - \bar{y}_1) \\ \vdots \\ y_{2n_2} - \bar{y}_1 - (\bar{y}_2 - \bar{y}_1) \\ \\ \vdots \\ \\ y_{p1} - \bar{y}_1 - (\bar{y}_p - \bar{y}_1) \\ \vdots \\ y_{pn_p} - \bar{y}_1 - (\bar{y}_p - \bar{y}_1) \\ \end{pmatrix*} \\ & = \frac{\sum_{j=1}^{n_1} (y_{1j} - \bar{y}_1)^2 + \sum_{j=1}^{n_2} (y_{2j} - \bar{y}_2)^2 + \ldots + \sum_{j=1}^{n_p} (y_{pj} - \bar{y}_p)^2}{n_1+\ldots+n_p - p} \\ & = \frac{\sum_{i=1}^p \sum_{j=1}^{n_i} (y_{ij} - \bar{y}_i)^2}{\sum_{i=1}^p n_i - p} \\ & =: s_{1 \ldots p}^2. \end{split} \end{align} $\hfill\Box$ # Modellschätzung \vspace{2mm} \tiny \setstretch{0.8} ```{r, echo = T} # Daten einlesen fname = "Daten/Einfaktorielle_Varianzanalyse_Daten.csv" D = read.table(fname, sep = ",", header = TRUE) # Datengruppen y_1 = D$dBDI[D$Setting == "F2F"] # BDI-Differenzwerte in der F2F-Gruppe y_2 = D$dBDI[D$Setting == "ONL"] # BDI-Differenzwerte in der ONL-Gruppe y_3 = D$dBDI[D$Setting == "WLC"] # BDI-Differenzwerte in der ONL-Gruppe # Modellformulierung p = 3 # drei Gruppen m = length(y_1) # balancierters Design mit n_i = 40 n = p*m # Datenvektordimension y = matrix(c(y_1, y_2, y_3), nrow = n) # Datenvektor Xt = cbind(matrix(1, nrow = n, ncol = 1), # Designmatrix kronecker(diag(p), matrix(1, nrow = m, ncol = 1))) X = Xt[,-2] # eliminiere zweite Spalte # Modellschätzung 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 s_sqr_123 = ((m-1)*var(y_1) + # gebündelte Stichprobenvarianz (m-1)*var(y_2) + (m-1)*var(y_3))/(m+m+m-p) # Ausgabe cat( "hat{beta} : ", beta_hat, "\nbar{y}_1,bar{y}_2,bar{y}_3 : ", c(mean(y_1),mean(y_2),mean(y_3)), "\nbar{y}_1,bar{y}_2-bar{y}_1,bar{y}_3-bar{y}_1 : ", c(mean(y_1),mean(y_2)-mean(y_1),mean(y_3)-mean(y_1)), "\nhat{sigsqr} : ", sigsqr_hat, "\ns_123^2 : ", s_sqr_123, "\ns_y^2 : ", var(y)) ``` # \large \setstretch{2.5} \vfill Anwendungsszenario Modellformulierung Modellschätzung **Modellevaluation** Anwendung/Praxis Selbstkontrollfragen \vfill # Modellevaluation \small \textcolor{darkblue}{Vorbemerkungen und Überblick} \footnotesize Prinzipiell sind alle Parameterschätzwerte in einem EVA-Modell von Interesse. Mit der T-Teststatistik können einzelne Komponenten oder lineare Kombinationen der Betaparameterschätzwerte im Sinne von Hypothesentests evaluiert werden. Nichtsdestotrotz steht im EVA-Szenario häufig die Evaluation der Nullhypothese, dass alle Effektparameter gleich Null sind, im Vordergrund. Intuitiv besagt diese Nullhypothese, dass die wahren, aber unbekannten Erwartungswertparameter aller Gruppen identisch sind, formal haben wir \begin{equation} H_0 : \alpha_i = 0 \quad \mbox{für alle} \quad i = 2,...,p \end{equation} und \begin{equation} H_1 : \alpha_i \neq 0 \quad \mbox{für mindestens ein} \quad i = 2,...,p \; . \end{equation} Zur Überprüfung von $H_0$ wird im Allgemeinen eine $F$-Teststatistik eingesetzt. Dabei wird das volle Modell mit dem reduzierten Modell verglichen, in der die Designmatrix nur die erste Spalte der Designmatrix des vollen Modells enthält. Im Folgenden entwickeln wir zunächst diese $F$-Teststatistik anhand einer *Quadratsummenzerlegung* der Datenvariabilität in einem EVA-Szenario und betrachten in diesem Zusammenhang auch das dem $R^2$ analoge Effektstärkemaß $\eta^2$. Ausgestattet mit der speziellen Form der $F$-Teststatistik in dem hier betrachteten Szenario diskutieren wir dann den traditionellen EVA-Hypothesentest. Wir verzichten auf eine Analyse der Testgütefunktion und eine Diskussion der Powerfunktion. # Modellevaluation \vspace{2mm} \footnotesize \begin{theorem}[Quadratsummenzerlegung bei einfaktorieller Varianzanalyse] \justifying \normalfont Für $i = 1,...,p$ und $j = 1,...,n_i$ sei $y_{ij}$ die $j$te Datenvariable in der $i$ten Gruppe eines EVA-Szenarios. Weiterhin seien mit $n := \sum_{i=1}^p n_i$ \center \vspace{1mm} \begin{tabular}{ll} $\bar{y} := \frac{1}{n}\sum_{i=1}^p \sum_{j=1}^{n_i} y_{ij}$ & das \textit{Gesamtstichprobenmittel} \\\\ $\bar{y}_i =\frac{1}{n_i}\sum_{j=1}^{n_i} y_{ij}$ & das \textit{$i$te Stichprobenmittel} \\\\ \end{tabular} \vspace{1mm} \flushleft sowie \center \vspace{1mm} \begin{tabular}{ll} $\mbox{SQT} := \sum_{i=1}^p \sum_{j=1}^{n_i} (y_{ij}-\bar{y})^{2}$ & die \textit{total sum of squares} \\\\ $\mbox{SQB} := \sum_{i=1}^p n_i(\bar{y}_i-\bar{y})^{2}$ & die \textit{``between'' sum of squares} \\\\ $\mbox{SQW} := \sum_{i=1}^p \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i)^{2}$ & die \textit{``within'' sum of squares} \; . \\\\ \end{tabular} \vspace{1mm} \flushleft Dann gilt \begin{equation} \mbox{SQT} = \mbox{SQB} + \mbox{SQW}. \end{equation} \end{theorem} \vspace{-1mm} Bemerkungen \vspace{-2mm} * Im Vergleich mit der Quadratsummenzerlegung mit Ausgleichsgerade (vgl. Einheit (2) in *Allgemeines Lineares Modell*) spielt $\mbox{SQB}$ die Rolle der erklärten Quadratsumme und $\mbox{SQW}$ die Rolle der residuellen Quadratsumme. * Die Ausdrücke "between" und "within" beziehen sich hierbei darauf, dass $\mbox{SQB}$ die Variabilität zwischen den Gruppen und $\mbox{SQW}$ die Variabilität innerhalb der Gruppen abbildet. # Modellevaluation \small \underline{Beweis} Es gilt \footnotesize \begin{align} \begin{split} \mbox{SQT} & = \sum_{i=1}^p \sum_{j=1}^{n_i}(y_{ij}-\bar{y})^{2} \\ & = \sum_{i=1}^p \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i+\bar{y}_i-\bar{y})^2 \\ & = \sum_{i=1}^p \sum_{j=1}^{n_i} \left((y_{ij}-\bar{y}_i)+(\bar{y}_i-\bar{y}) \right)^2 \\ & = \sum_{i=1}^p \sum_{j=1}^{n_i} \left( (y_{ij}-\bar{y}_i)^2 +2(y_{ij}-\bar{y}_i)(\bar{y}_i-\bar{y}) +(\bar{y}_i-\bar{y})^2 \right) \\ & = \sum_{i=1}^p \left( \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i)^2 +\sum_{j=1}^{n_i} 2(y_{ij}-\bar{y}_i)(\bar{y}_i-\bar{y}) +\sum_{j=1}^{n_i} (\bar{y}_i-\bar{y})^2 \right) \\ & = \sum_{i=1}^p \left( \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i)^2 +2(\bar{y}_i-\bar{y})\sum_{j=1}^{n_i}(y_{ij}-\bar{y}_i) +n_i(\bar{y}_i-\bar{y})^2 \right) \\ & = \sum_{i=1}^p \left( \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i)^2 +2(\bar{y}_i-\bar{y})\sum_{j=1}^{n_i} \left(y_{ij}-\frac{1}{n_i}\sum_{j=1}^{n_i}{y_{ij}}\right) +n_i(\bar{y}_i-\bar{y})^2 \right) \end{split} \end{align} # Modellevaluation \small \underline{Beweis (fortgeführt)} und weiter \footnotesize \begin{align*} \begin{split} \mbox{SQT} & = \sum_{i=1}^p \left( \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i)^2 +2(\bar{y}_i-\bar{y})\left(\sum_{j=1}^{n_i} y_{ij}-\sum_{j=1}^{n_i}\left(\frac{1}{n_i}\sum_{j=1}^{n_i}{y_{ij}}\right)\right) +n_i(\bar{y}_i-\bar{y})^2 \right) \\ & = \sum_{i=1}^p \left( \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i)^2 +2(\bar{y}_i-\bar{y})\left(\sum_{j=1}^{n_i} y_{ij}-\frac{n_i}{n_i}\sum_{j=1}^{n_i}{y_{ij}}\right) +n_i(\bar{y}_i-\bar{y})^2 \right) \\ & = \sum_{i=1}^p \left( \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i)^2 +2(\bar{y}_i-\bar{y})\left(\sum_{j=1}^{n_i} y_{ij}-\sum_{j=1}^{n_i}{y_{ij}}\right) +n_i(\bar{y}_i-\bar{y})^2 \right) \\ & = \sum_{i=1}^p \left( \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i)^2 +n_i(\bar{y}_i-\bar{y})^2 \right) \\ & = \sum_{i=1}^p \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_i)^2 + \sum_{i=1}^p n_i(\bar{y}_i-\bar{y})^2 \end{split} \end{align*} \small und damit \begin{equation} \mbox{SQT} = \mbox{SQB} + \mbox{SQW}. \end{equation} $\hfill\Box$ # Modellevaluation \footnotesize \begin{definition}[Effektstärkenmaß $\eta^2$] \justifying Für ein EVA-Szenario seien die \textit{between sum of squares} $\mbox{SQB}$ und die \textit{total sum of squares} $\mbox{SQT}$ definiert wie oben. Dann ist das \textit{Effektstärkenmaß} $\eta^2$ definiert als \begin{equation} \eta^2 := \frac{\mbox{SQB}}{\mbox{SQT}} \end{equation} \end{definition} Bemerkungen * $\eta^2$ ist analog zum Bestimmtheitsmaß $R^2$ der Regression definiert. * $\eta^2$ gibt den Anteil der Varianz zwischen den Gruppen an der Gesamtvarianz der Daten an. * Mit dem Theorem zur Quadratsummenzerlegung bei EVA folgt sofort $0 \le \eta^2 \le 1$, da \begin{align} \begin{split} \mbox{SQB} &= 0 \Rightarrow \mbox{SQT} = \mbox{SQW} \quad \mbox{und} \quad \eta^2 = 0 \\ \mbox{SQW} &= 0 \Rightarrow \mbox{SQT} = \mbox{SQB} \quad \mbox{und} \quad \eta^2 = 1 \; . \end{split} \end{align} # Modellevaluation \footnotesize \begin{theorem}[F-Teststatistik der einfaktoriellen Varianzanalyse] \justifying \normalfont Es sei \begin{equation} y = X\beta + \varepsilon \quad \mbox{mit} \quad \varepsilon \sim N(0_n,\sigma^2I_n) \end{equation} die Designmatrixform der Effektdarstellung mit Referenzgruppe des EVA-Modells und im Sinne der Definition der F-Statistik sei dieses Modell partioniert mit $p_0 := 1$ und $p_1 := p - 1$. Weiterhin seien \center \vspace{3mm} \begin{tabular}{ll} $\mbox{MSB} := \frac{\mbox{SQB}}{p-1}$ & die \textit{mean between sum of squares} und \\\\ $\mbox{MSW} := \frac{\mbox{SQW}}{n-p}$ & die \textit{mean within sum of squares} \\\\ \end{tabular} \flushleft Dann gilt \begin{equation} F = \frac{\mbox{MSB}}{\mbox{MSW}}. \end{equation} \end{theorem} Bemerkungen * $p_0 := 1$ impliziert, dass das reduzierte Modell die Designmatrix $X_0 := 1_n$ hat. * $p_0 := 1$ impliziert zudem, dass das reduzierte Modell den Betaparameter $\beta_0 := \mu_0$ hat. * $p_0 := 1$ impliziert damit auch, dass das reduzierte Modell keine Effektparameter $\alpha_i$ für $i = 2,...,p$ hat. * Die Zahl $p-1$ wird auch als "Zähler-Freiheitsgrade" bezeichnet. * Die Zahl $n-p$ wird auch als "Nenner-Freiheitsgrade" bezeichnet. # Modellevaluation \small \underline{Beweis} Wir halten zunächst fest, dass für den Betaparameterschätzer des reduzierten Modells gilt, dass \begin{align} \begin{split} \hat{\beta}_0 = (X_0^\mathrm{T}X_0)^{-1}X_0^\mathrm{T}y = (1_n^\mathrm{T} 1_n)^{-1}1_n^\mathrm{T}y = \frac{1}{n} \sum_{i=1}^p \sum_{j = 1}^{n_i} y_{ij} = \bar{y}. \end{split} \end{align} Weiterhin ergibt sich \begin{align} \begin{split} \hat{\varepsilon}_0^\mathrm{T} \hat{\varepsilon}_0 = (y - X_0\hat{\beta}_0)^\mathrm{T}(y - X_0\hat{\beta}_0) = (y - 1_n\bar{y})^\mathrm{T}(y - 1_n\bar{y}) = \sum_{i=1}^p \sum_{j=1}^{n_i} (y_{ij} - \bar{y})^2 = \mbox{SQT}. \end{split} \end{align} Der Betaparameterschätzer des vollständigen Modells ergibt sich wie oben gesehen zu \tiny \begin{equation} \renewcommand{\arraystretch}{1.4} \hat{\beta} = \begin{pmatrix} \hat{\mu}_0 \\ \hat{\alpha}_2 \\ \vdots \\ \hat{\alpha}_m \end{pmatrix} = \begin{pmatrix} \frac{1}{n_1}\sum_{j=1}^{n_1} y_{1j} \\ \frac{1}{n_2}\sum_{j=1}^{n_2} y_{2j} - \frac{1}{n_1}\sum_{j=1}^{n_1} y_{1j} \\ \vdots \\ \frac{1}{n_m}\sum_{j=1}^{n_m} y_{mj} - \frac{1}{n_1}\sum_{j=1}^{n_1} y_{1j} \\ \end{pmatrix} = \begin{pmatrix} \bar{y}_1 \\ \bar{y}_2 - \bar{y}_1 \\ \vdots \\ \bar{y}_p - \bar{y}_1 \\ \end{pmatrix}, \end{equation} \small sodass # Modellevaluation \small \underline{Beweis (fortgeführt)} \begin{align*} \renewcommand{\arraystretch}{1} \begin{split} \hat{\varepsilon}^\mathrm{T}\hat{\varepsilon} & = (y - X\hat{\beta})^\mathrm{T}(y - X\hat{\beta}) \\ & = \left( \begin{pmatrix} y_{11} \\ \vdots \\ y_{1n_1} \\ y_{21} \\ \vdots \\ y_{2n_2} \\ \\ \vdots \\ \\ y_{p1} \\ \vdots \\ y_{pn_p} \\ \end{pmatrix} - \begin{pmatrix} 1 & 0 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 0 \\ 1 & 1 & & 0 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 1 & & 0 \\ & & & \\ \vdots & \vdots & \cdots & \vdots \\ & & & \\ 1 & 0 & & 1 \\ \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 1 \\ \end{pmatrix} \begin{pmatrix} \bar{y}_1 \\ \bar{y}_2 - \bar{y}_1 \\ \vdots \\ \bar{y}_p - \bar{y}_1 \\ \end{pmatrix} \right)^\mathrm{T} (y - X\hat{\beta}) \end{split} \end{align*} # Modellevaluation \small \underline{Beweis (fortgeführt)} und weiter \begin{align*} \begin{split} \hat{\varepsilon}^\mathrm{T}\hat{\varepsilon} & = \begin{pmatrix} y_{11} - \bar{y}_1 \\ \vdots \\ y_{1n_1} - \bar{y}_1 \\ y_{21} - \bar{y}_1 - \bar{y}_2 + \bar{y}_1 \\ \vdots \\ y_{2n_2} - \bar{y}_1 - \bar{y}_2 + \bar{y}_1 \\ \\ \vdots \\ \\ y_{p1} - \bar{y}_1 - \bar{y}_p + \bar{y}_1 \\ \vdots \\ y_{pn_p} - \bar{y}_1 - \bar{y}_p + \bar{y}_1 \\ \end{pmatrix}^\mathrm{T} \begin{pmatrix} y_{11} - \bar{y}_1 \\ \vdots \\ y_{1n_1} - \bar{y}_1 \\ y_{21} - \bar{y}_1 - \bar{y}_2 + \bar{y}_1 \\ \vdots \\ y_{2n_2} - \bar{y}_1 - \bar{y}_2 + \bar{y}_1 \\ \\ \vdots \\ \\ y_{p1} - \bar{y}_1 - \bar{y}_p + \bar{y}_1 \\ \vdots \\ y_{pn_p} - \bar{y}_1 - \bar{y}_p + \bar{y}_1 \\ \end{pmatrix} = \begin{pmatrix} y_{11} - \bar{y}_1 \\ \vdots \\ y_{1n_1} - \bar{y}_1 \\ y_{21} - \bar{y}_2 \\ \vdots \\ y_{2n_2} - \bar{y}_2 \\ \\ \vdots \\ \\ y_{p1} - \bar{y}_p \\ \vdots \\ y_{pn_p} - \bar{y}_p \\ \end{pmatrix}^\mathrm{T} \begin{pmatrix} y_{11} - \bar{y}_1 \\ \vdots \\ y_{1n_1} - \bar{y}_1 \\ y_{21} - \bar{y}_2 \\ \vdots \\ y_{2n_2} - \bar{y}_2 \\ \\ \vdots \\ \\ y_{p1} - \bar{y}_p \\ \vdots \\ y_{pn_p} - \bar{y}_p \\ \end{pmatrix} \\ & = \sum_{i=1}^p \sum_{j = 1}^{n_i} (y_{ij} - \bar{y}_i)^2 \\ & = \mbox{SQW}. \end{split} \end{align*} # Modellevaluation \small \underline{Beweis (fortgeführt)} Mit dem Theorem zur Quadratsummenzerlegung bei einfaktorieller Varianzanalyse \begin{equation} \mbox{SQT} = \mbox{SQB} + \mbox{SQW} \quad \Leftrightarrow \quad \mbox{SQB} = \mbox{SQT} - \mbox{SQW} \end{equation} folgt sofort, dass \begin{align} \begin{split} \mbox{SQB} &= \mbox{SQT} - \mbox{SQW} \\ &= \hat{\varepsilon}_0^\mathrm{T}\hat{\varepsilon}_0 - \hat{\varepsilon}^\mathrm{T}\hat{\varepsilon}. \end{split} \end{align} Dann aber folgt auch direkt, dass \begin{equation} \frac{\mbox{MSB}}{\mbox{MSW}} = \frac{\frac{\mbox{SQB}}{p-1}}{\frac{\mbox{SQW}}{n-p}} = \frac{\frac{\hat{\varepsilon}_0^\mathrm{T}\hat{\varepsilon}_0 - \hat{\varepsilon}^\mathrm{T}\hat{\varepsilon}}{p-1}}{\frac{\hat{\varepsilon}^\mathrm{T}\hat{\varepsilon}}{n-p}} = F \end{equation} # Modellevaluation \footnotesize \begin{theorem}[Effektstärkenmaß $\eta^2$ und F-Teststatistik] \justifying \normalfont Für ein EVA-Szenario mit $p$ Gruppen und Gesamtdatenpunktanzahl $n$ seien das Effektstärkenmaß $\eta^2$ und die $F$-Teststatistik wie oben definiert. Dann gilt \begin{equation} \eta^2 = \frac{F(p-1)}{F(p-1) + (n-p)} \end{equation} \end{theorem} Bemerkungen * Das Verhältnis von $F$ und $\eta^2$ ist Analog zum Verhältnis von $T$ und Cohen's $d$. * Die gleichzeitige Angabe von $F$ und $\eta^2$ ist redundant. # Modellevaluation \small \setstretch{1} \underline{Beweis} Wir halten zunächst fest, dass \begin{equation} F = \frac{\mbox{SQB}}{\mbox{SQW}}\cdot \frac{n-p}{p-1} \quad \Leftrightarrow \quad \mbox{SQB} = \frac{p-1}{n-p}\cdot \mbox{SQW} \cdot F \; . \end{equation} Damit folgt dann \begin{align} \begin{split} \eta^2 & = \frac{\mbox{SQB}}{\mbox{SQT}} = \frac{\mbox{SQB}}{\mbox{SQB}+\mbox{SQW}} = \frac{\frac{p-1}{n-p}\cdot \mbox{SQW} \cdot F}{\frac{p-1}{n-p}\cdot \mbox{SQW} \cdot F +\mbox{SQW}} \\ & = \frac{\frac{F(p-1)}{n-p}\cdot \mbox{SQW} }{\frac{F(p-1)}{n-p}\cdot \mbox{SQW} +\mbox{SQW}} = \frac{\frac{F(p-1)}{n-p}\cdot \mbox{SQW} }{\left(\frac{F(p-1)}{n-p} + 1 \right)\cdot \mbox{SQW}} \\ & = \frac{\frac{F(p-1)}{n-p}}{\frac{F(p-1)}{n-p} + \frac{n-p}{n-p}} = \frac{\frac{F(p-1)}{n-p}}{\frac{F(p-1) + (n-p)}{n-p}} = \frac{F(p-1)}{F(p-1)+(n-p)} \; . \end{split} \end{align} $\hfill\Box$ # Modellevaluation \small \setstretch{2} \textcolor{darkblue}{Gliederung} (vgl. Einheit (11), Folie 29 in \textit{Wahrscheinlichkeitstheorie und Frequentistische Inferenz}) (1) Statistisches Modell $\checkmark$ (2) Testhypothesen $\checkmark$ (3) Teststatistik $\checkmark$ (4) Test (5) Testumfangkontrolle (6) p-Wert # Modellevaluation (3) Teststatistik \setstretch{3} \textcolor{darkblue}{Bezeichnungskonventionen für $f$-Zufallsvariablen} \small * Für eine (nichtzentrale) $f$-Zufallsvariable $\xi$ schreiben wir $\xi \sim f(n_1, n_2)$ (oder $\xi \sim f(\delta, n_1, n_2)$). * Die WDF einer (nichtzentralen) $f$-Zufallsvariable ist $f(\cdot; n_1, n_2)$ (oder $f(\cdot;\delta, n_1, n_2)$). * Die KVF einer (nichtzentralen) $f$-Zufallsvariable ist $\varphi(\cdot; n_1, n_2)$ (oder $\varphi(\cdot; \delta, n_1, n_2)$). * Die Teststatistik eines F-Test-Designs/Hypothesentests bezeichnen wir mit $F$. * Die Realisierung der F-Teststatistik für einen Datensatz bezeichnen wir mit $f$. \vfill # Modellevaluation (4) Test \footnotesize \begin{definition}[F-Test für einfaktorielle Varianzanalyse] \justifying Gegeben sie das Modell der einfaktoriellen Varianzanalyse sowie die zusammengesetzten Null- und Alternativhypothesen \begin{equation} H_0 : \alpha_i = 0 \quad \mbox{für alle} \quad i = 2,...,p \end{equation} bzw. \begin{equation} H_1 : \alpha_i \neq 0 \quad \mbox{für mindestens ein} \quad i = 2,...,p \; . \end{equation} Weiterhin sei die F-Teststatistik definiert durch \begin{equation} F = \frac{\mbox{MSB}}{\mbox{MSW}} \end{equation} mit der mean between sum of squares $\mbox{MSB}$ und der mean within sum of squares $\mbox{MSW}$. Dann ist der \textit{F-Test für die einfaktoriellen Varianzanalyse (EVA-F-Test)} definiert als der kritische Wert-basierte Test \begin{equation} \phi(y) := 1_{\{F \ge k\}} := \begin{cases} 1 & F \ge k \\ 0 & F < k \end{cases}. \end{equation} \end{definition} # Modellevaluation (5) Testumfangkontrolle \footnotesize \begin{theorem}[Testumfangkontrolle] \justifying \normalfont $\phi$ sei der F-Test zur einfaktoriellen Varianzanalyse. Dann ist $\phi$ ein Level-$\alpha_0$-Test mit Testumfang $\alpha_0$, wenn der kritische Wert definiert ist durch \begin{equation} k_{\alpha_0} := \varphi^{-1}(1-\alpha_0; p-1, n-p), \end{equation} wobei $\varphi^{-1}(\cdot; p-1, n-p)$ die inverse KVF der $f$-Verteilung mit Freiheitsgradparametern $p-1$ und $n-p$ ist. \end{theorem} \vfill # Modellevaluation (5) Testumfangkontrolle \small \center Wahl von $k_{\alpha_0} := \varphi^{-1}(1 - \alpha_0; p-1, n-p)$ mit $p = 3, m = 12, \alpha_0 := 0.05$ und Ablehnungsbereich ```{r, echo = F, eval = F} par( family = "sans", mfcol = c(1,2), pty = "m", bty = "l", lwd = 1, las = 1, mgp = c(2,1,0), xaxs = "i", yaxs = "i", font.main = 1, cex = 1, cex.main = 1.2) # Parameter p = 3 # Anzahl Gruppen m = 12 # Anzahl experimenteller Einheiten pro Gruppe n = p*m # Gesamtstichprobengröße alpha_0 = 0.05 # Signifikanzlevel k_alpha_0 = qf(1 - alpha_0, p-1,n-p) # kritischer Wert eff = seq(0,6,length=1e4 ) # F-Statistikwerte Peff = pf(eff,p-1,n-p) # F-Statistik KVF für H_0 peff = df(eff,p-1,n-p) # F-Statistik WDF für H_0 # KVF-Perspektive plot(eff, Peff, type = "l", xlab = " ", ylab = " ", ylim = c(0,1), main = TeX("$\\varphi(F;2,33)$")) lines(k_alpha_0, 0, type = "p", pch = 16, xpd = TRUE) lines(min(eff), 1 - alpha_0, type = "p", pch = 16, xpd = TRUE) arrows( x0 = min(eff), y0 = 1 - alpha_0, x1 = k_alpha_0, y1 = 1 - alpha_0, col = "darkorange", angle = 45, length = .1) arrows( x0 = k_alpha_0, y0 = 1-alpha_0, x1 = k_alpha_0, y1 = 0, col = "darkorange", angle = 45, length = .1) text(k_alpha_0-0.5, 0.1, TeX("$\\k_{\\alpha_0}$"), xpd = TRUE) text(1, 1, TeX("$1 - \\alpha_0$"), xpd = TRUE) # WDF-Perspektive plot(eff, peff, type = "l", ylab = " ", xlab = " ", ylim = c(0,.4), main = TeX("$f(F;2,33)$")) polygon(c(eff[eff >= k_alpha_0], max(eff), k_alpha_0), c(peff[eff >= k_alpha_0], 0, 0 ), col = "gray90", border = NA) lines(seq(k_alpha_0, max(eff), len = 1e2), rep(0,1e2), type = "l", lwd = 5, col = "darkorange") lines(k_alpha_0, 0, type = "p", pch = 16, xpd = TRUE) text(k_alpha_0-0.5, 0.03, TeX("$\\k_{\\alpha_0}$") , xpd = TRUE) text(5, 0.05, TeX("$P(F >= k_{\\alpha_0}) = \\alpha_0$"), xpd = TRUE, cex = 1, col = "gray50") # PDF-Speicherung dev.copy2pdf( file = "Abbildungen/eva_testumfang.pdf", width = 8, height = 4) ``` ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/eva_testumfang.pdf") ``` \vfill # Modellevaluation (5) Testumfangkontrolle \footnotesize \underline{Beweis} Die Testgütefunktion des betrachteten Tests im vorliegenden Testszenario ist definiert als \begin{equation} q : \mathbb{R}^p \to [0,1], \beta \mapsto q_\phi(\beta) := \mathbb{P}_\beta(\phi = 1). \end{equation} Wir haben bereits gesehen (siehe Einheit (8) in *Allgemeines Lineares Modell*), dass die F-Teststatistik für $p_1 = p - 1$ nach einer nichtzentralen $f$-Verteilung verteilt ist, \begin{equation} F \sim f(\delta,p-1, n-p). \end{equation} Weiterhin ist der Ablehnungsbereich des hier betrachteten Tests gegeben als $[k,\infty[$. Für die funktionale Form der Testgütefunktion ergibt sich also \begin{align} \begin{split} \mathbb{P}_\beta(\phi = 1) &= \mathbb{P}_\beta(F \in [k,\infty[) \\ &= \mathbb{P}_\beta(F \ge k) \\ &= 1 - \mathbb{P}_\beta(F \le k) \\ &= 1 - \varphi(k;\delta,p-1,n-p), \\ \end{split} \end{align} wobei $\varphi(k;\delta,p-1,n-p)$ den KVF-Wert der nichtzentralen $f$-Verteilung an Stelle $k$ mit Nichtzentralitätsparameter $\delta$ sowie Freiheitsgradparametern $p-1$ und $n-p$ bezeichnet (vgl. Einheit (8) in *Allgemeines Lineares Modell*). Damit der betrachtete Test ein Level-$\alpha_0$-Test ist, muss bekanntlich gelten, dass \begin{equation} q_\phi(\beta) \le \alpha_0 \quad \mbox{für alle} \quad \beta \in \Theta_0 \quad \mbox{mit} \quad \Theta_0 = \mathbb{R} \times \{0_{p-1}\}. \end{equation} # Modellevaluation (5) Testumfangkontrolle \vspace{2mm} \footnotesize \underline{Beweis} Der Nichtzentralitätsparameter ist gegeben durch (siehe Einheit (8) in *Allgemeines Lineares Modell*) \begin{equation} \delta = \frac{1}{\sigma^2} c^\mathrm{T}\beta \left(c^\mathrm{T}(X^\mathrm{T}X)^{-1}c\right)^{-1}c^\mathrm{T} \beta \; . \end{equation} Mit dem dazugehörigen $c$ und unter der Nullhypothese, dass $\beta \in \Theta_0$, also \begin{equation} c = \begin{pmatrix} 0 \\ 1_{p-1} \end{pmatrix} \in \mathbb{R}^p \quad \mbox{und} \quad \beta = \begin{pmatrix} \mu_0 \\ 0_{p-1} \end{pmatrix} \in \mathbb{R}^p \; , \end{equation} folgt dann aber $\delta = 0$ und somit \begin{equation} q_\phi(\beta) = 1 - \varphi(k; p-1, n-p) \quad \mbox{für alle} \quad \beta \in \Theta_0 \; , \end{equation} wobei $\varphi(k;p-1,n-p)$ den Wert der KVF der $f$-Verteilung an der Stelle $k$ mit Freiheitsgradparametern $p-1$ und $n-p$ bezeichnet (vgl. Einheit (8) in *Allgemeines Lineares Modell*). Der Testumfang des betrachteten Tests ergibt sich per Definition (siehe Einheit (11) in *Wahrscheinlichkeitstheorie und Frequentistische Inferenz*) als \begin{equation} \alpha = \max_{{\beta \in \Theta_0}} q_{\phi}(\beta) = 1 - \varphi(k;p-1,n-p), \end{equation} da $q_{\phi}(\beta)$ für $\beta \in \Theta_0$ nicht von $\mu_0$ abhängt. Wir müssen also lediglich zeigen, dass die Wahl von $k_{\alpha_0}$ wie im Theorem garantiert, dass $\phi$ den Testumfang $\alpha_0$ hat. Sei also $k := k_{\alpha_0}$. Dann gilt für alle $\beta \in \Theta_0$ \begin{equation} q_\phi(\beta) = 1-\varphi( \varphi^{-1}(1-\alpha_0; p-1, n-p); p-1, n-p) = 1-(1-\alpha_0) = \alpha_0 \end{equation} und damit ist alles gezeigt. $\hfill\Box$ # Modellevaluation (6) p-Wert \footnotesize Nach Definition ist der p-Wert das kleinste Signifikanzlevel $\alpha_0$, bei welchem man die Nullhypothese basierend auf einem vorliegenden Wert der Teststatistik ablehnen würde. Wir wollen einen vorliegenden Wert der $F$-Teststatistik hier mit $f$ bezeichnen. Bei $F = f$ würde $H_0$ für jedes $\alpha_0$ mit $f \ge \psi^{-1}(1-\alpha_0;p-1,n-p)$ abgelehnt werden. Für ein solches $\alpha_0$ gilt aber \begin{equation} \alpha_0 \ge \mathbb{P}(F \ge f), \end{equation} denn \begin{align} \begin{split} f & \ge \psi^{-1}(1-\alpha_0;p-1,n-p) \\ \Leftrightarrow \psi(f;p-1,n-p) & \ge \psi(\psi^{-1}(1-\alpha_0;p-1,n-p), p-1,n-p) \\ \Leftrightarrow \psi(f;p-1,n-p) & \ge 1-\alpha_0 \\ \Leftrightarrow \mathbb{P}(F \le f) & \ge 1-\alpha_0 \\ \Leftrightarrow \alpha_0 & \ge 1 - \mathbb{P}(F \le f) \\ \Leftrightarrow \alpha_0 & \ge \mathbb{P}(F \ge f). \end{split} \end{align} Das kleinste $\alpha_0 \in [0,1]$ mit $\alpha_0 \ge \mathbb{P}(F \ge f)$ ist dann $\alpha_0 = \mathbb{P}(F \ge f)$, also folgt \begin{equation} \mbox{p-Wert} = \mathbb{P}(F \ge f) = 1 - \varphi(f; p-1, n-p). \end{equation} # Modellevaluation (6) p-Wert Praktisches Vorgehen \footnotesize * \justifying Man nimmt an, dass ein vorliegender Datensatz aus $p$ Gruppendatensätzen besteht, wobei $y_{11}, ...,y_{1n_1}$, $y_{21}, ...,y_{2n_2}$, ..., $y_{p1}, ...,y_{pn_p}$ Realisationen von $y_{1j} \sim N(\mu_0,\sigma^2)$ und $y_{ij} \sim N(\mu_0 + \alpha_i,\sigma^2)$ für $i = 2,...,p$ mit wahren, aber unbekannten Parametern $\mu_0, \alpha_i, i = 2,...,p$ und $\sigma^2>0$ sind. * Man möchte entscheiden ob $H_0 : \alpha_i = 0$ für alle $i = 2,...,p$ eher zutrifft oder eher nicht. * Man wählt ein Signifikanzniveau $\alpha_0$ und bestimmt den zugehörigen Freiheitsgradparameter-abhängigen kritischen Wert $k_{\alpha_0}$. Zum Beispiel gilt bei Wahl von $\alpha_0 := 0.05, p = 3, m = 12, i = 1,2,3$ und somit $n = 36$, dass $k_{\alpha_0} = \varphi^{-1}(1 - 0.05; 3-1, 36-3) \approx 3.28$ ist. * Anhand der MSB und MSW berechnet man den realisierten Wert der F-Teststatistik, den wir hier mit $f$ bezeichnen. * Wenn $f$ größer-gleich $k_{\alpha_0}$ ist, lehnt man die Nullhypothese ab, andernfalls nicht. * Die oben entwickelte Theorie garantiert dann, dass man im Mittel in höchstens $\alpha_0 \cdot 100$ von 100 Fällen die Nullhypothese fälschlicherweise ablehnt. * Schließlich ergibt sich der assoziierte p-Wert der realisiertern F-Teststatistik $f$ zu \begin{equation} \mbox{p-Wert} = \mathbb{P}(F \ge f) = 1 - \varphi(f; p-1, n-p) \end{equation} # \large \setstretch{2.5} \vfill Anwendungsszenario Modellformulierung Modellschätzung Modellevaluation **Anwendung/Praxis** Selbstkontrollfragen \vfill # Anwendung/Praxis \small \vspace{2mm} \textcolor{darkblue}{Daten einlesen:} $j = 1,...,10$ für jede Gruppe \tiny \setstretch{0.6} \vspace{1mm} ```{r, echo = T} fname = "Daten/Einfaktorielle_Varianzanalyse_Daten.csv" D = read.table(fname, sep = ",", header = TRUE) ``` \vspace{-1mm} ```{r, echo = F} knitr::kable(D[c(1:10, 41:50, 81:90),2:6], align = "cccrr", "pipe") ``` # Anwendung/Praxis \textcolor{darkblue}{Boxplot} \vspace{2mm} \tiny \setstretch{1} ```{r, echo = T, eval = F} # Daten einlesen fname = "Daten/Einfaktorielle_Varianzanalyse_Daten.csv" D = read.table(fname, sep = ",", header = TRUE) # Abbildungsparameter 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, # Textvergrößerungsfaktor cex.main = 1.2) # Titeltextvergrößerungsfaktor # Boxplot boxplot(dBDI ~ Setting, # BDI ~ Setting enkodiert die Datengruppierung data = D) # PDF-Speicherung dev.copy2pdf( file = "Abbildungen/eva_boxplot.pdf", width = 7, height = 4) ``` # Anwendung/Praxis \textcolor{darkblue}{Boxplot} \vspace{-4mm} ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/eva_boxplot.pdf") ``` # Anwendung/Praxis \textcolor{darkblue}{Balkendiagramm} \vspace{2mm} \tiny \setstretch{0.8} ```{r, echo = T, eval = F} # Daten einlesen fname = "Daten/Einfaktorielle_Varianzanalyse_Daten.csv" D = read.table(fname, sep = ",", header = TRUE) # Datenselektion A = data.frame(F2F = D$dBDI[D$Setting == "F2F"], # BDI-Daten für F2F ONL = D$dBDI[D$Setting == "ONL"], # BDI-Daten für ONL WLC = D$dBDI[D$Setting == "WLC"]) # BDI-Daten für WLC # deskriptive Statistiken groupmeans = colMeans(A) # Gruppenmittelwerte groupstds = apply(A,2,sd) # Gruppenstandardabweichungen # Balkendiagramm 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 x = barplot(groupmeans, # Ausgabe der x-Ordinaten ylim = c(-5,15), col = "gray90", ylab = "dBDI", xlab = "Condition") 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) # Linielänge # PDF-Speicherung dev.copy2pdf( file = "Abbildungen/eva_barplot.pdf", width = 7, height = 4) ``` # Anwendung/Praxis \textcolor{darkblue}{Balkendiagramm} \vspace{2mm} \center Gruppenmittelwerte $\pm$ Gruppenstandardabweichungen \vspace{-4mm} ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/eva_barplot.pdf") ``` # Anwendung/Praxis Anwendungsbeispiel \tiny \setstretch{0.9} \vspace{2mm} ```{r, echo = T} # Daten einlesen fname = "Daten/Einfaktorielle_Varianzanalyse_Daten.csv" D = read.table(fname, sep = ",", header = TRUE) # Datengruppen y_1 = D$dBDI[D$Setting == "F2F"] # BDI-Differenzwerte in der F2F-Gruppe y_2 = D$dBDI[D$Setting == "ONL"] # BDI-Differenzwerte in der ONL-Gruppe y_3 = D$dBDI[D$Setting == "WLC"] # BDI-Differenzwerte in der WLC-Gruppe # Modellformulierung p = 3 # p = 3; drei Gruppen m = length(y_1) # m = 40; balancierters Design n = p*m # n = 120; Datenvektordimension y = matrix(c(y_1, y_2, y_3), nrow = n) # Datenvektor Xt = cbind(matrix(1, nrow = n, ncol = 1), # Designmatrix vollständiges Modell kronecker(diag(p), matrix(1, nrow = m, ncol = 1))) X = Xt[,-2] # eliminiere zweite Spalte X_0 = X[,1] # Designmatrix reduziertes Modell # Evaluation der F-Teststatistik beta_hat = solve(t(X) %*% X) %*% t(X) %*% y # Betaparameterschätzer vollständiges Modell beta_hat_0 = solve(t(X_0) %*% X_0) %*% t(X_0) %*% y # Betaparameterschätzer reduziertes Modell eps_hat = y - X %*% beta_hat # Residuenvektor vollständiges Modell eps_hat_0 = y - X_0 %*% beta_hat_0 # Residuenvektor reduziertes Modell SQT = t(eps_hat_0) %*% eps_hat_0 # total sum of squares SQW = t(eps_hat) %*% eps_hat # within sum of squares SQB = SQT - SQW # between sum of squares DFB = p - 1 # between degrees of freedom DFW = n - p # within degrees of freedom MSB = SQB/DFB # mean between sum of squares MSW = SQW/DFW # mean within sum of squares Eff = MSB/MSW # F-Teststatistik pval = 1 - pf(Eff, p-1, n-p) # p-Wert ``` # Anwendung/Praxis Anwendungsbeispiel \small \vspace{2mm} ```{r, echo = T} # Ausgabe cat( "DFB :", DFB, "\nDFW :", DFW, "\nSQB :", SQB, "\nSQW :", SQW, "\nMSB :", MSB, "\nMSW :", MSW, "\nF :", Eff, "\np :", paste(pval)) ``` # Anwendung/Praxis Anwendungsbeispiel \small \vspace{2mm} ```{r, echo = T} # Daten einlesen fname = "Daten/Einfaktorielle_Varianzanalyse_Daten.csv" D = read.table(fname, sep = ",", header = TRUE) # Benutzung von R's aov Funktion und Ausgabe res.aov = aov(D$dBDI ~ D$Setting, data = D) summary(res.aov) ``` \vfill # Anwendung/Praxis \small Datensatz: Zielzeiten beim Great 10k (08.10.2017), getrennt nach Altersklasse ($n = 4883$) ```{r, echo = F, eval = T} # Daten einlesen filename = "Daten/great_10k.csv" D = read.csv(filename) # Daten reduzieren excl = c("FU14", "FU16", "FU18", "FU20", "FU23", "F70", "F75", "F80", "MU14", "MU16", "MU18", "MU20", "MU23", "M70", "M75", "M80") incl = c("FE", "FH", "F30", "F35", "F40", "F45", "F50", "F55", "F60", "F65", "ME", "MH", "M30", "M35", "M40", "M45", "M50", "M55", "M60", "M65") D = D[!(D$AgeGroup %in% excl),] D = D[!(D$X5kTime == "NaN:NaN:NaN"),] D = D[!(D$X10kTime == "NaN:NaN:NaN"),] n = nrow(D) # Daten umsortieren D$AgeGroup = factor(D$AgeGroup, levels = incl) # Daten umrechnen t5_hms = D$X5kTime t5_min = rep(0, n) t10_hms = D$X10kTime t10_min = rep(0, n) for(i in 1:n){ t5 = t5_hms[i] t10 = t10_hms[i] t5_min[i] = as.numeric(substr(t5,1,2))*60 + as.numeric(substr(t5,4,5)) + as.numeric(substr(t5,7,8))/60 t10_min[i] = as.numeric(substr(t10,1,2))*60 + as.numeric(substr(t10,4,5)) + as.numeric(substr(t10,7,8))/60 } D$T5k1 = t5_min D$T5k2 = t10_min-t5_min D$T10k = t10_min D$Tdiff = D$T5k2 - D$T5k1 ``` ```{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$T10k ~ D$AgeGroup, D, col = "gray80", rectCol = "black", lineCol = "white", colMed = "black", border = "black", pchMed = 16, plotCentre = "lines", xlab = "Altersklasse", ylab = "Zielzeit [min]", ylim = c(25,85), drawRect = FALSE) # Datenpunkte stripchart(D$T10k ~ D$AgeGroup, D, method = "jitter", xaxt = "n", vertical = TRUE, pch = 19, col = "black", add = TRUE, cex = 0.4) # Mittelwerte dx = 0.3 for (i in 1:length(incl)) { segments( x0 = i - dx, x1 = i + dx, y0 = median(D$T10k[D$AgeGroup==incl[i]]), y1 = median(D$T10k[D$AgeGroup==incl[i]]), col = "black", lwd = 2) } # Speichern dev.copy2pdf( file = "Abbildungen/great_10k.pdf", width = 16, height = 9) ``` \vspace{2mm} ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/great_10k.pdf") ``` # Anwendung/Praxis \small Einfaktorielle Varianzanalyse: 10-km-Zeit, alle Altersgruppen ($n = 4883$) \vspace{2mm} \footnotesize ```{r, echo = T, eval = T} # Daten extrahieren D$T10k = t10_min # 10-km-Zeit [min] # Einfaktorielle Varianzanalyse (falsch) grps = c("FE", "FH", "F30", "F35", "F40", "F45", "F50", "F55", "F60", "F65", "ME", "MH", "M30", "M35", "M40", "M45", "M50", "M55", "M60", "M65") p = length(grps) # Anzahl Gruppen n_i = rep(0,p) # Anzahl Datenpunkte pro Gruppe y = matrix(rep(1,0)) # Initialisierung Datenvektor for (i in 1:p) { y_i = matrix(D$T10k[D$AgeGroup == grps[i]], ncol = 1) n_i[i] = length(y_i) # Anzahl Datenpunkte Gruppe i y = rbind(y, y_i) # Ergänzung Datenvektor } n = sum(n_i) # Anzahl Datenpunkte ``` \scriptsize ```{r, echo = F, eval = T} # Einfaktorielle Varianzanalyse X = matrix(rep(1,n), nrow = n) # Initialisierung Designmatrix for (i in 2:p) { n_i1 = sum(n_i[1:(i-1)]) n_i2 = n_i[i] # Anzahl Datenpunkte Gruppe i n_i3 = sum(n_i[(i+1):p]) if (i == p) { n_i3 = 0 } # Ergänzung Designmatrix X = cbind(X, matrix(c(rep(0,n_i1), rep(1,n_i2), rep(0,n_i3)), nrow = n)) } X_0 = X[,1] # Designmatrix reduziertes Modell # Evaluation der F-Teststatistik beta_hat = solve(t(X) %*% X) %*% t(X) %*% y # Betaparameterschätzer vollständiges Modell beta_hat_0 = solve(t(X_0) %*% X_0) %*% t(X_0) %*% y # Betaparameterschätzer reduziertes Modell eps_hat = y - X %*% beta_hat # Residuenvektor vollständiges Modell eps_hat_0 = y - X_0 %*% beta_hat_0 # Residuenvektor reduziertes Modell SQT = t(eps_hat_0) %*% eps_hat_0 # total sum of squares SQW = t(eps_hat) %*% eps_hat # within sum of squares SQB = SQT - SQW # between sum of squares DFB = p - 1 # between degrees of freedom DFW = n - p # within degrees of freedom MSB = SQB/DFB # mean between sum of squares MSW = SQW/DFW # mean within sum of squares Eff = MSB/MSW # F-Teststatistik pval = 1 - pf(Eff, p-1, n-p) # p-Wert # Ausgabe der Ergebnisse cat( "Datenpunkte pro Gruppe : ", n_i, "\nBetaparameterschätzer : ", round(beta_hat, digits = 1), "\nBetaparameterschätzer red. Mod. : ", round(beta_hat_0, digits = 3), "\nresiduelle Quadratsumme : ", round(SQW, digits = 3), "\nresiduelle Quadratsumme red. Mod. : ", round(SQT, digits = 3), "\nF-Teststatistik : ", round(Eff, digits = 3), "\np-Wert : ", round(pval, digits = 3), "\n\n") ``` # Anwendung/Praxis \small Einfaktorielle Varianzanalyse: 10-km-Zeit, Altersgruppen F30/F35/F40 ($n = 936$) \vspace{2mm} ```{r, echo = T, eval = T} # Daten extrahieren D$T10k = t10_min # 10-km-Zeit [min] # Einfaktorielle Varianzanalyse y1 = D$T10k[D$AgeGroup == "F30"] # Datenpunkte Gruppe 1 y2 = D$T10k[D$AgeGroup == "F35"] # Datenpunkte Gruppe 2 y3 = D$T10k[D$AgeGroup == "F40"] # Datenpunkte Gruppe 3 n = length(y1) + length(y2) + length(y3) # Anzahl Datenpunkte p = 3 # Anzahl Gruppen y = matrix(c(y1, y2, y3), nrow = n) # Datenvektor ``` ```{r, echo = F, eval = T} # Einfaktorielle Varianzanalyse n1 = length(y1) # Anzahl Datenpunkte Gr. 1 n2 = length(y2) # Anzahl Datenpunkte Gr. 2 n3 = length(y3) # Anzahl Datenpunkte Gr. 3 X = matrix(c(rep(1,n), # Designmatrix vollständiges Modell c(rep(0,n1), rep(1,n2), rep(0,n3)), c(rep(0,n1), rep(0,n2), rep(1,n3))), nrow = n) X_0 = X[,1] # Designmatrix reduziertes Modell # Evaluation der F-Teststatistik beta_hat = solve(t(X) %*% X) %*% t(X) %*% y # Betaparameterschätzer vollständiges Modell beta_hat_0 = solve(t(X_0) %*% X_0) %*% t(X_0) %*% y # Betaparameterschätzer reduziertes Modell eps_hat = y - X %*% beta_hat # Residuenvektor vollständiges Modell eps_hat_0 = y - X_0 %*% beta_hat_0 # Residuenvektor reduziertes Modell SQT = t(eps_hat_0) %*% eps_hat_0 # total sum of squares SQW = t(eps_hat) %*% eps_hat # within sum of squares SQB = SQT - SQW # between sum of squares DFB = p - 1 # between degrees of freedom DFW = n - p # within degrees of freedom MSB = SQB/DFB # mean between sum of squares MSW = SQW/DFW # mean within sum of squares Eff = MSB/MSW # F-Teststatistik pval = 1 - pf(Eff, p-1, n-p) # p-Wert # Ausgabe der Ergebnisse cat( "Datenpunkte pro Gruppe : ", c(n1,n2,n3), "\nBetaparameterschätzer : ", round(beta_hat, digits = 3), "\nBetaparameterschätzer red. Mod. : ", round(beta_hat_0, digits = 3), "\nresiduelle Quadratsumme : ", round(SQW, digits = 3), "\nresiduelle Quadratsumme red. Mod. : ", round(SQT, digits = 3), "\nF-Teststatistik : ", round(Eff, digits = 3), "\np-Wert : ", round(pval, digits = 3), "\n\n") ``` # Anwendung/Praxis \small Einfaktorielle Varianzanalyse: 10-km-Zeit, Altersgruppen M30/M35/M40 ($n = 1323$) \vspace{2mm} ```{r, echo = T, eval = T} # Daten extrahieren D$T10k = t10_min # 10-km-Zeit [min] # Einfaktorielle Varianzanalyse y1 = D$T10k[D$AgeGroup == "M30"] # Datenpunkte Gruppe 1 y2 = D$T10k[D$AgeGroup == "M35"] # Datenpunkte Gruppe 2 y3 = D$T10k[D$AgeGroup == "M40"] # Datenpunkte Gruppe 3 n = length(y1) + length(y2) + length(y3) # Anzahl Datenpunkte p = 3 # Anzahl Gruppen y = matrix(c(y1, y2, y3), nrow = n) # Datenvektor ``` ```{r, echo = F, eval = T} # Einfaktorielle Varianzanalyse n1 = length(y1) # Anzahl Datenpunkte Gr. 1 n2 = length(y2) # Anzahl Datenpunkte Gr. 2 n3 = length(y3) # Anzahl Datenpunkte Gr. 3 X = matrix(c(rep(1,n), # Designmatrix vollständiges Modell c(rep(0,n1), rep(1,n2), rep(0,n3)), c(rep(0,n1), rep(0,n2), rep(1,n3))), nrow = n) X_0 = X[,1] # Designmatrix reduziertes Modell # Evaluation der F-Teststatistik beta_hat = solve(t(X) %*% X) %*% t(X) %*% y # Betaparameterschätzer vollständiges Modell beta_hat_0 = solve(t(X_0) %*% X_0) %*% t(X_0) %*% y # Betaparameterschätzer reduziertes Modell eps_hat = y - X %*% beta_hat # Residuenvektor vollständiges Modell eps_hat_0 = y - X_0 %*% beta_hat_0 # Residuenvektor reduziertes Modell SQT = t(eps_hat_0) %*% eps_hat_0 # total sum of squares SQW = t(eps_hat) %*% eps_hat # within sum of squares SQB = SQT - SQW # between sum of squares DFB = p - 1 # between degrees of freedom DFW = n - p # within degrees of freedom MSB = SQB/DFB # mean between sum of squares MSW = SQW/DFW # mean within sum of squares Eff = MSB/MSW # F-Teststatistik pval = 1 - pf(Eff, p-1, n-p) # p-Wert # Ausgabe der Ergebnisse cat( "Datenpunkte pro Gruppe : ", c(n1,n2,n3), "\nBetaparameterschätzer : ", round(beta_hat, digits = 3), "\nBetaparameterschätzer red. Mod. : ", round(beta_hat_0, digits = 3), "\nresiduelle Quadratsumme : ", round(SQW, digits = 3), "\nresiduelle Quadratsumme red. Mod. : ", round(SQT, digits = 3), "\nF-Teststatistik : ", round(Eff, digits = 3), "\np-Wert : ", round(pval, digits = 3), "\n\n") ``` # \large \setstretch{2.5} \vfill Anwendungsszenario Modellformulierung Modellschätzung Modellevaluation Anwendung/Praxis **Selbstkontrollfragen** \vfill # Selbstkontrollfragen \vspace{2mm} \footnotesize \setstretch{1.5} \begin{enumerate} \item Erläutern Sie das Anwendungsszenario einer einfaktoriellen Varianzanalyse (EVA). \item Geben Sie die Definition des EVA-Modells in Erwartungswertparameterdarstellung wieder. \item Geben Sie die strukturelle Form des EVA-Modells in Effektdarstellung wieder. \item Erläutern Sie die Motivation für die Reparameterisierung des EVA-Modells. \item Welche Bedeutung haben die Parameter $\mu_0,\alpha_2,...,\alpha_p$ in der Effektparameterdarstellung des EVA-Modells? \item Warum gibt es bei $p$ Gruppen eines EVA-Szenarios nur die $p-1$ Effektparameter $\alpha_2,...,\alpha_p$? \item Geben Sie die Designmatrixform des EVA-Modells in Effektdarstellung wieder. \item Formulieren Sie die Designmatrix eines EVA-Modells mit $n_i = 3$ und $p = 2$. \item Formulieren Sie die Designmatrix eines EVA-Modells mit $n_i = 2$ und $p = 5$. \item Geben Sie das Theorem zur Betaparameterschätzung im EVA-Modell wieder. \item Mit welchen deskriptiven Statistiken werden die Parameter $\mu_0,\alpha_2,...,\alpha_p$ eines EVA-Modells geschätzt? \end{enumerate} \vfill # Selbstkontrollfragen \vspace{2mm} \footnotesize \setstretch{1.5} \begin{enumerate} \setcounter{enumi}{11} \item Geben Sie das Theorem zur Quadratsummenzerlegung bei einfaktorieller Varianzanalyse wieder. \item Erläutern Sie die Begriffe total sum of squares, between sum of squares und within sum of squares. \item Geben Sie die Definition des Effektstärkenmaßes $\eta^2$ an. \item Wann nimmt das Effektstärkenmaß $\eta^2$ der EVA seinen Minimalwert an und wie lautet dieser? \item Wann nimmt das Effektstärkenmaß $\eta^2$ der EVA seinen Maximalwert an und wie lautet dieser? \item Geben Sie das Theorem zur F-Teststatistik der EVA wieder. \item Erläutern Sie die Begriffe mean between sum of squares und mean within sum of squares. \item Geben Sie das Theorem zum Zusammenhang von Effektstärkenmaß $\eta^2$ und F-Teststatistik der EVA wieder. \item Geben Sie die Definition des EVA-F-Test wieder. \item Erläutern sie die Null- und Alternativhypothesen des EVA-F-Tests. \item Geben Sie das Theorem zur Testumfangkontrolle der EVA wieder. \item Skizzieren Sie den Beweis zur Testumfangkontrolle der EVA. \item Geben Sie den p-Wert zum F-Test der EVA wieder. \end{enumerate} \vfill # Referenzen