Multiple lineare Regression

Lernziele: Am Ende des Kapitels können Sie

• das multiple lineare Regressionsmodell erläutern und Anwendungsbeispiele nennen.
• die Methode der kleinsten Quadrate vorstellen und motivieren.
• die Parameter im multiplen linearen Regressionsmodell schätzen und Konfidenzintervalle bestimmen.
• Daten mithilfe des multiplen linearen Regressionsmodells in R analysieren und die R-Ausgabe interpretieren.

1 Allgemeines lineares Modell: Einleitung

Wir wollen in diesem Abschnitt das allgemeine lineare Modell vorstellen, das sich vom einfachen linearen Regressionsmodell zunächst einmal dadurch unterscheidet, dass es nicht nur eine erklärende Variable gibt, sondern gleich mehrere. Wir werden das Modell zunächst sehr allgemein formulieren und im weiteren Verlauf dieses Abschnitts zeigen, dass sich viele Beispiele als Spezialfälle dieses allgemeinen Modells auffassen lassen.

Wir betrachten ein Experiment, dessen Ergebnis y vom Wert von p reellwertigen erklärenden Variablen x_1,\ldots, x_p sowie vom Zufall abhängt. Wir fassen y als Realisierung einer Zufallsvariablen Y auf, für die wir die Modellannahme machen, dass Y=\beta_1\, x_1 +\beta_2\, x_2 +\ldots +\beta_p\, x_p+\epsilon, \tag{1} wobei \beta_1,\ldots,\beta_p unbekannte Parameter sind und \epsilon jegliche Abweichung von der idealen Beziehung y=\sum_{j=1}^p \beta_j \, x_j beinhaltet. Formal fassen wir \epsilon hier als Zufallsvariable auf. Klassisch denkt man dabei an Messfehler bei einem natur- oder ingenieurwissenschaftlichen Experiment; daher kommt auch die Notation \epsilon (für error). In den meisten Beispielen geht es aber weniger um Messfehler, sondern um zufällige Fluktuationen des Wertes von Y, die nicht durch die Variation der erklärenden Variablen hervorgerufen werden, sondern durch andere Ursachen, die nicht Teil des Modells sind. Wir machen dabei stets die Annahme, dass \operatorname{E}(\epsilon)=0.

Anmerkung.

1. Oft wird das lineare Modell mit einem konstanten Term \beta_0 formuliert. Die Modellgleichung lautet dann Y= \beta_0+\beta_1\, x_1 +\beta_2\, x_2 +\ldots +\beta_p\, x_p+\epsilon. \tag{2} In unserer Darstellung Gleichung 1 des linearen Modells kann man den konstanten Term dadurch integrieren, dass man x_1\equiv 1 setzt. Unsere Formulierung ist daher allgemeiner, weil sie auch Modellgleichungen ohne einen konstanten Term zulässt. Bei der Analyse eines linearen Modells mit dem R-Befehl lm wird immer ein konstanter Term aufgenommen, d.h. es wird die Modellgleichung Gleichung 2 verwendet.

2. Aus der Modellgleichung Gleichung 1 und der Annahme \operatorname{E}(\epsilon)=0 folgt direkt \operatorname{E}(Y)=\beta_1\, x_1 +\beta_2\, x_2 +\ldots +\beta_p\, x_p. Wir können diese Gleichung so interpretieren, dass wir bei Durchführung einer Vielzahl an Messungen im Mittel der y-Werte das Ergebnis \beta_1\, x_1 +\beta_2\, x_2 +\ldots +\beta_p\, x_p erhalten, was der idealen linearen Beziehung zwischen den erklärenden Variablen und y entspricht. Da der Erwartungswert die beste Vorhersage des Ergebnisses ist, erhalten wir bei bekannten Werten der erklärenden Variablen x_1,\ldots,x_p als beste Vorhersage für y \widehat{Y}=\widehat{Y}(x_1,\ldots,x_p) = \beta_1\, x_1 +\beta_2\, x_2 +\ldots +\beta_p\, x_p.

Das Modell Gleichung 1 ist das einfachste Modell zur Beschreibung des Zusammenhangs zwischen den erklärenden Variablen x_1,\ldots,x_p und der abhängigen Variablen y, in dem vom einem linearen Zusammenhang ausgegangen wird. Zur Rechtfertigung für die Verwendung eines derart einfachen Modells kann man vorbringen, dass viele Funktionen in erster Näherung linear sind und entsprechend das lineare Modell die Wirklichkeit in einem begrenzten Bereich von Werten der erklärenden Variablen gut abbildet.

Das Modell Gleichung 1 beschreibt eine einzelne Beobachtung. Im Allgemeinen haben wir n Beobachtungen bei unterschiedlichen Werten der erklärenden Variablen. Den Wert der j-ten erklärenden Variablen bei der i-ten Beobachtung bezeichnen wir mit x_{ij}, 1\leq i\leq n, 1\leq j\leq p. Mit y_i bezeichnen wir den Wert der abhängigen Variablen bei der i-ten Beobachtung, und fassen y_i als Realisierung einer Zufallsvariablen Y_i auf. Entsprechend erhalten wir die Modellgleichungen Y_i=\beta_1\, x_{i1}+\beta_2\, x_{i2} +\ldots +\beta_p\, x_{ip}+\epsilon_i, \; 1\leq i\leq n. \tag{3} Diese n Gleichungen können wir in kompakter Form in Vektor-Matrix Schreibweise darstellen. Wir führen dazu die Bezeichnungen Y=(Y_1,\ldots,Y_n)^t, \beta=(\beta_1,\ldots,\beta_p)^t, X=(x_{ij})_{1\leq i\leq n, 1\leq j\leq p} sowie \epsilon =(\epsilon_1,\ldots,\epsilon_n)^t ein und erhalten die zu Gleichung 3 äquivalente Gleichung Y=X\, \beta +\epsilon. Die Matrix X gibt die Werte der p erklärenden Variablen bei den n Beobachtungen an. X wird oft als Designmatrix bezeichnet, weil die Werte der erklärenden Variablen in vielen Anwendungen zu Beginn der Untersuchungen festgelegt werden und das Design des Experiments beschreiben.

Anmerkung.

1. Wenn man im linearen Modell wie in Gleichung 2 einen konstanten Term integrieren will, so erhält man als erste Spalte der Designmatrix einen Vektor mit n identischen Einträgen, die alle gleich 1 sind.
2. In unserer Darstellung gehen wir davon aus, dass die erklärenden Variablen deterministisch sind und im Idealfall von den Untersuchenden selber frei gewählt werden konnten. Dennoch wenden wir die Methoden auch auf Daten an, bei denen die erklärenden Variablen selber zufällig sind. Man muss dann bei der Interpretation der Ergebnisse etwas vorsichtig sein; sie gelten generell nur, wenn man auf die Werte der erklärenden Variablen bedingt.

2 Allgemeines lineares Modell: Beispiele

Wir wollen im Folgenden einige Beispiele linearer Modelle vorstellen und dabei illustrieren, dass eine überraschende Vielzahl konkreter Modelle als Spezialfall des allgemeinen linearen Modells aufgefasst werden können.

2.1 Einfaches lineares Regressionsmodell

Das einfache lineare Regressionsmodell, das wir bereits im ersten Abschnitt dieses Kapitels vorgestellt hatten, lautet Y_i = \alpha + \beta \, x_i + \epsilon_i. Zunächst verwirrt hier vielleicht der konstante Term \alpha, den man im allgemeinen linearen Modell so direkt nicht findet. Diese Herausforderung lösen wir dadurch, dass wir eine künstliche Variable x_{i1}\equiv 1 einführen und weiter x_{i2}=x_i setzen. Damit erhalten wir Y_i = \alpha \, x_{i 1} + \beta \, x_{i 2}+\epsilon_i, \; 1\leq i\leq n. Dieses Modell hat entsprechend die Designmatrix X = \left( \begin{array}{cc} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_n \end{array} \right) \in \mathbb{R}^{n\times 2}. Die erste Spalte besteht nur aus Einsen, während die zweite Spalte die Werte der erklärenden Variablen bei den n Experimenten enthält.

