6  Model Fit en Modellen Vergelijken

Tot nu lag de focus op individuele coëfficiënten. Nu verschuiven we de focus naar model ‘fit’, oftewel de mate waarin het model bij de data past. We bekijken ook hoe we de fit van meerdere modellen kunnen vergelijken.

#Packages
library(rio)             #laden van data
library(tidyverse)       #datamanipulatie en grafieken
library(broom)           #samenvattingen regressiemodellen

##Importeer data en data management
demdata <- import("demdata.rds") |> 
  as_tibble()

demdata <- demdata |>  
  mutate(TYPEDEMO1984_factor = factorize(TYPEDEMO1984), 
         Typeregime2006_factor = factorize(Typeregime2006))

6.1 R2, Adjusted R2 en de F-Test

Ons voorbeeld hier is een regressiemodel waarin we de electorale democratiescore van een land in 2020 (v2x_polyarchy) voorspellen aan de hand van gepercipieerde corruptie in dat land (cpi), politiek geweld (v2caviol), en regimestatus in 1984 (TYPEDEMO1984_factor).

model_multiple <- lm(v2x_polyarchy ~ cpi + v2caviol + TYPEDEMO1984_factor, data=demdata)

De meeste model fit statistieken verkrijgen we simpelweg via de summary() functie:

summary(model_multiple)

