4 Meervoudige Lineaire Regressie

In dit hoofdstuk ligt de focus op multiple (meervoudige) lineaire regressie, waarbij meerdere onafhankelijke variabelen gebruikt worden. We bespreken ook hoe gestandaardiseerde regressiecoëfficiënten te verkrijgen.

We beginnen weer met het laden van relevante R packages. Deze packages zijn reeds geïnstalleerd op de universitaire computers, maar moeten eerst geladen worden. We laden ook onze dataset.

#Packages
library(rio)             #laden van data
library(tidyverse)       #data manipulatie en grafieken
library(broom)           #samenvattingen regressiemodellen
library(parameters)      #berekenen gestandaardiseerde coëfficiënten

##Import data
demdata <- import("demdata.rds") |> 
  as_tibble()

4.1 Uitvoeren van de meervoudige lineaire regressie

In dit voorbeeld voorspellen we het niveau van electorale democratie in een land (v2x_polyarchy) aan de hand van 3 onafhankelijke variabelen (2 continue en 1 binair):

cpi: CPI staat voor “corruption perception index” en meet de mate van corruptie in de publieke sector van een land. Hogere waarden staan voor minder corruptie.
v2caviol: De variabele meet de mate van politiek geweld uitgevoerd door niet-statelijke actoren. Hogere waarden betekenen meer geweld.
TYPEDEMO1984: Binaire variabele die meet of een land in 1984 een democratie of autocratie was.

Voor we de regressie kunnen uitvoeren, moeten we eerst de binaire variabele transformeren naar een factor:

#omzetten naar factor variabele
demdata <- demdata |> 
  mutate(TYPEDEMO1984_factor = factorize(TYPEDEMO1984))

Voor meervoudige regressie gebruiken we ook de lm() functie. We kunnen meerdere onafhankelijke variabelen toevoegen met een ‘+’ teken:

#Model schatten en opslaan in data-object 
model_multiple <- lm(v2x_polyarchy ~ cpi + v2caviol + TYPEDEMO1984_factor, 
                     data=demdata)

multiple <-: We slaan de resultaten op in een data object dat we ‘multiple’ noemen. Deze naam kun je zelf bepalen.
lm(v2x_polyarchy ~): We voeren een lineaire regressie uit met de afhankelijke variabele “v2x_polyarchy”. Deze plaatsen we links van de tilde (~).
cpi + v2caviol + TYPEDEMO1984,: De onafhankelijke variabelen worden rechts van de tilde toegevoegd, van elkaar gescheiden door een ‘+’ teken. De volgorde maakt geen verschil voor de resultaten (wel de volgorde van de coëfficiënten in de output).
data = demdata): De naam van de dataset komt aan het einde van de syntax.

De resultaten bekijken we via summary():

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

Interpretatie

De interpretatie van de coëfficiënten is gelijkaardig aan die van bivariate modellen, maar we moeten wel de inclusie van meerdere predictoren in rekening brengen.

De “(Intercept)” waarde geeft weer welke waarde we kunnen verwachten voor de afhankelijke variabele als alle onafhankelijke variabelen de waarde 0 aannemen. We verwachten op basis van het model dat een land met score 0 op zowel cpi, v2caviol, als TYPEDEMO1984 (de referentiecategorie, namelijk een autocratie in 1984) gemiddeld een score op electorale democratie in 2020 van 0.19 zal hebben.

De coëfficiënten van de onafhankelijke variabelen vertellen ons nog steeds welke verandering we verwachten in de afhankelijke variabele als de predictor met 1 eenheid stijgt. Nu wordt dit effect echter “gecontroleerd op” de andere predictoren in het model. Het effect geldt als de andere variabelen constant worden gehouden (‘ceteris paribus’). Bijvoorbeeld:

v2caviol: Op basis van het model verwachten we dat electorale democratiescores dalen met -0.01 eenheden als politiek geweld met 1 eenheid stijgt, met de effecten van regimestatus in 1984 en corruptie constant gehouden.
TYPEDEMO1984_factor: Als we landen met dezelfde corruptie en politieke geweldscores vergelijken, verwachten we dat de electorale democratiescore in 2020 0.15 eenheden hoger is voor landen die in 1984 democratieën waren dan landen die autocratieën waren.

4.2 Gestandaardiseerde coëfficiënten

We kunnen in plaats van de ongestandaardiseerde coëfficiënten ook de gestandaardiseerde coëfficiënten berekenen. We kunnen hiervoor de standardize_parameters() functie gebruiken uit het parameters package.

multiple_std <- standardize_parameters(model_multiple, 
                       method = "refit")

De syntax lees je zo:

multiple_std <-: We slaan de resultaten op in een nieuw data object “multiple_std”.
standardize_parameters(multiple,: We passen de functie toe op het model tussen haakjes
method = 'refit'): We gebruiken de refit methode, de standaardmethode. Met deze methode worden de afhankelijke en onafhankelijke variabelen gestandaardiseerd en dan wordt het model opnieuw geschat met deze gestandaardiseerde versies.

We kunnen de resultaten vergelijken:

standardize_parameters() creëert een data frame met volgende kolommen:

glimpse(multiple_std)

Rows: 4
Columns: 5
$ Parameter       <chr> "(Intercept)", "cpi", "v2caviol", "TYPEDEMO1984_factor…
$ Std_Coefficient <dbl> -0.23661987, 0.49272847, -0.05393177, 0.60000039
$ CI              <dbl> 0.95, 0.95, 0.95, 0.95
$ CI_low          <dbl> -0.3957424, 0.3306681, -0.2037814, 0.3287296
$ CI_high         <dbl> -0.07749731, 0.65478881, 0.09591783, 0.87127121