2.2 k-Stichproben Modell

Aus k verschiedenen Populationen werden Stichproben der Größe n_1, \ldots, n_k genommen. Wir nehmen an, dass die Messwerte Realisierungen normalverteilter Zufallsvariablen mit möglicherweise unterschiedlichen Mittelwerten \mu_1, \ldots, \mu_k und gleicher Varianz \sigma^2 sind, d.h. \begin{align*} X_{1 1},\ldots, X_{1n_{1}} & \sim N(\mu_1,\sigma^2) \\ X_{2 1},\ldots, X_{2n_{k}} & \sim N(\mu_2, \sigma^2) \\ & \; \; \vdots \\ X_{k 1},\ldots, X_{kn_{k}} & \sim N(\mu_k, \sigma^2). \end{align*} Dies ist eine direkte Verallgemeinerung des Zwei-Stichproben Problems, bei dem man Daten x_1,\ldots,x_m und y_1,\ldots,y_n hat, die als Realisierungen unabhängiger normalverteilter Zufallsvariablen X_1,\ldots,X_m und Y_1,\ldots,Y_n mit Erwartungswerten \mu_X und \mu_Y und identischen Varianzen \sigma^2 aufgefasst werden.

Um zu erkennen, dass es sich bei dem k-Stichproben Modell um ein lineares Modell handelt, schreiben wir jede der Zufallsvariablen als Summe von Erwartungswert plus Zufallsterm, d.h. X_{ij}=\mu+\epsilon_{ij}, wobei alle \epsilon_{ij} eine N(0,\sigma^2)-Verteilung haben. So erhalten wir \begin{align*} X_{11} & = \mu_1 + \epsilon_{11}\\ &\; \; \vdots \\ X_{1 n_1} & = \mu_1 + \epsilon_{1n_1}\\ X_{21} & = \mu_2 + \epsilon_{21}\\ &\; \; \vdots \\ X_{2n_2} & = \mu_2 + \epsilon_{2 n_2}\\ &\; \; \vdots \\ X_{k1}& =\mu_k +\epsilon_{k 1} \\ &\; \; \vdots \\ X_{k n_k} & =\mu_k +\epsilon_{kn_k} \end {align*} Wenn wir die Zufallsvariablen zu einem Vektor Y=(X_{11},\ldots,X_{1n_1},\ldots,X_{k1},\ldots, X_{kn_k})^t \in \mathbb{R}^{n} zusammenfassen, wobei n=n_1+\ldots +n_k, so können wir die obigen Gleichungen in Vektor-Matrix Form schreiben als Y =X\, \left(\begin{array}{c}\mu_1\\ \vdots \\ \mu_k \end{array}\right) +\epsilon. Dabei ist die Designmatrix X\in \mathbb{R}^{n\times k} gegeben durch X = \left( \begin{array}{cccc} 1 & 0 & \cdots & 0\\ \vdots & \vdots & \cdots &\vdots \\ 1&0& \cdots & 0 \\ 0& 1 &\cdots & 0 \\ \vdots & \vdots & &0 \\ 0&1& \cdots & 0 \\ & & \vdots & \\ 0 & 0 & \cdots & 1 \\ \vdots & \vdots & &\vdots \\ 0 & 0 & \cdots & 1 \end {array} \right) \in \{0,1\}^{n \times k}

2.3 Ein-Stichproben Modell

Auch das einfachste Modell für wiederholte fehlerbehaftete Messungen einer physikalischen Größe \mu, gegeben durch Y_i = \mu + \epsilon, kann als lineares Modell aufgefasst werden. In diesem Fall ist die Designmatrix gegeben durch X=(1,\ldots,1)^t \in \mathbb{R}^{n\times 1}.

2.4 Polynomiale Regression

Das polynomiale Regressionsmodell stellt eine wesentliche Erweiterung des einfachen linearen Regressionsmodells dar, indem es zur Modellierung der Abhängigkeit zwischen der erklärenden Variablen und dem Ergebnis des Experiments nicht nur lineare Funktionen zulässt, sondern beliebige Polynome. Wir bezeichnen wieder mir x_i den Wert der erklärenden Variablen beim i-ten Experiment und betrachten dann das Modell Y_i = \beta_0 + \beta_1 \, x_i + \ldots + \beta_p \, x^p_i + \epsilon_i Dieses polynomialre Regressionsmodell hat eine Desginmatrix, deren (p+1) Spalten die j-ten Potenzen der Werte der erklärenden Variablen beim i-ten Experiment sind, das heißt X = \left( \begin{array}{ccccc} 1 & x_1 & x_1^2& \cdots & x_1^p \\ 1 & x_2 & x_2^2 & \cdots & x_2^p\\ \vdots & \vdots & \vdots && \vdots \\ 1 & x_n & x_n^2 &\cdots &x_n^p \end{array} \right) \in \mathbb{R}^{n\times (p+1)}. Im ersten Moment mag es überraschen, dass ein lineares Modell auch nichtlineare Funktionen enthalten kann. Bei genauer Betrachtung stellt sich aber heraus, dass es bei einem linearen Modell um lineare Abhängigkeit von den Parametern geht. Entsprechend kann man die Monome x^j, 0\leq j\leq p, durch beliebige Funktionen f_j(x) ersetzen und als Modell beliebige Linearkombinationen dieser Funktionen betrachten. Hingegen ist etwa das Modell Y_i=\exp(\alpha + \beta\, x_i) +\epsilon_i kein lineares Modell mehr.

3 Kleinste Quadrate Schätzer

Wir wollen uns jetzt mit der Schätzung der Parameter im allgemeinen linearen Modell befassen. Gegeben sind also die n Gleichungen \begin{align*} Y_1 & = \beta_1 \, x_{11}+\beta_2 \, x_{12}+\ldots + \beta_p \, x_{1p}+\epsilon_1\\ Y_2 & = \beta_2\, x_{21}+ \beta_{2}\, x_{22} +\ldots + \beta_p\, x_{2p}+ \epsilon_2\\ &\; \; \vdots \\ Y_n & = \beta_1\, x_{n1}+ \beta_2 \,x_{n2} + \ldots + \beta_p \, x_{n p} + \epsilon_n \end {align*} Nach Ablauf des Experiments kennen wir die x_{ij} sowie die Realisierungen der Zufallsvariablen Y_i und wollen daraus die Modellparameter, also \beta_1, \ldots, \beta_p sowie die Varianz von \epsilon_1, schätzen. Die Methode für das Schätzen von \beta_1,\ldots,\beta_p ist die Kleinste-Quadrate-Methode, die darin besteht, folgenden Ausdruck zu minimieren: \sum^{n}_{i=1}\big(y_i - (\beta_1\, x_{i 1}+ \ldots +\beta_p\, x_{i p})\big)^2\\ = \sum^n_{i=1} \big(y_i - \sum^p_{j=i} \beta_j\, x_{ij}\big)^2. Es gibt viele Motivationen der Kleinste-Quadrate-Methode. Eine davon besteht in der Überlegung, dass \widehat{Y} :=\operatorname{E}(Y)=\sum_{j=1}^p x_{j}\, \beta_j die beste Vorhersage für das Ergebnis des Experiments ist, wenn die erklärenden Variablen die Werte x_{1},\ldots, x_{p} haben. Dann ist \sum^n_{i=1} \big(y_i - \sum^p_{j=i} \beta_j\, x_{ij}\big)^2 die Summe der quadratischen Vorhersagefehler bei den vorliegenden Daten, und entsprechend suchen wir mit der Kleinste-Quadrate-Methode unter allen Modellen dasjenige, das diese Summe minimiert.

