Bisher waren wir damit beschäftigt unsere Modelle Stück für Stück immer komplizierter zu machen. Angefangen haben wir mit dem einfachen linearen Modell, das dann zum additiven multiplen linearen Modell mit mehreren \(X\)-Variablen wurde. Die \(X\)-Variablen konnten im nächsten Schritt miteinander interagieren, während im letzten Schritt die Anforderung aufgehoben wurde, dass die \(x\)-Variablen kontinuierlich sein mussten. Mit Hilfe von Indexvariablen können wir nun auch nominale Variablen modellieren. Im Kern wurde aber letztendlich immer das einfache Modell der Punkt-Steigungsform aus der Schule beibehalten. Im folgenden Abschnitt werden wir nun eine direkte Verbindung zwischen dem Regressionsmodell und der Varianzanalyse erarbeiten.
Wir beginnen wieder mit dem Modell der einfachen linearen Regression, mit nur einer \(X\)- und einer \(Y\)-Variable.
mod0 <- lm(y ~ x, simple)
summary(mod0)
Call:
lm(formula = y ~ x, data = simple)
Residuals:
1 2 3 4
-0.5817 0.9898 -0.2345 -0.1736
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.8414 0.7008 2.628 0.119
x 0.4574 0.3746 1.221 0.346
Residual standard error: 0.8376 on 2 degrees of freedom
Multiple R-squared: 0.4271, Adjusted R-squared: 0.1406
F-statistic: 1.491 on 1 and 2 DF, p-value: 0.3465Schauen wir uns noch einmal genauer die Residuen an, d.h. die Abweichungen der beobachteten Daten von der Regressionsgeraden. Bei der Besprechung des Determinationskoeffizienten \(R^2\) haben wir schon Quadratsummen und deren Unterteilung kennengelernt. Dort hatten wird die Aufteilung der Varianz von \(Y\), bezeichnet als \(SSTO\), die totale Varianz in die beiden Komponenten Regressionsvarianz \(SSR\) und Fehlervarianz \(SSE\) besprochen.
\[\begin{equation} SSTO = SSR + SSE \end{equation}\]
Die Fehlerquadratsumme SSE, die Summe der quadrierten Abweichungen zwischen dem beobachteten Wert \(y_i\) und dem vorhergesagten Wert \(\hat{y}_i\) ist definiert mittels:
\[\begin{equation} SSE = \sum_{i=1}^n (y_i - \hat{y}_i)^2 \end{equation}\]
Also letztendlich wieder nur die Summe der quadrierten Residuen.
Um für gegebene Daten die vorhergesagten Werte \(\hat{y}_i\) berechnen zu können, benötigen wir das Modell und die dazugehörigen Modellkoeffizienten \(\beta_0\) und \(\beta_1\). Das einfache lineare Modell hat in diesem Fall zwei Parameter. Streng genommen besitzt das Modell noch einen weiteren Parameter nämlich \(\sigma^2\). Diesen lassen wir aber erst wieder links liegen, da \(\sigma^2\) auch für die Berechnung der \(\hat{y}_i\) nicht unbedingt notwendig ist. Formalisierung wir nun die Parameteranzahl indem ihr eine eigene Variable spendieren. Per Konvention wird die Parameteranzahl mit \(p\) bezeichnet. In unserem einfachen Fall mit den beiden Parametern \(\beta_0\) und \(\beta_1\) gilt daher \(p = 2\).
Die Anzahl der Parameter \(p\) spielt eine Rolle wenn die sogenannten Freiheitsgrade \(df\) (degrees of freedom) bestimmt werden müssen. Umgangssprach bezeichnen die Freiheitsgerade die Anzahl der Daten die frei variiert werden können. Die tatsächliche Definition ist komplizierter und auch nicht immer eindeutig. Die Freiheitsgrade von \(SSE\) berechnen sich mittels der Formel \(N-p\). \(N\) bezeichnet wie immer die Anzahl der Beobachtungen, die Anzahl der Datenpunkte.
\[\begin{equation} df_E := N - p \end{equation}\]
Die Freiheitsgerade \(df\) bestimmen die effektive Anzahl der Beobachtungen die zur Verfügung stehen um die Varianz \(\sigma^2\) des Modells abzuschätzen. Dadurch, dass zwei Parameter anhand der Daten für das Modell bestimt werden, fallen zwei Datenpunkt als unabhängige Informationsquellen weg. Anders ausgedrückt, wenn ich die beiden Modellparameter \(\beta_0\) und \(\beta_1\) kenne, dann sind nur noch \(N-2\) Datenpunkt frei variierbar. Sobald ich die Werte von \(N-2\) Datenpunkten und die beiden Parameter kenne, kann ich die verbleibenden beiden Werte berechnen. Daher der Begriff der Freiheitsgrade.
Wenn \(SSE\) durch die Anzahl der Freiheitsgerade geteilt wird, dann lässt sich zeigen, das dieser Wert ein erwartungstreuer Schätzer für die Residualvarianz \(\sigma^2\) unter der Verteilungsannahme \(\epsilon_i \sim \mathcal{N}(0,\sigma^2)\) der Daten ist. Das Verhältnis von \(SSE\) zu \(df\) wird als Mean squared error (\(MSE\)) bezeichnet.
\[\begin{equation} MSE = \frac{SSE}{df_{E}} = \frac{\sum_{i=1}^n (y_i - \hat{y}_i)^2}{N-p} = \hat{\sigma}^2 \label{eq-mlm-hier-MSE} \end{equation}\]
Im Allgemeinen bedeutet groß \(SS\) immer, das es sich um die Summe von quadrierten werden handelt (\(S\)um of \(S\)squares) während ein \(M\) bedeutet das die Quadratsumme durch ihre Freiheitsgrade beteilt wurde (\(M\)ean).
Tatsächlich ist der Ansatzu in Formel \(\eqref{eq-mlm-hier-MSE}\) schon bekannt uns ist uns bei der Stichprobenvarianz \(s^2\) schon begegnet. Wenn wir eine Stichprobe der Größe \(N\) mit Werten \(y_i\) haben, dann haben wir die Stichprobenvarianz mittels der uns bekannten Formel berechnet:
\[\begin{equation} \hat{\sigma}^2 = s^2 = \frac{1}{N-1}\sum_{i=1}^2(y_i - \bar{y})^2 \end{equation}\]
Was bei dieser Formel vielleicht schon immer etwas undurchsichtig gewesen ist, ist der Nenner mit \(N-1\) anstatt einfach \(N\) wie wir es vom Mittelwert kennen. Hier kommen wieder die Freiheitsgerade zum Einsatz. Um die Varianz überhaupt berechnen zu können benötigen wir einen Parameter, den Mittelwert \(\hat{y}\) der \(y_i\). Den Mittelwert berechnen wir aber anhand der Daten. Dies führt dazu, dass nach Berechnung von \(\hat{y}\) wiederum nur noch \(N-1\) Werte frei variiert werden können. Sobald wir, neben dem Mittelwert \(\bar{y}\), \(N-1\) Werte kennen, können wir den verbleibenden \(N\)-ten Wert berechnen. Durch die Berechnung des Mittelwerts ist also ein Freiheitsgerad verloren gegangen. Der Zähler bei der Varianz ist wieder eine Summe quadrierter Abweichungen also ein \(SS\) und damit wird die Varianz \(s^2\) insgesamt zu einer \(MS\).
Nach dieser Wiederholung kommen wie wieder zur Regression zurück. Das Ziel ist eine Metrik zu entwickeln die es erlaubt abschätzen, ob die Hinzunahme von zusätzlichen Modellparametern zu einer Verbesserung der Modellvorhersage führt. D.h. wir starten mit einem Modell und fügen dann weitere Variablen dazu um zu überprüfen ob dadurch die Modellvorhersage verbessert wird. Eine Verbesserung versuchen wir über die Veränderung in \(SSE\) zu bestimmen. Wenn zusätzliche Modellparameter zu einer Verbesserung führen, dann sollte eine relevante Reduktion der Fehlervarianz \(SSE\) beobachtet werden. Die Residuen werden kleiner und wir können präzisere Vorhersagen treffen. Als Randbedinungen fordern wir, dass die zu vergleichenden Modelle in einer Hierarchie zueinander stehen. Dies bedeutet, dass einfachere Modelle als Teilmodelle von komplexeren Modellen interpretiert werden.
Wir müssen dazu zunächst die Unterscheidung in ein volles Modell (\(F\)ull model) und ein reduziertes Modell (\(R\)educed model) verstehen. Als das volle Modell bezeichnen wir bei einem Vergleich immer das komplizierte Modell. Dementsprechend ist das reduzierte Modell das weniger komplizierte Modell welches einem Spezialfall des vollen Modell entspricht. Dabei ist die Bezeichnung voll bzw. reduziertes Modell immer nur temporär, ein volles Modell kann bei einem weiteren Vergleich zum reduzierten Modell werden. Im weiteren wird bei der Bezeichnung der Modell die englische Bezeichung (full vs. reduced) verwendet. Für das Beispiel der einfachen lineare Regression ist das full model das uns schon bekannte:
\[ Y_i = \beta_0 + \beta_1 X_i + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0,\sigma^2) \]
Um die Residualvarianzen \(SSE\) für die beiden Modell voneinander unterscheiden zu können bezeichnen wir diese mit \(SSE(F)\) für das full model und \(SSE(R)\) für das reduced model. Die Residualvarianz \(SSE(F)\) berechnet sich wie oben wiederholt mittels:
\[\begin{equation} \textrm{SSE(F)} = \sum_{i=1}^N(y_i - \hat{y}_i)^2 = \sum_{i=1}^N [y_i - (\beta_0 + \beta_1 x_i)]^2 \label{eq-mlm-hier-ssef} \end{equation}\]
Da wir \(p = 2\) Modellparameter haben, hat das Modell \(df_{F} = N - p = N - 2\) Freiheitsgrade. Wir könnten uns jetzt die Frage stellen, ob wir wirklich den Modellparameter \(\beta_1\) benötigen. Vielleicht zeigt die \(X\)-Variable gar keinen Zusammenhang mit der \(Y\)-Variable und wir fitten nur Datenrauschen mit dem Modell. Aus dieser Überlegung heraus, können wir jetzt ein reduziertes Modell formulieren bei dem der Parameter \(\beta_1\) fehlt.
\[ Y_i = \beta_0 + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0,\sigma^2) \]
Die Residualvarianz \(SSE(R)\) für dieses reduzierte Modell berechnet sich nun mittels.
\[ \textrm{SSE(R)} = \sum_{i=1}^N (y_i - \beta_0)^2 \]
Im reduzierten Modell ist nur noch ein Paramter \(\beta_0\) vorhanden, somit gilt \(p = 1\) und entsprechend \(df_{R} = N - 1\).
Mit etwas Algebra lässt sich zeigen, dass im Allgemeinen \(SSE(F) \leq SSE(R)\) gilt. Dieser Zusammenhang lässt sich auch heuristisch herleiten. Wenn es keinen Zusammenhang zwischen \(X\) und \(Y\) gibt, dann wird der \(\beta_1\) im full model nahezu \(0\) sein und Formel \(\eqref{eq-mlm-hier-ssef}\) wird zu \(SSE(R)\). Im realistischen Fall wird aber, selbst wenn kein Zusammenhang besteht zwischen \(X\) und \(Y\), ein Teil des immer vorhandenen Rauschens in den Daten mittels \(\beta_1\) gefittet. Das führt dazu, das praktisch immer \(SSE(F)\) etwas kleiner ist als \(SSE(R)\). Somit folgt der beschriebene Zusammenhang zwischen \(SSE(F)\) und \(SSE(R)\).
Zwischen der Residualvarianz im reduzierten Modell SSE(R), dem optimalen Modellparameter \(\beta_0\) und der Stichprobenvarianz besteht ein enger Zusammenhang bzw. Identität wie sich anhand der folgenden Herleitung sehen lässt. Wir wollen den Modellparameter unter der Minimierung der Summer der Quadrate der Abweichung, sprich SSE, ermitteln.
\[\begin{align*} SSE &= \sum_{i=1}^n(y_i - \beta_0)^2 \\ &= \sum_{i=1}^n (y_i^2 - 2y_i\beta_0 + \beta_0^2) \end{align*}\]
Wir bestimmen \(min[SSE]\), wie wir das schon vorher beim einfachen Regressionmodell gemacht haben. Wir setzen den Term \(=0\) und leiten nach dem Modellparameter \(\beta_0\) ab. Ein bisschen Algebra führt zu:
\[\begin{align*} 0 &= \frac{\mathrm{d}}{\mathrm{d} \beta_0}\sum_{i=1}^n (y_i^2 - 2y_i\beta_0 + \beta_0^2) \\ 0 &= \sum_{i=1}^n (-2y_i + 2 \beta_0) = -2\sum_{i=1}^n y_i + 2\sum_{i=1}^n \beta_0\\ n\beta_0 &= \sum_{i=1}^n y_i \\ \beta_0 &= \frac{\sum_{i=1}^n y_i}{n} \end{align*}\]
Somit ist Wert von \(\beta_0\) der Abweichungen minimiert, der Mittelwert \(\bar{y}\). Daraus folgt, das unsere Schätzer für \(\sigma^2\),
\[\begin{equation*} \hat{\sigma}^2 = \frac{SSE}{n-p} = \frac{SSE(R)}{n-1} = \sum_{i=1}^n (y_i - \bar{y})^2 = s^2 \end{equation*}\]
einfach nur unser bekannter Schätzer der Stichprobenvarianz \(s^2\) ist. Somit gilt für das reduced model:
\[ \textrm{SSE(R)} = \sum_{i=1}^N (y_i - \beta_0)^2 = \sum_{i=1}^N(y_i - \bar{y})^2 = \textrm{SSTO} \] Kommen wir zurück zur Entwicklung unsere Metrik um das volle und das reduzierte Modell miteinander zu vergleichen. Gehen wir davon aus, das das reduzierte Modell korrekt ist. D.h. die Hinzunahme von \(X\) sollte keine Verbesserung des Modells, keine Verbesserung in der Modellvorhersage, nach sich ziehen. Konkret bedeutet dies, dass \(SSE(R)\) und \(SSE(F)\) in etwas gleich groß sein sollten, bzw. \(SSE(F)\) nur wenig besser (also kleiner) ist als \(SSE(R)\) ist. Wenn wir die Differenz der beiden Quadratsummen berechnen, dann sollte die Differenz entsprechend eher klein sein.
\[\begin{equation} \textrm{SSE(R)} - \textrm{SSE(F)} \label{eq-mlm-hier-div} \end{equation}\]
Beide Modelle sind in etwas gleich gut, fitten die Daten also in etwa gleich gut. Das volle Modell etwas besser, das es durch den zusätzlichen Parameter etwas flexibler als das reduzierte Modell ist.
Gehen wir von nun von der entgegengesetzen Annahme aus. Das reduzierte Modell ist falsch und wir benötigen die Variable \(X\) um die Varianz in \(Y\) aufzuklären. In diesem Fall sollte die Differenz \(\eqref{eq-mlm-hier-div}\) einen deutlich größeren Wert annehmen. Das reduzierte Modell kann denjenigen Teil der Varianz von \(Y\) nicht aufklären der durch \(X\) entsteht. Dadurch wandert die durch \(X\) verursachte Varianz in die Residuen, was dazu führt das \(SSE(R)\) größer wird. Im vollem Modell dagegen, kann die durch \(X\) entstehende Varianz durch den zusätzlichen Modellparameter \(\beta_1\) erklärt werden und entsprechend sind die Residuen geringer was wiederum zu einem kleineren \(SSE(F)\) führt.
Zusammengefasst haben wir heuristisch eine Metrik hergeleitet, die uns erlaubt verschiedene Modell miteinander zu vergleichen. Wenn das reduzierte, das einfachere, Modell ausreicht um die Daten zu erklären, dann wird die Differenz \(\eqref{eq-mlm-hier-div}\) eher klein ausfallen. Wenn dagegen die zusätzlichen Parameter im vollem Modell benötigt werden und die Varianz in \(Y\) zu erklären, dann wird der Unterschied \(\eqref{eq-mlm-hier-div}\) eher groß werden.
Wir bringen nun noch einen zusätzlichen Parameter in unseren Modellvergleich ein. Die Bedeutsamkeit des Unterschieds zwischen den beiden Modellen ist auch noch abhängig davon in wie vielen Parametern die beiden Modelle sich voneinander unterscheiden. Wenn im vollen Modell \(p = 10\) Parameter sind und im reduzierten Modell eben nur \(p = 1\) Parameter ist, dann ist ein gegebener Unterschied in den \(SSE\)s zwischen den Modellen anders zu bewerten, als wenn im vollem Modell \(p = 2\) Parameter geschätzt werden. Bei gleichem Unterschied zwischen den Modellen ist der Unterschied im ersteren Fall weniger bedeutsam im Vergleich zum letzteren Fall. Daher kalibrieren wir den Unterschied zwischen den Modellen noch anhand des Unterschieds in der Anzahl der Parameter. Anders formuliert, wir betrachten die Veränderung in der Varianzverkleinerung pro Freiheitsgrad an. In unserem einfachen Fall passiert da nichts, da der Unterschied in der Parameteranzahl \(= 1\) ist. Wir werden aber später sehen, dass auch Modelle mit größeren Unterschieden in der Parameteranzahl \(p\) miteinander verglichen werden können.
Mit \(p_F\) = Anzahl der Parameter im vollen Modell, \(p_R\) = Anzahl der Parameter im reduzierten Modell gilt:
\[ p_{F} - p_{R} = p_{F} - p_{R} + N - N = N - p_{R} - (N - p_{F}) = df_{R} - df_{F} \]
Nach dem ersten \(=\) haben wir \(0\) in Form von \(+N-N\) addiert umd die Schreibweise mit den Freiheitsgraden zu erhalten. Achtung, die Reihenfolge der Modellparameter ändert sich von \(p\) zu \(df\). Als Merkhilfe, der Unterschied muss positiv sein, d.h. es wird immer der größere Wert vom kleineren Wert abgezogen. Somit schreiben wir den Unterschied zwischen den beiden Modellen folgendermaßen auf und geben dem Term auch noch eine eigene Bezeichnung \(MS_{\textrm{test}}\) für mean squared test.
\[\begin{equation} MS_{\textrm{test}} = \frac{\textrm{SSE(R)} - \textrm{SSE(F)}}{df_{E(R)} - df_{E(F)}} \label{eq-mlm-hier-mstest} \end{equation}\]
Unter der Annahme, das das reduzierte Modell korrekt und den üblichen Modellannahmen im Regressionsmodell zur Normalverteilung der Residuen \(\epsilon_i \sim \mathcal{N}(0,\sigma^2)\), lässt sich zeigen, dass \(MS_{\textrm{test}}\) ein Schätzer für die Varianz \(\sigma^2\) ist. D.h. es gilt:
\[\begin{equation*} MS_{\textrm{test}} = \hat{\sigma}^2 \end{equation*}\]
D.h. wir haben neben dem uns schon bekannten Schätzer \(MSE\) für \(\sigma^2\) noch einen weiteren Schätzer für \(\sigma^2\) erhalten. Aber, wenn das reduzierte Modell korrekt ist, dann ist auch das volle Modell korrekt, da das volle Modell das reduzierte Modell als einen Spezialfall beinhaltet. Der zusätzlich Parameter sollte dann normalerweise in der Nähe von \(0\) sein, da kein Zusammenhang zwischen \(X\) und \(Y\) besteht wenn das reduzierte Modell korrekt ist. Daher ist \(MSE(F)\) ebenfalls ein Schätzer für die Varianz \(\sigma^2\) unter den Modellannahmen \(\epsilon_i \sim \mathcal{N}(0,\sigma^2)\). Letztendlich haben wir diese Eigenschaft schon die ganze Zeit verwendet wenn wir Regressionen berechnet haben.
\[\begin{equation} MS_{E(F)} = \frac{\textrm{SSE(F)}}{df_{E(F)}} = \hat{\sigma}^2 \end{equation}\]
Zusammengenommen haben wir nun zwei Ansätze um \(\sigma^2\) zu schätzen. Einmal über das uns schon bekannte \(MSE\) und nun neu dazu gekommen über \(MS_{\text{test}}\).
Die Beziehung zwischen dem vollen und dem reduzierten Modell kann auch dahingehend interpretiert werden, das das reduzierte Modell das volle Modell mit einer zusätzlichen Randbedigung ist. Wenn das volle Modell die folgende Form hat:
\[\begin{equation*} Y = \beta_0 + \beta_1 \cdot X_1 + \epsilon \end{equation*}\]
Und das reduzierte Modell die Form:
\[\begin{equation*} Y = \beta_0 + \epsilon \end{equation*}\]
Dann kann das reduzierte Modell als das volle Modell mit der Randbedingung das \(\beta_1 = 0\) interpretiert werden.
\[\begin{equation*} Y = \beta_0 + \beta_1 \cdot X_1 + \epsilon, \quad \beta_1 = 0 \end{equation*}\]
Daher läuft ein Vergleich der beiden Modell darauf hinaus zu testen ob der Parameter \(\beta_1 = 0\) ist. Dies erklärt auch noch mal den Satz, das das volle Modell korrekt ist, wenn das reduzierte Modell korrekt ist, da das reduzierte Modell nur ein Spezialfall des vollen Modells ist.
Eine Sache fehlt uns noch um die Größe von \(MS_{\textrm{test}}\) einordnen zu können. Da wir es mit Varianzen zu tun haben, können wir die Größe der Quadratsummen verändern indem wir die Einheiten der abhängigen Variablen \(Y\) verändern. Würden wir z.B. von \([m]\) auf \([cm]\) gehen, da würde sich der Unterschied in Formel \(\eqref{eq-mlm-hier-mstest}\) um den Faktor \(10\times 10=100\) vergrößern ohne das wirklich eine Veränderung in den Modellen stattgefunden hat. Daher kalibrieren wir \(MS_{\textrm{test}}\) indem wir den Term durch \(MS_{E(F)}\) teilen. Damit fallen alle Problem mit Veränderungen durch Änderungen in den Einheiten weg (siehe auch Maxwell, Delaney, und Kelley 2004, p.75).
\[\begin{equation} \frac{MS_{\textrm{test}}}{MS_{E(F)}} = \frac{\frac{\textrm{SSE(R)} - \textrm{SSE(F)}}{df_{E(R)} - df_{E(F)}}}{ \frac{\textrm{SSE(F)}}{df_{E(F)}}} \label{eq-mlm-hier-Ftest} \end{equation}\]
Die dadurch entstehende Metrik hat unter der \(H_0\), das das reduzierte Modell korrekt ist, eine uns bekannte theoretische Verteilung, nämlich die \(F\)-Verteilung mit \(df_{E(R)} - df_{E(F)}\) und \(df_{E(F)}\) Freiheitsgeraden (Warum dies so ist, können wir im Rahmen des Skripts leider nicht herleiten).
\[\begin{equation*} F = \frac{MS_{\textrm{test}}}{MS_{E(F)}} \sim F(df_{E(R)}-df_{E(F)},df_{E(F)}) \end{equation*}\]
Zur Erinnerung sind in Abbildung 20.1 nochmal ein paar Beispiele für \(F\)-Verteilung mit verschiedenen Freiheitsgeraden abgebildet.
Da beide Terme in Formel \(\eqref{eq-mlm-hier-Ftest}\) die Varianz abschätzen deutet ein Wert in der Nähe von \(1\) daraufhin, das das reduzierte Modell adäquat ist um die Daten zu beschreiben und die Hinzunahme des Parameters im vollen Modell keine Verbesserung liefert.
Sobald wir eine bekannte theoretische Verteilung unter einer \(H_0\) haben, können wir unser Hypothestinstrumentarium auf die Verteilung los lassen und entsprechend einen kritischen Bereich für ein gegebenes \(\alpha\) bestimmen (siehe Abbildung 20.2 für ein Beispiel). Wenn der beobachtete \(F\)-Wert in den kritischen Bereich fällt, interpretieren wir dies wie immer als Evidenz gegen die \(H_0\). Wir sind überrascht diesen Wert unter der \(H_0\) zu beobachten und lehen die \(H_0\) das das einfachere Modell korrekt ist ab und werten dies als Evidenz dafür, das das komplexere, volle Modell die Daten besser abbildet und statistisch signifikant mehr Varianz der abhängigen Variable modellieren kann.
Schauen wir uns diesen letzten Schritt noch einmal mit einer Simulation genauer an. Gegeben sein ein DGP der folgenden Form:
\[\begin{equation} y_i = 3 + 2 \cdot x_i + \epsilon_i \quad \epsilon_i \sim \mathcal{N}(0,1) \label{eq-mlm-hier-mstest-sim} \end{equation}\]
D.h. wir haben ein einfache lineares Regressionsmodell mit \(\beta_0 = 3\), \(\beta_1 = 2\) und normalverteilten Residuen mit \(\sigma = 1\). Wir simulieren den DGP \(1000\)mal mit \(N = 30\) Datenpunkten und fitten drei verschiedene Modelle an die Daten. Einmal ein reduziertes Modell ohne \(\beta_1\), ein korrektes Modell mit \(\beta_0\) und \(\beta_1\) und ein überparameterisiertes Modell mit \(\beta_0\), \(\beta_1\) und \(\beta_2\). Als \(X_2\) generieren wir eine Zufallsvariable die der Standardnormalverteilung folgt (D.h. \(X_2 \sim \mathcal{N}(0,1)\). Zwischen \(X_2\) und \(Y\) besteht somit kein Zusammenhang. Daher sollte die Hinzunahme von \(X_2\) zu keiner Verbesserung des Modell führen.
In Abbildung 20.3 ist das Ergebnis der Simulation zu sehen. Für jede der \(1000\) Simulation generieren wir \(n = 30\) Datenpaare \((y_i,x_i)\) nach der Formel \(\eqref{eq-mlm-hier-mstest-sim}\). Wir fitten dann an jeden dieser \(1000\) Datensätze die drei verschiedenen Modelle.
\[\begin{align*} \text{under} &: y_i = \beta_0 + \epsilon_i \\ \text{correct} &: y_i = \beta_0 + \beta_1 \cdot x_i + \epsilon_i \\ \text{over} &: y_i = \beta_0 + \beta_1 \cdot x_{1i} + \beta_2 \cdot x_{2i} + \epsilon_i \end{align*}\]
Für jeden Modellfit ist \(MSE = \hat{\sigma}^2\) berechnet worden. Da wir die Daten simuliert haben, kennen wir den wahren Wert \(\sigma^2 = 1\).
In Abbildung 20.3 ist in rot die tatsächliche Varianz bzw. Standardabweichung des DGP mit \(\sigma = 1\) eingezeichnet. Im reduzierten Modell (under) wird die Residualvarianz klar in allen Durchgängen überschätzt. Da der Modellparameter für \(X\) fehlt, wir diejenige Varianz von \(Y\) die auf \(X\) zurückgeht in die Residualvarianz mit aufgenommen. Im mittleren Graphen, beim korrekten Modell, sind die abgeschätzten Residualvarianzen schön um den tatsächlichn Wert herum verteilt. Im einzelnen Fall kommt es natürlich auf Grund der Stichprobenvariabilität zur Überschätzung bzw. Unterschätzung von \(\sigma^2\). Im Mittel sind die Werte trotzdem korrekt und der Schätzer \(MSE\) ist Erwarungstreu für \(\sigma^2\). Im letzten Fall, für das überparameterisierte Modell (over) wird die Residualvarianz ebenfalls korrekt abgeschätzt. Die Hinzunahme der zufälligen Variable \(X_2\) führt wie erwartet zu keiner Verschlechterung von \(\hat{\sigma}^2\) aber eben auch zu keiner relevanten Verbesserung.
Schauen wir uns als Nächstes an, wie sich der \(F\)-Wert aus Formel\(\eqref{eq-mlm-hier-Ftest}\) verhält.
In Abbildung 20.4 sind die Verteilungen der \(F\)-Werte für die \(1000\) Simulationen einmal für den Vergleich correct vs. under (a) und over vs. correct (b) abgetragen. Schauen wir uns zunächst Abbildung 20.4 (b) an. Wenn wir die beobachtete Verteilung (Balken) mit der theoretischen Verteilung unter der \(H_0\) (rot) vergleichen, ist zu erkennen, dass die beobachteten Werte sehr gut der unter der \(H_0\) erwarteten Verteilung folgen. Der Großteil der Wert liegt in der Umgebung von \(1\). Dies entspricht auch unserer Erwartung, die Hinzunahme der Variable \(X_2\), die keinen Zusammenhang mit \(Y\) hat, führt zu keiner Verbesserung des Modells. Das heißt, dass wir in \(\alpha\)-Prozent die \(H_0\) ablehnen würden und uns dementsprechend irren würden. Im Gegensatz dazu folgt die Verteilung der beobachteten \(F\)-Wert in Abbildung 20.4 (a) nicht einmal annährend der erwarteten unter der \(H_0\). Wir haben in der Umgebung von \(1\) überhaupt keine Werte beobachtet. Die \(H_0\) in diesem Fall ist, das die Hinzunahme von \(X_1\) zu keiner Verbesserung führt. Dementsprechend würde in praktisch allen Fällen die \(H_0\) abgelehnt werden, da der beobachtete \(F\)-Wert deutlich größer als der kritische Wert der \(F\)-Verteilung ist. Dementsprechend würden wir bei der Modellauswahl das korrekte Modell identizieren.
Durch den Vergleich von hierarchisch zueinnander in Beziehung stehenden Modellen, sind wir in der Lage, die Verbesserung/Verschlechterung der Modellvorhersage bei Hinzunahme von Modellvariablen statistisch zu überprüfen. Über die Statistik \(\eqref{eq-mlm-hier-Ftest}\) erhalten wir Teststatistik die einer bekannten theoretischen Verteilung folgt. Wenn dieser \(F\)-Test statistisch signifikant wird, dann werten wir dies als Evidenz dafür, das das volle Modell die Daten so viel besser modelliert, das wir dieses Modell dem reduzierten Modell vorziehen sollten.
Schauen wir uns hergeleitete Metrik in Aktion an einem konkreten Beispiel an. Dabei führen wir nun den Begriff der Modellhierarchien ein. In Abbildung 20.5 ist ein exemplarischer Datensatz abgebildet.
In einer Studie wurde der Zusammenhang zwischen wie gut ein Bonbon bewertet wurde (like, Skala 0-100) und dem Süßegrad und dem Saftanteil untersucht. Wir sehen, dass Bonbons umso besser bewertet werden umso höher der Süßegrad war, aber das dieser Effekt durch den Saftanteil beeinflusst wird und umso stärker ist, umso höher der Saftanteil ist. Daher spricht dies für ein Interaktionsmodell.
Das volle Modell kann dementsprechend mit \(X_1\) = Süßegrad und \(X_2\) = Saftanteil folgendermaßen modelliert werden.
\[ y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{1i}x_{2i} + \epsilon_i \]
Das volle Modell besitzt daher \(p = 4\) Modellparameter und es können drei verschiedene reduzierte Modelle betrachtet werden. Das einfachste Modell bezeichnen wir mit \(m_0\) und dementsprechend ansteigend bis zum vollen Modell \(m_3\).
\[\begin{align*} m_0&: y_i = \beta_0 + \epsilon_i \\ m_1&: y_i = \beta_0 + \beta_1 x_{1i} + \epsilon_i \\ m_2&: y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \epsilon_i \\ m_3&: y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{1i}x_{2i} + \epsilon_i \end{align*}\]
Das Modell \(m_0\) besitzt nur einen \(y\)-Achsenabschnitt \(\beta_0\), der, wie wir oben gesehen haben, zu \(\bar{y}\) wird. Modell \(m_1\) hat einen zusätzlichen Parameter mit einem Steigungskoeffizienten \(\beta_1\) für den Süßegrad, \(m_1\) hat zusätzlich noch einen Parameter \(\beta_2\) für den Saftanteil. Wir können nun diese Modelle in folgender Hierarchie anordnen.
\[\begin{equation*} m_0 \subseteq m_1 \subseteq m_2 \subseteq m_3 \end{equation*}\]
Es gilt ebenfalls
\[\begin{equation*} p_{m_0} = 1 < p_{m_1} = 2< p_{m_2}=3< p_{m_3}=4 \end{equation*}\]
In R können wir die entsprechenden Modelle wie folgt fitten.
mod_0 <- lm(like ~ 1, candy)
mod_1 <- lm(like ~ sweetness, candy)
mod_2 <- lm(like ~ sweetness + moisture, candy)
mod_3 <- lm(like ~ sweetness * moisture, candy)Betrachten wir zunächst den Vergleich zwischen \(m_0\) gegen \(m_1\). D.h. wir wollen \(MS_{\text{test}}\) berechnen bei dem \(m_0\) das reduzierte Modell und \(m_1\) das volle Modell ist.
\[\begin{align*} m_0: y_i &= \beta_0 + \epsilon_i \\ m_1: y_i &= \beta_0 + \beta_1 x_{1i} + \epsilon_i \end{align*}\]
Mit R können wir diesen Vergleich mit der anova()-Funktion machen. Dazu übergeben wird die beiden gefitteten Modelle (tatsächlich können auch mehr Modelle übergeben werden) als Parameter an die Funktion.
anova(mod_0, mod_1)Analysis of Variance Table
Model 1: like ~ 1
Model 2: like ~ sweetness
Res.Df RSS Df Sum of Sq F Pr(>F)
1 77 54237
2 76 34695 1 19542 42.805 6.303e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wir erhalten eine Tabelle. Die Einträge in der ersten Spalte unter Res.Df sind die jeweiligen Freiheitsgrade der beiden Modelle \(df_{E}\). \(m_0\) hat einen Parameter bei \(N=78\) und somit \(78-1=77\) Freiheitsgrade. Entsprechnend hat \(m_0\) mit zwei Parametern \(76\) Freiheitsgrade. Unter RSS sind die jeweiligen Fehlerquadratsummen \(SSE\) der beiden Modelle dokumentiert. Einmal per Hand zur Kontrolle:
SSE_R <- sum(resid(mod_0)**2)
SSE_F <- sum(resid(mod_1)**2)
c(SSE_R, SSE_F)[1] 54236.83 34695.38In der Spalte Df ist der Unterschied in der Anzahl der Modellparameter zwischen den beiden Modellen angeben. In der Spalte Sum of Sq wird die \(MS_{\text{test}}\)-Statistik dokumentiert berechnet als der Unterschied zwischen \(SSE(R) - SSE(F)\).
SSE_R - SSE_F[1] 19541.46Die resultierenden \(F\)-Wert nach Formel \(\eqref{eq-mlm-hier-Ftest}\) folgt in der nächsten Spalte zusammen mit dem p-Wert unter der \(H_0\).
F <- (SSE_R - SSE_F)/(77-76)/(SSE_F/76)
F[1] 42.805431 - pf(F, 77-76, 76)[1] 6.302609e-09In Abbildung 20.6 ist die Verteilung der F-Werte unter der \(H_0\) eingezeichnet zusammen mit dem beobachteten Werte von \(F = 42.805\). Es ist zu erkennen das der beobachtete Werte so weit rechts in der Verteilung ist, dass unter der \(H_0\) es extrem unwahrscheinlich ist einen solchen Wert zu beobachten.
Dementsprechend beobachten wir ein statistisch signifikantes Ergebnis und lehnen daher die \(H_0\) das das einfachere Modell \(m_0\) ausreicht zugunsten des komplexeren Modells ab.
Diesen Ansatz können wir nun einfach so weiter treiben und auf den Vergleich zwischen \(m_1\) gegen \(m_2\) anwenden.
\[\begin{align*} m_1: y_i &= \beta_0 + \beta_1 x_{1i} + \epsilon_i \\ m_2: y_i &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \epsilon_i \end{align*}\]
anova(mod_1, mod_2)Analysis of Variance Table
Model 1: like ~ sweetness
Model 2: like ~ sweetness + moisture
Res.Df RSS Df Sum of Sq F Pr(>F)
1 76 34695
2 75 6459 1 28236 327.86 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wiederum deutet der Test einen darauf, das das komplexere Modell das beide Variablen enthält zu einer Verbesserung der Modellvorhersage führt. Kommen wir zum letzten Vergleich zwischen von \(m_3\) gegen \(m_2\).
\[\begin{align*} m_2: y_i &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \epsilon_i \\ m_3: y_i &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{1i} x_{2i} + \epsilon_i \end{align*}\]
anova(mod_2, mod_3)Analysis of Variance Table
Model 1: like ~ sweetness + moisture
Model 2: like ~ sweetness * moisture
Res.Df RSS Df Sum of Sq F Pr(>F)
1 75 6459.2
2 74 273.5 1 6185.7 1673.8 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wieder beobachten wir ein statistisch signifikantes Ergebnis. Wenn wir Modell miteinander vergleichen, die sich nur um einen Parameter voneinander unterschieden, dann ist der beobachtete \(F\)-Test äquivalent zum \(t\)-Wert den wir unter summary() angezeigt bekommen. Der \(F\)-Wert ist gleich dem quadrierten \(t\)-Wert.
summary(mod_3)
Call:
lm(formula = like ~ sweetness * moisture, data = candy)
Residuals:
Min 1Q Median 3Q Max
-4.9943 -1.1231 -0.1122 1.0661 5.8159
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.029961 0.871354 -1.182 0.240981
sweetness 0.189563 0.047908 3.957 0.000173 ***
moisture 0.168627 0.070880 2.379 0.019937 *
sweetness:moisture 0.155395 0.003798 40.912 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.922 on 74 degrees of freedom
Multiple R-squared: 0.995, Adjusted R-squared: 0.9948
F-statistic: 4867 on 3 and 74 DF, p-value: < 2.2e-16t_mod_3[1] 40.91157t_mod_3**2[1] 1673.756Der Vergleich zwischen hierarchisch miteinander in Beziehung stehenden Modell ist nicht darauf beschränkt das die Modelle sich immer nur um einen zusätzlichen Parameter unterscheiden. Wir können genauso das Modell \(m_3\) mit \(p = 4\) Parametern gegen das \(m_0\) Modell mit \(p = 1\) Parametern vergleichen.
\[\begin{align*} m_0: y_i &= \beta_0 + \epsilon_i \\ m_3: y_i &= \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{1i} x_{2i} + \epsilon_i \end{align*}\]
anova(mod_0, mod_3)Analysis of Variance Table
Model 1: like ~ 1
Model 2: like ~ sweetness * moisture
Res.Df RSS Df Sum of Sq F Pr(>F)
1 77 54237
2 74 273 3 53963 4867.2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In diesem Fall wird getestet ob die Hinzunahme der Parameter \(\beta_1, \beta_2\) und \(\beta_3\) eine statistisch signifikante Verbesserung im Modellfit bedeutet. Entsprechend sind in der Spalte DF jetzt drei Parameterunterschied dokumentiert. Das Prinzip ist aber immer noch das Gleiche. Die Fehlerquadratsummen werden voneinander abgezogen durch den Unterschied in der Parameteranzahl geteilt und mittels der \(MSE(F)\) kalibriert.
Mit diesem Test verstehen wir auch endlich die letzte Zeile im Output von summary(). Tatsächlich wird dieser Test in der Ausgabe von summary() dokumentiert.
summary(mod_3)
Call:
lm(formula = like ~ sweetness * moisture, data = candy)
Residuals:
Min 1Q Median 3Q Max
-4.9943 -1.1231 -0.1122 1.0661 5.8159
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.029961 0.871354 -1.182 0.240981
sweetness 0.189563 0.047908 3.957 0.000173 ***
moisture 0.168627 0.070880 2.379 0.019937 *
sweetness:moisture 0.155395 0.003798 40.912 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.922 on 74 degrees of freedom
Multiple R-squared: 0.995, Adjusted R-squared: 0.9948
F-statistic: 4867 on 3 and 74 DF, p-value: < 2.2e-16D.h. dieser Test prüft immer ob das spezifizierte Modell im Gegensatz zum einfachsten Modell mit nur einem Parameter eine Verbesserung im Modellfit ergibt.
Insgesamt haben wir jetzt einen flexiblen Ansatz gewonnen um Modellparameter zu testen. Bisher haben wir immer nur einen einzelnen Parameter anhand des zum jeweiligen Koeffizienten \(\beta_i\) gehörenden \(t\)-Werts getestet. Die Interpretation war hier die resultierende Steigung der zum Koeffizienten gehörenden Grade. Nun haben wir eine zweite Interpretation in Form der Reduktion der Fehlerquadrate. Der Modellvergleich erlaubt darüber hinaus dasjenige Modell zu bestimmen, welches die Daten am Besten bzw. am einfachsten (mit den wenigstens Parametern) modelliert.
Schauen wir uns als Nächstes an, was passiert wenn wir eine nominale Variable in unser Modell haben und diese durch Indexvariablen darstellen. In Abbildung 20.7 haben wir wieder unser hypothetishces Beispiel mit den Reaktionszeiten unter vier verschiedenen Konditionen A, B, C und D.
Klassischerweise würde wir diese Daten mit einer Varianzanalyse (ANOVA) untersuchen, aber wir haben ja schon gesehen das wir diese Analyse auch mit einem linearen Modell durchführen können. In der Varianzanalyse unterteilen wir die Varianz in die drei Komponenten \(SS_{total}, SS_{between}\) und \(SS_{within}\).
\[\begin{equation*} SS_{total} = SS_{between} + SS_{within} \end{equation*}\]
Mit etwas motiviertem auf die Gleichungstarren fällt uns natürlich auf, das das ziemlich genauso aussieht wie die Regressionszerlegung in \(SSR\) und \(SSE\). Bei der Herleitung der Varianzanalyse sind euch möglicherweise die folgenden Formeln zur Berechnung der \(F\)-Statistik untergekommen (\(K\) = Anzahl der Faktorstufen).
\[\begin{align*} \hat{\sigma}^2_{within} = s_{within}^2 &= \frac{1}{N-K}\sum_{j=1}^K\sum_{i=1}^{N_j}(y_{ji}-\bar{y}_{j.})^2 \\ &= \frac{1}{N-K}\sum_{j=1}^K(N_j-1)s_j^2 \\ \hat{\sigma}^2_{between} = s_{between}^2 &= \frac{1}{K-1}\sum_{j=1}^K N_j (\bar{y}_{j.}-\bar{y})^2 \\ F &= \frac{\hat{\sigma}_{between}^2} {\hat{\sigma}_{within}^2} \sim F(K-1,N-K) \end{align*}\]
Schauen wir uns \(s_{between}^2\) etwas genauer an, dann wird vom beobachteten Wert \(y_{ij}\) der jeweiligen Gruppenmittelwert \(\bar{y}_{j.}\) abgezogen, quadriert, aufsummiert und durch \(N-K\) geteilt. Wenn wir die Daten mittels eines linearen Modells modellieren, dann bildet der \(y\)-Achsenabschnitt den Mittelwert der Referenzgruppe ab und die anderen Koeffizienten \(\beta_i\) den Unterschied zu den anderen Mittelwerten. Wenn die Residuen berechnet werden, dann sind das letztendlich nichts anderes als die Abweichungen von den Gruppenmittelwerten. Im linearen Modell hätten wir dann vier Modellparameter \(p = 4\) dementsprechend ist \(N-K = N-p\) und \(s_{within}^2\) ist nichts anderes als \(SSE(F)\). Wenig überraschend lässt sich nun auch noch zeigen, dass der Term \(s_{between}^2\) auch nichts anderes als \(MS(test)\) ist. Daher ist die Varianzanalyse nichts anderes als ein Vergleich von Modellen miteinandern und geht daher vollkommen in unserem Regressionansatz und letztendlich in der Punkt-Steigungsform aus der 7. Klasse auf.
In R können wir eine ANOVA mit der aov()-Funktion berechnen. Die Modellformulierung ist gleich derjenigen wie sie bei lm() verwendet wird.
mod_aov <- aov(rt ~ group, rt_tbl)Anwendung von summary() auf das gefittete aov-Modell liefert die übliche ANOVA-Tabelle.
Df Sum Sq Mean Sq F value Pr(>F)
group 3 988935 329645 157.3 <2e-16 ***
Residuals 76 159221 2095
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wie sieht das Ganz aus, wenn wir den Ansatz mit Modellhierarchien anwenden? Das vollständige Modell ist dasjenige, bei dem wir für beispielsweise Stufe A als Referenz nehmen und mittels dreier Dummy-Variablen (\(K-1\)) die Abweichungen der Stufen B-D von A modellieren, woraus folgt \(p = 4\).
\[\begin{equation*} y_i = \beta_0 + \beta_{\Delta_{B-A}} x_1 + \beta_{\Delta_{C-A}} x_2 + \beta_{\Delta_{D-A}} x_3 + \epsilon_i \end{equation*}\]
Als reduziertes Modell wählen wir das Modell mit nur einem \(y\)-Achsenabschnitt und daher nur Parameter \(p = 1\).
\[\begin{equation*} y_i = \beta_0 + \epsilon_i \end{equation*}\]
Wenn das reduzierte Modell die Daten gleich gut fittet wie das vollständige Modell dann bedeutet dass, das die zusätzliche Information über die verschiedenen Konditionstufen meine Vorhersage von \(y_i\) nicht verbessert.
mod <- lm(rt ~ group, rt_tbl)
summary(mod)
Call:
lm(formula = rt ~ group, data = rt_tbl)
Residuals:
Min 1Q Median 3Q Max
-120.261 -25.789 -0.857 27.520 114.994
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 509.53 10.23 49.784 < 2e-16 ***
groupB 90.15 14.47 6.228 2.4e-08 ***
groupC 197.41 14.47 13.639 < 2e-16 ***
groupD 295.56 14.47 20.420 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 45.77 on 76 degrees of freedom
Multiple R-squared: 0.8613, Adjusted R-squared: 0.8559
F-statistic: 157.3 on 3 and 76 DF, p-value: < 2.2e-16Wenn wir uns die letzte Zeile im summary() Ausdruck anschauen, dann sehen wir dort die gleichen Werte die wir mit aov() erhalten haben. Wenn wir anova() das gefittete lm Modell übergeben erhalten wir sogar genau die gleiche ANOVA Tabelle.
anova(mod)Analysis of Variance Table
Response: rt
Df Sum Sq Mean Sq F value Pr(>F)
group 3 988935 329645 157.35 < 2.2e-16 ***
Residuals 76 159221 2095
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Tatsächlich verwendet aov() im Hintergrund auch nichts anderes als die lm() Funktion. Insgesamt bedeutet das für uns, das die ANOVA und Regression letztendlich auf das gleiche Modell zurück gehen, nämlich auf das Allgemeine Lineare Modell.
Viele der Ideen die hier diskutiert werden sind aus Christensen (2018, p.57–64). Daher wer noch mal etwas mehr Hintergrund dazu haben möchte, sei diese Quelle wärmstens empfohlen.