library(rio) #laden van data
library(tidyverse) #data manipulatie en grafieken
library(broom) #samenvattingen regressiemodellen
library(marginaleffects) #voorspelde waarden
# inladen data en data management
demdata <- import("demdata.rds") |>
as_tibble()
demdata <- demdata |>
mutate(TYPEDEMO1984_factor = factorize(TYPEDEMO1984),
Typeregime2006_factor = factorize(Typeregime2006))
# Modellen schatten
model_continuous <- lm(v2x_polyarchy ~ gini_2019, data = demdata)
model_binary <- lm(v2x_polyarchy ~ TYPEDEMO1984_factor, data=demdata)
model_categorical <- lm(v2x_polyarchy ~ Typeregime2006_factor, data=demdata)
model_multiple <- lm(v2x_polyarchy ~ cpi + v2caviol + TYPEDEMO1984_factor, data=demdata)5 Voorspellingen en Fouten
Een regressiemodel geeft voorspellingen van de waarde op de afhankelijke variabele op basis van de waarden op de onafhankelijke variabele(n). We kunnen deze voorspelde waarden onderzoeken met R om onze resultaten beter te begrijpen.
We laden de relevante R packages eerst. Deze packages zijn reeds geïnstalleerd op de universitaire computers, maar moeten eerst geladen worden. We laden ook onze dataset en schatten enkele modellen die we gebruiken in dit hoofdstuk.
5.1 Voorspellingen en fouten voor de observaties in het model
Op basis van het lineaire regressiemodel kunnen we voor elke observatie gebruikt in het model een voorspelling maken van de waarde voor de afhankelijke waarde. Het verschil tussen deze voorspelling en de echte waarde die we vinden in de dataset is de fout (‘error’) of ‘residual’.
De predictions() functie uit het marginaleffects package kan gebruikt worden om voorspellingen te maken voor elke observatie gebruikt in het model. 1
model_binary_predictions <- predictions(model_binary, newdata = demdata) |>
as_tibble() #as_tibble() niet strikt nodig, zie waarschuwingsvak hieronderZo lees je de syntax:
model_binary_predictions-
We slaan de output hier op in een nieuw data object “model_binary_predictions”. Deze naam kun je zelf bepalen.
predictions(model_binary,-
We gebruiken de
predictionsfunctie op het model tussen haakjes. newdata = demdata)-
Hier verduidelijken we de originele dataset voor deze voorspellingen. Deze syntax vertelt R dat we in ons nieuwe data object de voorspellingen willen, maar ook alle variabelen uit de originele dataset, niet enkel de variabelen gebruikt in het model. Dit is nuttig als we specifieke observaties willen identificeren (bv. door te kijken naar de naam van het land). Als je dit niet specificeert krijg je een dataset zonder de overige variabelen in de originele dataset.
De output kunnen we printen met behulp van de volgende code:
model_binary_predictions# A tibble: 179 × 51
rowid estimate std.error statistic p.value s.value conf.low conf.high
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0.690 0.0287 24.1 3.16e-128 424. 0.634 0.746
2 2 0.418 0.0233 17.9 1.16e- 71 236. 0.372 0.463
3 3 0.690 0.0287 24.1 3.16e-128 424. 0.634 0.746
4 4 0.690 0.0287 24.1 3.16e-128 424. 0.634 0.746
5 5 0.418 0.0233 17.9 1.16e- 71 236. 0.372 0.463
6 6 0.418 0.0233 17.9 1.16e- 71 236. 0.372 0.463
7 7 0.690 0.0287 24.1 3.16e-128 424. 0.634 0.746
8 8 0.418 0.0233 17.9 1.16e- 71 236. 0.372 0.463
9 9 NA NA NA NA NA NA NA
10 10 0.418 0.0233 17.9 1.16e- 71 236. 0.372 0.463
# ℹ 169 more rows
# ℹ 43 more variables: country_name <chr>, year <dbl>, v2x_polyarchy <dbl>,
# v2x_libdem <dbl>, v2x_egaldem <dbl>, v2cacamps <dbl>, v2caviol <dbl>,
# e_peaveduc <dbl>, cpi <dbl>, e_regiongeo <dbl>, e_regionpol_6C <dbl>,
# v2elcomvot <dbl>, compulsory_voting <dbl>, bicameral <dbl>, dem_diff <dbl>,
# dem_increase <dbl>, dem_decrease <dbl>, TypeSoc2005 <dbl>,
# TypeEcon2006 <dbl>, HDI2005 <dbl>, GDP2006 <dbl>, TYPEDEMO1984 <dbl>, …
Output uitleg
- estimate: Dit is de voorspelde waarde op de afhankelijke variabele voor elke observatie in het model. Observaties die niet in het model werden opgenomen (omwille van ontbrekende data) krijgen hier ‘NA’.
- std.error, statistic, p.value, conf.low, en conf.high: de standaardfout van de voorspelling, t-statistiek, p-waarde en het 95% betrouwbaarheidsinterval. De s-waarde is een andere manier om onzekerheid weer te geven maar behoort niet tot de leerstof. 2
- De overige kolommen bevatten de variabelen uit de originele dataset.
We kunnen het nieuwe dataobject gebruiken om ook de residuals te berekenen. Dit doen we door het verschil tussen echte en voorspelde waarde in een variabele op te nemen. 3
model_binary_predictions <- model_binary_predictions |>
mutate(residual_value = v2x_polyarchy - estimate) #residual = echte waarde - voorspelde waardeDeze variabele kunnen we gebruiken om na te gaan welke observaties goed of slecht worden voorspeld. Dit kan nuttig zijn bij het nagaan of aan assumpties voldaan is, zie Hoofdstuk 7 .
model_binary_predictions |>
select(country_name, v2x_polyarchy, estimate, residual_value) # A tibble: 179 × 4
country_name v2x_polyarchy estimate residual_value
<chr> <dbl> <dbl> <dbl>
1 Mexico 0.647 0.690 -0.0432
2 Suriname 0.761 0.418 0.343
3 Sweden 0.908 0.690 0.218
4 Switzerland 0.894 0.690 0.204
5 Ghana 0.72 0.418 0.302
6 South Africa 0.703 0.418 0.285
7 Japan 0.832 0.690 0.142
8 Myanmar 0.436 0.418 0.0184
9 Russia 0.262 NA NA
10 Albania 0.485 0.418 0.0674
# ℹ 169 more rows5.2 Voorspellingen voor bepaalde waarden van de onafhankelijke variabele (Bivariaat)
We kunnen ook nagaan welke waarde op de afhankelijke we kunnen verwachten volgens het model als de onafhankelijke variabele bepaalde waarden aanneemt. Bijvoorbeeld: welke democratiescore kunnen we gemiddeld verwachten voor landen die in 1984 een autocratie waren? Of voor landen die een lage of hoge economische ongelijkheid kennen? We kunnen hier ook de predictions() functie voor gebruiken.
Eerst voorspellen we de verwachte democratiescore in 2020 voor landen die in 1984 een autocratie versus democratie waren op basis van ons bivariaat model (model_binary).
predictions(model_binary,
by = 'TYPEDEMO1984_factor') |>
as_tibble() # A tibble: 2 × 10
rowid TYPEDEMO1984_factor estimate std.error statistic p.value s.value
<int> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 Autocracies 0.418 0.0233 17.9 1.16e- 71 236.
2 2 Democracies 0.690 0.0287 24.1 3.16e-128 424.
# ℹ 3 more variables: conf.low <dbl>, conf.high <dbl>, rowid_dedup <int>predictions(model_binary,-
We passen de functie toe op het model tussen haakjes.
by = "TYPEDEMO1984_factor")-
Hier vragen we de voorspelling voor elk niveau (level) van de factor “TYPEDEMO1984_factor”. De “by=” syntax wordt enkel gebruikt met factor variabelen. We maken geen gebruik van “newdata=” omdat we hier geen voorspellingen vragen voor alle observaties.
We kunnen ook voorspellingen maken op basis van een continue onafhankelijke variabele. We kunnen bijvoorbeeld de score voor electorale democratie voorspellen aan de hand van economische ongelijkheid (gini_2019). Hier gaan we na welke democratiescore we kunnen verwachten als ongelijkheid laag (25) versus hoog is (45).
predictions(model_continuous,
newdata = datagrid(gini_2019 = c(25,45))) |>
as_tibble()# A tibble: 2 × 10
rowid estimate std.error statistic p.value s.value conf.low conf.high
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0.764 0.0451 16.9 3.19e-64 211. 0.675 0.852
2 2 0.527 0.0493 10.7 1.10e-26 86.2 0.430 0.623
# ℹ 2 more variables: gini_2019 <dbl>, v2x_polyarchy <dbl>newdata = datagrid(gini_2019 = c(25,45))-
Hier bepalen we de waarden van de continue onafhankelijke variabele waar we voorspellingen voor willen maken. Je kunt de naam van de variabele veranderen, alsook de waarden waarvoor je voorspellingen maakt. De rest van de syntax blijft gelijk.
We kunnen eventueel meerdere waarden toevoegen om voorspellingen voor te maken door de code op de volgende manier uit te breiden in het c() gedeelte van de syntax:
predictions(model_continuous,
newdata = datagrid(gini_2019 = c(25,30,35,40,45))) |>
as_tibble()# A tibble: 5 × 10
rowid estimate std.error statistic p.value s.value conf.low conf.high
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0.764 0.0451 16.9 3.19e- 64 211. 0.675 0.852
2 2 0.705 0.0316 22.3 2.20e-110 364. 0.643 0.766
3 3 0.645 0.0267 24.2 6.30e-129 426. 0.593 0.698
4 4 0.586 0.0345 17.0 1.04e- 64 213. 0.518 0.654
5 5 0.527 0.0493 10.7 1.10e- 26 86.2 0.430 0.623
# ℹ 2 more variables: gini_2019 <dbl>, v2x_polyarchy <dbl>
Waarschuwing!
We eindigden de predictions() functie met as_tibble(). Deze stap is niet strikt noodzakelijk. Dit is het resultaat zonder de toevoeging:
predictions(model_binary, by = 'TYPEDEMO1984_factor')
TYPEDEMO1984_factor Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
Autocracies 0.418 0.0233 17.9 <0.001 235.6 0.372 0.463
Democracies 0.690 0.0287 24.1 <0.001 423.5 0.634 0.746
Type: response Het verschil zit hem in de weergave van de output in R: standaard geeft predictions() andere namen aan de kolommen (bv., Estimate i.p.v. estimate, 2.5% i.p.v. conf.low) om de zaken netter te maken, maar dit bemoeilijkt de zaken eigenlijk vaak voor ons omdat dit niet de echte variabelenamen zijn zoals ze opgeslagen worden in het object. Later in het vak gebruiken we deze variabelen om verdere bewerkingen te doen. Daarvoor moeten we de juiste variabelenamen opgeven: estimate dus en niet Estimate.
5.3 Voorspelde waarden (Meervoudige Lineaire Regressie)
Voorspelde waarden en fouten kunnen we ook voor meervoudige regressie bekijken via de predictions() functie. De procedure om voorspelde waarden te vinden voor alle observaties in het model is dezelfde als hierboven dus herhalen we deze niet. De procedure voor voorspellingen op basis van waarden van een onafhankelijke variabele is gelijkaardig, met 1 belangrijk verschil voor factor variabelen.
5.3.1 Voorspellingen voor een continue predictor
Dit waren de resultaten van ons meervoudig lineair regressiemodel:
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 cpi meet gepercipieerde corruptie in een land op een schaal van 0 tot 100 (hogere waarden staan voor minder corruptie). In de praktijk is het bereik van de variabele in ons model 12 tot 88. Voor we voorspellingen doen gaan we het werkelijke bereik eerst na:
predictions(model_multiple) |>
select(cpi) |>
summary()- 1
-
We gebruiken de
predictions()functie hier om enkel observaties te selecteren die gebruikt werden in het model (observaties met ontbrekende waarden op ‘NA’ worden weggefilterd). - 2
- We selecteren de cpi variabele
- 3
- En vragen de beschrijvende statistieken voor de variabele.
cpi
Min. :12.00
1st Qu.:28.00
Median :39.50
Mean :43.37
3rd Qu.:56.75
Max. :88.00 We kunnen voorspelde waarden gebruiken om een inschatting te maken over het verwachte niveau van democratie bij lage en hoge corruptie. Een regressiecoëfficiënt zegt ons wat er gebeurt als corruptie met 1 eenheid stijgt, maar voorspelde waarden kunnen vaak een intuïtiever beeld geven over de sterkte van een effect. Hier gebruiken we predictions() om verwachte democratiescores te berekenen voor corruptiescores (cpi) van 20 tot 80 met verhogingen van telkens 10 eenheden.
preds1 <- predictions(model_multiple,
newdata = datagrid(cpi = c(20,30,40,50,60,70,80))) |>
as_tibble()preds1 <--
We slaan de resultaten op in een data object omdat we ze ook voor andere doeleinden zullen gebruiken. De naam bepaal je zelf.
predictions(multiple,-
We passen de functie toe op het model tussen haakjes.
newdata = datagrid(cpi = c(20,30,40,50,60,70,80))-
Hier bepalen we voor welke onafhankelijke variabele we voorspellingen willen (cpi) en voor welke waarden (20…80). De waarden zijn numeriek en gaan niet tussen aanhalingstekens.
We printen de voorspellingen:
preds1# A tibble: 7 × 12
rowid estimate std.error statistic p.value s.value conf.low conf.high
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0.318 0.0278 11.4 2.41e- 30 98.4 0.264 0.373
2 2 0.382 0.0217 17.6 3.30e- 69 227. 0.339 0.424
3 3 0.445 0.0199 22.4 2.59e-111 367. 0.406 0.484
4 4 0.509 0.0233 21.9 6.10e-106 350. 0.463 0.555
5 5 0.573 0.0302 18.9 4.92e- 80 263. 0.513 0.632
6 6 0.636 0.0389 16.4 2.78e- 60 198. 0.560 0.713
7 7 0.700 0.0483 14.5 1.17e- 47 156. 0.605 0.795
# ℹ 4 more variables: v2caviol <dbl>, TYPEDEMO1984_factor <fct>, cpi <dbl>,
# v2x_polyarchy <dbl>
Output uitleg
- estimate: De voorspelde waarde
- kolommen “std.error” tot “conf.high”: informatie met betrekking tot onzekerheid van de schatting
- “4 more variables”: Dit zegt dat ons tidied dataframe nog 4 variabelen heeft (dit verschilt naargelang het model dat je gebruikt). De kolommen zijn genoemd naar de variabelen gebruikt in het model. Voor de onafhankelijke variabelen (hier:
v2caviol,TYPEDEMO1984_factor, encpi) tonen ze de waarden die gebruikt worden voor deze variabelen om de voorspellingen te maken.
In bovenstaand voorbeeld houdt predictions() automatisch de 2 overige onafhankelijke variabelen (v2caviol en TYPEDEMO1984_factor) constant op dezelfde waarde bij de berekening van elke voorspelde waarde. Continue predictoren worden constant gehouden op hun gemiddelde, voor factor variabelen wordt de modus (de meest voorkomende categorie) gebruikt. Dit kunnen we nagaan:
preds1 |>
select(estimate, cpi, v2caviol, TYPEDEMO1984_factor)# A tibble: 7 × 4
estimate cpi v2caviol TYPEDEMO1984_factor
<dbl> <dbl> <dbl> <fct>
1 0.318 20 -0.394 Autocracies
2 0.382 30 -0.394 Autocracies
3 0.445 40 -0.394 Autocracies
4 0.509 50 -0.394 Autocracies
5 0.573 60 -0.394 Autocracies
6 0.636 70 -0.394 Autocracies
7 0.700 80 -0.394 Autocracies 5.3.2 Voorspellingen voor een factor predictor
We kunnen een gelijkaardige procedure gebruiken om voorspellingen te maken voor de verschillende niveaus van factor variabelen. Dit kunnen we doen met behulp van de by= optie i.p.v newdata = datagrid().4 Om ervoor te zorgen dat we voor de overige onafhankelijke variabelen het gemiddelde of de modus nemen, moeten we hier wel nog syntax toevoegen via newdata:
preds2 <- predictions(model_multiple, by= "TYPEDEMO1984_factor",
newdata = "mean") |>
as_tibble()by="TYPEDEMO1984_factor"-
Hier verduidelijken we dat we voorspellingen willen voor elk niveau van de factor variabele.
newdata = "mean")-
Hier zeggen we dat voor de overige onafhankelijke variabelen het gemiddelde of de modus aangehouden moet worden. Dit gebeurde automatisch in vorig voorbeeld, maar moet toegevoegd worden als we “by=” gebruiken.
De resultaten zijn als volgt:
preds2# A tibble: 2 × 12
rowid TYPEDEMO1984_factor estimate std.error statistic p.value s.value
<int> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 Autocracies 0.467 0.0205 22.8 4.94e-115 380.
2 2 Democracies 0.620 0.0260 23.8 2.47e-125 414.
# ℹ 5 more variables: conf.low <dbl>, conf.high <dbl>, cpi <dbl>,
# v2caviol <dbl>, rowid_dedup <int>Opnieuw kunnen we zien dat predictions() de andere onafhankelijke variabelen constant houdt:
preds2 |>
select(estimate, TYPEDEMO1984_factor, cpi, v2caviol)# A tibble: 2 × 4
estimate TYPEDEMO1984_factor cpi v2caviol
<dbl> <fct> <dbl> <dbl>
1 0.467 Autocracies 43.4 -0.394
2 0.620 Democracies 43.4 -0.3945.3.3 Voorspellingen voor specifieke waarden van de onafhankelijke variabelen
We kunnen predictions() ook gebruiken om voorspellingen te maken voor specifieke, hypothetische casussen. Bijvoorbeeld, hier vragen we de voorspelde waarde voor de afhankelijke variabele electorale democratie voor een land dat: een democratie was in 1984, 88 scoort op de corruptieperceptie-index (de maximumwaarde in de dataset) en -3.429 voor politiek geweld (de minimumwaarde in de dataset).
We bepalen deze waarden in het newdata = datagrid() gedeelte van de syntax. Indien we een variabele niet zouden specificeren zou deze constant gehouden worden op het gemiddelde of de modus.
predictions(model_multiple,
newdata = datagrid(cpi=c(88),
v2caviol=c(-3.429),
TYPEDEMO1984_factor=c("Democracies"))) |>
as_tibble()# A tibble: 1 × 12
rowid estimate std.error statistic p.value s.value conf.low conf.high cpi
<int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0.930 0.0392 23.7 1.89e-124 411. 0.853 1.01 88
# ℹ 3 more variables: v2caviol <dbl>, TYPEDEMO1984_factor <fct>,
# v2x_polyarchy <dbl>De
augment()functie uit hetbroompackage kunnen we ook gebruiken om de residuals te bestuderen. Dit gebruiken we in een ander hoofdstuk. Hier richten we ons oppredictions()omdat deze functie gemakkelijker een dataframe produceert met de voorspellingen, fouten en de overige data in de originele dataset.↩︎De s-waarde is een poging om de p-waarde te vertalen naar een maat die volgens sommigen gemakkelijker te interpreteren is. In het bijzonder vertelt het ons: “Hoeveel opeenvolgende”kop”-worpen zouden dezelfde hoeveelheid bewijs (of “verrassingen”) leveren tegen de nulhypothese dat de munt eerlijk is?” Een p-waarde van 0,05 zou bijvoorbeeld een overeenkomstige s-waarde van 4,3 of zo hebben. We zouden dan kunnen zeggen dat een p-waarde van 0,05 ongeveer net zo verrassend is als vier keer een eerlijke munt opgooien en de munt alle vier de keren op kop zien landen. Zou je je gerust voelen om een verklaring af te leggen dat de munt vals is in plaats van eerlijk op basis van die reeks muntworpen? In de context van de output van
predictions()(en van deslopes()-functie die we in latere hoofdstukken zien), zouden hogere s-waarden aangeven dat we steeds verraster zouden moeten zijn om onze resultaten te zien als de waarde van het ding dat we schatten eigenlijk 0 is. Deze statistiek is niet zo nuttig voor onze voorspelde waarden, maar zou nuttiger kunnen zijn om te begrijpen hoe verrassend een schatting van een coëfficiënt of “marginaal effect” is. Als je wilt, kun meer lezen over wat p-waarden zijn, enkele van de complicaties die onderzoekers tegenkomen bij het interpreteren ervan, en een discussie over wat s-waarden zijn en hoe ze kunnen helpen in deze blogpost. De s-waarde is geen onderdeel van de leerstof.↩︎Rusland heeft hier een ‘NA’ waarde voor
estimateenresidual_valueomdat het omwille van ontbrekende waarden niet is opgenomen in het regressiemodel.↩︎We zouden technisch gezien wel
newdata = datagrid()kunnen gebruiken maar dan moeten we de niveaus van de factor variabele manueel typen (bv.newdata = datagrid(TYPEDEMO1984_factor = c("Autocracies', "Democracies")). Deby =functie is dus gemakkelijker.↩︎