Es gibt an dieser Stelle (mindestens) zwei Möglichkeiten, das Minimalisierungsproblem zu lösen. Zum einen kann man alle partiellen Ableitungen nach den \beta_j gleich 0 setzen. Eine andere Möglichkeit besteht in einer geometrischen Interpretation des obigen Problems und einer Lösung mit Methoden der linearen Algebra. Beide Verfahren liefern am Ende dieselbe Lösung. Wir beginnen mit dem analytischen Zugang und setzen die partiellen Ableitungen des totalen quadratischen Fehlers nach jedem der Parameter \beta_1,\ldots,\beta_p gleich 0. Auf diesem Wege erhalten wir das lineare Gleichungssystem 0 = \frac{\partial}{\partial \beta_k} \sum_{i=1}^n\big(y_i-\sum_{j=1}^p \beta_j\, x_{ij} \big)^2 = -2 \sum^{n}_{i=1} \big(y_i - \sum^n_{j=1}\,\beta_j \, x_{ij} \big) \, x_{ik}, \; k=1, \ldots, p, das wir wie folgt umformen können: \sum ^n_{i=1}y_i \, x_{ik}= \sum^n_{i=1}\sum^p_{j=1}\beta_j \, x_{ij}\,x_{ik}, \, k=1, \ldots, p. Mit etwas genauem Hinschauen entdeckt man, dass dieses System wie folgt in einer geschlossenen Form als Vektor-Matrix Gleichung geschrieben werden kann: X^t\, y= (X^t X)\, \beta. wobei y=(y_1,\ldots,y_n)^t der Vektor der Werte der abhängigen Variablen ist. Ist jetzt zusätzlich \text{Rang}(X) = p, so ist X^t \, X invertierbar, und damit kann man diese Gleichung nach \beta auflösen und erhält die Formel \hat{\beta} = (X^t X)^{-1}\, X^t\, y.

Wir wollen jetzt noch einen geometrischen Zugang zur Bestimmung des Kleinste-Quadrate-Schätzers vorstellen, der nicht nur viel eleganter ist, sondern darüber hinaus weitere Einsichten verschafft, die uns im Zusammenhang mit Tests von Hypothesen über die Parameter des linearen Modells noch weiterhelfen werden. Der geometrische Zugang beginnt mit der Beobachtung, dass die zu minimierende Summe der Quadrate sich als quadrierter euklidischer Abstand von zwei Vektoren auffassen lässt, nämlich \sum^{n}_{i=1}(y_i - \sum^p_{j=1}\beta_j \, x_{ij})^2= ||y-X\, \beta||^2. Wir suchen jetzt dasjenige \beta, für das der Abstand von X\, \beta zu dem gegebenen y minimal wird. Wir betrachten dazu den Raum \omega=\{ X \beta: \beta \in \mathbb{R}^p\} \subset\mathbb{R}^n. Der Raum \omega ist der Spaltenraum der Designmatrix, d.h. der Unterraum des \mathbb{R}^n, der von den Spalten X_1,\ldots , X_p der Designmatrix X aufgespannt wird. So betrachtet, suchen wir also den Vektor {\widehat \xi} \in \omega, der y am nächsten liegt. Aus der linearen Algebra ist bekannt, dass man diesen Vektor durch orthogonale Projektion von y auf den Unterraum \omega finden kann. Diese Projektion \widehat {\xi} ist eindeutig charakterisiert durch die beiden Eigenschaften \widehat {\xi} \in \omega und y - \widehat {\xi} \perp \omega. Um hieraus den Kleinste-Quadrate-Schätzer \widehat{\beta} zu bestimmen, bemerken wir, dass \widehat{\xi} dargestellt werden kann als {\widehat \xi} = X{\widehat \beta}. Da \omega = \, \text{span} \, (X_1, \ldots, X_p), ist die zweite Eigenschaft der Projektion, also y-\widehat{\xi}\perp \omega, äquivalent dazu, dass y-X{\hat \beta} \perp X_j für alle j\in \{1,\ldots, p\}. Letzteres ist äquivalent dazu, dass X_j^t(y-X\, \widehat{\beta})=0, j=1,\ldots,p, was sich in einer einzigen Gleichung schreiben lässt als X^t y = (X^t X) \widehat {\beta}. \tag{4} Wir sind also am Ende bei demselben linearen Gleichungssystem für \widehat{\xi} angekommen wie beim analytischen Zugang. Auch hier lässt sich das Gleichungssystem Gleichung 4 genau dann auflösen, wenn \text{Rang}(X) = p. Wir werden ab jetzt stets diese Annahme machen. Wir fassen unsere Ergebnisse jetzt in einem Theorem zusammen

Theorem 1

Der Kleinste-Quadrate-Schätzer für \beta im allgemeinen linearen Modell y = X \beta + \varepsilon mit \text{Rang}(X) = p ist gegeben durch {\widehat \beta} = (X^t X)^{-1} X^t y.

Neben den Regressionskoeffizienten \beta_1,\ldots,\beta_p ist auch die Varianz \sigma^2 des Zufallsterms \epsilon unbekannt. Diese Varianz spiegelt sich im Vorhersagefehler wider, sodass der mittlere quadratische Vorhersagefehler \frac{1}{n} \sum_{i=1}^n (y_i-\sum_{j=1}^p \widehat{\beta}_j x_{ij})^2 ein naheliegender Schätzer ist. Dieser Schätzer ist jedoch nicht erwartungstreu, sodass wir stattdessen eine Variante nehmen, von der wir im weiteren Verlauf dieses Kapitels zeigen werden, dass sie erwartungstreu ist: s_{y|x}^2:=\frac{1}{n-p} \sum_{i=1}^n (y_i-\sum_{j=1}^p \widehat{\beta}_j x_{ij})^2.

Beispiel 1 In R findet man den Datensatz USairpollution, der Daten zur Luftverschmutzung in n=41 Großstädten in den USA enthält. Die Luftverschmutzung wird charakterisiert durch den Gehalt an Schwefeldioxid (y=SO2). Als erklärende Variablen werden die durchschnittliche Jahrestemperatur (x_1=temp), die Anzahl der Industrieunternehmen mit mindestens 20 Beschäftigten (x_2=manu), die Einwohnerzahl (x_3=popul), die durchschnittliche jährliche Windgeschwindigkeit (x_4=wind), die jährliche Niederschlagsmenge (x_5=precip) sowie die durchschnittliche Anzahl an Regentagen (x_6=predays) verwendet.

data("USairpollution", package = "HSAUR3")
USairpollution

Wir modellieren die Daten durch ein multiples lineares Regressionsmodell, indem wir den Schwefeldioxidgehalt als Linearkombination der erklärenden Variablen und eines Zufallsterms auffassen, d.h. wir fassen die y-Werte als Realisierung einer Zufallsvariablen Y auf, die gegeben ist durch Y= \beta_0+ \beta_1\, x_1 + \beta_2\, x_2 + \beta_3\, x_3 + \beta_4\, x_4+ \beta_5\, x_5 + \beta_6\, x_6 + \epsilon. Zur Analyse dieses Modells inRverwenden wir den Befehl

attach(USairpollution)
summary(lm(SO2 ~ temp + manu + popul + wind + precip + predays))

Hier haben wir durch SO2 ~ temp + manu + popul + wind + precip + predays angegeben, dass der Schwefeldioxidgehalt durch eine Linearkombination der sechs erklärenden Variablen modelliert werden soll. Wenn wir nur eine Auswahl der erkärenden Variablen betrachten wollen, müssen wir den R-Befehl entsprechend anpassen. So würden wir zum Beispiel mit dem Befehl lm(SO2~temp+precip) ein lineares Modell betrachten, dass ausschließlich die durchschnittliche Jahrestemperatur und die jährlichen Niederschlagsmenge als erklärende Variablen hat.

