12  Model Fit en Modellen Vergelijken

In vorige hoofdstukken lag de focus op de effecten van de onafhaneklijke variabelen en hoe ze te begrijpen op basis van coëfficiënten, odds ratios, marginale effecten en voorspelde kansen. In dit hoofdstuk richten we ons op het model als geheel en hoe goed het bij de data past (‘fit’).

We laden de packages en data:

#Packages
#Packages
library(rio)             #laden van data
library(tidyverse)       #datamanipulatie en grafieken
library(performance)     #goodness-of-fit statistieken en tests

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

12.1 Fit statistieken met summary()

Laten we teruggaan naar het logistisch model waar we ook in vorige hoofdstukken mee werkten: wat is de kans dat iemand gaat stemmen op basis van informatie over gender, leeftijd, vertouwen in politici en links-rechtsideologie?

#Data management
ESS9NL <- ESS9NL |>
  #Factor maken van categorische variabelen
  mutate(gndr_F = factorize(gndr), 
         vote_F = factorize(vote))  |> 
  #Not Eligible op missing zetten
  mutate(vote_F = na_if(vote_F,"Not eligible to vote")) |> 
  #Relevel van variabelen
  mutate(vote_F = relevel(vote_F, "No"), 
         gndr_F = relevel(gndr_F, "Female"))

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

#Resultaten printen
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
Output uitleg

Zoals bij een lineair model (lm), zal de summary() functie in het onderste gedeelte van de output informatie bevatten over model fit. We krijgen informatie over de “Null” en “Residual Deviance” statistieken. De Residual Deviance statistiek duidt het verschil (“deviance”) aan tussen het geschatte model en een “perfect’ model dat precies bij de data past. De Null Deviance statistiek doet dezelfde vergelijking, maar ten opzichte van een nulmodel dat enkel het intercept bevat.

Kleinere Residual Deviance waarden duiden beter passende modellen aan. Echter is het niet gewenst om de deviance statistiek op zich te interpreteren gezien de schaal onduidelijk is en er geen maximumwaarde is. Daarom maken we gebruik van andere statistieken en een test gebaseerd op de deviance statistiek (zie onder).

12.2 Modellen vergelijken: Likelihood Ratio Test

We kunnen de likelihood ratio test gebruiken om verschillende logistische regressiemodellen te vergelijken met elkaar en na te gaan welke beter past. De LRT berekent de ratio tussen de deviance statistieken van de modellen en gaat na of er een significant verschil is.

Als we modellen willen verglijken moeten we net zoals bij lineaire regressie (zie Paragraaf 6.2) zorgen dat de modellen een gelijke N hebben en dat complexere modellen alle predictors bevatten van simpelere modellen (nested). We zorgen dat we eerst een dataset maken met complete observaties voor het meest complexe model (alle predictors).

ESS9NL_glm <- ESS9NL |>
  filter(complete.cases(vote_F,  gndr_F,  agea,  trstplt,  lrscale))

Dan schatten we een reeks modellen waaraan we telkens 1 van de onafhankelijke varaibelen toevoegen. We beginnen met een nulmodel dat enkel het intercept bevat. Het model met 1 onafhankelijke variabele kan daar dan mee vergeleken worden.

#Nulmodel
Vote_model0 <- glm(vote_F ~ 1,
                      data = ESS9NL_glm, family = "binomial")
# + gndr
Vote_model1 <- glm(vote_F ~ gndr_F , 
                data = ESS9NL_glm, family = "binomial")
# + agea
Vote_model2 <- glm(vote_F ~ gndr_F + agea , 
                data = ESS9NL_glm, family = "binomial")

# + trst
Vote_model3 <- glm(vote_F ~ gndr_F + agea + trstplt, 
                data = ESS9NL_glm, family = "binomial")

# + lrscale
Vote_model4 <- glm(vote_F ~ gndr_F + agea + trstplt + lrscale, 
                data = ESS9NL_glm, family = "binomial")

Nu kunnen we de likelihood ratios van deze modellen met elkaar vergelijken en significantietoetsen uitvoeren. We gebruiken het performance package met de test_likelihoodratio functie. De test vergelijkt de deviance statistiek (-2LL) van elk model, de verandering in vrijheidsgraden (df= degrees of freedom) per model, en gebruikt een Chi2 (\(\chi^2\)) toets.

test_likelihoodratio(Vote_model0,
                     Vote_model1,
                     Vote_model2,
                     Vote_model3,
                     Vote_model4)
# Likelihood-Ratio-Test (LRT) for Model Comparison (ML-estimator)

