Lektion 5 – lineare Regression

Nach der Behandlung von einfachen Hypothesentest in Lektion 4, wird in dieser Lektion erklärt, wie man eine lineare Regression rechnen kann. Um eine lineare Regression rechnen zu können, muss die abhängige Variable intervallskaliert sein.

Univariate Modelle

In folgendem Video wird erklärt, wie ein Modell mit einer unabhängigen Variable geschätzt werden kann, welche kontinuierlich oder diskrete sein kann. Im Anschluss folgt eine Zusammenfassung der wichtigsten Punkte.

Ein lineares Modell kann in R mit der Funktion lm() geschätzt werden. Der erste Parameter beinhaltet die Modell-Formel in der Form abhängige_Variable ~ unabhängige_Variable. Anschliessend folgt der Datensatz. Optional können noch gewichtige angegeben werden, welche bei repräsentativen Umfragen oft enthalten sind, um die Umfrage repräsentative zu machen. Ein Beispiel:

model = lm(lr_self ~ age, data = df_selects,
        weights = weight_total)
summary(model)

R liefert folgendes Ergebnis:

Die erste Zeile erinnert uns nochmals, welches Modell wir geschätzt haben. Anschliessend zeigt es uns die Verteilung der gewichten Residuen an anhand der Quartile. Residuen sind die Differenz zwischen dem geschätzten Wert und dem tatsächlichen Wert.

Unter Coefficients werden die geschätzten Werte angezeigt. Die Kolonne Estimate beinhaltet die Koeffizienten: im obigen Beispiel β0 (Grundwert, Achsenabschnitt) und β1 (Alter, die Steigung der Kurve). Alter kann folgendermassen interpretiert werden: Eine Person, welche ein Jahr älter ist, schätzt sich im Durchschnitt um 0.01 Einheiten rechter ein auf der Links-Rechts-Achse. β0 alleine zu interpretieren macht wenig Sinn. Es wäre der Wert eines frischgeborenen Babys. Ob dieses schon eine politische Meinung hat? Wohl eher nicht. β0 ist aber notwendig, um den vorausgesagte Wert für eine Person zu berechnen. Eine 20-Jährige hätte beispielsweise einen vorausgesagten Wert bei der Links-Rechts-Selbsteinschätzung von 5.13 (ŷ = β0 + β1age = 4.914 + 0.011 * 20). Hat eine konkrete 20-Jährige im Datensatz hingegen einen tatsächlichen Wert von 4 beträgt das Residuum für diesen Fall 1.13. Die Summe aller quadrierten Residuen (grössere Werte werden also stärker bestraft) ist mit der Methode OLS am kleinsten (OLS wird von lm() verwendet). OLS steht Ordinary Least Squares, also die Methode der kleinsten Quadrate. Wenn die Annahmen zutreffen, liefert diese Methode immer das optimale Ergebnis.

Die Kolonne Std. Error gibt die Standardfehler zurück. Diese sind wichtig für die Signifikanz. Der t-Wert ist Beta durch den Standardfehler dividiert. Ist sein absolute Wert 1.96 oder grösser, ist der Koeffizient auf einem α-Wert von 0.05 signifikant. Eine andere Möglichkeit ist der p-Wert [«Pr(>|t|)»]. Ist dieser kleiner als 0.05 (oder was immer man als α-Wert festgelegt hat, bei kleiner Fallzahl wird manchmal 0.1 genommen, es gibt auch Forscher, welcher immer 0.005 nehmen), so ist der Koeffizient signifikant. Der p-Wert ist die Wahrscheinlichkeit, dass wir die Nullhypothese fälschlicherweise verwerfen, obwohl sie in Wirklichkeit wahr ist. Die Nullhypothese sagt, dass die Koeffizient gleich 0 ist. Können wir sie verwerfen, wissen wir, dass sich der Wert signifikant von Null unterscheidet und somit ein Effekt vorhanden ist. Im obigen Beispiel ist das der Fall. Alter hat einen signifikanten Einfluss wo sich jemand auf der Links-Rechts-Achse positioniert.

