9  Logistische Regressie & Odds Ratios

Logistische regressiemodellen worden gebruikt voor binaire afhankelijke variabelen. We maken hier gebruik van de Nederlandse survey voor ronde 9 van de European Social Survey (ESS). Deze dataset is op de ESS website beschikbaar in SPSS format (.sav). We kunnen de view_df functie in sjplot gebruiken om de variabelen en hun labels in de dataset te inspecteren. Zie Paragraaf 1.1

#Packages
library(rio)             #laden van data
library(sjPlot)          #overzichten van data objecten
library(tidyverse)       #datamanipulatie en grafieken
library(broom)           #resultaten van regressiemodellen

#Data
ESS9NL <- import("ESS9e03, Netherlands.sav")

#view_df example op subset van de dataset
ESS9NL |> 
  select(polintr, ppltrst) |> 
  sjPlot::view_df()

Data frame: select(ESS9NL, polintr, ppltrst)

ID	Name	Label	Values	Value Labels
1	polintr	How interested in politics	1 2 3 4 7 8 9	Very interested Quite interested Hardly interested Not at all interested Refusal Don't know No answer
2	ppltrst	Most people can be trusted or you can't be too careful	0 1 2 3 4 5 6 7 8 9 10 77 88 99	You can't be too careful 1 2 3 4 5 6 7 8 9 Most people can be trusted Refusal Don't know No answer

Deview_df() output toont dat er numerieke waarden zijn die geassocieerd zijn met ontbrekende waarden (bv. respondenten die “Don’t Know” antwoorden krijgen een score van 8 op polintr). Deze waarden werden al naar ‘NA’ omgezet zoals je kunt zien in onderstaande tabel en dus is er geen verdere data managament nodig. Zie Section 4.2 van het Statistiek 1 boek om deze leerstof op te frissen.

table(ESS9NL$polintr)

  1   2   3   4 
213 832 478 149 

Waarschuwing!

Kijk altijd goed naar instructies van opdrachten alsook verdere informatie over de dataset (bv. het codeboek, of informatie over de dataset te verkrijgen met functies zoals view_df() of attributes()) om te weten welke datamanagement stappen je moet zetten vooraleer je je analyse kan uitvoeren.

9.1 Logistische regressieanalyse

9.1.1 Data Management

In dit voorbeeld zullen we eerst onderzoeken of gender (gndr) opkomst bij verkiezingen bepaalt (vote).

We kijken eerst of data management nodig is:

#Variabele attributen
ESS9NL |> 
  select(gndr, vote) |> 
  view_df()

Data frame: select(ESS9NL, gndr, vote)

ID	Name	Label	Values	Value Labels
1	gndr	Gender	1 2 9	Male Female No answer
2	vote	Voted last national election	1 2 3 7 8 9	Yes No Not eligible to vote Refusal Don't know No answer

#Tabellen
table(ESS9NL$gndr)

  1   2 
833 840 

table(ESS9NL$vote)

   1    2    3 
1291  247  130 

De onafhankelijke variabele heeft 2 categorieën en moet in een factor variabele worden omgezet. De afhankelijke variabele heeft 3 categorieën (Yes, No, Not Eligible). We maken hier een binaire factor variabele van door eerst de “Not Eligible” categorie op NA te zetten:

#Factor maken
ESS9NL <- ESS9NL |>
  mutate(gndr_F = factorize(gndr),
         vote_F = factorize(vote))

#not eligible op NA
ESS9NL <- ESS9NL |>
  mutate(vote_F = na_if(vote_F,"Not eligible to vote"))

1

De categorie met de laagste numerieke waarde wordt hierbij de referentiecategorie. We gebruiken factorize gezien datawaarden labels hebben.

2

We maken hier geen nieuwe variabele aan voor de gehercodeerde variabelen maar overschrijven de originele variabelen. Dit is doorgaans niet aangeraden voor studenten, gelukkig weten wij meestal wel waar we mee bezig zijn.

We checken de niveaus van de variabelen, in het bijzonder van de vote variabele:

levels(ESS9NL$vote_F)

[1] "Yes"                  "No"                   "Not eligible to vote"
[4] "Refusal"              "Don't know"           "No answer"           

table(ESS9NL$vote_F)

                 Yes                   No Not eligible to vote 
                1291                  247                    0 
             Refusal           Don't know            No answer 
                   0                    0                    0 

De vote variabele is nu een factor variabele met “Yes” als de referentiecategorie. Dit is niet wat we willen gezien we stemmen willen voorspellen. Als we de variabele zo laten voorspelt het model of een persoon niet heeft gestemd. We veranderen dit met de relevel functie (zie Paragraaf 2.1.1).¹

#Relevel de variabele
ESS9NL <- ESS9NL |> 
  mutate(vote_F = relevel(vote_F, "No"))

#en controleer het resultaat
levels(ESS9NL$vote_F)

[1] "No"                   "Yes"                  "Not eligible to vote"
[4] "Refusal"              "Don't know"           "No answer"           

mutate(vote_F = relevel(vote_F, "No"))

We gebruiken de relevel functie op de vote variabele. We creëren hier geen nieuwe variabele, maar overschrijven de oude. Je zou er ook voor kunnen kiezen een nieuwe variabele te maken. De categorie tussen dubbele aanhalingstekens zal de eerste categorie worden en dus de referentiecategorie in de regressie. We gebruiken het label “No” omdat de variabele reeds tot factor is getransformeerd (en dus niet de originele numerieke waarde ‘2’).

Laten we nu kijken naar de gndr_F variabele:

table(ESS9NL$gndr_F)

     Male    Female No answer 
      833       840         0 

levels(ESS9NL$gndr_F)

[1] "Male"      "Female"    "No answer"

Als we de niveaus bekijken zien we dat ‘Male’ als referentie zal worden genomen. Dit is prima, maar bij wijze van voorbeeld veranderen we dit hieronder naar ‘Female’. We zien ook een derde categorie ‘No answer’. Dit label werd gegeven aan de waarde in de variabele maar is leeg. R zal deze dus verwijderen in de analyses.

ESS9NL <- ESS9NL |> 
  mutate(gndr_F = relevel(gndr_F, "Female"))

#controleer je codering
levels(ESS9NL$gndr_F)

[1] "Female"    "Male"      "No answer"

Waarschuwing!

De afhankelijke variabele voor een logistische regressie is een (factor) binaire variabele. Zorg ervoor dat de referentiecategorie de uitkomst is die je niet wil voorspellen en de hoogste categorie net die is die je wil voorspellen. Anders zal je interpretatie foutief zijn.

9.1.2 Logistische regressie uitvoeren

Het uitvoeren van logistische regressie in R is gelijkaardig aan lineare regressie. In plaats van de ingebouwde functie lm(), gebruiken we de eveneens ingebouwde glm() functie. De afkorting staat voor ‘generalized linear model’.

#Schat het model
Vote_model <- glm(vote_F ~ gndr_F, 
                data = ESS9NL, family = "binomial")

Vote_model <-

We slaan de resultaten op in een data object met naam naar keuze.

glm(vote_F ~ gndr_F,

We voeren de glm functie uit met gefactoriseerde vote als afhankelijke variabele, voorspeld (~) door onze enige onafhankelijke variabele: gndr (gefactoriseerd). We kunnen meerdere onafhankelijke variabelen toevoegen, gescheiden van elkaar met een ‘+’ teken.

data = ESS9NL,

We geven aan welke dataset gebruikt wordt.

family = "binomial")

We verduidelijken de familie van modellen voor ons generalized linear model. Voor logistische regressie is dit “binomial”. Dit gedeelte van de code blijft onveranderd. Zie de Veelvoorkomende Fouten appendix ( Paragraaf A.3) voor een error die je kunt tegenkomen als je dit gedeelte vergeet.

De resultaten bekijken we weer met de summary() functie:

summary(Vote_model)