Output uitleg

Parameter: Naam van de term of variabele in het model
Std_Coefficient: De waarde van de gestandaardiseerde coëfficiënt voor elke variabele
CI: Niveau van het betrouwbaarheidsinterval voor de gestandaardiseerde coëfficiënt.
CI_low en CI_high: De onder -en bovengrenzen van het betrouwbaarheidsinterval. Deze waarden worden gecombineerd in 1 cel als we de waarden straks printen.

We kunnen de resultaten vergelijken met het ongestandaardiseerde model. We gebruiken tidy() hier om de output te vereenvoudigen.

#Oorspronkelijk model
tidy(model_multiple)

# A tibble: 4 × 5
  term                           estimate std.error statistic      p.value
  <chr>                             <dbl>     <dbl>     <dbl>        <dbl>
1 (Intercept)                     0.187     0.0426      4.40  0.0000219   
2 cpi                             0.00636   0.00106     6.01  0.0000000155
3 v2caviol                       -0.00872   0.0123     -0.712 0.478       
4 TYPEDEMO1984_factorDemocracies  0.153     0.0349      4.37  0.0000239

#gestandaardiseerd model
multiple_std

# Standardization method: refit

Parameter                         | Std. Coef. |         95% CI
---------------------------------------------------------------
(Intercept)                       |      -0.24 | [-0.40, -0.08]
cpi                               |       0.49 | [ 0.33,  0.65]
v2caviol                          |      -0.05 | [-0.20,  0.10]
TYPEDEMO1984 factor [Democracies] |       0.60 | [ 0.33,  0.87]

Voor de continue variabelen geven de gestandaardiseerde coëfficiënten weer hoeveel standaardafwijkingen de afhankelijke variabele gaat veranderen als de onafhankelijke variabele met 1 standaardafwijking stijgt.¹

Voor factor variabelen ligt de interpretatie anders. De gestandaardiseerde coëfficiënt die we krijgen is de ongestandaardiseerde coëfficiënt gedeeld door de standaardafwijking van de afhankelijke variabele. De gestandaardiseerde coëfficiënten van continue en factor variabelen kunnen niet direct vergeleken worden.²

Interpretatie

We verwachten dat democratiescores met-0.05 standaardafwijkingen dalen als politiek geweld met 1 standaardafwijking stijgt (en met de effecten van corruptie en regimestatus in het verleden constant gehouden).

Als we landen met dezelfde corruptie en politieke geweldscores vergelijken, verwachten we dat de electorale democratiescore in 2020 0.6 standaardafwijkingen hoger is voor landen die in 1984 democratieën waren dan landen die autocratieën waren.

Waarschuwing!

Je zult opgemerkt hebben dat we noch summary() noch tidy() gebruikt hebben om de gestandaardiseerde coëfficiënten te printen in R. Deze functies zijn niet nodig omdat de output van standardize_parameters() reeds opgeslagen is in een dataframe.

Indien je summary() zou gebruiken zou je samenvattende statistieken vinden voor elke kolom in het dataframe:

summary(multiple_std)

  Parameter         Std_Coefficient         CI           CI_low        
 Length:4           Min.   :-0.2366   Min.   :0.95   Min.   :-0.39574  
 Class :character   1st Qu.:-0.0996   1st Qu.:0.95   1st Qu.:-0.25177  
 Mode  :character   Median : 0.2194   Median :0.95   Median : 0.06247  
                    Mean   : 0.2005   Mean   :0.95   Mean   : 0.01497  
                    3rd Qu.: 0.5195   3rd Qu.:0.95   3rd Qu.: 0.32921  
                    Max.   : 0.6000   Max.   :0.95   Max.   : 0.33067  
    CI_high        
 Min.   :-0.07750  
 1st Qu.: 0.05256  
 Median : 0.37535  
 Mean   : 0.38612  
 3rd Qu.: 0.70891  
 Max.   : 0.87127

Met tidy() krijg je een foutmelding gezien tidy() bedoeld is voor objecten die afkomstig zijn van statistische modellen:

tidy(multiple_std)

Warning in tidy.data.frame(multiple_std): Data frame tidiers are deprecated and
will be removed in an upcoming release of broom.

Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
returning NA

Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm =
na.rm): NAs introduced by coercion

Warning in mean.default(sort(x, partial = half + 0L:1L)[half + 0L:1L]):
argument is not numeric or logical: returning NA

Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
returning NA

Warning in mean.default(sort(x, partial = half + 0L:1L)[half + 0L:1L]):
argument is not numeric or logical: returning NA

Error in x - stats::median(x, na.rm = na.rm): non-numeric argument to binary operator

We zouden ook kunnen vragen enkel de onafhankelijke variabelen te standaardiseren en de schaal van de afhankelijke variabele te behouden met de optie “include_response = F” (F=False). Dit zou ons zeggen hoeveel Y verwacht wordt te veranderen op de originele schaal als de onafhankelijke variabele met 1 standaardafwijking stijgt. We kunnen dit doen als de schaal van de afhankelijke variabele zeer intuïtief is, bijvoorbeeld percentage stemmen voor een bepaalde partij.↩︎
De gestandaardiseerde coëfficiënten van continue en factor variabelen kunnen meer direct vergeleken worden als we de optie “two_sd = TRUE” toevoegen. De coëfficiënt van de continue onafhankelijke variabele geeft dan weer wat er gebeurt met Y als de onafhankelijke met 2 standaardafwijkingen stijgt, ongeveer het volledige bereik van de onafhankelijke variabele.↩︎