fml <- y ~ 0 + x1 + x2
class(fml)[1] "formula"7 Lineare Modelle
Prognosen sind schwierig, besonders wenn sie die Zukunft betreffen.
Unbekannt
Lineare Modelle haben zum Ziel, eine Zielvariable (auch abhängige Variable genannt) mit Hilfe von Prädiktoren (unabhängige Variablen) zu erklären und vorherzusagen.
Die einfachste Form ist die Einfache Lineare Regression. Sie findet in Situationen Anwendung, bei denen ein linearer Zusammenhang zwischen einem unabhängigen und einem abhängigen Merkmal unterstellt werden kann. Genauer wird ein Modell der Form \[ Y = \beta_0 + \beta_1 X + \varepsilon \] zugrunde gelegt, wobei \(x\) die unabhängige und \(y\) die abhängige Variable ist, \(\beta_0\) und \(\beta_1\) die (unbekannten) Koeffizienten des linearen Zusammenhangs und \(\varepsilon\) ist eine Zufallsvariable. Das Merkmal \(Y\) kann also nicht vollständig über den linearen Zusammenhang durch \(X\) erklärt werden sondern nur ‘bis auf einen zufälligen Rest’. In der Regel wird außerdem angenommen, das \(\varepsilon\) normalverteilt ist. Wir folgen dieser Annahmen im weiteren Verlauf dieses Abschnitts. Ziel der linearen Regression ist nun insbesondere, die Koeffizienten \(\beta_0\) und \(\beta_1\) aus erhobenen Daten \((x_i,y_i), i=1,\dots,n\) für die beiden Merkmale zu schätzen.
Das Einfache Lineare Regressionsmodell lässt sich zu einem Multiplen Regressionsmodell verallgemeinern, in dem mehr als nur eine Kovariable erlaubt sind (vgl. Kapitel 7.2.2). Sind die Kovariablen (ausschließlich) kategoriell spricht man auch von Analysis of Variance (vgl Kapitel 7.2.3). Für nicht-metrische oder nicht normalverteilte abhängige Variablen lässt sich das Konzept der Linearen Modelle weiter verallgemeinern – man spricht dann auch von Verallgemeinerten Linearen Modellen (vgl. Kapitel 7.3), für die die Logistische Regression ein Spezialfall darstellt.
Zur Beschreibung dieses Modellzusammenhangs nutzt R eine Formelnotation, mit der man die abhängige und die unabhängigen Variablen festlegen kann. Kapitel 7.1 befasst sich zunächst mit der Formelnotation, bevor im Anschluss die unterschiedlichen Modelltypen beleuchtet werden.
7.1 Formeln
Einfache Formeln
Eine Formel ist ein Objekttyp in R, mit dem sich Modellgleichungen spezifizieren lassen. Die Tilde ~ ist der Kernbestandteil solch einer Formel, welche die abhängige von den unabhängigen Variablen trennt: die abhängige Variable steht immer links der Tilde, während die unabhängigen Variablen rechts stehen. Die Formel für ein einfaches Modell mit y als Zielvariable und x als Prädiktor schreibt man also als y ~ x. Mehrere einfache abhängige Variablen gibt man mit einem + getrennt an, also z. B. mit Prädiktoren x1 und x2 ergibt sich y ~ x1 + x2. Dabei wird in dieser Formulierung davon ausgegangen, dass auch eine Konstante (im linearen Regressionsmodell also der y-Achsenabschnitt) berücksichtigt wird. Die dazugehörige Regressionsgerade hat also die Form \[
Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \varepsilon.
\]
Soll der y-Achsenabschnitt nicht berücksichtigt werden, erreicht man dies durch das Addieren einer 0 in der Formel, welche dann y ~ 0 + x_1 + x_2 lautet und der Regressionsgleichung \[
Y = \beta_1 X_1 + \beta_2 X_2 + \varepsilon
\] entspricht. Rechenoperationen für die Zielvariable oder Prädiktoren können direkt in der Formel verwendet werden, so entspricht beispielsweise y ~ x + x^2 + log(x) der Regressionsgleichung \[
Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \beta_3 \ln({X})+ \varepsilon.
\]
Möchte man komplexere Berechnungen machen, z. B. Kovariablen miteinander addieren, nutzt man in der Formel dazu die Identitätsfunktion I() und schreibt y ~ x1 + I(x2 + x3), was der Regressionsgleichung \[
Y = \beta_0 + \beta_1 X_1 + \beta_2 (X_2 + X_3) + \varepsilon
\] entspricht.
Formeln werden als Formelobjekt abgespeichert:
Man kann Formeln auch mit Hilfe von as.formula aus Text erzeugen, was besonders praktisch ist, wenn man sehr viele Variablen miteinander verknüpfen möchte. Möchten wir beispielsweise die Ankunftsverspätung der New Yorker Flüge aus 2013 in Abhängigkeit von allen anderen Variablen des flights Datensatzes modellieren, machen wir dies mit
library(nycflights13)
target <- "arr_delay"
covariables <- names(flights)[-which(names(flights) == target)]
fml_txt <- paste0(target, " ~ ", paste(covariables, collapse = " + "))
fml <- as.formula(fml_txt)
fmlarr_delay ~ year + month + day + dep_time + sched_dep_time +
dep_delay + arr_time + sched_arr_time + carrier + flight +
tailnum + origin + dest + air_time + distance + hour + minute +
time_hourAls Platzhalter lässt sich in diesem Fall auch ein Punkt . verwenden. Wie wir später sehen werden, werden Formeln in Linearen Modellen immer im Zusammenhang mit einem Data Frame genutzt. Die Formel arr_delay ~ . wird im Zusammenhang mit dem Data Frame flights alle Spalten außer arr_delay als Kovariablen verwenden. Somit erzeugt
fml <- arr_delay ~ .
fmlarr_delay ~ .ebenfalls die Formel arr_delay ~ year + month + day + dep_time + sched_dep_time + dep_delay + arr_time + sched_arr_time + carrier + flight + tailnum + origin + dest + air_time + distance + hour + minute + time_hour + arr_delay_cat.
Interaktionseffekte
Interaktionseffekte (auch Wechselwirkungseffekte genannt) bezeichnen Effekte von Kombinationen zweier oder mehrerer Variablen. Eine Interaktion liegt vor, wenn der Effekt einer Kovariablen davon abhängt, welche Ausprägung eine andere Kovariable hat. Grafisch erkennt man dies häufig daran, dass die Linien in einem Gruppenvergleich nicht parallel verlaufen.
Abbildung 7.1 zeigt drei typische Fälle. Im gezeigten Beispiel untersuchen wir, wie sich ein Rabatt auf die Kaufwahrscheinlichkeit auswirkt und ob sich dieser Effekt zwischen Neu- und Bestandskunden unterscheidet.
Im ersten Panel verlaufen die Linien parallel. Das bedeutet: Der Rabatt wirkt in beiden Gruppen gleich stark. Zwar unterscheiden sich die beiden Kundentypen insgesamt in ihrer Kaufwahrscheinlichkeit, aber der Abstand zwischen den Linien bleibt konstant. Man könnte also sagen, dass es hier zwar einen Unterschied zwischen Neu- und Bestandskunden gibt und zugleich einen allgemeinen Rabatteffekt, dass sich diese beiden Effekte aber nicht gegenseitig verändern. Genau deshalb liegt in diesem Fall kein Interaktionseffekt vor.
Im zweiten Panel kreuzen sich die Linien. Hier sehen wir einen besonders starken Fall einer Interaktion: Der Effekt des Rabatts kehrt sich je nach Kundentyp sogar um. Bei den Neukunden steigt die Kaufwahrscheinlichkeit mit Rabatt deutlich an, bei den Bestandskunden sinkt sie dagegen. Deshalb lässt sich der Einfluss des Rabatts nicht mehr sinnvoll mit einer einzigen Zahl für alle Kunden gemeinsam beschreiben. Man muss den Rabatteffekt getrennt nach Kundentyp betrachten.
Im dritten Panel sind die Linien ebenfalls nicht parallel, sie kreuzen sich aber nicht. Auch hier hängt der Effekt des Rabatts also vom Kundentyp ab, nur weniger drastisch als im zweiten Fall. Der Rabatt wirkt in beiden Gruppen in dieselbe Richtung, aber nicht gleich stark: Bei Neukunden ist der Anstieg deutlich größer als bei Bestandskunden. Auch das ist bereits ein Interaktionseffekt, denn die Wirkung einer Kovariablen verändert sich abhängig von der Ausprägung der anderen.
Möchte man bei der Bestimmung von Regressionsmodellen in R Interaktionseffekte berücksichtigen, müssen diese in die Formel eingefügt werden. Für das Beispiel aus Abbildung 7.1 könnte man etwa kaufwahrscheinlichkeit ~ rabatt + kundentyp + rabatt:kundentyp schreiben; die Kurzform dafür ist kaufwahrscheinlichkeit ~ rabatt * kundentyp. Allgemein wird die Interaktion zwischen zwei Kovariablen mit Hilfe des Doppelpunkts : ausgedrückt. Ein Modell mit Kovariablen x1 und x2 sowie deren Interaktion schreiben wir also als y ~ x1 + x2 + x1:x2. In R bedeutet der Interaktionsterm dabei, dass sich der Effekt von x1 zwischen den Ausprägungen von x2 unterscheiden darf. Eine Kurzschreibweise für eine Formel mit den Haupteffekten und Interaktionseffekten von zwei Variablen wäre y ~ x1*x2. Die entsprechende Modellgleichung lautet dann \[
Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 \left( X_1 X_2 \right) + \varepsilon.
\] Auch Wechselwirkungseffekte zwischen mehr als zwei Kovariablen sind möglich, beispielsweise ergeben y ~ x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3 + x1:x2:x3 und y ~ (x1 + x2 + x3)^3 dasselbe Modell mit 2-fach und 3-fach Interaktionen. Die entsprechende Gleichung ist dann \[
\begin{split}
Y = \beta_0 & + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 \\
& + \beta_4 (X_1 X_2) + \beta_5 (X_1 X_3) + \beta_6 (X_2 X_3) \\
& + \beta_7 (X_1 X_2 X_3) + \varepsilon.
\end{split}
\]
Möchte man nur die 2-fach Interaktionen berücksichtigen, nutzt man y ~ x1 + x2 + x3 + x1:x2 + x1:x3 + x2:x3 oder in kurz y ~ (x1 + x2 + x3)^2, was dann \[
\begin{split}
Y = \beta_0 & + \beta_1 X_1 + \beta_2 X_2 + \beta_3 X_3 \\
& + \beta_4 (X_1 X_2) + \beta_5 (X_1 X_3) + \beta_6 (X_2 X_3) \\
& + \varepsilon
\end{split}
\] entspricht.
TippTipp
Die Formelnotation kann auch bei vielen anderen Funktionen in R außerhalb von statistischen Modellen verwendet werden. So haben wir in Kapitel 5.3.3 schon kennengelernt, wie wir mehrere Boxplots mit Hilfe einer Formel darstellen können.
7.2 Lineare Modelle
Lineare Modelle sind statistische Verfahren, mit denen eine (oder mehrere) beobachtete abhängige Variable(n) durch eine oder mehrere unabhängige Variablen erklärt wird. Diese Modelle sind in dem Sinne linear, als dass die abhängige Variable eine Linearkombination der Regressionskoeffizienten darstellt (aber nicht notwendigerweise der unabhängigen Variablen). Das bedeutet, die Modellgleichung stellt eine Linearkombination in den Koeffizienten dar, nicht zwangsweise aber in den Kovariablen. Wir haben im vorangegangenen Kapitel schon gesehen, dass auch z. B. quadratische Terme der Kovariablen erlaubt sind.
7.2.1 Einfaches Lineares Regressionsmodell
Berechnung des Modells
Das einfache lineare Regressionsmodell geht von lediglich zwei metrischen Größen aus: einer Einflussgröße \(X\) und einer Zielgröße \(Y\). Durch die einfache lineare Regression wird mithilfe von zwei Parametern eine Gerade so durch eine Punktwolke gelegt, dass der lineare Zusammenhang zwischen \(x\) und \(y\) möglichst gut beschrieben wird. Die Gleichung der linearen Einfachregression ist gegeben durch \[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \quad i = 1,\dots, n, \] wobei \(\varepsilon_i\) Fehlerterm zum Beobachtungspaar \((x_i, y_i)\) entspricht.
Für die Schätzung der unbekannten Parameter \(\beta_0\) und \(\beta_1\) erfolgt mit Hilfe der Methode der Kleinsten Quadrate. Dabei werden die Koeffizienten derart bestimmt, dass die Summe der quadratischen Residuen \[ \sum_{i=1}^n e_i = \sum_{i=1}^n (y_i-\hat{y}_i)^2 = \sum_{i=1}^n (y_i - \hat{\beta}_0 - \hat{\beta}_1 x_i)^2 \] minimiert wird. Dann bezeichnen \(\hat{\beta}_0\) und \(\hat{\beta}_1\) die Koeffizientenschätzer für \(\beta_0\) bzw. \(\beta_1\) und \(\hat{y}_i\) ist die Vorhersage der Regressionsgleichung für einen Wert \(x_i\) der Kovariable: \[ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i. \]
Um in R die Regressionsgerade nach der Methode der Kleinsten Quadrate zu berechnen, kann die Funktion lm() genutzt werden. Die wichtigsten Argumente der Funktion zeigt Tabelle 7.1.
lm Funktion
| Argument | Beschreibung |
|---|---|
formula |
Ein Formelobjekt, das die symbolische Beschreibung des Linearen Models angibt. |
data |
Ein Data Frame, der die Daten zu den in der Formel angegebenen Variablen enthält; die Spaltennamen des Data Frames müssen dazu den Namen der Variablen aus der Formel entsprechen. |
subset |
Ein Vektor, der die für die Modellschätzung zu verwendenden Zeilen angibt; dies kann zum Beispiel für die Auswahl eines Trainings- und Testdatensatzes genutzt werden. |
weights |
Ein Vektor mit Gewichten für die Beobachtungen des Data Frames. |
na.action |
Spezifizierung des Umgangs mit fehlenden Werten; standardmäßig werden Zeilen mit fehlenden Werten nicht für die Bestimmung des Modells genutzt. |
Betrachten wir nun als Beispiel die Daten des diamonds Datensatzes aus dem ggplot2 Paket. Dieser Datensatz enthält Preise und weitere Eigenschaften von etwa 54.000 Diamanten (für weitere Informationen zu den Variablen des Datensatzes vgl. ?ggplot2::diamonds). Zum Verständnis der Merkmale hilft auch diese Abbildung:
Wir wollen zunächst den Zusammenhang zwischen dem Preis eines Diamanten und der Masse des Edelsteins (gemessen in Karat) betrachten. Dazu stellen wir zunächst das Ergebnis in einem Scatterplot dar (vgl. Abbildung 7.2).
library(ggplot2)
ggplot(diamonds, aes(x = carat, y = price)) +
geom_point(alpha = 0.25, size = 0.5) +
labs(x = "Gewicht [in Karat]", y = "Preis [in US Dollar]")
Da es nur wenige Ausreißer mit mehr als 2,5 Karat gibt, wollen wir diese zunächst vernachlässigen. Außerdem können wir anhand der Abbildung erkennen, dass es sehr viele Diamanten mit kleiner Anzahl Karat und niedrigem Preis gibt. Um die Daten zu entzerren, transformieren wir beide Variablen mit dem natürlichen Logarithmus und verwenden diese transformierten Variablen anschließend auch für die Modellschätzung.
library(dplyr)
diamonds_new <- diamonds %>% filter(carat <= 2.5) %>%
mutate(lcarat = log(carat), lprice = log(price))
ggplot(diamonds_new, aes(x = lcarat, y = lprice)) +
geom_point(alpha = 0.25, size = 0.6) +
labs(x = "Log(Gewicht [in Karat])", y = "Log(Preis [in US Dollar])")
Die Log-Transformation hat dazu geführt, dass der Zusammenhang zwischen Gewicht und Preis linear und damit durch eine Regression gut abbildbar ist. Wir wollen nun ein lineares Modell bestimmen. Dazu nutzen wir die lm Funktion
lr_diamonds <- lm(lprice ~ lcarat, data = diamonds_new)
lr_diamonds
Call:
lm(formula = lprice ~ lcarat, data = diamonds_new)
Coefficients:
(Intercept) lcarat
8.452 1.681 Die Ausgabe des neuen Objekts lr_diamonds enthält eine knappe Zusammenfassung des Modells, namentlich die Koeffizientenschätzer für den y-Achsenabschnitt (engl. Intercept) und den Steigungsparameter des logarithmierten Gewichts. Diesem konkreten Ergebnis können wir entnehmen, dass der Preis mit wachsendem Gewicht des Diamanten zunimmt. Konkret nimmt der (logarithmierte) Preis durchschnittlich je zusätzlichem (logarithmierten) Karat um 1,681 USD zu.
WichtigAchtung
Durch das Logarithmieren der Zielgröße und der Kovariable lässt sich das geschätzte Regressionsmodell auch schreiben als \(Y = \alpha_0 \cdot X^{\alpha_1}\), denn
\[ \begin{split} & \log(Y) = \beta_0 + \beta_1 \log(X) \\ \Leftrightarrow \quad & Y = \exp\left(\beta_0 + \beta_1 \log(X)\right) \\ \Leftrightarrow \quad & Y = \underbrace{\exp(\beta_0)}_{\alpha_0} \cdot \exp\left(\beta_1 \log(X)\right) \\ \Leftrightarrow \quad & Y = \alpha_0 \cdot \exp\left(\log(X)\right)^{\beta_1} \\ \Leftrightarrow \quad & Y = \alpha_0 \cdot X^{\underbrace{{\beta_1}}_{\alpha_1}} \\ \Leftrightarrow \quad & Y = \alpha_0 \cdot X^{\alpha_1} \end{split} \]
Dieses Modell ist eigentlich kein Lineares Modell (es ist nicht Linear in den Koeffizienten \(\alpha_0\) und \(\alpha_1\)), allerdings mit Hilfe der Transformation dennoch als Lineares Modell darstellbar.
Die Koeffzienten erhalten damit auch eine neue Interpretation: der Steigungsparamter \(\beta_1\) bzw. \(\alpha_1\) ist somit eine Elastizität, d. h. wenn \(X\) um 10% steigt, so steigt \(Y\) um \(\beta_1\) * 10%.
Neben dieser Ausgabe der Modellzusammenfassung in der Konsole kann man auf die einzelnen Elemente direkt zugreifen. Die wesentlichen Elemente sind in Tabelle 7.2 aufgelistet.
lm Funktion
| Element | Beschreibung |
|---|---|
coefficients |
Ein Vektor mit den geschätzten Koeffizienten |
residuals |
Ein Vektor mit den Residuen |
fitted.values |
Die durch das Modell vorhergesagten Werte der Zielvariable |
rank |
Der Rang des angepassten Modells |
weights |
Die für die Modellierung genutzten Gewichte |
df.residual |
Die Freiheitsgrade der Residuen |
xlevels |
Die Faktorlevels der Faktoren |
y |
Die Zielvaraible |
x |
Die Modellmatrix |
model |
Der Model Frame |
na.action |
Der Umgang mit NAs |
Um beispielsweise auf die Koeffizienten zuzugreifen, nutzen wir das Element coefficients
lr_diamonds$coefficients(Intercept) lcarat
8.452142 1.681371 Ausführliche Zusammenfassung der Modellergebnisse
Eine ausführlichere Beschreibung des Modells liefert die summary Funktion:
summary(lr_diamonds)
Call:
lm(formula = lprice ~ lcarat, data = diamonds_new)
Residuals:
Min 1Q Median 3Q Max
-1.36139 -0.17016 -0.00585 0.16587 1.34114
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.452142 0.001364 6194.5 <2e-16 ***
lcarat 1.681371 0.001936 868.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2611 on 53812 degrees of freedom
Multiple R-squared: 0.9334, Adjusted R-squared: 0.9334
F-statistic: 7.542e+05 on 1 and 53812 DF, p-value: < 2.2e-16Zunächst wird im Abschnitt Call der Funktionsaufruf angegeben. Der Residuals Bereich gibt wesentliche Kennzahlen der Residuen an, mit denen z. B. die Normalverteilungsannahme geprüft werden kann. Anschließend wird unter Coefficients eine Tabelle mit den Details zu den Koeffizientschätzern angegeben. In der ersten Spalte Estimate werden die geschätzten Koeffizienten angegeben, mit Std. Error die Standardabweichung der Schätzung, außerdem der Wert der \(t\)-Teststatistik (t value) für die Hypothese \[
\textrm{H}_0: \beta_i = 0
\] versus \[
\textrm{H}_1: \beta_i \neq 0
\] für den Koeffizienten \(\beta_i\) und der dazugehörige p-Wert (Pr(>|t|)).
Die Signifikanzlevels werden je Koeffizienten mit Sternen hervorgehoben (siehe Legende unter der Tabelle der Koeffizienten): drei Sterne (***) stehen für einen p-Wert zwischen 0 und 0,001, also für hoch-signifikante Zusammenhänge. Zwei Sterne (**) für einen p-Wert zwischen 0,001 und 0,01, ein Stern für einen p-Wert zwischen 0,01 und 0,05. Das heißt, alle Koeffizienten mit mindestens einem Stern * haben einen signifikanten Effekt zum Niveau von 5% auf die Zielgröße. Ein Punkt . wird noch vergeben, wenn der p-Wert zwischen 0,05 und 0,1 liegt, ansonsten ist der Platz am Ende der entsprechenden Zeile leer.
Nach der Tabelle mit den Regressionskoeffizienten folgen in der summary noch Informationen zu den Residuen und Freiheitsgraden, zum Bestimmtheitsmaß (R-squared) und adjustierten Bestimmtheitsmaß (Adjusted R-squared).
Mit Hilfe der Bestimmtheitsmaße lässt sich erklären, inwiefern die abhängige Variable \(Y\) durch die unabhängige Variable \(X\) erklärt werden kann. Um zu quantifizieren, wie groß der durch den linearen Zusammenhang erklärte Anteil der Varianz von \(Y\) ist, werden zwei Größen miteinander verglichen:
- Summe der quadratischen Residuen (residual sum of squares): \[ \textrm{RSS} = \sum_{i=1}^n{(y_i - (\hat{\beta}_0 + \hat{\beta}_1 x_i))^2} \]
- Summe der quadratischen Abweichungen vom Mittelwert (total sum of squares): \[ \textrm{TSS} = \sum_{i=1}^n{(y_i - \bar{y})^2} \]
Die Größe \(\textrm{TSS}\) ist ein Maß für die gesamte Variabilität des Merkmals \(Y\), losgelöst von der Annahme der Abhängigkeit von \(X\). Die Größe \(\textrm{RSS}\) dagegen misst die Variabilität, die nicht durch den linearen Zusammenhang \(Y \approx \beta_0 + \beta_1 X\) erklärt werden kann. Um nun ein Maß für den Anteil zu erhalten, der durch die lineare Regression erklärt wird, wird der folgende Quotient gebildet: \[ R^2 = \frac{\text{TSS} - \text{RSS}}{\text{TSS}}. \] Die \(R^2\) Statistik ist also gerade derjenige Anteil an der gesamten Variabilität von \(Y\), welcher durch die Regression erklärt wird. Dieser Wert liegt immer zwischen \(0\) und \(1\). Dabei ist ein hoher Wert (nahe bei \(1\)) ein Indiz dafür, dass die Variable \(Y\) größtenteils durch die Variable \(X\) erklärt werden kann. Das Bestimmtheitsmaß allerdings hat die Eigenschaft, dass es umso größer wird, je größer die Zahl der unabhängigen Variablen ist. Hierfür ist unerheblich, ob die zusätzlichen unabhängigen Variablen einen Beitrag zur Erklärungskraft liefern. Daher ist es ratsam, das adjustierte Bestimmtheitsmaß (auch bereinigtes, korrigiertes oder angepasstes Bestimmtheitsmaß genannt) zu Rate zu ziehen. Das adjustierte Bestimmtheitsmaß wird berechnet als \[ R^2_{\textrm{adj}} = 1-(1-R^2) \frac{n-1}{n-p-1}. \] Definitionsgemäß ist das adjustierte Bestimmtheitsmaß für mehr als eine erklärende Variable stets kleiner als das unadjustierte. Beim adjustierten Bestimmtheitsmaß wird die Erklärungskraft des Modells – repräsentiert durch \(R^2\) – ausbalanciert mit der Komplexität des Modells, repräsentiert durch die Anzahl der Parameter (also allgemein \(p+1\)). Je komplexer das Modell ist, desto mehr „bestraft“ das adjustierte Bestimmtheitsmaß jede neu hinzugenommene erklärende Variable.
Neben den beiden Bestimmheitsmaßen wird eine weitere Kennzahl zur Beurteilung der Modellgüte insgesamt ausgegeben, nämlich die \(F\)-Statistik und deren p-Wert (F-statistic). Mit Hilfe des \(F\)-Tests wird an dieser Stelle Die Nullhypothese \[
\textrm{H}_0: \beta_0 = \beta_1 = \dots = \beta_p = 0
\] gegen die Alternative \[
\textrm{H}_1: \exists i: \beta_i \neq 0
\] geprüft. Das bedeutet, es wird getestet, ob überhaupt irgendeiner der Koeffizienten des Modells von Null verschieden ist.
In dem vorliegenden Beispiel können wir erkennen, dass sowohl der y-Achsenabschnitt als auch der Effekt des Gewichts hoch-signifikant sind und einen p-Wert kleiner als 0,001 aufweisen. Das Bestimmtheitsmaß \(R^2\) weist einen sehr hohen Wert aus, es liegt bei 0,9334 – die Anpassung des Modells ist also sehr gut gelungen. Abbildung 7.4 zeigt den Scatterplot der (logarithmierten) Karatzahl und (logarithmierten) Preise mit eingezeichneter Regressionsgeraden (hier über geom_smooth(method = "lm")) und macht deutlich, dass der Zusammenhang durch die Regressionsgerade sehr gut angepasst wird. Das bedeutet, dass der Preis eines Diamanten ganz wesentlich von dessen Gewicht abhängt.
ggplot(diamonds_new, aes(x = lcarat, y = lprice)) +
geom_point(alpha = 0.25, size = 0.6) +
geom_smooth(method = "lm", se = FALSE, color = "red", linewidth = 1) +
labs(x = "Log(Gewicht [in Karat])", y = "Log(Preis [in US Dollar])")
Auch die summary des Modells enthält mehrere Elemente, auf die wir zugreifen können. Diese sind in Tabelle 7.3 aufgelistet.
summary.lm Funktion
| Element | Beschreibung |
|---|---|
residuals |
Die Residuen des Modells |
coefficients |
Die Koeffizientenmatrix mit Koeffizienten, Standardfehlern, \(t\)-Statistik und p-Wert |
fstatistic |
Die F-Statistik des Modells |
r.squared |
Das Bestimmtheitsmaß |
adj.r.squared |
Das adjustierte Bestimmtheitsmaß |
So erhalten wir beispielsweise die Koeffizientenmatrix mit den Befehlen
lr_diamonds_summary <- summary(lr_diamonds)
lr_diamonds_summary$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.452142 0.001364466 6194.4689 0
lcarat 1.681371 0.001936057 868.4514 0str(lr_diamonds_summary$coefficients) num [1:2, 1:4] 8.45 1.68 1.36e-03 1.94e-03 6.19e+03 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:2] "(Intercept)" "lcarat"
..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"Modelldiagnostik
Für lm Objekte existiert auch eine plot Funktion (also plot.lm). Mit dieser lassen sich diagnostische Plots der Regressionsanalyse abbilden.
Details zu den Plots findet man in der Hilfe zu plot.lm. Mit den Plots lässt sich feststellen, ob Homoskedastizität der Residuen vorliegt (Residuals vs Fitted, Scale-Location), ob sie normalverteilt sind (Normal Q-Q) und ob es Ausreißer in den Daten gibt, die das Modell stark beeinflussen (Cook’s distance).
HinweisLink
Weitere Informationen zur Diagnostik von Regressionsmodellen finden Sie bei Belsley, Kuh, und Welsch (1980) und Cook und Weisberg (1982). Eine Zusammenfassung gibt es unter https://data.library.virginia.edu/diagnostic-plots/.
Konfidenzintervalle und Vorhersagen
Konfidenzintervalle für die Koeffizienten erhält man mit der confint Funktion, also
confint(lr_diamonds, level = 0.95) 2.5 % 97.5 %
(Intercept) 8.449467 8.454816
lcarat 1.677576 1.685166Mit Hilfe der predict Methode können wir Vorhersagen eines Modells erstellen. Die Funktion erwartet dazu als Input neben dem Regressionsmodell außerdem noch mit dem Argument newdata einen Data Frame, der die je Kovariable des Modells eine Spalte enthält. Möchten wir beispielsweise wissen, welchen (logarithmierten) Preis wir bei einem (logarithmierten) Gewicht von 0.5 erwarten können, erhalten wir dies mit Hilfe von
new_data <- data.frame(lcarat = 0.5)
predict(lr_diamonds, newdata = new_data) 1
9.292827 Auch für die Vorhersagen lassen sich Konfidenzintervalle ausgeben, so zum Beispiel für das Konfidenzlevel von 95% mit
new_data <- data.frame(lcarat = 0.5)
predict(lr_diamonds, newdata = new_data, interval = "predict", level = 0.95) fit lwr upr
1 9.292827 8.780985 9.80467Wir nutzen diese Funktion nun, um den Zusammenhang zwischen Preis und Karat (ohne Logarithmus) darzustellen (vgl. Abbildung 7.6).
new_data <- data.frame(lcarat = log(seq(0,2.5,0.01)))
new_data[c("lprice", "lprice_lwr", "lprice_upr")] <- predict(lr_diamonds, newdata = new_data, interval = "predict", level = 0.95)
new_data <- new_data %>% mutate(carat = exp(lcarat), price = exp(lprice),
price_lwr = exp(lprice_lwr), price_upr = exp(lprice_upr))
ggplot(diamonds_new, aes(x = carat, y = price)) +
geom_ribbon(
data = new_data,
aes(x = carat, ymin = price_lwr, ymax = price_upr),
inherit.aes = FALSE,
fill = "grey80",
alpha = 0.8
) +
geom_point(alpha = 0.2, size = 0.5) +
geom_line(data = new_data, aes(x = carat, y = price), color = "red", linewidth = 1) +
labs(x = "Gewicht [in Karat]", y = "Preis [in US Dollar]")Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_ribbon()`).
7.2.2 Multiples Lineares Regressionsmodell
Ein Multiples Lineares Regressionsmodell verwendet mehrere (allgemein \(p\)) unabhängige Variablen und ein abhängiges Merkmal. Die dazugehörige Regressionsgleichung lautet
\[ \begin{split} y_i & = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip} + \varepsilon_i \\ & = \mathbf{x}^{\top}{\beta} + \varepsilon_i, \quad\quad i = 1,\dots, n \end{split} \] bzw. in Matrixschreibweise (alle \(i\) Beobachtungen zusammenfassend) \[ \mathbf{y} = \mathbf{X} {\beta} + {\varepsilon} \] mit \[ \begin{split} \mathbf{y} & = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}, \quad \mathbf{X} = \begin{pmatrix} \mathbf{x}_1^{\top} \\ \mathbf{x}_2^{\top} \\ \vdots \\ \mathbf{x}_n^{\top} \end{pmatrix} = \begin{pmatrix} x_{11} & x_{12} & \dots & x_{1p} \\ x_{11} & x_{12} & \dots & x_{1p} \\ \vdots & \vdots & \ddots & \vdots \\ x_{n1} & x_{n2} & \dots & x_{np} \end{pmatrix},\\ {\beta} & = \begin{pmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_p \end{pmatrix}, \quad {\varepsilon} = \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{pmatrix}. \end{split} \]
WichtigAchtung
Häufig wird im Zusammenhang von Regressionsanalysen mit mehreren unabhängigen Variablen allgemein von Multivariater Regression gesprochen. Es wird allerdings zwischen zwei Fällen unterschieden:
- Ein Regressionsmodell mit einer abhängigen und mehreren unabhängigen Variablen nennt man Multiples Regressionsmodell (engl. multivariable regression).
- Ein Regressionsmodell mit mehreren abhängigen und mehreren unabhängigen Variablen nennt man Multivariates Regressionsmodell (engl. multivariate regression).
Im Beispiel der Diamanten können wir einfach zusätzliche Kovariablen (depth, table) hinzufügen und erhalten das entsprechende Modell und die Ergebnisse über die summary
lr_diamonds <- lm(lprice ~ lcarat + depth + table, data = diamonds_new)
summary(lr_diamonds)
Call:
lm(formula = lprice ~ lcarat + depth + table, data = diamonds_new)
Residuals:
Min 1Q Median 3Q Max
-1.24492 -0.16698 -0.00545 0.16280 1.35321
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.2843453 0.0664152 169.91 <2e-16 ***
lcarat 1.6969244 0.0019469 871.61 <2e-16 ***
depth -0.0285355 0.0008135 -35.08 <2e-16 ***
table -0.0185189 0.0005304 -34.91 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2567 on 53810 degrees of freedom
Multiple R-squared: 0.9356, Adjusted R-squared: 0.9356
F-statistic: 2.608e+05 on 3 and 53810 DF, p-value: < 2.2e-16Durch die Hinzunahme der Variablen erhöht sich das Bestimmtheitsmaß leicht von 0,9334 auf 0,9356.
In diesem Beispiel wurden alle Variablen als stetig behandelt. Das heißt, der Koeffizient für depth besagt, das je Prozentpunkt an Tiefe der Preis um 0.03 USD abnimmt.
Möchten wir kategorielle Variablen im Regressionsmodell verwenden, können wir Dummy-Variablen durch Dichotomisierung erzeugen. Dies wird von R automatisch vorgenommen, wenn man Variablen als factor im Modell übergibt1.
Wir verwenden nun also zusätzlich cut (Qualität des Schliffs), color (Qualität der Farbe) und clarity (Klarheit des Diamanten) als kategorielle Kovariablen und vernachlässigen depth und table. Zunächst betrachten wir die Ausprägungen der drei kategoriellen Variablen:
table(diamonds_new$cut)
Fair Good Very Good Premium Ideal
1589 4889 12063 13745 21528 table(diamonds_new$color)
D E F G H I J
6771 9792 9538 11280 8269 5389 2775 table(diamonds_new$clarity)
I1 SI2 SI1 VS2 VS1 VVS2 VVS1 IF
707 9116 13055 12256 8169 5066 3655 1790 Es gibt also für den cut 5 Level (von ‘Fair’ bis ‘Ideal’), für die color die 7 Ausprägungen D-E und für clarity die 8 Stufen I1 bis IF. Alle kategoriellen Variablen sind dabei als geordnete Kategorien (also ordinal skaliert) hinterlegt. Zur Vereinfachung entfernen wir diese Ordnung:
diamonds_new$cut <- factor(diamonds_new$cut, ordered = FALSE)
diamonds_new$color <- factor(diamonds_new$color, ordered = FALSE)
diamonds_new$clarity <- factor(diamonds_new$clarity, ordered = FALSE)Dann können wir das Regressionsmodell bestimmen.
lr_diamonds <- lm(lprice ~ lcarat + cut + color + clarity, data = diamonds_new)
summary(lr_diamonds)
Call:
lm(formula = lprice ~ lcarat + cut + color + clarity, data = diamonds_new)
Residuals:
Min 1Q Median 3Q Max
-0.81367 -0.08621 -0.00065 0.08263 1.92918
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.878793 0.005824 1352.76 <2e-16 ***
lcarat 1.886239 0.001124 1677.81 <2e-16 ***
cutGood 0.079667 0.003879 20.54 <2e-16 ***
cutVery Good 0.116721 0.003609 32.34 <2e-16 ***
cutPremium 0.139116 0.003570 38.97 <2e-16 ***
cutIdeal 0.160826 0.003539 45.45 <2e-16 ***
colorE -0.054404 0.002103 -25.87 <2e-16 ***
colorF -0.095218 0.002126 -44.78 <2e-16 ***
colorG -0.160878 0.002082 -77.29 <2e-16 ***
colorH -0.250963 0.002210 -113.56 <2e-16 ***
colorI -0.372398 0.002477 -150.36 <2e-16 ***
colorJ -0.511213 0.003061 -166.99 <2e-16 ***
claritySI2 0.408518 0.005248 77.84 <2e-16 ***
claritySI1 0.572443 0.005216 109.75 <2e-16 ***
clarityVS2 0.721966 0.005243 137.71 <2e-16 ***
clarityVS1 0.792239 0.005319 148.94 <2e-16 ***
clarityVVS2 0.927682 0.005475 169.45 <2e-16 ***
clarityVVS1 0.999500 0.005626 177.65 <2e-16 ***
clarityIF 1.094550 0.006071 180.30 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1328 on 53795 degrees of freedom
Multiple R-squared: 0.9828, Adjusted R-squared: 0.9828
F-statistic: 1.706e+05 on 18 and 53795 DF, p-value: < 2.2e-16Wir können erkennen, dass sich das Modell insgesamt nochmals verbessert und jetzt ein Bestimmtheitsmaß von 0,9828 aufweist.
Wir erhalten allerdings nicht für jeden Kategorie der neuen factor Variablen einen Koeffizientenschätzer, da diese nicht identifizierbar wären. Da sämtliche Dummyvariablen eines kategoriellen Merkmals zusammen linear abhängig wären, können die Koeffizienten im Regressionsmodell dann nicht mehr geschätzt werden. R wählt standardmäßig das erste Level jedes Factors als entsprechende Referenzkategorie. Die Koeffizientenschätzer einer Dummy-Variable einer Regression entsprechen damit den Abweichungen der Gruppenmittelwerte von der Referenzgruppe. Dummy-Variablen lassen sich auch anders Definieren, über sogenannte Kontraste. Diese lassen sich mit Hilfe der contrasts Funktion spezifizieren oder über das contrasts Argument der lm Funktion. So kann man beispielsweise festlegen, dass es nicht eine Referenzkategorie gibt, sondern dass sich die Effekte der kategoriellen Variable zu 0 aufsummieren.
lr_diamonds <- lm(lprice ~ lcarat + cut + color + clarity, data = diamonds_new,
contrasts = list(cut = "contr.sum", color = "contr.sum", clarity = "contr.sum"))Konfidenzintervalle mit confint und Vorhersagen mit predict funktionieren in der Multiplen Linearen Regression übrigens genau wie im einfachen Fall. So können wir beispielsweise für alle Diamanten den erwarteten Preis als Vorhersage durch das Regressionsmodell bestimmen und herausfinden, welche Diamanten weiter unter Wert verkauft wurden (vgl. Tabelle 7.4).
pricediff_diamonds <- exp(predict(lr_diamonds)) - diamonds_new$price| carat | cut | color | clarity | depth | price | price_expected | price_diff |
|---|---|---|---|---|---|---|---|
| 2.46 | Premium | E | SI2 | 59.7 | 10470 | 23622.21 | 13152.206 |
| 2.38 | Premium | H | VS1 | 60.9 | 18522 | 26762.16 | 8240.157 |
| 2.50 | Premium | H | SI2 | 62.7 | 12071 | 20006.41 | 7935.406 |
| 2.31 | Premium | J | VVS1 | 62.1 | 16100 | 23991.26 | 7891.258 |
| 2.50 | Ideal | I | VS2 | 61.9 | 16955 | 24773.52 | 7818.520 |
| 2.29 | Premium | J | IF | 61.4 | 18426 | 25954.29 | 7528.288 |
| 2.29 | Premium | J | IF | 61.4 | 18594 | 25954.29 | 7360.288 |
| 2.42 | Premium | I | VS1 | 59.8 | 17262 | 24458.85 | 7196.846 |
| 2.37 | Premium | I | VS2 | 61.6 | 14837 | 21918.66 | 7081.660 |
| 2.48 | Premium | I | VS2 | 60.1 | 16820 | 23876.97 | 7056.970 |
7.2.3 Varianzanalyse (ANOVA)
Bei der (einfaktoriellen) Varianzanalyse geht es darum, zu entscheiden, ob die Erwartungswerte mehrerer Verteilungen übereinstimmen oder nicht. Üblicherweise liegen hierbei also mehrere Stichproben vor, von denen angenommen wird, dass sie aus Normalverteilungen mit identischer Standardabweichung stammen. Die zu beantwortende Frage ist dann, ob auch die Erwartungswerte übereinstimmen.
Bezeichnen wir die Erwartungswerte der zugrunde liegenden theoretischen Verteilungen mit \(\mu_1,\mu_2, \dots, \mu_k\) , so wollen wir die Hypothese \[ \textrm{H}_0 :\mu_1 =\mu_2 =\dots=\mu_k \] gegen die Alternative \[ \textrm{H}_1 : \exists i, j \in \{ 1,2,\dots,k\}: \mu_i \neq \mu_j \] testen. Die Teststatistik dieses Tests folgt der \(F\)-Verteilung und wird deshalb auch \(F\)-Test genannt2. Eine typische Anwendungssituation für solch einen Mittelwertvergleich stellt das Diamantenbeispiel dar. Möchten wir herausfinden, ob die Farbe einen signifikanten Effekt auf den Preis eines Diamanten hat, interessiert uns also, ob der durchschnittliche Preis aller Diamanten einer bestimmten Farbe sich von den anderen unterscheidet. Dies entspricht also der Nullhypothese des \(F\)-Tests, wenn wir die Diamanten in Gruppen je nach Farbe sortieren und deren durchschnittliche Preise betrachten.
Bei einer allgemeinen (mehrfaktoriellen) ANalysis Of VAriance (ANOVA) werden solche Vergleiche für mehrere kategoriellen gleichzeitig angestellt. Wir erhalten dann für jede kategorielle Variable das Ergebnis eines \(F\)-Tests, der besagt, ob eines der Faktorlevels einen signifikant unterschiedlichen Effekt von den anderen Levels hat. Das ist gleichbedeutend mit der Frage, ob die kategorielle Variable insgesamt einen signifikanten Einfluss auf die Zielgröße hat.
In R können wir uns eine solche Varianzanalysetafel aus einem lm Objekt mit der Funktion anova erstellen
lr_diamonds <- lm(lprice ~ lcarat + cut + color + clarity, data = diamonds_new)
anova(lr_diamonds)Analysis of Variance Table
Response: lprice
Df Sum Sq Mean Sq F value Pr(>F)
lcarat 1 51430 51430 2917226 < 2.2e-16 ***
cut 4 223 56 3169 < 2.2e-16 ***
color 6 684 114 6463 < 2.2e-16 ***
clarity 7 1814 259 14699 < 2.2e-16 ***
Residuals 53795 948 0
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1oder anstelle der lm Funktion die Funktion aov und anschließend summary nutzen
lr_diamonds <- aov(lprice ~ lcarat + cut + color + clarity, data = diamonds_new)
summary(lr_diamonds) Df Sum Sq Mean Sq F value Pr(>F)
lcarat 1 51430 51430 2917226 <2e-16 ***
cut 4 223 56 3169 <2e-16 ***
color 6 684 114 6463 <2e-16 ***
clarity 7 1814 259 14699 <2e-16 ***
Residuals 53795 948 0
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Das Ergebnis ähnelt der Ausgabe eines Linearen Regressionsmodells. Ein Unterschied ist, dass hier nun nur noch eine Zeile je Variable ausgegeben wird (auch für kategorielle Variablen werden nicht alle Dummy-Variablen angezeigt). In jeder Zeile wird das Ergebnis des \(F\)-Tests angegeben (Teststatistik F value und p-Wert Pr(>F)) und mit der ‘Sternchen-Notation’ versehen.
Anstelle der Koeffizientenschätzer wird eine Quadratsumme ausgegeben (Sum Sq). Diese lässt sich ähnlich wie das Bestimmtheitsmaß interpretieren: Die Gesamtquadratsumme \(\textrm{TSS}\) (Summe über die Spalte Sum Sq) entspricht \[
\textrm{TSS} = \sum_{i=1}^n (y_i-\bar{y})^2
\] und die Quadratsumme der Residuen (Residuals) ist \[
\textrm{RSS} = \sum_{i=1}^n (y_i-\hat{y})^2,
\] wobei \(\hat{y_i}\) der Vorhersage des Modells für \(x_i\) entspricht. Wie beim Bestimmtheitsmaß ist also der Teil der Quadratsumme, der durch das Modell erklärt wird \[\textrm{TSS}-\textrm{RSS}.\] Diese Quadratsumme lässt sich nun weiter aufteilen, nämlich in den Teil der durch jede einzelne Kovariable erklärt wird. Das bedeutet, dass im obigen Beispiel die Quadratsumme der Klarheit (clarity) eines Diamanten also \[
\textrm{SS}_{\texttt{clarity}} = 1814
\] entspricht. Da \[
\begin{split}
\textrm{TSS} & = \textrm{SS}_{\texttt{lcarat}} + \textrm{SS}_{\texttt{cut}} + \textrm{SS}_{\texttt{color}} + \textrm{SS}_{\texttt{clarity}} + \textrm{RSS} \\
& = 51430 + 223 + 684 + 1814 + 984 \\
& = 55135,
\end{split}
\] ist der Anteil der durch clarity erklärten Quadratsumme/Varianz also ca. 3,3%, denn \[
\frac{\textrm{SS}_{\texttt{clarity}}}{\textrm{TSS}} = \frac{1814}{55135} = 0,0329.
\] Für nicht-kategorielle Variablen kann ebenso die Quadratsumme berechnet werden; für lcarat beträgt sie hier 51430, was einem Anteil an erklärter Varianz von 93,33% entspricht, denn \[
\frac{\textrm{SS}_{\texttt{lcarat}}}{\textrm{TSS}} = \frac{51460}{55135} = 0,9333.
\]
7.3 Verallgemeinerte Lineare Modelle
Mit der linearen Regression haben wir in den letzten Abschnitten eine Methode diskutiert, die zur Vorhersage eines abhängigen quantitativen Merkmals auf Basis von Beobachtungen unabhängiger Merkmale genutzt werden kann. Dabei wird in der Regel angenommen, dass die Zielgröße normalverteilt ist. Beide Annahmen – das abhängige Merkmal ist metrisch und es ist normalverteilt – sind in vielen Situationen nicht gegeben. Für diese Fälle haben McCullagh und Nelder (1989) Verallgemeinerte Lineare Modelle entwickelt (engl. Generalized Linear Models, kurz GLM). Wir wollen hierzu zunächst den Spezialfall der Logistischen Regression betrachten, bevor wir uns dem allgemeinen Fall nähern.
7.3.1 Logistische Regression
Oftmals ist man in Anwendungen daran interessiert ein qualitatives Merkmal auf Grundlage von unabhängigen Variablen vorherzusagen. Ein qualitatives Merkmal nimmt im Gegensatz zu quantitativen Merkmalen Werte an, die nicht numerisch sind, sondern aus einer endlichen Anzahl an Kategorien stammen. Zum Beispiel nimmt das Merkmal Geschlecht eine der beiden Kategorien weiblich oder männlich an. Im Zusammenhang mit Schrifterkennung könnte man das Merkmal Buchstabe betrachten, welches einen der 26 Buchstaben aus unserem Alphabet annehmen kann. Die Vorhersage eines abhängigen qualitativen Merkmals auf Basis anderer unabhängiger Merkmale wird als Klassifizierung oder auch Klassifikation bezeichnet. Beispiele für derartige Problemstellungen sind in ganz unterschiedlichen Bereichen zu finden. Um nur einige wenige zu nennen:
- In einer Notfallambulanz möchte man auf Basis von Symptomen entscheiden, welche von drei in Frage kommenden Diagnosen die richtige ist.
- Beim Online Banking möchte man auf Basis von bspw. IP-Adresse, Transaktions-Historie, etc. vorhersagen, ob bei einer durchgeführten Transaktion ein Betrug vorliegt.
- Auf Grundlage von DNA Sequenzen möchte man ermitteln, ob eine bestimmte Krankheit wahrscheinlich ausbrechen wird oder nicht.
- Auf Basis diverser Kundendaten möchte man vorhersagen, ob ein bestimmter Kunde sich in näherer Zukunft vermutlich vom Unternehmen abwenden wird oder nicht.
Eine Methode, die zur Klassifikation genutzt werden kann, ist die logistische Regression, welche wir im Folgenden diskutieren werden. Weitere Klassifikationsverfahren werden Sie in der Vorlesung ‘Data Mining’ kennen lernen.
Bei der logistischen Regression wird, wie bei vielen anderen Klassifizierungsmethoden auch, das interessierende qualitative Merkmal nicht direkt modelliert. Anstatt dessen wird ein Modell zur Prognose der Wahrscheinlichkeit dafür betrachtet, dass das Merkmal einer bestimmten Kategorie angehört.
Wir betrachten hier die Situation, dass das interessierende Merkmal \(Y\) (die Zielvariable) nur zwei mögliche Merkmalsausprägungen annehmen kann, welche mit \(1\) und \(0\) bezeichnet werden (z. B. ‘ja’ oder ‘nein’; ‘weiblich’ oder ‘männlich’; ‘wahr’ oder `falsch’ etc.). Es wird nun ein Modell betrachtet, welches die Wahrscheinlichkeit dafür, dass \(Y = 1\) für einen betrachteten Merkmalsträger gilt, als Funktion weiterer unabhängiger Merkmale \(X_1,\ldots,X_q\) (der Inputvariablen) des Merkmalsträgers beschreibt. Bezeichnen wir mit \(\mathbf{X}\) nun den Vektor der Merkmale \(X_1,\ldots,X_q\), und mit \(p(\mathbf{X})\) die von \(\mathbf{X}\) abhängende Wahrscheinlichkeit für \(Y=1\), dann wird \(p(\mathbf{X})\) wie folgt modelliert: \[ p(\mathbf{X}) = \text{sigm}(\beta_0+\beta_1X_1+\ldots+\beta_q X_q), \tag{7.1}\] wobei \[ \text{sigm}(t) = \frac{\exp(t)}{1 + \exp(t)} \tag{7.2}\] die sogenannte Sigmoid-Funktion ist. Mit Gleichung 7.1 und Gleichung 7.2 ergibt sich also \[ p(\mathbf{X}) = \frac{\exp(\beta_0+\beta_1X_1+\ldots+\beta_q X_q)}{1 + \exp(\beta_0+\beta_1X_1+\ldots+\beta_q X_q)}. \tag{7.3}\] Die Sigmoid Funktion nimmt ausschließlich Werte zwischen \(0\) und \(1\) an und ist monoton steigend (vgl. Abbildung 7.7); somit liegt \(p(\mathbf{X})\) also auch zwischen \(0\) und \(1\) und entspricht damit der Wahrscheinlichkeit für \(Y=1\).
Gleichung Gleichung 7.1 ist äquivalent zu \[ \log\left(\frac{p(\mathbf{X})}{1-p(\mathbf{X})}\right)= \beta_0+\beta_1X_1+\ldots+\beta_q X_q. \tag{7.4}\]
Die Linke Seite von Gleichung 7.4 nennt man auch die Logit Funktion (vgl. auch Abbildung 7.8).
Das Ziel ist nun, die unbekannten Koeffizienten \(\beta_0,\ldots, \beta_q\) aus einer vorliegenden Stichprobe der Merkmale \(X_1,\ldots, X_q, Y\) zu schätzen. Wir bezeichnen mit \((x_{1i}, \ldots, x_{qi}, y_i)\) die Realisierungen der Merkmale des \(i\)-ten Untersuchungsobjekts, wobei \(i=1,\ldots,n\) und \(n\) ist der Stichprobenumfang.
Mithilfe der Schätzungen für \(\beta_0,\ldots, \beta_q\) kann dann die modellierte Wahrscheinlichkeit \(p(\mathbf{X})\) prognostiziert werden, um auf Grundlage dieser Prognose zu entscheiden, ob die Outputvariable vermutlich den Wert \(1\) oder doch den Wert \(0\) annimmt. Zum Beispiel könnte man ab einem prognostizierten Wert von \(0.5\) für die Wahrscheinlichkeit \(p(\mathbf{X})\) zur Vorhersage der Output Variable den Wert \(1\) festlegen. Bei einer prognostizierten Wahrscheinlichkeit unter \(0.5\) würde man entsprechend den Wert \(0\) für die Variable \(Y\) vorhersagen.
Die bei der logistischen Regression verwendete Methode zum Schätzen der Parameter wird als Maximum Likelihood Methode (auch kurz ML für Maximum Likelihood) bezeichnet. Hierbei werden die Schätzwerte \(\hat{\beta_0}, \ldots, \hat{\beta_q}\) für die unbekannten Modell-Parameter derart gewählt, dass das Produkt \[ \prod_{y_i=1}p(x_{1i},\ldots,x_{qi}) \prod_{y_i = 0} (1-p(x_{1i},\ldots,x_{qi})) \] maximiert wird. Dieses Produkt entspricht gerade der Wahrscheinlichkeit dafür, dass die vorliegende Stichprobe für das Merkmal \(Y\) bei der Erhebung beobachtet wird. Der Maximum Likelihood Ansatz bestimmt die unbekannten Parameter so, dass diese die wahrscheinlichsten Parameter sind, welche die realisierte Stichprobe unter der Modellannahme erzeugen. Das bedeutet, dass die vorliegende Stichprobe beobachtet wurde, weil die Wahrscheinlichkeit dafür groß ist, dass genau diese Stichprobe realisiert wird. Daher werden die Schätzwerte für die Parameter gerade so gewählt, dass das obige Produkt maximal wird.
Aus dem ML Ansatz lassen sich neben den Koeffizientenschätzern \(\hat{\beta_0}, \ldots, \hat{\beta_q}\) auch Teststatistiken und dazugehörige p-Werte für die Koeffizienten bestimmen.
Anwendung auf simulierte Daten
Bevor wir im nächsten Abschnitt einen konkreten Anwendungsfall für die logistische Regression betrachten, werden wir nun zunächst anhand eines abstrakten Beispiels das Verfahren in R illustrieren. Hierzu betrachten wir einen simulierten Datensatz, der zwei unabhängige quantitative Merkmale \(X_1\) und \(X_2\) umfasst sowie ein abhängiges qualitatives Merkmal \(Y\), welches die Werte \(1\) oder \(0\) annimmt. So könnten \(X_1\) und \(X_2\) zwei metrische Merkmale von Kunden sein (Alter und Einkommen) und die Zielgröße \(Y\) die Information enthalten, ob ein Kunde gekündigt hat oder nicht. Dann können wir die beiden metrischen Merkmale in einem Scatterplot darstellen und das Kündigungsverhalten durch eine Einfärbung der Punkte kenntlich machen (\(Y=1\) als rot, \(Y = 0\) als blau, vgl. Abbildung 7.9).
Was wir nun mit Hilfe der Logistischen Regression suchen ist eine Möglichkeit, die blauen Punkte möglichst trennscharf von den roten Punkten zu trennen.
Schätzen der Parameter und Prognose
Die logistische Regression wird in R mit der Funktion glm durchgeführt. Diese Funktion erwartet als erstes Argument eine Formel, die das zu verwendende Modell beschreibt. Da die Funktion glm eine sehr generische Funktion ist, welche auch zur Anpassung vieler anderer linearer Modelle verwendet werden kann, muss zudem für die logistische Regression das Argument family auf binomial(link = "logit") gesetzt werden.
fit1 <- glm(Y ~ X1 + X2, family = binomial(link = "logit"), data = df_logit)
fit1
Call: glm(formula = Y ~ X1 + X2, family = binomial(link = "logit"),
data = df_logit)
Coefficients:
(Intercept) X1 X2
-35.475 2.947 1.048
Degrees of Freedom: 149 Total (i.e. Null); 147 Residual
Null Deviance: 207.9
Residual Deviance: 61.35 AIC: 67.35Das von glm erzeugte Objekt enthält nun Berechnungsergebnisse, die im Zusammenhang mit der logistischen Regression relevant sind. Die Ausgabe auf der Konsole liefert insbesondere die Schätzwerte für die Parameter \(\beta_0\) ((Intercept)), \(\beta_1\) und \(\beta_2\). Die summary Funktion liefert wieder zusätzliche Details zum Modell:
summary(fit1)
Call:
glm(formula = Y ~ X1 + X2, family = binomial(link = "logit"),
data = df_logit)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -35.4748 7.0709 -5.017 5.25e-07 ***
X1 2.9469 0.5929 4.970 6.68e-07 ***
X2 1.0481 0.2224 4.714 2.43e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 207.917 on 149 degrees of freedom
Residual deviance: 61.349 on 147 degrees of freedom
AIC: 67.349
Number of Fisher Scoring iterations: 7Die Teststatistik zur Überprüfung der Nullhypothese, ob sich ein Koeffizient signifikant von 0 unterscheidet heißt hier \(z\)-Statistik (z value) und der dazugehörige p-Wert entsprechend Pr(>|z|).
Für die Interpretation ist wichtig: Die Koeffizienten wirken linear auf die Log-Odds, also auf \[ \log\left(\frac{p(\mathbf{X})}{1-p(\mathbf{X})}\right). \] Wird ein Koeffizient exponentiert, erhält man das zugehörige Odds Ratio. Der Einfluss auf die Wahrscheinlichkeit selbst ist dagegen nicht konstant, sondern hängt vom Ausgangsniveau der Wahrscheinlichkeit ab.
exp(coef(fit1)) (Intercept) X1 X2
3.922037e-16 1.904762e+01 2.852239e+00 Prädiktiver Workflow in R
Für die Prognose verwenden wir wie schon bei lm die Funktion predict. Bei der logistischen Regression sind dabei zwei Ausgaben besonders wichtig:
type = "link"liefert den linearen Prädiktor, also den geschätzten Logit.type = "response"liefert die dazugehörige vorhergesagte Wahrscheinlichkeit.
new <- data.frame(X1 = c(3, 5, 7), X2 = c(18, 22, 25))
workflow_new <- new
workflow_new$logit_hat <- predict(fit1, newdata = new, type = "link")
workflow_new$p_hat <- predict(fit1, newdata = new, type = "response")
workflow_new$Y_hat <- ifelse(workflow_new$p_hat > 0.5, 1, 0)
workflow_new$logit_hat <- round(workflow_new$logit_hat, 3)
workflow_new$p_hat <- round(workflow_new$p_hat, 3)
workflow_new X1 X2 logit_hat p_hat Y_hat
1 3 18 -7.768 0.00 0
2 5 22 2.318 0.91 1
3 7 25 11.356 1.00 1Damit ist die typische Arbeit mit einem logistischen Modell in R im Kern dieselbe wie bei der linearen Regression: Modell fitten, Vorhersagen erzeugen, die Vorhersagen an Daten anhängen und anschließend interpretieren. Der häufig verwendete Schwellenwert \(0.5\) ist dabei nur eine Konvention und nicht naturgegeben.
df_logit_pred <- df_logit
df_logit_pred$logit_hat <- predict(fit1, type = "link")
df_logit_pred$p_hat <- predict(fit1, type = "response")
df_logit_pred$Y_hat <- ifelse(df_logit_pred$p_hat > 0.5, 1, 0)
workflow_view <- head(df_logit_pred[c("Y", "logit_hat", "p_hat", "Y_hat")], 8)
workflow_view$logit_hat <- round(workflow_view$logit_hat, 3)
workflow_view$p_hat <- round(workflow_view$p_hat, 3)
workflow_view Y logit_hat p_hat Y_hat
1 1 1.149 0.759 1
2 0 -1.028 0.264 0
3 1 -0.849 0.300 0
4 0 -0.457 0.388 0
5 0 -0.164 0.459 0
6 0 -1.279 0.218 0
7 0 -1.738 0.150 0
8 0 -0.814 0.307 0Decision Boundary
Die Funktion \(p(\mathbf{X})\) ist hier nun gerade eine Funktion von der Ebene \(R^2\) in das Intervall \([0,1]\), da der Vektor \(X\) die zwei Variablen \(X_1\) und \(X_2\) umfasst. Durch die Bedingung \(p(\mathbf{X})>0.5\) wird nun der Bereich in der Ebene beschrieben, der zu Merkmalsausprägungen von \(X_1\) und \(X_2\) korrespondiert, für die der Wert \(1\) für das Merkmal \(Y\) vorhergesagt wird. Analog ist durch die Bedingung \(p(\mathbf{X})\leq 0.5\) der Bereich charakterisiert, in dem \(Y=0\) prognostiziert wird. Diese beiden Bereiche werden durch eine Kurve bzw. eine Höhenlinie voneinander getrennt, die durch die Bedingung \(p(\mathbf{X})=0.5\) bestimmt ist. Diese Kurve wird als Decision Boundary bezeichnet. Zur weiteren Veranschaulichung wollen wir nun einmal die Decision Boundary in das obige Diagramm einzeichnen.
Hierzu verwenden wir die Funktion contour. Für Details zu dieser Funktion schauen Sie sich die Hilfe an. Wir können die Decision Boundary mit der Funktion contour in unser Diagramm einzeichnen (vgl. Abbildung 7.10).
plot(df_logit$X1, df_logit$X2,
xlab = expression(X[1]), ylab = expression(X[2]),
col = ifelse(df_logit$Y==1, "red", "blue"), pch = 19)
x1 <- seq(0,10,0.1)
x2 <- seq(14,26,0.1)
fun_pred <- function(x1, x2) predict(fit1,
newdata = data.frame(X1 = x1, X2 = x2),
type = "response")
py <- outer(X = x1, Y = x2, FUN = "fun_pred")
contour(x1, x2, py, levels = 0.5, add = TRUE)
Diese Decision Boundary zeigt uns nun an, wie wir die Entscheidung zwischen rot und blau treffen können, indem wir prüfen, ob eine Beobachtung rechts oder links dieser Linie liegt.
Fehlerrate (Error Rate)
Ein nahe liegender und verbreiteter Ansatz zur Beurteilung der Präzision, mit der eine Klassifizierungsmethode das interessierende qualitative Merkmal vorhersagen kann, ist die Betrachtung der Error Rate. Die Fehlerrate ist definiert als der Anteil aller falsch klassifizierten Untersuchungsobjekte. Für unseren simulierten Datensatz ergibt sich:
c(`Anteil richtiger Vorhersagen` = mean(df_logit_pred$Y_hat == df_logit_pred$Y),
`Fehlerrate auf Trainingsdaten` = mean(df_logit_pred$Y_hat != df_logit_pred$Y)) Anteil richtiger Vorhersagen Fehlerrate auf Trainingsdaten
0.91333333 0.08666667 Diese Fehlerrate wird hier auf denselben Daten berechnet, auf denen das Modell geschätzt wurde. Sie ist daher nur eine grobe deskriptive Kennzahl und keine faire Bewertung der Generalisierungsfähigkeit des Modells.
Nicht lineare Erweiterung
Die Decision Boundary im obigen Beispiel ist eine Gerade. In anderen Situationen kann ein solches Modell jedoch zu unflexibel sein. Betrachten wir hierzu einen Datensatz, bei dem die Trennung der Klassen nicht mehr linear erscheint (vgl. Abbildung 7.11).
Wir können die Flexibilität des Modells erhöhen, indem wir transformierte Terme direkt in der Formel ergänzen. Die Formel enthält dann neben \(X_1\) und \(X_2\) beispielsweise auch \(X_1^2\), \(X_2^2\) und das Produkt \(X_1X_2\). Diese Terme müssen dafür nicht dauerhaft als neue Spalten im Data Frame gespeichert werden.
fit2 <- glm(Y2 ~ X1 + X2 + I(X1^2) + I(X2^2) + I(X1 * X2),
family = binomial(link = "logit"),
data = df_logit)Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredprob2 <- predict(fit2, type = "response")
pred2 <- ifelse(prob2 > 0.5, 1, 0)
c(`Fehlerrate auf Trainingsdaten` = mean(pred2 != df_logit$Y2))Fehlerrate auf Trainingsdaten
0.01333333 
Abbildung Abbildung 7.12 zeigt, dass diese erweiterte Decision Boundary die Daten wesentlich besser beschreibt. Die Idee ist dieselbe wie bei der linearen Regression mit transformierten Kovariablen: Nicht die Methode wird geändert, sondern die Formel wird um geeignete Terme ergänzt.
Anwendung zur Churn-Analyse
Eine typische Anwendung für Klassifizierungsmethoden ist die Abwanderungsanalyse. Hierbei möchte man auf Grundlage diverser Merkmale eines Kunden prognostizieren, ob dieser Kunde voraussichtlich in Kürze die Kundenbeziehung mit dem Unternehmen beenden wird. Wir werden nun eine Abwanderungsanalyse mithilfe der logistischen Regression durchführen. Hierzu betrachten wir einen Beispiel-Datensatz, der Kundendaten eines Telekommunikations-Unternehmens für einen bestimmten Monat beinhaltet. Dieser liegt im Vorlesungsverzeichnis in der Datei telecom.RData. Der darin enthaltene Data Frame heißt telecom. Tabelle 7.5 zeigt einen Auszug aus diesem Datensatz, der insgesamt 3.000 Beobachtungen enthält.
str(telecom)'data.frame': 3000 obs. of 7 variables:
$ laufzeit : Factor w/ 2 levels "1 Jahr","1 Monat": 2 2 2 2 2 2 2 2 2 2 ...
$ support : Factor w/ 2 levels "Ja","Nein": 2 2 1 2 2 1 1 2 1 1 ...
$ internet : Factor w/ 3 levels "DSL","Glasfaser",..: 2 2 2 2 2 2 1 2 1 2 ...
$ telefon : Factor w/ 2 levels "Ja","Nein": 1 1 1 1 1 1 1 1 1 2 ...
$ ue50 : Factor w/ 2 levels "Ja","Nein": 1 2 2 2 1 2 2 2 2 1 ...
$ vertragsmonate: num 18 3 2 31 2 5 8 16 10 39 ...
$ abwanderung : Factor w/ 2 levels "Ja","Nein": 1 1 1 2 1 2 2 2 1 2 ...summary(telecom) laufzeit support internet telefon ue50
1 Jahr : 820 Ja : 883 DSL :1264 Ja :2663 Ja : 725
1 Monat:2180 Nein:2117 Glasfaser :1635 Nein: 337 Nein:2275
Kein Internet: 101
vertragsmonate abwanderung
Min. : 1.00 Ja : 641
1st Qu.: 6.00 Nein:2359
Median :23.00
Mean :26.52
3rd Qu.:43.00
Max. :70.00 telecom Datensatzes.
| laufzeit | support | internet | telefon | ue50 | vertragsmonate | abwanderung |
|---|---|---|---|---|---|---|
| 1 Monat | Nein | Glasfaser | Ja | Ja | 18 | Ja |
| 1 Monat | Nein | Glasfaser | Ja | Nein | 3 | Ja |
| 1 Monat | Ja | Glasfaser | Ja | Nein | 2 | Ja |
| 1 Monat | Nein | Glasfaser | Ja | Nein | 31 | Nein |
| 1 Monat | Nein | Glasfaser | Ja | Ja | 2 | Ja |
| 1 Monat | Ja | Glasfaser | Ja | Nein | 5 | Nein |
| 1 Monat | Ja | DSL | Ja | Nein | 8 | Nein |
| 1 Monat | Nein | Glasfaser | Ja | Nein | 16 | Nein |
| 1 Monat | Ja | DSL | Ja | Nein | 10 | Ja |
| 1 Monat | Ja | Glasfaser | Nein | Ja | 39 | Nein |
Das Merkmal abwanderung gibt nun für jeden Kunden an, ob dieser im betrachteten Monat seinen Vertrag gekündigt hat oder nicht. Außerdem gibt es Informationen zur Vertragslaufzeit (laufzeit), ob der Kundensupport genutzt wird (support), die Art des Internetanschlusses (internet), ob ein Telefonanschluss besteht (telefon), ob der Kunde alter als 50 Jahre ist (ue50) und die Anzahl der Monate des aktuellen Vertrags (vertragsmonate). Mit Hilfe einer logistischen Regression soll nun ein Modell angepasst werden, das zur Prognose der möglichen Abwanderung von Kunden verwenden werden kann.
Für die Modellierung codieren wir die Zielvariable zunächst so um, dass type = "response" direkt die Wahrscheinlichkeit für eine Kündigung liefert. Die kategorialen Kovariablen werden von glm() dabei analog zu lm() automatisch über Dummy-Variablen verarbeitet. Zusätzlich berücksichtigen wir einen Wechselwirkungseffekt zwischen Internetanschluss und Alter.
telecom_glm <- telecom
telecom_glm$abwanderung_ja <- factor(telecom_glm$abwanderung,
levels = c("Nein", "Ja"))
fit <- glm(
abwanderung_ja ~ laufzeit + support + internet + telefon + ue50 +
vertragsmonate + internet:ue50,
data = telecom_glm,
family = binomial(link = "logit")
)
summary(fit)
Call:
glm(formula = abwanderung_ja ~ laufzeit + support + internet +
telefon + ue50 + vertragsmonate + internet:ue50, family = binomial(link = "logit"),
data = telecom_glm)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.822617 0.435415 -8.779 < 2e-16 ***
laufzeit1 Monat 1.268706 0.218797 5.799 6.69e-09 ***
supportNein 0.961142 0.168737 5.696 1.23e-08 ***
internetGlasfaser 4.476357 0.385447 11.613 < 2e-16 ***
internetKein Internet 1.580427 0.691339 2.286 0.0223 *
telefonNein 0.917686 0.228164 4.022 5.77e-05 ***
ue50Nein -0.274118 0.381835 -0.718 0.4728
vertragsmonate -0.111012 0.005761 -19.271 < 2e-16 ***
internetGlasfaser:ue50Nein -0.930951 0.417012 -2.232 0.0256 *
internetKein Internet:ue50Nein -14.888448 421.775786 -0.035 0.9718
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 3112.6 on 2999 degrees of freedom
Residual deviance: 1531.5 on 2990 degrees of freedom
AIC: 1551.5
Number of Fisher Scoring iterations: 16Die Koeffizienten werden wie oben auf der Log-Odds-Skala interpretiert. Mit predict(..., type = "response") erhalten wir nun direkt die geschätzte Wahrscheinlichkeit für abwanderung = "Ja".
telecom_glm$p_abwanderung <- predict(fit, type = "response")
telecom_glm$prognose <- factor(ifelse(telecom_glm$p_abwanderung > 0.5,
"Ja", "Nein"),
levels = levels(telecom_glm$abwanderung))
mean(telecom_glm$abwanderung != telecom_glm$prognose)[1] 0.05533333Auch hier ist diese Fehlerrate nur eine grobe Kennzahl auf den Trainingsdaten. Für eine echte Bewertung der Prognosegüte würden wir später separate Evaluationsverfahren verwenden.
telecom_view <- head(telecom_glm[c("abwanderung",
"p_abwanderung",
"prognose")], 10)
telecom_view$p_abwanderung <- round(telecom_view$p_abwanderung, 3)
knitr::kable(telecom_view,
booktabs = TRUE,
caption = "Tatsächliche Abwanderung, vorhergesagte Wahrscheinlichkeit und Prognose für die ersten 10 Kunden.")| abwanderung | p_abwanderung | prognose |
|---|---|---|
| Ja | 0.708 | Ja |
| Ja | 0.793 | Ja |
| Ja | 0.621 | Ja |
| Nein | 0.146 | Nein |
| Ja | 0.935 | Ja |
| Nein | 0.540 | Ja |
| Nein | 0.024 | Nein |
| Nein | 0.476 | Nein |
| Ja | 0.019 | Nein |
| Nein | 0.184 | Nein |
7.3.2 Weitere GLMs
Während man in linearen Modellen annimmt, dass die Zielgröße normalverteilt ist, kann sie in GLMs auch anderen Verteilungen angehören, etwa der Binomial- oder Poisson-Verteilung. Praktisch werden GLMs in R über drei Bausteine festgelegt:
- Formel: Sie beschreibt den linearen Prädiktor, also welche Kovariablen und welche Interaktionen oder Transformationen verwendet werden.
- Verteilungsfamilie: Sie legt fest, welche Verteilung für die Zielgröße angenommen wird.
- Linkfunktion: Sie beschreibt, wie der Erwartungswert der Zielgröße mit dem linearen Prädiktor verknüpft wird.
Bei der logistischen Regression wird dafür family = binomial(link = "logit") verwendet. Eine Alternative wäre etwa der Probit-Link:
glm(y ~ x1 + x2, family = binomial(link = "probit"), data = daten)Ist die Zielvariable keine binäre Größe, sondern eine Ereignisanzahl, verwendet man häufig eine Poisson-Regression. In R wird diese typischerweise mit logarithmischer Linkfunktion angepasst:
glm(y ~ x1 + x2, family = poisson(link = "log"), data = daten)Das Grundprinzip bleibt also gleich: Die Formel beschreibt den systematischen Teil des Modells, family beschreibt Verteilung und Linkfunktion, und predict() liefert anschließend die gewünschten Vorhersagen.