Politikwissenschaft

 

Webseite durchsuchen

Logistische Regression in R

Benjamin Schlegel | 18. April 2016

PDF

Eine logistische Regression kann in R mit der Funktion glm() gerechnet werden. Wichtig dabei ist, dass als Familie binomial angegeben wird. Doch vor dem rechnen einen Regression muss zuerst der Datensatz eingelesen und rekodiert werden. Der Artikel setzt die Artikel logistische Regression und R Grundlagen voraus.

Es wird der Datensatz London Survey gelesen, welcher unter Data heruntergeladen werden kann. Anschliessend wird die Variable illness rekodiert, welche angibt, ob jemand eine schwere Krankheit oder sonst ein Gebrechen hat. Die Variable livingLong wird aus der Variable Q7 gewonnen und gibt an, ob jemand schon länger als zehn Jahr in London lebt. Die Variable neighbourhood gibt an, ob sich eine Person in der Nacht in ihrer Umgebung draussen sicher fühlt und kann vier Ausprägungen annehmen.


data = read.csv(file.choose(),stringsAsFactors = FALSE) # london_survey_2010.csv
library(car)
data$illness = factor(recode(data$Q61,"'No'=0;'Yes'=1;else=NA"))
data$livingLong = factor(recode(data$Q7, "c('All my life/born in London',
  'Over 10 years up to 15 years','Over 20 years or more',
  'Over 15 years up to 20 years')=1; c('Up to 1 year','Over 1 year up to 2 years',
  'Over 2 years up to 5 years', 'Over 5 years up to 10 years')=0;else=NA"))
data$neighbourhood = factor(recode(data$Q13, "'1 - Very unsafe'='very unsave';
  '2 - A bit unsafe'='a bit unsafe';'3 - Fairly safe'='fairly safe';
  '4 - Very safe'='very safe';else=NA"))

Nun wollen wir nur diese drei Variablen aus dem Datensatz und verwenden dazu die Funktion select() aus der Bibliothek dplyr. Mit dem Befehl na.omit() schmeissen wir alle Fälle raus, welche nicht vollständig sind.

.

library(dplyr)
data = select(data,illness,livingLong,neighbourhood)
data = na.omit(data)

Regression

Nun sind die Daten vorbereiten. Wir wollen den Einfluss von schweren Krankheiten und die Nachbarschaft auf die Zeit in London überprüfen. Dazu stellen wir folgendes Modell auf:


model1 = glm(livingLong ~ illness + neighbourhood,data=data,
    family=binomial(link= "logit"))
summary(model1)

Call:
glm(formula = livingLong ~ illness + neighbourhood, 
    family = binomial(link = "logit"), data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.6717   0.2500   0.6893   0.7549   0.7549  

Coefficients:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)               1.22542    0.21108   5.805 6.42e-09 ***
illness1                  2.22425    0.59060   3.766 0.000166 ***
neighbourhoodfairly safe -0.11583    0.23087  -0.502 0.615875    
neighbourhoodvery safe    0.09064    0.24290   0.373 0.709029    
neighbourhoodvery unsave  0.68362    0.39875   1.714 0.086449 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1348.5  on 1308  degrees of freedom
Residual deviance: 1311.2  on 1304  degrees of freedom
AIC: 1321.2

Number of Fisher Scoring iterations: 6

Mit der Funktion glm() kann ein generalisiertes Modell gerechnet werden. Mit family=binomial wird angegeben, dass es sich um eine dichotome Variable handelt (Bei einer Poisson Regresson wäre es z.B. family=poisson.). Der Standardwert für den Link bei binomial ist logit, die Angabe wäre also nicht nötig gewesen. Will man ein Probit rechnet, muss dort probit angegeben werden. Unter data geben wir den Datensatz an mit dem gerechnet werden soll.

Anders als bei der linearen Regression können die Koeffizienten nicht direkt interpretiert werden. Die Koeffizienten geben uns einzig die Richtung des Zusammenhangs an. Zudem sieht man an den p-Werten, ob die Koeffizienten signifikant sind (Kleiner als 0.05 bedeutet in der Wissenschaft signifikant.). In diesem Beispiel sieht man, dass schwer kranke Menschen er lang in London leben als gesunde Menschen, der Einfluss ist signifikant. Anders das Sicherheitsgefühl am Abend im Quartier: Es hat keinen signifikanten Einfluss ob man lange in London lebt.

Interpretation

vorausgesagte Wahrscheinlichkeiten und diskrete Änderungen

Eine Möglichkeit der Interpretation von logistischen Modellen ist das berechnen von vorausgesagten Wahrscheinlichkeiten und diskreten Änderungen. Das geht ganz einfach mit dem Paket glm.predict (Unter dem Link befindet sich ein Installationsanleitung.).


library(glm.predict)
predicts(model1, "0,1;F", position = 1)

  val1_mean val1_lower val1_upper val2_mean val2_lower val2_upper    dc_mean   dc_lower    dc_upper illness_val1 illness_val2 neighbourhood
1 0.7711344  0.6935842  0.8417224 0.9632906  0.9010585  0.9901658 -0.1921562 -0.2658895 -0.11865659            0            1  a bit unsafe
2 0.7510950  0.7158122  0.7852164 0.9594590  0.8966024  0.9893704 -0.2083640 -0.2568494 -0.14459021            0            1   fairly safe
3 0.7884394  0.7492191  0.8271662 0.9671972  0.9174076  0.9907619 -0.1787578 -0.2263789 -0.12312399            0            1     very safe
4 0.8665730  0.7830854  0.9290180 0.9804139  0.9414450  0.9959582 -0.1138410 -0.1835975 -0.05657875            0            1   very unsave

Zuerst braucht die Funktion predicts() das Modell mit dem gerechnet wurde. Als zweites Argument einen Character mit den Werten, welche der Funktion übergeben werden sollen. Für illness wollen wir beide Werte (0 und 1) und für neighbourhood alle Werte des Faktors (F). Die diskrete Änderung wird für die erste Variable durchgeführt, weshalb es dort immer min. zwei Werte braucht (position = 1). Die Funktion gibt ein data.frame zurück. Dort sieht man, dass Personen, welche sich ein wenig unsicher fühlen und nicht krank sind eine Wahrscheinlichkeit von 77.2 Prozent haben lange in London zu wohnen (mean1). Wäre die gleiche Person jedoch schwer krank, hätte sie eine Wahrscheinlichkeit von 96.3 Prozent lange in London zu wohnen (mean2). Eine sich etwas unsicher fühlende gesunde Person hat eine um 19.1 Prozentpunkte geringere Wahrscheinlichkeit über zehn Jahre in London zu leben als eine sich unsicher fühlende schwer kranke Person (mean.diff). Der Unterschied ist signifikant (Die Bandbreite von lower.diff und upper.diff (95%-Konfidenzintervall) geht nicht durch 0 hindurch.). Bei den anderen Ausprägungen von neighbourhood ist es äquivalent.

Odds Ratio

Die obige Interpretation funktioniert bei allen generalisierten linearen Modellen. Bei logistischen Modellen kommt zusätzlich die Interpretationsmöglichkeit der Odds Ratio (Quotenverhältnis) hinzu.


exp(coef(model1))

            (Intercept)                 illness1 neighbourhoodfairly safe 
               3.4056128                9.2465892                0.8906291 
  neighbourhoodvery safe neighbourhoodvery unsave 
               1.0948743                1.9810421 

Eine schwer kranke Person hat eine 9.24-Fache Chance mehr als zehn Jahr in London zu leben als eine gesunde Person.

Geschachtelte Modelle vergleichen

Zuerst wird ein Modell gerechnet, welches nur die Variable illness enthält.


model2 = glm(livingLong ~ illness,data=data,family=binomial(link= "logit"))

Das Modell 1 ist das unbeschränkte (unrestricted) Modell und Modell zwei das beschränkte (restricted) Modell. Mit dem BIC (Bayessches Informationskriterium) und dem AIC (Akaikes Informationskriterium) können die Modelle verglichen werden. Ein kleinerer Wert bedeutet ein besseres Modell. Das BIC ist bei grösseren Stichproben zu bevorzugen. Man kann auch immer nur das BIC berechnen.


BIC(model1,model2)

       df      BIC
model1  5 1347.107
model2  2 1332.448

AIC(model1,model2)

       df      AIC
model1  5 1321.222
model2  2 1322.094

Beim BIC ist das Modell 2 ohne die Variable neighbourhood besser, beim AIC Modell 1. Da es sich um eine doch eher grosse Stichprobe von 1309 handelt, wir das BIC bevorzugt. Die Variable neighbourhood bringt keinen Mehrwert und kann weggelassen werden.

Eine weitere Möglichkeit Modelle zu vergleichen ist der Waldtest. Beim Waldtest wird nur das unbeschränkte Modell angegeben und zusätzlich die Variable die geprüft werden soll.


library(lmtest)
waldtest(model1,2,test="Chisq")

Wald test

Model 1: livingLong ~ illness + neighbourhood
Model 2: livingLong ~ illness
  Res.Df Df  Chisq Pr(>Chisq)
1   1304                     
2   1307 -3 6.1579     0.1042

Die Funktion waldtest() gibt uns einen p-Wert von 0.1 zurück. Dieser Wert ist grösser als 0.05 und deshalb bringt die zusätzliche Variable neighbourhood nichts. Die Information des BIC wird bestätigt. Es können auch mehr als eine Variable angegeben werden.

Grafik

Will man die vorausgesagten Wahrscheinlichkeiten, resp. diskreten Änderungen zeichnen, kann mit ggplot2 eine Grafik mit dem data.frame von glm.predic, resp. discrete.changes gezeichnet werden. Ein Beispiel mit einer kontinuierlichen Variable findet sich hier, eines mit einer diskreten Variable hier.