Ordinal Logit in R
Eine ordinale logistische Regression kann in R mit der Funktion polr() aus der Bibliothek MASS gerechnet werden. Doch bevor ein Modell gerechnet werden kann, müssen die Daten vorbereitet werden. Zuerst wird der Datensatz eingelesen und anschliessend die abhängige Variable recodiert. Der London Data Datensatz kann unter Data heruntergeladen werden.
data = read.csv("london_survey_2010.csv")
data$local_influence =
car::recode(data$Q54,"'Don\\'t know'=NA")
Anschliessend wird die abhängige Variable als ordinal skaliert definiert. Die Variable misst die Einstellung, ob man lokale Entscheide beeinflussen kann.
data$local_influence = ordered(data$local_influence,
levels = c('Definitely disagree','Tend to disagree','Tend to agree','Definitely agree'))
Nun sind alle Vorbereitungen abgeschlossen und das Modell kann gerechnet werden. Als unabhängige Variable wird die Variable Q59E genommen. Sie misst, ob man tagtäglich Angestellte leitet. Als Kontrollvariable wird das Geschlecht genommen (Q1). Es wird folgende Hypothese überprüft: Wer einem unterstellte Mitarbeiter hat, hat eher das Gefühl, lokale Entscheidungen beeinflussen zu können. Dies, weil die Person bereits bei der Arbeit mehr Verantwortung hat und mitentscheiden kann.
library(MASS)
model1 = polr(local_influence ~ Q59E + Q1, data=data, Hess=TRUE)
summary(model1)
Call:
polr(formula = local_influence ~ Q59E + Q1, data = data, Hess = TRUE)
Coefficients:
Value Std. Error t value
Q59EYes 0.3226 0.1219 2.647
Q1Male 0.1348 0.1011 1.334
Intercepts:
Value Std. Error t value
Definitely disagree|Tend to disagree -1.1703 0.0838 -13.9623
Tend to disagree|Tend to agree -0.1332 0.0767 -1.7373
Tend to agree|Definitely agree 2.9061 0.1276 22.7686
Residual Deviance: 3263.704
AIC: 3273.704
(97 observations deleted due to missingness)
Im Output sieht man, dass Q59E signifikant ist (|t|>1.96) und Q1 nicht (|t|<1.96). Zudem ist der Effekt von Q59E auf local_influance positiv. Die Nullhypothese kann verworfen werden. Wer Angestellte hat, hat eher das Gefühl, lokale Entscheidungen beeinflussen zu können.
Effektstärke
In einem nächsten Schritt wird die Stärke des Effekts genauer angeschaut. Dazu werden die vorausgesagten Wahrscheinlichkeiten und die diskreten Änderungen berechnet. Beides kann mit der Bibliothek glm.predict berechnet werden.
library(glm.predict)
output = predicts(model1,"F;F(2)", position = 1)
output
val1_mean val1_lower val1_upper val2_mean val2_lower val2_upper dc_mean dc_lower dc_upper Q59E_val1 Q59E_val2 Q1
1 0.21331163 0.18511333 0.24322306 0.1643911 0.13288511 0.2018225 0.04892057 0.015684712 0.082162051 No Yes Male
2 0.21976527 0.19708472 0.24212440 0.1918018 0.16630695 0.2185673 0.02796350 0.007989933 0.049680668 No Yes Male
3 0.50713093 0.47395550 0.54054475 0.5627642 0.51935344 0.6022262 -0.05563325 -0.093022030 -0.017453762 No Yes Male
4 0.05979217 0.04696452 0.07448802 0.0810430 0.06155168 0.1080690 -0.02125083 -0.039681574 -0.005920574 No Yes Maleee
Bei den negativen Einstellungen ist die Wahrscheinlichkeit für Personen mit Angestellten tiefer als bei Arbeitern ohne leitende Funktion. Bei den positiven Einstellungen ist es genau umgekehrt. Die Unterschiede sind überall signifikant, da im Konfidenzintervall von dc_mean die 0 nicht enthalten ist.
In einem nächsten Schritt kann das Ganze noch grafisch dargestellt werden. Dazu wird ggplot verwendet.
library(ggplot2)
ggplot(output, aes(x=level,y=dc_mean,ymin=dc_lower,ymax=dc_upper)) +
geom_point() + geom_errorbar() +
geom_hline(yintercept=0,color="red",linetype="dashed") +
theme_bw() + ylab("discrete change of gender") +
xlab("Can you influence local decisions?")
Modellgüte
Als dritter Schritt wird die Modellgüte überprüft. Ist das Modell besser als ein Modell mit weniger Variablen? Dazu gibt es verschiedene Tests. Einer davon ist der Loglikelihood Ratio Test. Er kann mit der Funktion lrtest() aus der Bibliothek lmtest durchgeführt werden. Um den Test durchzuführen, muss zuerst ein Modell mit nur einem Teil der Variablen (restricted) aus dem ganzen Modell (unrestricted) gebildet werden. Anschliessend können die beiden Modelle verglichen werden.
library(lmtest)
model2 = polr(local_influence ~ Q1, data=data, Hess=TRUE)
lrtest(model1,model2)
Likelihood ratio test
Model 1: local_influence ~ Q59E + Q1
Model 2: local_influence ~ Q1
#Df LogLik Df Chisq Pr(>Chisq)
1 5 -1631.8
2 4 -1635.4 -1 7.0777 0.007805 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Der p-Wert ist mit 0.0078 kleiner als 0.05. Damit ist das Modell mit Q59E besser als das Modell ohne diese Variable. Oder anders ausgedrückt: Die Variable Q59E bietet einen Mehrwert.
model3 = polr(local_influence ~ Q59E, data=data, Hess=TRUE)
lrtest(model1,model3)
Likelihood ratio test
Model 1: local_influence ~ Q59E + Q1
Model 2: local_influence ~ Q59E
#Df LogLik Df Chisq Pr(>Chisq)
1 5 -1631.8
2 4 -1632.7 -1 1.7826 0.1818
Beim Vergleich mit oder ohne Q1 bekommt man einen p-Wert von 0.1818, was grösser ist als 0.05. Damit bringt die Variable Q1 keinem Mehrwert und kann weggelassen werden.
Die Resultate der Tests zeigen das gleiche wie die t-Werte. Bei diesen Tests ist es jedoch auch möglich mehrere Variablen zu vergleichen. So kann zum Beispiel getestet werden, ob das Modell besser ist als das Null-Modell.
model3 = polr(local_influence ~ Q59E, data=data, Hess=TRUE)
lrtest(model1,model3)
Likelihood ratio test
Model 1: local_influence ~ Q59E + Q1
Model 2: local_influence ~ Q59E
#Df LogLik Df Chisq Pr(>Chisq)
1 5 -1631.8
2 4 -1632.7 -1 1.7826 0.1818
Der p-Wert ist 0.01 und damit kann die Null-Hypothese verworfen werden. Das Modell ist signifikant besser als das Null-Modell. Eine weitere Möglichkeit sind BIC/AIC, welche bei der logistischen Regression erklärt werden.