Call:
glm(formula = vote_F ~ gndr_F, family = "binomial", data = ESS9NL)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.56485    0.09535  16.412   <2e-16 ***
gndr_FMale   0.18359    0.13925   1.318    0.187    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1355.5  on 1537  degrees of freedom
Residual deviance: 1353.7  on 1536  degrees of freedom
  (135 observations deleted due to missingness)
AIC: 1357.7

Number of Fisher Scoring iterations: 4

Output uitleg

De output is gelijkaardig aan die van de lineaire regressie met lm().

• Call: Het model dat geschat werd
• Deviance Residuals: beschrijvende statistieken over de residuals van het model.
• Coefficients: De logistische regressiecoëfficiënten (Estimate), hun standaardfouten (Std. Error), en de teststatistiek (z-waarde; de Z-statistiek is gelijk aan \(\frac{\textrm{Coefficient}}{\textrm{Std. Error}}\)), en de p-waarde voor de z-statistiek (Pr(>|z|). Symbolen m.b.t. statistische significantie vind je rechts van de p-waarde waar van toepassing. De interpretatie van de symbolen wordt uitgelegd in de legende onder de coëfficiënten (“Signif. Codes:”).
• (Dispersion parameter…): Te negeren.
• Area that begins with Null deviance: Informatie over de fit van het model, besproken in een volgend hoofdstuk.
• Number of Fisher Scoring Iterations: Aantal iteraties van het algoritme.

We kunnen meerdere predictoren toevoegen aan het model op een gelijkaardige manier als bij lineaire regressie: door ze te scheiden met een + teken. Hier voegen we leeftijd (agea), vertrouwen in politici (trstplt), en links-rechts positie (lrscale) toe. Datamanagement voor deze variabelen is niet nodig: ze zijn continue en missing waarden zijn reeds als NA aangeduid.

#Schat het model
Vote_model_mp <- glm(vote_F ~ gndr_F + agea + trstplt + lrscale, 
                data = ESS9NL, family = "binomial")

#Bekijk de output
summary(Vote_model_mp)

Call:
glm(formula = vote_F ~ gndr_F + agea + trstplt + lrscale, family = "binomial", 
    data = ESS9NL)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.284194   0.380455  -0.747    0.455    
gndr_FMale   0.043281   0.154201   0.281    0.779    
agea         0.018349   0.004503   4.075 4.61e-05 ***
trstplt      0.195020   0.038706   5.039 4.69e-07 ***
lrscale      0.029257   0.039306   0.744    0.457    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1173.9  on 1424  degrees of freedom
Residual deviance: 1135.3  on 1420  degrees of freedom
  (248 observations deleted due to missingness)
AIC: 1145.3

Number of Fisher Scoring iterations: 4

Interpretatie

Logistische regressiecoëfficiënten geven een schatting van de verandering in de log van de odds dat Y=1 als X met 1 eenheid stijgt. Ze zijn dus niet makkelijk direct te interpreteren. We kunnen ze wel gebruiken om iets over de richting van de relatie te zeggen. Een positieve coëfficiënt toont dat de kans dat Y=1 stijgt als de onafhankelijke variabele stijgt. Een negatieve coëfficiënt toont dat de kans dat Y=1 daalt als de onafhankelijke variabele stijgt. Voor verdere interpretatie maak je best gebruik van de odds ratios (in minder mate), de gemiddelde marginale effecten (average marginal effects) (zie Hoofdstuk 10) of de voorspelde kansen (zie Hoofdstuk 11) .

Voor dit voorbeeld:

• Mannen hebben een grotere kans om te stemmen dan vrouwen, maar het verschil is niet statistisch significant (p = 0.28).
• Oudere respondenten hebben een grotere kans om te gaan stemmen en dit verband is statistisch significant (p < 0.001).
• Respondenten met meer vertrouwen in politici hebben een grotere kans om te gaan stemmen. deze relatie is statistisch significant (p < 0.001).
• Stemmen is meer waarschijnlijk naarmate respondenten zich rechtser positioneren op de ideologieschaal, maar dit effect is niet statistisch significant (p = 0.74).

9.2 Odds Ratios

Logistische regressiecoëfficiënten kunnen omgezet worden in odds ratios die (iets) intuïtiever zijn om te interpereteren.

We kunnen de odds ratios en 95% betrouwbaarheidsintervallen verkrijgen met de tidy functie uit het broom package:

# logistische regressiecoëfficiënten en hun betrouwbaarheidsintervallen
tidy(Vote_model_mp, conf.int = TRUE)

# A tibble: 5 × 7
  term        estimate std.error statistic     p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>       <dbl>    <dbl>     <dbl>
1 (Intercept)  -0.284    0.380      -0.747 0.455       -1.03       0.463 
2 gndr_FMale    0.0433   0.154       0.281 0.779       -0.259      0.346 
3 agea          0.0183   0.00450     4.07  0.0000461    0.00958    0.0272
4 trstplt       0.195    0.0387      5.04  0.000000469  0.119      0.271 
5 lrscale       0.0293   0.0393      0.744 0.457       -0.0480     0.106 

# odds ratios (i.e. 'exponentiële coëfficiënten') en hun betrouwbaarheidsintervallen
tidy(Vote_model_mp, conf.int = TRUE, exp = TRUE)

# A tibble: 5 × 7
  term        estimate std.error statistic     p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>       <dbl>    <dbl>     <dbl>
1 (Intercept)    0.753   0.380      -0.747 0.455          0.357      1.59
2 gndr_FMale     1.04    0.154       0.281 0.779          0.772      1.41
3 agea           1.02    0.00450     4.07  0.0000461      1.01       1.03
4 trstplt        1.22    0.0387      5.04  0.000000469    1.13       1.31
5 lrscale        1.03    0.0393      0.744 0.457          0.953      1.11

Zo lees je de syntax:

tidy(Vote_model_mp

We gebruiken de tidy functie op het model tussen haakjes.

conf.int = TRUE

We vragen R om de betrouwbaarheidsintervallen weer te geven. We kunnen ‘FALSE’ schrijven of deze code weglaten als we de betrouwbaarheidsintervallen niet willen.

exp = TRUE)

We vragen hier om de exponentiële (exponentiated) logistische regressiecoëfficiënten, oftewel de odds ratios. We kunnen ‘FALSE’ schrijven of deze code weglaten als we de logistische regressiecoëfficiënten willen.

Interpretatie

Voor de interpretatie van odds ratios zijn er 3 zaken waar je op moet letten.

Ten eerste, odds ratios vertellen ons iets over de relatieve odds dat Y = 1 (bv. iemand gaat stemmen). Dit is verschillend van de coëfficiënten. de coefficiënten zijn de gelogde versies van de relatieve odds..

Ten tweede, odds ratios zijn multiplicatief en worden geïnterpreteerd in termen van 1 in plaats van 0. Een odds ratio groter dan 1 betekent een hogere kans dat Y=1. Een odds ratio lager dan 1 betekent een lagere kans dat Y=1. Een odds ratio van 1 betekent dat er geen effect is. Een betrouwbaarheidsinterval voor een odds ratio waar 1 niet in voorkomt duidt een statistisch significant effect aan. Een odds ratio kan niet negatief zijn.

Ten derde interpreteren we ook de odds ratios met een multiplicatieve logica. In het voorbeeld vinden we dat de odds om te stemmen 1.04 keer groter zijn voor mannelijke respondenten dan vrouwelijke respondenten, ceteris paribus. Het effect is wel niet significant. De odds om te stemmen vermenigvuldigen met 1.02 telkens leeftijd met 1 eenheid omhoog gaat (de andere onafhankelijke variabelen constant gehouden).

---

1. Met factor() zouden we direct in de syntax de volgorde van de niveaus aan kunnen duiden: mutate(vote_binary = factor(vote, levels = c(2, 1), labels = c("Did not vote", "Voted")). Het gebruik van factor vermijdt een veelvoorkomnde fout besproken in volgend hoofdstuk.