Das Multiple R-Square sagt aus, welcher Anteil der totalen Varianz der abhängigen Variable von der unabhängigen Variable erklärt wird. In Beispielmodell beträgt der Wert 0.6% und ist damit sehr gering. Mit Alter alleine die Links-Rechts-Selbsteinschätzung von jemandem zu erklären. In Sozialwissenschaften gilt ein Wert von rund 26% als sehr gut.

Nun noch ein Beispiel einer diskreten Variable (=factor in R).

model = lm(lr_self ~ opinion_eu_membership, 
           data = df_selects, weights = weight_total)
summary(model)

R liefert folgendes Ergebnis:

Anders als beim Alter bekommen wir nun mehrere βs zurück. Die Variable opinion_eu_membership hat fünf Ausprägungen. Im Modell fehlt jedoch die erste, sie ist im Achsenabschnitt [«(Intercept)»] abgebildet. Alle Werte von opinion_eu_membership sind eine Veränderung des Achsenabschnitts und ein Vergleich mit der Grundwert. Im Beispiel ist der Grundwert eine Person, welche einen Beitritt der Schweiz zur EU wünscht. Sie hat einen vorausgesagten Wert von 3.71 bei der Links-Rechts-Selbsteinschätzung. Eine Person, welche eher für den EU-Beitritt ist unterscheidet sich nicht signifikant von einer Person die klar für den EU Beitritt ist. Eine unentschlossene Person hingegen ist 0.77 Einheiten rechter und damit einen vorausgesagten Wert von 4.48. Der Unterschied ist signifikant, was auch für die anderen beiden Werte gilt.

Alle diese Werte sind ein Vergleich zur Basiskategorie. Wollen wir hingegen die Werte untereinander vergleichen, müssen wir einen eigenen Hypothesentest durchführen. Dies können wir mit der Funktion linearHypothesis() aus der Bibliothek car tun. Das folgende Beispiel vergleicht Personen, welche sicher draussen bleiben wollen mit Personen, die eher draussen bleiben sollen. Dafür testen wir die Differenz der beiden Koeffizienten auf Null.