Der obige R-Befehl liefert uns die folgende Ausgabe:

Call:
lm(formula = SO2 ~ temp + manu + popul + wind + precip + predays)

Residuals:
    Min      1Q  Median      3Q     Max 
-23.004  -8.542  -0.991   5.758  48.758 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 111.72848   47.31810   2.361 0.024087 *  
temp         -1.26794    0.62118  -2.041 0.049056 *  
manu          0.06492    0.01575   4.122 0.000228 ***
popul        -0.03928    0.01513  -2.595 0.013846 *  
wind         -3.18137    1.81502  -1.753 0.088650 .  
precip        0.51236    0.36276   1.412 0.166918    
predays      -0.05205    0.16201  -0.321 0.749972    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.64 on 34 degrees of freedom
Multiple R-squared:  0.6695,    Adjusted R-squared:  0.6112 
F-statistic: 11.48 on 6 and 34 DF,  p-value: 5.419e-07

Die erste Spalte enthält die Schätzwerte der Regressionsparameter \beta_0,\ldots,\beta_6, berechnet mithilfe der Kleinste-Quadrate-Methode. Den Schätzwert der Standardabweichung \sigma findet man als Residual standard error. Wir erhalten somit das lineare Modell Y=111.73 - 1.27 \, x_1 +0.06\ x_2 - 0.04\, x_3 -3.18 \, x_4 +0.51\, x_5 -0.05\, x_6 +\epsilon. Für die Standardabweichung von \epsilon erhalten wir den Schätzwert s_{y|x} = 14.64.

Einige der Schätzwerte sind auf den ersten Blick überraschend. So sollte man erwarten, dass die Luftverschmutzung in Städten mit vielen Einwohnern größer ist als in Städten mit wenigen Einwohnern. Entgegen dieser Intuition ist der entsprechende Regressionskoeffizient \beta_3=-0.04 negativ. Die Erklärung liegt darin, dass der Regressionskoeffizient \beta_3 angibt, um wie viel die erwartete Luftverschmutzung zunimmt, wenn die Einwohnerzahl um 1.000 wächst, während alle anderen erklärenden Variablen festgehalten werden. Damit vergleicht man Städte mit derselben Anzahl großer Industrieunternehmen, aber unterschiedlicher Einwohnerzahl, und dann ist das Ergebnis nicht mehr so überraschend, weil die größere Stadt weniger industrialisiert ist.

Aufgabe 1 Welchen Schwefeldioxidgehalt hätten Sie aufgrund des linearen Regressionsmodells mit den oben berechneten Schätzwerten für die Parameter in Albany erwartet?

4 Kleinste Quadrate Schätzer: Beispiele

Beispiel 2

Als erstes Beispiel betrachten wir das einfache lineare Regressionsmodell. In diesem Fall hat die Designmatrix X\in \mathbb{R}^{n\times 2} in der ersten Spalten nur Einsen und in der zweiten Spalte den Vektor der erklärenden Variablen x_1,\ldots,x_n. Wir berechnen jetzt den Kleinste-Quadrate-Schätzer \widehat{\beta}=(X^tX)^{-1}X^ty und beginnen mit der Matrix X^tX\in \mathbb{R}^{2\times 2}. Es gilt X^tX =\left( \begin{array}{ccc} 1 & \cdots & 1 \\ x_1 & \cdots & x_n \end{array}\right) \left(\begin{array}{cc} 1 &x_1 \\ \vdots & \vdots \\ 1 &x_n \end{array} \right) = \left( \begin{array}{cc}n &\sum_{i=1}^n x_i \\[2mm] \sum_{i=1}^n x_i &\sum_{i=1}^n x^2_i \end{array} \right). Im nächsten Schritt berechnen wir die Inverse dieser Matrix. Mithilfe der Regeln für das Invertieren von 2\times 2-Matrizen erhalten wir \begin{align*} (X^t X)^{-1} & = \frac {1}{n \, \sum^n_{i=1}x_i^2 -(\sum^n_{i=1}x_i)^2} \left(\begin{array}{cc} \sum^n_{i=1} x_i^2 & -\sum^n_{i=1} x_i\\[2mm] - \sum^n_{i=1} x_i & n \end{array}\right) \\[2mm] & = \frac{1}{\sum_{i=1}^n x_i^2 -\frac{1}{n}(\sum_{i=1}^n x_i)^2} \left(\begin{array}{cc} \frac{1}{n} \sum^n_{i=1} x_i^2 & -\bar{x}\\[2mm] - \bar{x} & 1\\ \end{array}\right) \\[2mm] &=\frac{1}{\sum_{i=1}^n (x_i-\bar{x})^2} \left(\begin{array}{cc} \frac{1}{n} \sum^n_{i=1} x_i^2 & -\bar{x}\\[2mm] - \bar{x} & 1\\ \end{array}\right). \end{align*}

Als letzten Zwischenschritt berechnen wir schließlich das Produkt X^t \, y:

X^t \, y= \left( \begin{array}{ccc} 1 & \ldots & 1 \\[1mm] x_1 & \ldots & x_n \end{array}\right) \left( \begin{array}{c} y_1 \\ \vdots \\ y_n \end{array} \right) =\left( \begin{array}{c} \sum^n_{i=1} y_1 \\[2mm] \sum^n_{i=1} x_i \, y_i \end{array} \right).

Mit diesen Zwischenergebnissen können wir jetzt mittels der Formel \widehat{\beta}=(X^t\,X)^{-1}X^t\, y den Kleinste-Quadrate-Schätzer für den Parametervektor (\alpha,\beta) berechnen. Wir erhalten \begin{align*} \left( \begin{array}{c} \hat{\alpha} \\[3mm] \hat{\beta} \end{array} \right) & = \frac{1}{\sum^n_{i=1}(x_i - \bar{x})^2} \left( \begin{array}{cc} \frac{1}{n} \sum^n_{i=1} x_i^2 & - \bar{x}\\[3mm] -\bar{x} & 1\\ \end{array}\right) \left( \begin{array}{c} \sum^n_{i=1} y_i\\[3mm] \sum^n_{i=1} x_i\, y_i \end{array}\right)\\[3mm] & = \frac{1}{\sum^n_{i=1}(x_i - \bar{x})^2} \left( \begin{array}{c} \bar{y} \, \sum^n_{i=1} x_i^2 - \bar{x}\sum_{i=1}^n x_i\, y_i\\[3mm] -n\, \bar{x}\, \bar{y} + \sum^n_{i=1} x_i \, y_i\end{array}\right). \end{align*} Wir betrachten jetzt die beiden Komponenten des Parametervektors einzeln und beginnen mit \widehat{\beta}, dem Schätzer für die Steigung der Regressionsgeraden. Unter Verwendung der Identität \sum_{i=1}^n (x_i-\bar{x})(y_i-\bar{y})=\sum_{i=1}^n x_i\, y_i - n\, \bar{x}\, \bar{y} erhalten wir \begin{align*} \hat{\beta} &= \frac{1}{\sum_{i=1}^n (x_i-\bar{x})^2} \Big(\sum_{i=1}^n x_i\, y_i - n\, \bar{x}\, \bar{y} \Big) \\ & = \frac{1}{\sum^n_{i=1}(x_i-\bar{x})^2} \sum_{i=1}^n (x_i-\bar{x})(y_i -\bar{y}), \end{align*} also genau die Formel, die wir bereits zu Beginn dieses Kapitels bei der Betrachtung des Spezialfalls der einfachen linearen Regression gefunden hatten. Für den geschätzten y-Achsenabschnitt der Regressionsgeraden kennen wir ebenfalls bereits eine Formel, die wir aus der obigen Darstellung mit etwas ergebnisoriertem Rechnen erneut herleiten können: \begin{align*} \hat{\alpha} & = \frac{1}{\sum^n_{i=1}(x_i-\bar{x})^2} \Big( \bar{y} \, \sum^n_{i=1}x_i^2 - \bar{x}\, \sum^n_{i=1}x_i y_i\Big) \\ & = \frac{1}{\sum_{i=1}^n (x_i -\bar {x})^2} \Big( \bar{y}\big (\sum_{i=1}^n x_i^2 -n\, (\bar{x})^2 \big) + n\, (\bar{x})^2 \, \bar{y} - \bar{x}\, \sum_{i=1}^n x_i y_i \Big) \\ & = \bar{y} - \frac{1}{\sum_{i=1}^n (x_i-\bar{x})^2 }\, \bar{x}\, \Big(\sum^n_{i=1}x_i y_i - n\, \bar{x}\, \bar{y}\Big) \\ & = \bar{y} - \bar{x}\, \hat{\beta}, \end{align*} also auch hier genau die Formel, die wir bereits bei der Betrachtung des Spezialfalls der einfachen linearen Regression gefunden hatten.