Name        | Model | df | df_diff |  Chi2 |      p
---------------------------------------------------
Vote_model0 |   glm |  1 |         |       |       
Vote_model1 |   glm |  2 |       1 |  0.11 |  0.744
Vote_model2 |   glm |  3 |       1 | 13.59 | < .001
Vote_model3 |   glm |  4 |       1 | 24.38 | < .001
Vote_model4 |   glm |  5 |       1 |  0.55 |  0.457
test_likelihoodratio(

: We voeren de likelihood ratio test uit op de modellen tussen haakjes. Er moeten minstens twee modellen aangeduid zijn en de volgorde bepaalt welke vergelijking gemaakt wordt.

Output uitleg

De output lees je als volgt:

  • Name: naam van het object waarin het model is opgeslagen
  • Model: informatie over het type model. Te negeren.
  • df: Geeft weer hoeveel termen gebruikt werden in het model. Vote_model0 heeft een df van 1 gezien er maar 1 term is: het intercept. Vote_model4 heeft 5 df omdat er 5 termen zijn: het intercept en de coëfficiënten voor de 4 onafhankelijke variabelen.
  • df_diff: Geeft weer hoeveel het model verschilt van het vorige model in termen van df. Dit is telkens 1 hier omdat we telkens maar 1 nieuwe predictor hebben toegevoegd.
  • Chi2 & p: Dit is de Chi2 statistiek en bijhorende p-waarde. De test gaat na of de fit van een model beter is dan de fit van het model in de rij erboven. De nulhypothese is dat er geen verschil is in fit. Een significante toets betekent dat het model beter past.

In dit voorbeeld:

  • Model 1 heeft geen significant betere fit dan een nulmodel (vote_model0)
  • Model 2 heeft een betere fit dan Model 1
  • Model 3 heeft een betere fit dan Model 1 Model 2
  • Model 4 heeft geen significant betere fit dan Model 3.

We kunnen concluderen dat Model 3 (vote_model3) het best passende model is zonder inclusie van nietszeggende variabelen (i.e. het model is het meest ‘parsimonious’).

Zoals het geval was voor de anova() functie bij lineaire regressie kun je ook specifieke groepen van modellen vergelijken:

#Past Model 4 beter dan Model 1?: Ja!
test_likelihoodratio(Vote_model1, Vote_model4)
# Likelihood-Ratio-Test (LRT) for Model Comparison (ML-estimator)

Name        | Model | df | df_diff |  Chi2 |      p
---------------------------------------------------
Vote_model1 |   glm |  2 |         |       |       
Vote_model4 |   glm |  5 |       3 | 38.53 | < .001
#Past Model 3 beter dan een nulmodel?: Ja!
test_likelihoodratio(Vote_model0, Vote_model3)
# Likelihood-Ratio-Test (LRT) for Model Comparison (ML-estimator)

Name        | Model | df | df_diff |  Chi2 |      p
---------------------------------------------------
Vote_model0 |   glm |  1 |         |       |       
Vote_model3 |   glm |  4 |       3 | 38.08 | < .001
Waarschuwing!

De volgorde waarin we onze modellen aanduiden in de test_likelihoodratio() syntax bepaalt welke modellen precies vergeleken worden net zoals met anova() bij lineaire regressie ( Paragraaf 6.2). Bij een verkeerde volgorde krijg je een error in R. Een juiste volgorde houdt in dat je van minder naar meer complex gaat:

test_likelihoodratio(Vote_model0,
                     Vote_model4,
                     Vote_model2,
                     Vote_model1,
                     Vote_model3)
Error: The models are not nested, which is a prerequisite for
  `test_likelihoodratio()`.
  See the 'Details' section.
  You may try `test_vuong()` instead.

Als we 2 modellen testen en de eerste in de syntax is de meest complexe, dan krijgen we dezelfde Chi2 en p-waarde vergeleken met een juiste volgorde, maar de df_diff zal negatief zijn (-3 ipv +3 in dit voorbeeld). Op zich is dit geen probleem, zolang we maar weten wat we precies aan het vergelijken zijn zodat we geen interpretatiefouten maken.

test_likelihoodratio(Vote_model4, Vote_model1)
# Likelihood-Ratio-Test (LRT) for Model Comparison (ML-estimator)

Name        | Model | df | df_diff |  Chi2 |      p
---------------------------------------------------
Vote_model4 |   glm |  5 |         |       |       
Vote_model1 |   glm |  2 |      -3 | 38.53 | < .001

12.3 Pseudo R2

Bij een lineair regressiemodel beoordelen we fit met de R2 waarde. Een logistisch model is anders geschat en dus hebben we deze waarde niet. Verschillende zogenaamde pseudo R2 statistieken werden ontwikkeld om meer intuïtief inzicht te verkrijgen in de verklarende kracht van een model. De pseudo R2 maatstaven zijn gebaseerd op de likelihood ratio test en kunnen niet als ‘proportie verklaarde variantie’ geïnterpreteerd worden.

Hier gebruiken we de Nagelkerke R². De waarde van deze maatstaf ligt tussen 0 en 1. Lage waarden wijzen op een lage verklarende kracht, hoge waarden op een hoge verklarende kracht.

We kunnen de Nagelkerke R² statistiek opvragen met de r2_nagelkerke() functie uit het performance package.

# Nagelkerke R2: Model 3
r2_nagelkerke(Vote_model3)
Nagelkerke's R2 
     0.04698189 
# Nagelkerke R2: Model 4
r2_nagelkerke(Vote_model4)
Nagelkerke's R2 
     0.04765513
r2_nagelkerke(

Deze functie berekent de Nagelkerke R2 voor het model tussen haakjes. Er kan slecht 1 model opgegeven worden.

De Nagelkerke R2 is hoger voor Model 4 dan Model 3. Echter is een likelihood ratio test nodig om te weten of dit verschil significant is. Zoals we hierboven zagen is dit niet het geval.

Waarschuwing!

Er bestaan verschillende pseudo R2 statistieken om de fit van logistische regressiemodellen te helpen interpreteren. Geen enkele van hen kan geïnterpreteerd worden in termen van ‘proprotie verklaarde variantie’.