library(car)
linearHypothesis(model, 
   "opinion_eu_membershipRather in favour to stay out - 
    opinion_eu_membershipStrongly in favour to stay out")

R gibt uns folgenden Output zurück:

Da der p-Wert verworfen werden kann (p [«Pr(>F)»] < 0.05), ist der Unterschied statistisch signifikant.

Multivariate Modelle

Das folgende Video behandelt lineare Modelle mit mehr als einer unabhängigen Variable. Unterhalb des Videos werden die wichtigsten Punkte zusammengefasst.

Um mehrere unabhängige Variablen dem Modell hinzuzufügen, kann das + -Zeichen verwendet werden. Im folgenden Beispiel verwenden wir die beiden Variablen von oben plus das Geschlecht.

model = lm(lr_self ~ age + gender + opinion_eu_membership,
           data = df_selects, weights = weight_total)
summary(model)

R liefert folgenden Output zurück:

Das Modell kann folgendermassen interpretiert werden: Ein Person, welche ein Jahr älter ist, schätzt sich im Durchschnitt um 0.019 Einheiten rechter ein, ceteris paribus. Der Effekt ist statistisch signifikant. Ceteris paribus bedeutet, dass alle anderen Variablen konstant gehalten werden. Die 0.019 Einheiten treffen zum Beispiel zu, wenn wir einen 42-jährigen Mann, welcher eher gegen einen EU-Beitritt ist, mit einem 41-jährigen Mann vergleichen, der ebenfalls eher gegen einen EU-Beitritt ist. Es trifft aber nicht zu, wenn wir den 42-jährigen Mann mit einer 41-jährigen Frau vergleichen, die eher für einen EU-Beitritt ist. Ceteris paribus müssen wir immer bei unserer Interpretation hinzufügen, wenn wir einen Koeffizienten aus einen multivariaten Modell erklären wollen. Der Modelloutput sagt uns des weiteren, dass Frauen sich im Durchschnitt um 0.55 Einheiten linker einschätzen als Männer und das EU-Befürworter linker sind als solche, welche einen EU-Beitritt der Schweiz ablehnen. Die Effekte sind statistisch signifikant. Wenn wir bei EU-Beitritt Kategorien vergleichen wollen, können wir das gleich tun wie beim univariaten Modell.

Das Multiple R-squared ist nun 17.9%, was ein guter Wert ist. Mit diesen Variablen die links-rechts-Positionierung abzuschätzen, kann durchaus Sinn machen. Beim Multiple R-squared gibt es jedoch ein Problem. Wenn wir weitere Variablen zum Modell hinzufügen, wird es grösser und grösser. Egal wie unsinnig die Variablen sind, es wird nie kleiner. Um zwischen mehreren Modellen das beste auszuwählen, gibt es das Adjusted R-squared. Dieses kann auch kleiner werden, da es die Anzahl Variablen mit einbezieht. Fügen wir ein gute Variable hinzu, wachsen bei R-squared, fügen wir hingegen eine unbrauchbare unabhängige Variable hinzu, wird das multiple R-squared etwa gleich bleiben (oder minim grösser), während das Adjusted R-squared kleiner wird.

Dann gibt es noch die F-Statistik am Ende des Outputs. Sie sagt uns, ob das ganze Modell statistisch signifikant ist oder nicht. Das ist hier der Fall, da der p-Wert unter 0.05 liegt.

Wenn man theoretisch davon ausgeht, dass eine Variable einen Effekt auf den Effekt einer anderen Variable auf die abhängige Variable hat, kann eine Interaktion geschätzt werden. Eine Interaktion zwischen zwei Variablen wird mit dem *-Zeichen verwendet werden. Im folgenden Beispiel testen wir die Hypothese, dass der Effekt von Alter auf die Links-Rechts-Selbsteinschätzung bei Männern stärker ist als bei Frauen. Als Kontrollvariable nehmen wir die Meinung zu einem EU-Beitritt. Dazu schätzen wir eine Interaktion.

model = lm(lr_self ~ age * gender + opinion_eu_membership,
           data = df_selects, weights = weight_total)
summary(model)

R liefert uns folgenden Output zurück:

Der Output gleicht dem obigen Modell ohne die Interaktion. Hinzugekommen ist der Koeffizient age:genderfemale am Ende der Tabelle. Dieser ist nicht signifikant, wir müssen die Hypothese also verwerfen. Alter hat den gleichen Effekt für Männer und Frauen. Das :-Zeichen kann als Multiplikationszeichen gesehen werden. Der Effekt von Alter für Männer kann bei Koeffizienten age abgelesen werden (0.017), bei Frauen muss zusätzlich der Interaktionskoeffizient age:genderfemale hinzu addiert werden, was 0.02 ergibt. Wir können nur auch testen, ob der Geschlechtsunterschied noch signifikant ist für eine 110-jährige Person. Dazu können wir wiederum die Funktion linearHypothesis() aus car verwenden.

car::linearHypothesis(model, "genderfemale + 
                      110 * age:genderfemale")

Wer stellen fest, dass der Effekt nun nicht mehr signifikant ist. Warum hat der im Modell statistisch nicht signifikante Interaktionskoeffizient nun doch einen Effekt, indem er aber einem gewissen Alter (wenn auch unrealistisch hohen) den starken Effekt des Geschlechts neutralisiert? Die Signifikanz des Interaktionskoeffizienten ist immer für einen Faktor von 1, was in unserem Beispiel ein einjähriges Mädchen ist. Mit anderen Termen zusammen und einem anderen Faktor, können dann trotzdem Effekte ändern.

Die Ergebnisse eines Modells will man normalerweise in eine schöne Tabelle packen. Dazu gibt es zwei verschieden Bibliotheken: texreg und stargazer. Das folgende Beispiel zeigt, wie man die Tabelle in texreg in der Konsole ausgibt.

texreg::screenreg(model)

Annahmen überprüfen

OLS ist dann unverzerrt, wenn gewisse Annahmen erfüllt sind. Das folgende Video geht diese Annahmen durch und erklärt Schritt für Schritt wie diese überprüft werden können. Im Anschluss folgt wie immer eine Zusammenfassung.

Die lineare Regression geht von einer lineare Beziehung aus. Das können wir mit der Funktion crPlots() aus der Bibliothek car überprüfen.

library(car)
model = lm(lr_self ~ age, data = df_selects,
           weights = weight_total)
crPlots(model)

An der Grafik erkennen wir, dass die Annahme erfüllt ist, da die pinke Linie mehr oder weniger der blau gestrichenen Ideallinie folgt. Hat die pinke Linie eher die Form einer Parabel, so deutet es auf einen quadratischen Zusammenhang hin. Andere Formen sind beispielsweise logistisch (log(Variable)), was beispielsweise beim BIP gemacht wird.

Eine weitere Annahme ist Homoskedastizität. Es besagt, dass die Varianz über das Modell gleichmässig verteilt ist. Ist dies nicht der Fall, weil beispielsweise die Varianz von älteren Person viel höher ist als die von jungen, spricht man von Heteroskedastizität. In diesem Fall ist BLUE verletzt. Überprüft werden kann das ganze mit dem Breusch-Pagan Tast, welche in R der Funktion bptest() aus der Bibliothek lmtest implementiert ist.

library(lmtest)
bptest(model)

Die Nullhypothese besagt, dass Homoskedastizität vorherrscht. Da der p-Wert unter 0.05 liegt (0.006), müssen wir die Nullhypothese verwerfen. Wir haben also Heteroskedastizität und BLUE ist verletzt und dadurch die Standardfehler verzerrt. Nun gibt es glücklicherweise eine Möglichkeit die durch Heteroskedastizität verzerrten Standardfehler zu korrigieren: der Huber-Sandwich-Schätzer. Die Bibliothek sandwich enthält eine Funktion vcovHC() um die Varianz-Kovarianz Matrix zu korrigieren. Aus der Varianz-Kovarianz Matrix können wir die Standard Fehler aus den Diagonalelementen (Varianzen) berechnen. Aus diesen dann die t- und p-Werte.

library(sandwich)
vcov = vcovHC(model)
sd = sqrt(diag(vcov))
t = coef(model) / sd
p = 2 * pnorm(-abs(t))

Als nächstes Prüfen wir die Annahme, on die Residuen normalverteilt sind. Das können wir mit der Funktion qqPlot() aus car tun.

library(car)
qqPlot(model)

Die Punkte sollten mehr oder weniger zwischen den gestrichelten Linien sein, was ausser bei einigen Outliers der Fall ist. Eine anderen Möglichkeit ist die studentisierten Residuen zu plotten. Die studentisierten Residuen erhält man, in dem man für jeden Fall ein neues Modell rechnet ohne diesen Fall und dann das Residuum dieses Fall zu diesem Model berechnet. Berechnet werden können diese mit der Funktion studres() aus der Bibliothek MASS. Anschliessend können wir sie mit der Normalverteilung graphisch vergleichen.

library(MASS)
stud_res = studres(model)

library(tidyverse)
ggplot() + geom_density(aes(x = stud_res), color = "red") +
  geom_line(aes(x = seq(min(stud_res), max(stud_res), length.out = 50),
                y = dnorm(seq(min(stud_res), max(stud_res), length.out = 50)))) +
  theme_minimal()

Auch hier ist ersichtlich, dass diese einigermassen normalverteilt sind.

Als letzte Annahme sollte die Multikollinearität überprüft werden. Perfekte Multikollinearität erkennt R automatisch und ersetzt einen Koeffizienten durch NA. Überprüfen können wir das ganze mit der Funktion imcdiag() aus mctest. Ein Beispiel, wo R einen Output liefert trotz Multikollinearität: Alter und zweimal Alter plus einen kleiner Fehler.

df_selects$age2 = 2 * df_selects$age + rnorm(5337)
model = lm(lr_self ~ age + age2,
           data = df_selects, weights = weight_total)
imcdiag(model)

An der Zeile «Klein» sehen wir, das age und age2 zu stark korrelieren. In diesem Fall sollte man eine der Variablen aus dem Model entfernen. Sie bietet keinen Mehrwert und verzerrt das Model (Alter ist in diesem Model nicht signifikant, alleine aber schon.)