Der Schätzer für die Varianz ist in diesem Fall s_{y|x}^2 = \frac{1}{n-2} \| y-X\, \beta \|^2 = \frac{1}{n-2}\sum^n_{i=1}(y_i - \hat{\alpha}- \hat{\beta}x_i)^2, in Übereinstimmung mit dem Varianzschätzer im einfachen linearen Regressionsmodell.

Beispiel 3
Als zweites Beispiel betrachten wir das k-Stichprobenproblem. Wir erinnern nochmals an dieses Modell: \begin{align*} X_{11}, \ldots, X_{1n_1} & \sim N(\mu_1, \sigma^2)\\ X_{21},\ldots, X_{2n_2} & \sim N(\mu_2, \sigma^2)\\ &\; \; \vdots\\ X_{k1}, \ldots, X_{kn_k} & \sim N(\mu_k, \sigma^2) \end{align*} Wir hatten bereits erkannt, dass man dieses Modell als ein lineares Modell auffassen kann mit dem Datenvektor Y=(X_{11}, \ldots, X_{kn_k}) und der Designmatrix X = \left( \begin{array}{cccc} 1 & 0 & \cdots & 0\\ \vdots & \vdots & & \vdots\\ 1 & 0 & \cdots & 0\\ 0 & 1 & & 0 \\ \vdots & \vdots & &\vdots \\ 0 & 1 & & 0\\ && \ddots & \\ 0 & 0 && 1\\ \vdots & \vdots & & \vdots\\ 0 & 0 & & 1 \end{array}\right). Zur Berechnung des Kleinste-Quadrate-Schätzers für den Parametervektor \beta=(\mu_1,\ldots,\mu_k)^t mittels der Formel \widehat{\beta}=(X^t\,X)^{-1}\, X^t \,Y berechnen wir zunächst die Matrix X^t\,X\in \mathbb{R}^{k\times k}. Es gilt X^t \, X =\left(\begin{array}{cccc} n_1 &0& \cdots &0\\ 0 & n_2 & \cdots & 0\\ \vdots & \vdots& \ddots & \vdots \\ 0 & 0& \cdots & n_k \end{array}\right). Diese Diagonalmatrix kann leicht invertiert werden; wir erhalten (X^t \, X)^{-1} =\left(\begin{array}{cccc} \frac{1}{n_1} &0& \cdots &0\\ 0 & \frac{1}{n_2} & \cdots & 0\\ \vdots & \vdots& \ddots & \vdots \\ 0 & 0& \cdots & \frac{1}{n_k} \end{array}\right). Weiter gilt X^tY=\left(\begin{array}{c} \sum^{n_1}_{j=1} X_{1j}\\ \vdots \\ \sum^{n_k}_{j=1} X_{kj} \end{array}\right), und somit erhalten wir schließlich den Kleinste-Quadrate-Schätzer \hat{\beta} = (X^t X)^{-1} X^t Y=\left(\begin{array}{c} \frac{1}{n_1}\sum^{n_1}_{j=1} X_{ij}\\ \vdots \\ \frac{1}{n_k} \sum^{n_k}_{j=1} X_{kj} \end{array}\right). Wenn wir die Einträge auf der rechten und linken Seite vergleichen, finden wir folgende Formel für den Kleinste-Quadrate-Schätzer von \mu_i: \hat{\mu}_i = X_{i\cdot}:= \frac{1}{n_i} \sum_{j=1}^{n_i} X_{ij}. Dieses Ergebnis ist nicht überraschend, denn das arithmetische Mittel der i-ten Stichprobe ist ein naheliegender Schätzer für den Erwartungswert \mu_i. Wir haben an dieser Stelle noch eine sehr praktische Notation eingeführt: Wenn wir in einer mehrfach indizierten Folge einen oder mehrere der Indizes durch einen Punkt ersetzen, so meinen wir damit das arithmetische Mittel genommen über alle Werte dieser Indizes.

Zum Schluss berechnen wir noch den Schätzer für die Varianz und erinnern zunächst an die allgemeine Formel s^2=\frac{1}{n-k}\|Y-X\hat{\beta}\|^2. In diesem konkreten Fall erhalten wir \begin{align*} s^2 : & = \frac{1}{n-k}||Y-X {\hat\beta}||^2\\ & = \frac{1}{n-k}\sum^k_{i=1} \sum_{j=1}^{n_i}\Big( X_j - \frac{1}{n_i}\sum_{j=1}^{n_i} X_{ij}\Big)^2\\ & = \frac{1}{n-k}\sum_{i=1}^k \sum _{j=1}^{n_i} (X_{ij}-X_{i\cdot})^2. \end{align*}

5 Eigenschaften des Kleinste Quadrate Schätzers

In diesem Abschnitt werden wir einige Eigenschaften des Kleinste Quadrate Schätzers für die Regressionsparameter \beta_1,\ldots, \beta_p sowie des Schätzers s_{y|x}^2 für die Varianz des Fehlerterms \epsilon vorstellen. An dieser Stelle müssen wir Annahme über die Verteilung von \epsilon treffen. Folgende beiden Annahmen sind üblich:

\epsilon_1, \ldots, \epsilon_n \text{ sind unkorreliert mit } \operatorname{E}(\epsilon_i)= 0, \operatorname{Var}(\epsilon_1)=\sigma^2. \tag{5} \epsilon_1, \ldots, \epsilon_n \epsilon_i \text{ sind unabhängig und } N(0,\sigma^2)\text{-verteilt}. \tag{6}

Von diesen beiden Annahmen ist Gleichung 6 die restriktivere, weil aus Unabhängigkeit die Unkorreliertheit folgt. Unter der Annahme Gleichung 6 kann man entsprechend stärkere Aussagen über die Schätzer beweisen, die aber natürlich nur gelten, wenn diese Annahme wirklich erfüllt ist. Wir beginnen mit Ergebnissen, die unter der allgemeineren Annahme Gleichung 5 gelten.

Theorem 2
Durch Y = X\, \beta +\epsilon sei ein lineares Modell mit unkorrelierten \epsilon _i und \operatorname{E} (\epsilon_i) = 0, \, \operatorname{Var}(\epsilon_i)= \sigma^2 gegeben, wobei X \in \mathbb{R}^{n \times p} mit \text{Rang}( X) = p. Dann gilt
(i) \operatorname{E}(\hat{\beta})= \beta, d.h. \hat{\beta} ist ein erwartungstreuer Schätzer für \beta.
(ii) Die Kovarianzmatrix von \hat{\beta} ist gegeben durch \Sigma_{\hat \beta}= \sigma^2 (X^t X)^{-1}. (iii) \operatorname{E}(s^2) = \sigma^2, d.h. s^2 ist ein erwartungstreuer Schätzer für \sigma^2.

Beweis.