Call:
lm(formula = v2x_polyarchy ~ cpi + v2caviol + TYPEDEMO1984_factor, 
    data = demdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.56402 -0.09376  0.01442  0.12926  0.34206 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     0.187394   0.042634   4.395 2.19e-05 ***
cpi                             0.006365   0.001059   6.012 1.55e-08 ***
v2caviol                       -0.008724   0.012258  -0.712    0.478    
TYPEDEMO1984_factorDemocracies  0.152698   0.034915   4.373 2.39e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1807 on 138 degrees of freedom
  (37 observations deleted due to missingness)
Multiple R-squared:  0.5068,    Adjusted R-squared:  0.4961 
F-statistic: 47.26 on 3 and 138 DF,  p-value: < 2.2e-16
Output uitleg

Model fit statistieken vinden we onderaan de output. “Multiple R-Squared” geeft de \(R^2\) (R kwadraat) statistiek. “Adjusted R-Squared” geeft de \(R^2\) gecorrigeerd voor het aantal predictoren in het model. De F-statistiek geeft informatie over de statistische significantie van het model.

  • Multiple R-squared: Dit toont de \(R^2\) (R kwadraat) statistiek, die meestal geïnterpreteerd wordt in termen van % van de variatie in Y verklaard door de predictoren in het model
  • Adjusted R-squared: Dit toont de \(R^2\) gecorrigeerd voor het aantal predictoren in het model.
  • F-statistic…: De F-statistiek geeft informatie over de statistische significantie van het model. Het eerste getal is de F-statistiek zelf (47.26). Het cijfer achter “p-value:” is de p-waarde voor de F-statistiek. De nulhypothese die hierbij getest wordt is dat geen enkele van de onafhankelijke variabelen (hier: cpi, v2caviol, TYPEDEMO1984) statistisch significant is. Een statistisch signifcante F-statistiek betekent dat op z’n minst 1 predictor significant is.

Deze output kunnen we ook verkrijgen via de glance() functie uit het broom package:

glance(model_multiple)
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.507         0.496 0.181      47.3 4.46e-21     3   43.5 -77.0 -62.3
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

De relevante statistieken vind je bij r.squared (R2), adj.r.squared (Adjusted R2), statistic (F-statistic), en p.value (p-waarde voor de F-statistiek) kolommen. nobs toont het aantal observaties gebruikt in het model.

6.2 Modellen vergelijken

De F-statistiek gaat na of het model een significant verbeterde voorspelling geeft dan een ‘nul model’ zonder predictoren, oftewel het gemiddelde van de afhankelijke variabele. We kunnen ook meerdere modellen vergelijken met elkaar. Hier vergelijken we een model met enkel cpi als onafhankelijke, dan een model met zowel cpi als v2caviol, en ten slotte een model met alle predictoren. Deze modellen zijn ‘nested’, dat wil zeggen dat meer uitgebreide modellen alle variabelen bevatten van de meer simpele modellen.

Om deze vergelijking te maken moeten we er wel voor zorgen dat onze modellen met dezelfde observaties werken en dus dezelfde N hebben. Dit kunnen we bereiken door een nieuwe dataset aan te maken met complete waarden (non-missing) voor alle variabelen die gebruikt worden in het meest complete model.

demdata_complete <- demdata |> 
  filter(complete.cases(v2x_polyarchy, cpi, v2caviol, TYPEDEMO1984_factor))

Deze dataset gebruiken we om onze modellen te schatten. Om een volledige vergelijking mogelijk te maken, schatten we ook een nulmodel zonder onafhankelijke variabelen met enkel een intercept (~ 1). Dit intercept bevat de gemiddelde waarde voor Y in de dataset (i.e. onze beste voorspelling zonder predictoren):

#Null model
model1 <- lm(v2x_polyarchy ~ 1, data = demdata_complete)

#Model with just cpi
model2 <- lm(v2x_polyarchy ~ cpi, data = demdata_complete)

#Model with cpi & v2caviol
model3 <- lm(v2x_polyarchy ~ cpi + v2caviol, data = demdata_complete)

#Model with all predictors
model4 <- lm(v2x_polyarchy ~ cpi + v2caviol + TYPEDEMO1984_factor, data = demdata_complete)

We kunnen de R2/Adj. R-Squared van de modellen vergelijken om te bekijken welk model het beste past. Dit geeft ons echter geen significantietoets:

Model R2 Adj. R2
Model 1 0 0
Model 2 0.437 0.433
Model 3 0.438 0.43
Model 4 0.507 0.496

Model 4 lijkt het beste te passen, maar om de significantietoets uit te voeren moeten we de anova() functie gebruiken. Deze is ingebouwd in R.

anova(model1, model2, model3, model4)
Analysis of Variance Table

Model 1: v2x_polyarchy ~ 1
Model 2: v2x_polyarchy ~ cpi
Model 3: v2x_polyarchy ~ cpi + v2caviol
Model 4: v2x_polyarchy ~ cpi + v2caviol + TYPEDEMO1984_factor
  Res.Df    RSS Df Sum of Sq        F    Pr(>F)    
1    141 9.1323                                    
2    140 5.1436  1    3.9887 122.2045 < 2.2e-16 ***
3    139 5.1285  1    0.0151   0.4614    0.4981    
4    138 4.5043  1    0.6243  19.1269 2.392e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova()

We voeren de functie uit op de modellen tussen haakjes.

Output uitleg

Het eerste deel van de output toont welke modellen vergeleken worden met elkaar. De onderste helft bevat het volgende:

  • Res.Df: De residual degrees of freedom (vrijheidsgraden) van het model
  • RSS: Dit staat voor “residual sum of squares”. RSS meet de variatie tussen de residuals in het model. RSS = \(\sum(y_{i} - \hat{y}_{i})^2\), waarbij \(\sum\) staat voor “sum up”, \(y_{i}\) is de geobserveerde Y voor een observatie in het model, and \(\hat{y}_{i}\) is de voorspelde waarde voor diezelfde observatie.1 RSS vertelt ons hoeveel van de variatie in Y het model niet kan verklaren of voorspellen. Een model met een lagere RSS voorspelt de Y beter, maar het verschil in RSS tussen modellen is niet altijd significant. We hebben dus nog een significantietoets nodig.
  • DF: Vrijheidsgraden. In de praktijk de hoeveelheid onafhankelijke variabelen toegevoegd in vergelijking met het voorgaande model. Dit getal is 1 als er 1 predictor werd toegevoegd in vergelijking met het vorige model in bovenstaande rij.2
  • Sum of Sq: De model of “regression” sum of squares is gebaseerd op de volgende formule: \(\sum(\hat{y}_{i} - \bar{y})^2\), waarbij \(\hat{y}_{i}\) de voorspelde waarde is voor een observatie in het model, en \(\bar{y}\) de gemiddelde waarde voor Y op basis van alle observaties in het model.3 De model sum of squares meet de variatie in Y die verklaart wordt door de predictoren in het model. De Sum of Sq in de anova() output toont de verandering in Sum of Sq ten opzichte van het voorgaande model. Hoe hoger de stijging hoe beter, maar hier moet ook een signifcantietest voor gebeuren.
  • F & Pr(>F): De F-statistiek en bijhorende p-waarde. De nulhypothese is dat het model in de desbetreffende rij niet beter past dan het model in de voorgaande rij. In feite test dit of tenminste 1 van de variabelen toegevoegd aan het model significant is. Indien de nulhypothese verworpen wordt, dan kunnen we zeggen dat het nieuwe model beter past.

We kunnen de output als volgt lezen: Model 2 past hier beter dan 1 (nulmodel), Model 3 past niet beter dan 2, en Model 4 past beter dan 3. We kunnen ook Modellen 1 en 2 direct met Model 4 vergelijken:

#Model 4 vs. Model 2
anova(model2, model4)
Analysis of Variance Table

Model 1: v2x_polyarchy ~ cpi
Model 2: v2x_polyarchy ~ cpi + v2caviol + TYPEDEMO1984_factor
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    140 5.1436                                  
2    138 4.5043  2   0.63935 9.7941 0.0001053 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Model 4 vs. Model 1
anova(model1, model4)
Analysis of Variance Table

Model 1: v2x_polyarchy ~ 1
Model 2: v2x_polyarchy ~ cpi + v2caviol + TYPEDEMO1984_factor
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    141 9.1323                                  
2    138 4.5043  3     4.628 47.264 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Waarschuwing!

anova() De volgorde van de modellen is van belang voor de uitkomst. Hierboven vergelijken we telkens een complexer model met een simpeler model. Indien we schrijven “(model4, model1, model2”, “model3”), dan vergelijkt R model 1 tegen 4, model 2 tegen 1 enz.

anova(model4, model1, model2, model3)
Analysis of Variance Table

Model 1: v2x_polyarchy ~ cpi + v2caviol + TYPEDEMO1984_factor
Model 2: v2x_polyarchy ~ 1
Model 3: v2x_polyarchy ~ cpi
Model 4: v2x_polyarchy ~ cpi + v2caviol
  Res.Df    RSS Df Sum of Sq        F Pr(>F)    
1    138 4.5043                                 
2    141 9.1323 -3   -4.6280  47.2642 <2e-16 ***
3    140 5.1436  1    3.9887 122.2045 <2e-16 ***
4    139 5.1285  1    0.0151   0.4614 0.4981    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

De tweede rij vergelijkt nu het nulmodel (model1) met het meest complexe model (model4). We krijgen een negatieve waarde voor “DF” en “Sum of Sq” omdat Model 1 minder predictoren heeft en ook minder goed past. Het verschil is statistisch significant. Dit interpreteren we nu als: model 4 is beter dan model 1.De derde rij vergelijkt model2 (enkel cpi als predictor) met model1 (het nulmodel). De resultaten zijn dezelfde als hierboven. De resultaten voor de laatste rij zijn ook dezelfde.

Let erop dat de namen die de anova() functie geeft aan de modellen niet noodzakelijk dezelfde zijn als de namen die je zelf geeft (model 1 voor anova is nu ons model 4).

Je kunt de fit van bepaalde modellen testen tegenover elkaar om zo stapsgewijs het beste model te vinden. Meestal zullen we meer complexe modellen vergelijken met meer simpele modellen. In de syntax gaan we dan van meest simpel naar meest complex model.

  1. Deze vergelijking behoort niet tot de leerstof.↩︎

  2. De DF kolom geeft op zich niet weer hoeveel extra onafhankelijke variabelen werden toegevoegd, maar wel hoeveel nieuwe coëfficiënten (of termen) werden toegevoegd. Dit is vooral van belang bij factor variabelen (zeker als ze 3 of meer categorieën hebben). Hoewel je misschien 1 factor variabele toevoegt, kun je meer dan 1 coëfficiënt (en dus DF) krijgen als je voor meerdere categorieën dummy variabelen moet toevoegen.↩︎

  3. Deze vergelijking behoort niet tot de leerstof.↩︎