Die beiden ersten Aussagen folgen fast unmittelbar aus der expliziten Darstellung des Kleinste-Quadrate-Schätzers \widehat{\beta}= (X^t X)^{-1}\, X^t \, Y zusammen mit den Regeln für das Verhalten von Erwartungswert und Kovarianzmatrix unter linearen Transformationen. Für den Erwartungswert erhalten wir \operatorname{E}(\widehat{\beta}) = (X^t X)^{-1} X^t \operatorname{E}(Y) = (X^t X)^{-1} X^t \operatorname{E} (X\beta +\epsilon) = (X^t X)^{-1} X^t X \beta = \beta und für die Kovarianzmatrix \Sigma_{\widehat{\beta}} = (X^t X)^{-1}X^t \, \Sigma_Y ((X^t X)^{-1}X^t)^t = \sigma^2 (X^t X)^{-1}. (iii) Der Beweis von Teil (iii) ist etwas aufwendiger und erfordert die Einführung geeigneter Koordinaten im \mathbb{R}^n. Wir betrachten dazu das lineare Modell Y=X\, \beta +\epsilon in einer koordinatenfreien Darstellung, indem wir \xi=X\beta einführen. Wir definieren den Unterraum \omega: = {\rm span} \{ X_1, \ldots, X_p\} und bemerken, dass die Elemente von \omega genau die Vektoren \xi sind, die eine Darstellung als \xi=X\, \beta, \beta\in \mathbb{R}^p haben. Entsprechend ist unser lineares Modell äquivalent dazu, dass Y=\xi + \epsilon ,\quad \xi \in \omega. Wie wir bereits bei der Berechnung des Kleinste-Quadrate-Schätzers gesehen hatten, ist der Kleinste-Quadrate-Schätzer für \xi gegeben durch die Projektion von Y auf den Unterraum \omega und der Schätzer für \sigma^2 ist der quadratische Abstand des Datenvektors von \omega geteilt durch n-p, d.h. \begin{align*} \hat{\xi} &=\text{proj}_{\omega}Y, \\ s^2 &= \frac{1}{n-p}||Y - {\hat \xi}||^2. \end{align*} Wir führen nun eine neue orthonormale Basis im \mathbb{R}^n ein, in der der Kleinste-Quadrate-Schätzer und somit auch s^2 ganz einfach berechnet werden können. Dazu sei u_1, \ldots, u_p eine Orthonormalbasis für \omega und u_1, \ldots u_p, u_{p+1}, \ldots,u_n insgesamt eine Orthonormalbasis für \mathbb{R}^n. Diese lässt sich stets finden, indem man die Spalten X_1, \ldots X_p der Designmatrix mit dem Gram-Schmidt Verfahren orthogonalisiert und anschließend die so gewonnene Orthonormalbasis u_1, \ldots, u_p für \omega zu einer Basis für den gesamten \mathbb{R}^n ergänzt. Aus der linearen Algebra ist bekannt, dass man den neuen Koordinatenvektor z\in \mathbb{R}^n eines Vektors durch Multiplikation des alten Koordinatenvektors y \in \mathbb{R}^n mit einer orthogonalen Transformationsmatrix U erhält, d.h. z = U y, wobei U die Matrix mit den Zeilen u_1^t, \ldots, u_n^t ist. In den neuen Koordinaten lautet das lineare Modell Y = \xi + \epsilon dann Z = \zeta + \epsilon^\prime, wobei \epsilon^\prime= U\, \epsilon und \zeta=U\, \xi die neuen Koordinaten des Zufallsterms bzw. des Parametervektors sind. Nach den Transformationsregeln für Erwartungswerte und Kovarianzmatrizen unter linearen Abbildungen gilt \operatorname{E}(\epsilon^\prime)=0 und \Sigma_{\epsilon^\prime}= U \Sigma_{\epsilon} U^t = \sigma^2 I_n, d.h. die Zufallsvariablen \epsilon_i^\prime sind unkorreliert mit Erwartungswert 0 und Varianz \sigma^2. Für \zeta, den neuen Koordinatenvektor von \xi, bedeutet die Bedingung \xi \in \omega, dass \xi = \zeta_1 u_1 + \ldots + \zeta_p u_p, und somit \zeta_{p+1}= \ldots = \zeta_n = 0. Insgesamt erhalten wir auf diesem Wege das lineare Modell in kanonischer Form Z = \zeta + \epsilon^\prime, \; \zeta_{p+1} = \ldots = \zeta_n = 0. In den neuen Koordinaten bestimmt man ganz einfach die Projektion auf \omega. Die Projektion ist dadurch charakterisiert, dass || Y - \xi ||^2 = \sum^n_{i=1}(Z_i - \zeta_i)^2 unter der Nebenbedingung \xi \in \omega minimiert wird, d.h. unter \zeta_{p+1}= \ldots = \zeta_n =0. In dieser Form ist das Minimierungsproblem ganz einfach zu lösen; das Ergebnis lautet \hat{\zeta}_{1}= Z_1, \ldots, \hat{\zeta}_p =Z_p, \hat{\zeta}_{p+1} = \ldots = \hat{\zeta}= 0. Somit folgt || Y - \hat{\xi}||^2 = \sum^n_{i=1}(Z_i - \hat{\zeta}_i)^2 = \sum^n_{i=p+1}Z_i^2. Nun ist \operatorname{E}(Z_i) = \zeta_i, \, \operatorname{Var} (Z_i) = \sigma^2 und somit gilt \operatorname{E}(Z_i^2) = \sigma^2 für i\geq\, p+1. Hieraus folgt schließlich \operatorname{E} (\| Y - \hat {\xi}\|) = (n-p)\sigma^2, und damit ist bewiesen, dass s^2=\frac{1}{n-p}\|Y-X\, \hat{\beta}\|^2 ein erwartungstreuer Schätzer für \sigma^2 ist.

Unter der stärkeren Annahme der Normalverteilung kann man wesentlich mehr Eigenschaften der Schätzer \hat{\beta} und s^2 beweisen und sogar die gemeinsame Verteilung berechnen. Alle Ergebnisse sind im folgenden Satz festgehalten.

Theorem 3
Wir betrachten das lineare Modell Y = X \beta + \epsilon mit Designmatrix X\in \mathbb{R}^{n \times p} und normalverteiltem Zufallsterm \epsilon \sim N(0,\sigma^2 I_n) unter der Annahme, dass \text{Rang}(X) = p. Dann gilt für die Schätzer \hat{\beta} und s^2:
(i) \hat{\beta} \sim N(\beta, \sigma^2(X^t X)^{-1})
(ii) (n-p) s^2/ \sigma^2 \sim \chi^2_{n-p}
(iii) \hat{\beta} und s^2 sind unabhängige Zufallsvariablen.

Beweis.
Die Aussage (i) folgt mit den Transformationsregeln für multivariate Normalverteilungen unter linearen Transformationen unmittelbar aus der Darstellung \hat{\beta}=(X^t X)^{-1}Y, wobei Y \sim N(X\beta, \sigma^2 I_n). Zum Beweis von (ii) und (iii) verwenden wir noch einmal die Darstellung in kanonischen Koordinaten, die wir im Beweis von Theorem 2 eingeführt haben. Wir bezeichnen wieder mit Z, \zeta und \epsilon^\prime die neuen Koordinatenvektoren von Y, \xi=X\, \beta bzw. \epsilon bezüglich der Orthonormalbasis u_1,\ldots,u_p für \omega, den Spaltenraum der Designmatrix X. In den neuen Koordinaten lautet das lineare Modell Z = \zeta + \epsilon^\prime. Wir erinnern daran, dass \xi \in \omega äquivalent dazu ist, dass \xi als Linearkombination von u_1,\ldots,u_p geschrieben werden kann, was wiederum äquivalent dazu ist, dass \zeta_{p+1}= \ldots =\zeta_n = 0. Nach den Transformationsregeln für multivariate Normalverteilungen sind sowohl \epsilon^\prime als auch Z normalverteilt und es gelten \epsilon^\prime \sim N(0, \sigma^2 I_n) sowie Z \sim N(\zeta, \sigma^2 I_n). Die neuen Koordinaten der Projektion von Y auf \omega sind \hat{\zeta}=(Z_1, \ldots, Z_p, 0, \ldots,0). Hieraus folgt direkt \begin{align*} \hat{\xi}&= Z_1 u_1 + \ldots + Z_p u_p,\\ ||Y-\hat{\xi}||^2 &= \sum^n_{i=p+1} Z_i^2. \end{align*} Aus dieser Darstellung von \hat{\xi} und || Y-\hat{\xi}||^2 kann man jetzt direkt Teil (ii) und (iii) des Satzes herleiten. Es gilt s^2 = \frac{1}{n-p} \|Y-\hat{\xi}\|^2= \frac{1}{n-p} \sum^n_{i=p+1}Z_i^2 und somit (n-p)\, s^2/ \sigma^2= \sum^n_{i=p+1}(Z_i/\sigma)^2 \sim \chi^2_{n-p}, da Z_i\sim N(0,\sigma^2) für i=p+1,\ldots, n. Weiter sind die Zufallsvariablen Z_1, \ldots, Z_p unabhängig von Z_{p+1},\ldots, Z_n und somit ist s^2 unabhängig von \hat{\xi} und somit auch von \hat{\beta}.

Anmerkung.
Im Fall normalverteilter Fehler ist der Kleinste-Quadrate-Schätzer zugleich der Maximum-Likelihood-Schätzer für den Parameter \beta, was uns eine tiefere Rechtfertigung für die Verwendung des Kleinste-Quadrate-Schätzers liefert, die bereits Gauß bekannt war. Um dies zu erkennen, schreiben wir zunächst die gemeinsame Dichte der Beobachtungen Y_1, \ldots, Y_n auf. Nach den Transformationsregeln für eindimensionale Normalverteilungen hat Y_i = \sum^{p}_{j=1}x_{ij} \beta_j + \epsilon_i eine N(\sum^p_{j=1}x_{ij} \beta_{j}, \sigma^2)-Verteilung und somit haben Y_1,\ldots,Y_n die gemeinsame Dichtefunktion

\begin{align*} f_{\beta, \sigma^2}(y_1,\ldots, y_n)& = \frac{1}{(2 \pi\sigma^2)^{n/2}} \exp\Big(-\frac{1}{2\sigma^2}\sum^n_{i=1}\big(y_i - \sum^p_{j=1} x_{ij}\beta_j \big)^2\Big)\\ & = \frac{1}{(2 \pi\sigma^2)^{n/2}} \exp\Big( - \frac{1}{2 \sigma^2}||y - X \beta||^2\Big). \end{align*}

Entsprechend ist die Likelihood-Funktion gegeben durch L(\beta,\sigma)= \frac{1}{(2 \pi\sigma^2)^{n/2}}\, \exp\Big( - \frac{1}{2 \sigma^2}\, ||y - X \beta||^2\Big). Auch hier ist es einfacherer, die Log-Likelihood-Funktion zu maximieren, also l(\beta,\sigma) = -\frac{n}{2}\log(2 \pi\sigma^2)-\frac{1}{2\sigma^2}||y -X \beta||^2. Zur Bestimmung der Maximum-Likelihood-Schätzer für \beta und \sigma^2 gehen wir so vor, dass wir erst \sigma^2 festhalten und \beta\mapsto l(\beta,\sigma) maximieren. Dies führt direkt zur Kleinste-Quadrate-Methode ||y -X \beta||^2 \rightarrow \min, deren Lösung \hat{\beta} wir bereits gefunden haben. Die Suche nach dem Maximum-Likelihood-Schätzer für den gesamten Parameter wird jetzt erheblich vereinfacht durch die Tatsache, dass die Stelle, an der l(\beta,\sigma) bei festgehaltenem \sigma ihr Minimum annimmt, nicht von \sigma abhängt.
Wir können also jetzt den bereits gefundenen Schätzer für \beta in die Log-Likelihood-Funktion einsetzen und haben dann noch \sigma \mapsto l(\hat{\beta},\sigma)= -\frac{n}{2} \log 2 \pi - n \log \sigma - \frac{1}{2 \sigma^2}|| \, y-X \hat {\beta}||^2 zu minimieren. Dies ist ein eindimensionales Minimierungsproblem, das wir durch Bestimmung der Nullstellen der Ableitung der obigen Funktion nach \sigma lösen können. Es gilt \frac{\partial}{\partial \sigma} l(\hat{\beta},\sigma)= - \frac{n}{\sigma} + \frac{1}{\sigma^3}||y-X\hat{\beta}||^2, und somit ist der Maximum-Likelihood-Schätzer gegeben durch \hat{\sigma}_{ML}^2= \frac{1}{n}\, ||y-X \hat{\beta}||^2. Wir weisen darauf hin, dass der Maximum-Likelihood-Schätzer von dem eingangs eingeführten Schätzer s^2 insofern abweicht, als wir durch n teilen und nicht durch n-p. Entsprechend ist der Maximum-Likelihood-Schätzer auch nicht erwartungstreu.

6 Konfidenzintervalle für die Regressionsparameter

Wir können die bislang erzielten Ergebnisse jetzt verwenden, um Konfidenzintervalle für die Parameter des linearen Modells anzugeben. Vorab formulieren und beweisen wir eine Folgerung aus dem obigen Satz.

Theorem 4 (Korollar)
Unter denselben Annahmen wie in Theorem 3 gilt für j = 1, \ldots, p T := \frac{\hat{\beta}_j - \beta_j}{\sqrt{s^2((X^t X)^{-1})_{jj}}}\sim t_{n-p}.

Beweis.
Als Folgerung aus Theorem 3 erhalten wir, dass \hat{\beta}_j eine Normalverteilung mit Erwartungswert \beta_j und Varianz \sigma^2 ((X^tX)^{-1})_{jj} hat, wobei ((X^tX)^{-1})_{jj} der (j,j)-Eintrag der Matrix (X^tX)^{-1} ist. Also hat \frac{\hat{\beta}_j-\beta_j}{ \sqrt{\sigma^2 ((X^tX)^{-1})_{jj} }} eine Standardnormalverteilung. Weiter hat (n-p)\, s^2/\sigma^2 eine \chi_{n-p}^2-Verteilung und ist unabhängig von \hat{\beta}_j. Wir schreiben T jetzt wie folgt um: T=\frac{(\hat{\beta}_j -\beta_j)/ \sqrt{\sigma^2((X^t X)^{-1})_{jj}}}{\sqrt{\big((n-p)\frac{s^2}{\sigma^2}\big)/(n-p)}} und erkennen, dass die rechte Seite t_{n-p}-verteilt ist.

Anmerkung.
Mithilfe des obigen Korollars können wir ein (1-\alpha)-Konfidenzintervall für \beta_j bestimmen. Es gilt mit Wahrscheinlichkeit 1-\alpha, dass -t_{n-p, 1- \alpha/2}\leq \frac{\hat{\beta}_j - \beta_j}{\sqrt{s^2((X^t X)^{-1})_{jj}}} \leq t_{n-p, 1- \alpha/2}, wobei t_{n-p,1-\alpha/2} das (1-\alpha/2)-Quantil der t_{n-p}-Verteilung bezeichnet. Die obige Ungleichungskette ist äquivalent zu \hat{\beta}_j - t_{n-p, 1-\alpha/2}\, \sqrt{s^2\, ((X^t X)^{-1})_{jj}} \leq \beta_j \leq \hat{\beta}_j + t_{n-p, 1-\alpha/2}\, \sqrt{s^2\, ((X^t X)^{-1})_{jj}}, und somit bildet \Big[\hat{\beta}_j - t_{n-p, 1-\alpha/2}\, \sqrt{s^2\, ((X^t X)^{-1})_{jj}}\, ,\, \hat{\beta}_j + t_{n-p, 1-\alpha/2}\, \sqrt{s^2\, ((X^t X)^{-1})_{jj}}\Big] ein (1-\alpha)-Konfidenzintervall für den Parameter \beta_j.

Beispiel 4
Wir wollen jetzt die bisherigen Erkenntnisse anwenden auf das einfache lineare Regressionmodell Y_i = \alpha + \beta x_i + \epsilon_i,\; 1\leq i\leq n. In diesem Fall ist (X^t X)^{-1} = \frac{1}{\sum_{i=1}^n (x_i -\bar{x})^2}\left( \begin{array}{cc} \frac {1}{n} \sum^n_{i=1} x_i^2 & -\bar{x}\\[2mm] -\bar{x} & 1 \end{array}\right), und somit folgt \operatorname{E}(\hat{\alpha})=\alpha und \operatorname{E}(\hat{\beta})=\beta sowie \begin{align*} \operatorname{Var}(\hat{\alpha})& = \sigma^2 \frac{\sum^n_{i=1}x_i^2} {n \sum^n_{i=1} (x_i - \bar{x})^2}, \\ \operatorname{Var}( \hat{\beta})& = \sigma^2 \frac{1}{\sum^n_{i=1}(x_i - \bar{x})^2}, \end{align*} in Übereinstimmung mit den Ergebnissen, die wir bereits mit elementaren Methoden bei der Untersuchung des einfachen linearen Regressionsmodells gefunden hatten. Interessant ist noch, dass die Varianz von \hat{\beta} umgekehrt proportional zu \sum_{i=1}^n(x_i - \bar{x})^2 ist, also zu der totalen Streuung der x-Werte. Je größer diese Streuung ist, umso präziser kann man die Steigung der Regressionsgeraden schätzen. Wenn möglich, sollte man also Experimente bei weit auseinanderliegenden x-Werten durchführen. Allerdings sollte man dabei in der Praxis unbedingt bedenken, dass lineare Modelle immer nur in einem begrenzten Intervall von x-Werten anwendbar sind, und dass man diesen Bereich nicht verlässt.

Schließlich können wir noch Konfidenzintervalle für die Parameter \alpha und \beta der Regressionsgerade angeben. Wir erhalten die (1-\alpha)-Konfidenzintervalle \begin{align*} \Big[\hat{\alpha}- s \, \sqrt{\frac{\sum^n_{i=1}x_i^2}{n \sum^n_{i=1}(x_i - \bar{x})^2}} \cdot t_{n-2, 1-\alpha/2} &\; , \, \hat{\alpha}+ s \, \sqrt{\frac{\sum^n_{i=1}x_i^2}{n \sum^n_{i=1}(x_i - \bar{x})^2}} \cdot t_{n-2, 1-\alpha/2} \Big], \\[2mm] \Big[ \hat{\beta} - s \, \sqrt \frac{1}{\sum^n_{i=1}(x_i-\bar{x})^2} \cdot t_{n-2, 1-\alpha/2} &\; ,\, \hat{\beta} + s\, \sqrt \frac{1}{\sum^n_{i=1}(x_i-\bar{x})^2} \cdot t_{n-2,1- \alpha/2} \Big]. \end{align*} Die oben gemachte Bemerkung über die Varianz von \hat{\beta} spiegelt sich hier in der Länge des Konfidenzintervalls wider, das umgekehrt proportional zu \sqrt{\sum_{i=1}^n (x_i-\bar{x})^2} ist.

Anmerkung (Konfidenzintervalle für die Regressionsparameter in R).

Ein (1-\alpha)-Konfidenzintervall für den Regressionsparameter \beta_j ist gegeben durch \hat{\beta}_j \pm t_{n-p, 1-\alpha/2} \, \sqrt{s_{y|x}^2\, ((X^tX)^{-1})_{jj}}, wobei der Faktor \sqrt{s_{y|x}^2\, ((X^tX)^{-1})_{jj}} genau die geschätzte Standardabweichung des Schätzers \hat{\beta}_j ist. Diese geschätzte Standardabweichung findet man in R in der Spalte neben dem Schätzwert für den Parameter \beta_j. So erhalten wir etwa für den Parameter \beta_2 das 95\%-Konfidenzintervall mit den Grenzen 0.065 \pm 0.032, also das Intervall [0.033, 0.097].

Aufgabe 2 Betrachten Sie erneut den Datensatz USairpllution, aber berücksichtigen Sie nur die erklärendenVariablen x_2 (Anzahl der Industrieunternehmen) und x_4 (mittlerejährliche Windgeschwindigkeit).

1. Bestimmen Sie die Schätzwerte für die Parameter \beta_2 und \beta_4 im linearen Modell Y_i=\beta_0+\beta_2\, x_{i2}+\beta_4 \, x_{i4} +\epsilon_i.

2. Bestimmen Sie ein 95\%-Konfidenzintervall für den Parameter \beta_2.

Lösungen der Aufgaben

Lösung zu Aufgabe 1

Den erwarteten Schwefeldioxidgehalt für Albany berechnen wir, indem wir in die Modellgleichung \widehat{Y} (x_1,x_2,x_3,x_4,x_5,x_6)= 111.73 -1.27\, x_1 +0.06\, x_2 -0.04 \, x_3 -3.18\, x_4 +0.51\, x_5 -0.05\, x_6 die Werte der sechs erklärenden Variablen für Albany einsetzen, d.h. x_{11}=47.6, x_{12}=44, x_{13}=116, x_{14}=8.8, x_{15}=33.37 und x_{16}=134. Auf dem Wege erhalten wir folgenden erwarteten Wert \hat{y}_1=31.61, den wir mit dem tatsächlich beobachteten Schwefeldioxidgehalt y_1=46 vergleichen können. Die Differenz der beiden Werte, also y_1-\hat{y}_1=14.39, wird durch das lineare Modell nicht erklärt.

zurück zu Aufgabe 1

Lösung zu Aufgabe 2

1. Wir bestimmen die Schätzwerte für die Parameter inRwie folgt:

attach(USairpollution)
summary(lm(SO2~manu+wind))

Call:
lm(formula = SO2 ~ manu + wind)

Residuals:
    Min      1Q  Median      3Q     Max 
-25.748 -13.224  -3.788   6.150  68.433 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 26.984360  19.521702   1.382    0.175    
manu         0.027476   0.005301   5.183 7.49e-06 ***
wind        -1.022835   2.090935  -0.489    0.628    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.35 on 38 degrees of freedom
Multiple R-squared:  0.4194,    Adjusted R-squared:  0.3888 
F-statistic: 13.72 on 2 and 38 DF,  p-value: 3.265e-05

So erhalten wir folgende (gerundeten) Schätzwerte für die Regressionsparameter: \widehat{\beta}_0=26.984,\; \widehat{\beta}_2=0.027,\; \widehat{\beta}_4=-1.023. 2. Ein 95\%-Konfidenzintervall für \beta_2 finden wir gemäß der Formel „Schätzwert \pm\, t_{n-p,0.975}\cdot Geschätzte Standardabweichung des Schätzers”. Es gilt hier n=41, p=3 und t_{38,0.975}=2.02, sodass das 95\% Konfidenzintervall gegeben ist durch 0.027476 \pm 2.02 \cdot 0.005301 =[0.017,0.038]. Eine Anwendung der Faustregel „Schätzwert \pm\, 2\cdot Geschätzte Standardabweichung des Schätzers ”, hätte dasselbe Resultat geliefert.

zurück zu Aufgabe 2

Autor:innen

Die Lerneinheit “Multiple lineare Regression” wurde von Herold Dehling und Daniel Meißner unter Mithilfe von Elias Kaiser an der Ruhr-Universität Bochum entwickelt. Es ist lizenziert unter der CC-BY-SA 4.0 Lizenz und ist verfügbar auf ORCA.nrw.