Binaire logistische regressie
Op deze pagina vind je een demonstratie van een statistische techniek aan de hand van een voorbeeld.
Meer informatie over hoe je deze pagina kan gebruiken vind je in deze handleiding.
De analyse gebeurt met behulp van R en RStudio. Een inleiding tot deze software vind je hier.
Het voorbeeld op deze pagina is afkomstig van een studie van Van Nieuwenhove & De Wever (2023). Er zijn kleine wijzigingen aangebracht om didactische redenen. De studie is ook uitgebreider dan wat hier wordt gedemonstreerd.
1 Doel
Logistische regressie is een statistische techniek die toelaat om een categorische uitkomstvariabele te modelleren in functie van een reeks predictoren. Dat wordt gedemonstreerd op deze pagina. De predictoren kunnen zowel continu als categorisch zijn.
Ook interacties tussen predictoren zijn mogelijk. Dat komt aan bod in een andere demonstratie op Statlas.
Binaire logistische regressie verwijst naar de specifieke situatie waarin de categorische uitkomstvariabele exact twee waarden (“niveaus”) kan aannemen. Er bestaat ook multinomiale logistische regressie, waarbij de uitkomst meer dan twee waarden kan aannemen.
Logistische regressie maakt deel uit van een familie van statistische technieken die gezamenlijk het “Veralgemeend Lineair Model” worden genoemd. (in het Engels: “Generalized Linear Model”)
Het onderzoek van Van Nieuwenhove & De Wever (2023)
In deze demonstratie belichten we een deel van de studie van Van Nieuwenhove & De Wever (2023). Zij proberen te verklaren waarom mensen wel of niet deelnemen aan opleidingen in het volwassenenonderwijs. Hun onderzoek is in belangrijke mate gebaseerd op een courante theorie in de gedragswetenschappen: de theorie van gepland gedrag (afgekort TGG; in het Engels: “theory of planned behavior”). Kort samengevat stelt deze theorie dat gedrag een gevolg is van de intentie om dat gedrag te stellen, waarbij die intentie op zijn beurt het gevolg is van drie grote groepen van oorzaken:
- Perceptie van controle over gedrag (in het Engels: “Perceived behavioral control”)
- Perceptie van sociale normen (in het Engels: “Perceived social norms”)
- Attitudes
Een van de onderzoeksvragen die Van Nieuwenhove & De Wever (2023) zich stelden was of de intentie om wel of niet een opleiding te volgen in het volwassenenonderwijs kan worden verklaard met behulp van deze drie variabelen1. Het is duidelijk dat de intentie om wel of niet deel te nemen aan een opleiding een categorische afhankelijke variabele is met twee niveaus. De bedoeling is om die te modelleren in functie van een aantal predictoren, meer bepaald de drie grote variabelen uit de TGG. Om die onderzoeksvraag met die variabelen te bestuderen gebruikten de onderzoekers binaire logistische regressie. Hoe dat werkt tonen we in de demonstratie hieronder.
Van Nieuwenhove & De Wever (2023) bouwen daarna voort op de TGG met als doel om tekortkomingen in de bestaande wetenschappelijke literatuur over volwassenenonderwijs te remediëren. Dat doen zij door interactie-effecten te toetsen. Een demonstratie van die analyse vind je terug op een andere pagina op Statlas.
Van Nieuwenhove & De Wever (2023) gaan in hun studie nog een stap verder en onderzoeken ook de rol van psychosociale barrières bij de intentie om een opleiding te volgen. Dat is het belangrijkste doel van hun artikel (dat kan je al afleiden uit de titel). Dergelijke barrières zijn onderbelicht in de wetenschappelijke literatuur en dus de moeite waard om te onderzoeken. Specifiek vermoeden ze dat eerdere ervaringen met onderwijs een psychosociale barrière zouden kunnen vormen om als volwassene opnieuw een opleiding te volgen. Voor dit alles gebruiken zij een model met mediatie.
2 Het model
Het doel van de demonstratie op deze pagina is om na te gaan of de
intentie om een opleiding te volgen (intentie) kan worden
verklaard door de drie variabelen van de TGG, namelijk perceptie van
controle over gedrag (pcg), perceptie van sociale normen
(psn) en attitudes (attitudes). Dit is een
eenvoudig logistisch regressiemodel.
In een diagram ziet dit model er zo uit (klik op de afbeelding om te vergroten):
×
De regressievergelijking ziet er zo uit:
\[log\left(\frac{\pi}{1-\pi}\right) = \beta_0 + \beta_1 \: pcg+ \beta_2 \: psn+ \beta_3 \: attitudes\]
\(\pi\) staat daarbij voor de kans op de intentie om een opleiding te volgen. \(\frac{\pi}{1-\pi}\) is de odds op de intentie om een opleiding te volgen. De logaritme van de odds wordt ook de logit genoemd.
Voor een uitgebreidere uitleg over logaritmes, odds en oddsratio’s, kan je in deze syllabus terecht.
3 Dataset & packages
De data die voor de demonstratie worden gebruikt zijn online beschikbaar. Gebruik de volgende code om de data in te laden in RStudio.
barrieres <- read.csv('https://statlas.ugent.be/datasets/barrieres.csv')
De dataset barrieres bevat de gegevens van 563
respondenten die een vragenlijst invulden. Dat resulteerde in 8
variabelen.
Relevante variabelen voor deze demonstratie zijn:
-
intentie“Ik ben van plan om tijdens de volgende 12 maand deel te nemen aan een vorm van levenslang leren.” De antwoordopties “nee” en “ja” werden gecodeerd als 0 en 1. -
pcgPerceptie van controle over het eigen gedrag -
psnPerceptie van sociale normen -
attitudesAttitudes tegenover het gedrag
Om een logistisch regressiemodel te fitten heb je genoeg aan base R.
Er zijn geen packages nodig. Om enkele specifieke zaken te bekomen
zullen we gebruik maken van deze packages: fmsb,
effects, pROC en generalhoslem.
Je kan die eventueel al installeren. Verderop zullen we uitdrukkelijk
vermelden wanneer we die packages gebruiken.
install.packages(c("fmsb", "effects", "pROC", "generalhoslem")) # vergeet niet dat R hoofdlettergevoelig is, dus "pROC" is niet hetzelfde als "proc"
4 Visuele exploratie
Je kan meteen een eerste indicatie krijgen van het antwoord op de
onderzoeksvraag door de associatie te bekijken tussen de
uitkomstvariabele intentie en telkens een andere predictor.
In de code hieronder vertel je aan R dat je boxplots wil van de
predictor pcg, opgesplitst volgens de
intentie. Het argument horizontal=TRUE doet de
horizontale en de verticale as van plaats wisselen. Dat is niet echt
nodig, maar het zorgt ervoor dat de predictor op de horizontale as komt
en de uitkomst op de verticale as, zoals je waarschijnlijk gewoon
bent.
boxplot(pcg ~ intentie, data=barrieres, horizontal=TRUE)
Wat leert deze grafiek je? Bij hele lage waarden voor
pcg zie je enkel maar intentie=0. Wanneer je
geleidelijk naar rechts opschuift, dan kom je eerst voornamelijk
gevallen tegen waar intentie=0. Pas aan de rechterkant,
verschijnen de meeste gevallen waar intentie=1. Hogere
waarden voor pcg lijken dus min of meer geassocieerd met
“hogere waarden” voor intentie (dat wil zeggen
1 in plaats van 0). Op het eerste zicht zou er
dus sprake kunnen zijn van een associatie tussen deze predictor en de
uitkomst.
Op dezelfde manier kan je even goed psn en
attitudes plotten tegenover intentie.
boxplot(psn ~ intentie, data=barrieres, horizontal=TRUE)
boxplot(attitudes ~ intentie, data=barrieres, horizontal=TRUE)
Deze visuele exploratie is natuurlijk maar een hele ruwe analyse. Die moet verfijnd worden door formele toetsen uit te voeren. Dat komt verderop aan bod.
5 Modelspecificatie in R
Je kan het model specifiëren in R door de uitkomstvariabele
intentie te laten volgen door een tilde ~, met
daarna de verschillende predictoren gescheiden door een plusteken
+. Dit is identiek dezelfde werkwijze als bij lineaire
regressiemodellen.
De formule met de modelspecificatie sla je op als een string genaamd
formule1.
formule1 <- "intentie ~ pcg + psn + attitudes"
Die string geef je vervolgens aan de functie glm. Je
vertelt erbij uit welk dataframe de gegevens moeten komen (argument
data) en welk soort analyse je wil uitvoeren (argument
family). Dat laatste is nodig omdat glm ook
andere soorten analyses dan logistische regressie aankan. Je kiest voor
link="logit" omdat je een logaritme van de odds wil
modelleren.2
model1 <- glm(formula=formule1, data=barrieres, family=binomial(link="logit"))
De opsplitsing van de modelspecificatie in twee stappen doen we louter om de code leesbaarder te maken.
6 De output opvragen
Wanneer je het object model1 geeft aan de functie
summary() krijg je uitgebreide informatie terug.
summary(model1)
Call:
glm(formula = formule1, family = binomial(link = "logit"), data = barrieres)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.1605 -0.8501 -0.3922 0.8907 2.6769
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.36288 0.76557 -9.617 < 2e-16 ***
pcg 0.33941 0.11529 2.944 0.00324 **
psn 0.54778 0.07876 6.955 3.53e-12 ***
attitudes 0.50774 0.12637 4.018 5.87e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 726.63 on 538 degrees of freedom
Residual deviance: 559.08 on 535 degrees of freedom
(24 observations deleted due to missingness)
AIC: 567.08
Number of Fisher Scoring iterations: 5
Wat leer je hier nu uit? Je kan de parameterschattingen interpreteren en je ziet voor elke predictor een significantietoets. Vooraleer je in die informatie duikt is het verstandig om eerst het model in zijn geheel te beoordelen.
7 Het model beoordelen
Wat is de verklarende kracht van het model met deze drie predictoren? Er zijn enkele manieren om dit te beoordelen.
Nagelkerke’s \(R^2\)
Een eerste mogelijkheid is om Nagelkerke’s \(R^2\) te berekenen. Voor een logistisch
regressiemodel kan je dit berekenen met het package fmsb.
Je zal het eerst (eenmalig) moeten installeren en vervolgens laden in
R.
install.packages("fmsb") # eenmalig
library(fmsb) # in elke R-sessie
Nu zijn de functies uit het package beschikbaar en kan je
Nagelkerke’s \(R^2\) berekenen. Je
geeft het object model1 aan de functie
NagelkerkeR2.
NagelkerkeR2(model1)$N
[1] 539
$R2
[1] 0.3609118
Het eerste deel van de output $N geeft gewoon de
steekproefgrootte.
Nagelkerke’s \(R^2\) is ongeveer gelijk aan 0.36. Net zoals de determinatiecoëfficiënt \(R^2\) in lineaire regressiemodellen geeft het een idee van de model fit. Het is ook een getal tussen 0 en 1. Maar in tegenstelling tot \(R^2\) mag je Nagelkerke’s \(R^2\) niet interpreteren als de proportie verklaarde variantie van de uitkomstvariabele. Het is een zogenaamde pseudo-\(R^2\) maat. De interpretatie is niet altijd eenduidig. Soms zal je een “verrassend” lage waarde krijgen voor Nagelkerke’s \(R^2\) ondanks een op het eerste zicht goede fit. Hierdoor is het af te raden om deze maat te gebruiken, zeker als dit de enige manier is die je gebruikt om het model te beoordelen.
Modelvergelijking met het nulmodel
Door het model met drie predictoren te vergelijken met het nulmodel krijg je een globale beoordeling van het model. Een modelvergelijking is een toets, dus begin met uitdrukkelijk de nulhypothese en alternatieve hypothese te formuleren.
-
De nulhypothese stelt dat geen enkele predictor een effect heeft op de
logit van
intentie. - De alternatieve hypothese stelt dat minstens één predictor een effect heeft.
Deze modelvergelijking doe je door zowel model1 als het
nulmodel aan de functie anova() te geven. Eerst bouwen we
het nulmodel.
formule0 <- "intentie ~ 1"
model0 <- glm(formula=formule0, data=barrieres, family=binomial(link="logit"))Met deze code is op zich niets mis, maar toch zal je hiermee geen
modelvergelijking kunnen uitvoeren met de functie anova().
Je zal volgende foutmelding krijgen.
anova(model0, model1, test="LRT")Error in anova.glmlist(c(list(object), dotargs), dispersion = dispersion, : models were not all fitted to the same size of dataset
Het probleem is dat het nulmodel en model1 niet geschat
zijn op basis van een even grote dataset. Hoe komt dat?
Het nulmodel maakt enkel gebruik van de kolom intentie.
Bij het schatten van dit model zal R geen rekening houden met de
eventuele rijen/observaties waarbij de waarde voor intentie
ontbreekt. Deze rijen zullen worden geschrapt uit de data.
model1 maakt gebruik van vier kolommen:
pcg, psn, attitudes en
intentie. Bij het schatten van dit model zal R geen
rekening houden met de eventuele rijen waarbij de waarde voor
intentie ontbreekt of waarbij de waarde voor
pcg ontbreekt of waarbij de waarde voor
psn ontbreekt of waarbij de waarde voor
attitudes ontbreekt. Al deze rijen worden geschrapt en dat
zijn er natuurlijk meer dan bij het nulmodel.
Hieronder zie je een deel van de dataset om dit te illustreren. R zal
al deze rijen gebruiken om het nulmodel te schatten. Er is immers overal
een waarde voor intentie. Om model1 te
schatten zal R de rijen 104, 105 en 113 schrappen want daar is geen
waarde voor psn.
X intentie opl.niveau attitudes pcg psn leeftijd geslacht
104 1 2 6.000000 6.25 NA 63 0
105 0 1 4.500000 5.00 NA 53 0
106 0 1 6.000000 5.00 3.333333 62 1
107 1 3 6.833333 5.25 5.000000 37 1
108 0 2 4.500000 5.75 4.000000 28 1
109 1 1 7.000000 6.50 4.333333 55 1
110 1 3 6.833333 6.25 6.000000 28 1
111 1 3 5.500000 4.75 4.000000 48 1
112 0 1 5.333333 6.25 1.000000 63 0
113 0 1 5.333333 4.75 NA 32 0
Om de modelvergelijking te kunnen uitvoeren heb je een nulmodel nodig
dat gebaseerd is op dezelfde rijen/observaties als model1.
Er zijn verschillende manieren om dat te doen. De makkelijkste werkwijze
is waarschijnlijk om eerst model1 te schatten (wat we
eerder op deze pagina al hebben gedaan). Het object model1
bevat een onderdeel model1$model. Als je dit opvraagt, zie
je de data die hiervoor werden gebruikt.
model1$model intentie pcg psn attitudes
1 1 5.75 5.666667 5.833333
2 1 5.25 6.000000 6.500000
3 0 5.50 5.333333 6.666667
4 0 6.25 3.000000 5.000000
5 0 5.25 4.000000 5.166667
6 1 5.50 4.000000 6.500000
7 0 7.00 4.666667 7.000000
9 1 5.25 4.333333 3.666667
10 0 4.75 5.666667 5.500000
11 1 6.00 5.000000 5.666667
12 0 6.00 1.000000 4.833333
13 0 6.00 6.000000 4.666667
14 0 7.00 6.333333 7.000000
15 1 5.75 5.000000 5.666667
16 0 3.50 3.333333 4.000000
17 1 7.00 5.666667 5.833333
18 1 6.25 6.000000 7.000000
19 0 6.50 5.000000 7.000000
20 1 6.25 6.333333 5.333333
21 0 5.50 2.333333 5.500000
22 0 4.25 1.666667 4.666667
23 0 6.75 5.000000 6.666667
24 1 6.25 4.000000 7.000000
25 0 6.25 4.000000 6.000000
26 0 5.25 2.333333 4.833333
28 1 5.75 5.000000 7.000000
29 1 4.75 3.333333 5.000000
30 0 4.50 2.000000 6.000000
31 0 7.00 6.333333 7.000000
32 0 5.00 5.333333 4.666667
33 0 5.75 6.000000 5.500000
34 0 5.50 3.333333 5.666667
35 1 6.75 4.666667 7.000000
36 0 5.25 4.000000 6.166667
37 0 4.50 1.000000 6.333333
38 1 6.50 5.000000 5.833333
39 1 5.75 4.666667 6.833333
40 1 6.00 6.666667 6.666667
41 0 5.25 4.333333 6.333333
42 0 5.00 1.000000 6.166667
43 1 7.00 5.000000 6.500000
44 0 6.50 4.000000 5.666667
45 0 3.50 2.333333 5.000000
46 0 6.25 7.000000 4.333333
47 0 3.75 1.333333 4.666667
48 1 5.00 4.333333 6.333333
49 1 6.50 5.000000 6.666667
50 0 6.50 6.333333 6.666667
51 0 5.00 4.000000 5.833333
52 0 3.25 1.000000 4.000000
53 0 7.00 3.333333 7.000000
54 0 2.75 1.666667 1.833333
55 0 4.00 3.333333 5.666667
56 0 4.00 4.000000 4.000000
57 1 5.75 4.333333 6.166667
58 1 7.00 5.666667 7.000000
59 0 5.25 4.000000 6.333333
60 1 3.75 6.000000 5.000000
61 1 4.50 5.000000 3.000000
62 0 4.00 4.000000 5.000000
63 0 6.75 5.333333 5.000000
64 1 6.75 6.000000 6.500000
65 0 4.25 4.000000 4.833333
66 0 7.00 5.333333 6.333333
67 0 5.75 5.333333 5.333333
68 0 4.50 4.666667 6.166667
69 1 6.00 6.000000 5.500000
70 0 5.25 4.000000 6.000000
71 1 7.00 6.333333 6.500000
72 1 5.50 5.666667 5.666667
73 0 3.00 4.000000 4.000000
74 1 4.25 4.333333 4.500000
75 1 3.75 7.000000 5.500000
76 1 6.75 4.666667 7.000000
77 0 4.75 4.000000 5.666667
78 1 5.00 5.333333 5.000000
79 1 6.50 3.333333 6.166667
80 0 5.50 5.666667 6.500000
81 0 4.25 2.333333 5.333333
82 0 4.00 2.333333 4.666667
83 0 5.25 5.000000 5.833333
84 0 5.00 4.000000 6.666667
85 0 6.25 5.000000 5.500000
86 0 6.50 4.333333 6.500000
87 0 6.50 5.000000 6.000000
88 0 4.75 1.000000 3.500000
89 0 4.25 2.000000 5.166667
90 0 6.00 2.666667 5.833333
91 0 6.75 2.333333 6.166667
92 0 3.75 3.666667 5.333333
93 0 5.75 2.333333 4.833333
94 0 6.25 1.000000 7.000000
95 1 6.00 1.333333 4.666667
96 1 5.50 5.333333 6.166667
97 0 5.00 4.333333 2.333333
98 0 5.25 4.333333 5.166667
99 1 4.25 5.333333 5.833333
100 0 5.00 3.333333 5.666667
101 0 5.50 1.000000 6.333333
102 0 1.00 1.000000 2.000000
103 0 3.50 5.000000 5.166667
106 0 5.00 3.333333 6.000000
107 1 5.25 5.000000 6.833333
108 0 5.75 4.000000 4.500000
109 1 6.50 4.333333 7.000000
110 1 6.25 6.000000 6.833333
111 1 4.75 4.000000 5.500000
112 0 6.25 1.000000 5.333333
114 0 3.00 2.666667 4.500000
115 1 5.50 3.666667 6.666667
116 0 6.00 2.666667 4.833333
117 1 7.00 5.666667 6.500000
118 1 4.50 6.333333 6.333333
119 0 6.50 2.333333 6.333333
120 1 4.75 4.333333 7.000000
121 0 4.25 1.333333 5.666667
122 1 5.00 4.666667 6.833333
123 0 6.50 2.000000 4.000000
124 0 5.25 5.000000 5.833333
125 1 5.75 1.500000 6.166667
126 0 2.00 1.000000 2.333333
127 0 4.00 1.000000 4.000000
128 1 4.00 5.666667 5.500000
129 0 4.00 4.000000 4.166667
130 0 3.50 4.000000 3.500000
131 0 4.25 4.000000 6.000000
132 0 5.75 3.000000 6.000000
133 0 4.75 3.000000 3.166667
134 0 4.00 3.333333 4.333333
135 0 5.50 1.000000 4.500000
136 1 6.50 7.000000 6.500000
137 0 6.50 5.000000 5.500000
138 1 4.75 5.333333 4.166667
139 1 6.00 6.333333 5.666667
140 0 6.00 4.333333 6.833333
142 0 4.00 4.666667 5.000000
143 1 5.75 5.333333 6.333333
144 0 4.25 3.333333 3.833333
146 1 7.00 6.666667 7.000000
147 0 4.50 4.000000 5.500000
148 1 7.00 4.000000 7.000000
149 0 5.75 4.666667 7.000000
150 1 7.00 5.000000 7.000000
151 0 5.00 5.000000 6.166667
152 0 6.00 1.000000 3.666667
153 0 6.50 7.000000 7.000000
154 1 6.75 2.333333 7.000000
155 1 5.50 1.333333 3.666667
156 1 4.25 3.000000 6.500000
157 0 2.50 1.000000 6.500000
158 0 5.50 1.666667 5.833333
159 1 7.00 7.000000 5.000000
160 0 6.50 5.666667 6.333333
161 1 6.00 6.000000 7.000000
162 0 4.25 2.333333 5.000000
163 0 4.50 3.000000 5.333333
165 0 4.00 4.333333 4.666667
166 1 5.75 4.666667 6.833333
168 0 2.50 7.000000 1.000000
169 1 6.00 6.666667 6.833333
170 0 4.00 4.000000 4.000000
171 0 4.00 1.000000 4.166667
172 1 6.25 5.333333 7.000000
173 1 6.25 5.000000 6.333333
174 0 5.00 3.000000 5.000000
175 0 5.25 5.000000 5.000000
176 1 6.25 2.333333 7.000000
177 1 4.75 4.666667 5.333333
178 1 5.75 5.000000 4.666667
179 0 4.75 4.000000 6.166667
180 0 4.00 5.000000 4.833333
181 1 3.75 4.666667 6.666667
182 1 6.00 4.333333 6.000000
183 1 7.00 6.000000 6.333333
184 0 4.50 4.000000 5.500000
185 0 6.75 4.666667 6.166667
187 1 6.75 6.666667 6.833333
188 0 4.75 4.333333 4.500000
189 1 6.00 6.666667 6.333333
190 0 4.25 5.000000 4.500000
191 0 5.25 5.666667 4.333333
192 1 4.75 6.000000 4.833333
193 1 6.00 3.666667 6.833333
194 0 5.00 5.666667 7.000000
195 1 6.25 5.000000 6.500000
196 1 6.75 5.000000 6.000000
197 1 6.50 3.666667 5.500000
198 0 6.50 3.000000 5.500000
199 0 6.25 6.000000 5.833333
200 0 4.75 5.333333 5.000000
201 0 6.00 6.666667 6.500000
202 0 6.50 5.333333 5.333333
203 0 5.00 4.000000 5.500000
205 0 4.00 4.000000 4.000000
206 1 6.50 4.000000 6.333333
207 1 6.00 4.000000 5.666667
208 0 5.25 2.666667 6.166667
209 1 5.50 5.000000 5.500000
210 1 5.75 7.000000 6.833333
211 0 6.50 2.666667 6.500000
212 1 6.75 6.000000 6.833333
213 1 4.75 4.333333 6.500000
214 1 4.50 2.666667 4.833333
215 0 4.50 3.000000 6.666667
216 1 6.00 5.000000 6.666667
217 0 6.25 4.000000 6.166667
218 1 5.75 5.333333 7.000000
219 1 7.00 4.333333 6.833333
220 1 6.00 3.000000 6.333333
221 1 6.75 5.333333 7.000000
222 0 5.25 3.333333 6.333333
223 0 6.00 4.666667 6.000000
224 1 5.75 6.000000 6.500000
225 1 6.75 5.666667 6.833333
226 1 5.75 6.000000 6.666667
228 1 5.75 4.666667 6.500000
229 0 5.50 4.000000 4.833333
230 1 6.25 6.333333 6.000000
231 1 6.75 4.000000 5.833333
232 0 6.50 2.333333 5.166667
233 1 7.00 4.333333 5.333333
234 0 6.25 3.666667 6.000000
235 1 6.75 7.000000 7.000000
236 1 7.00 3.666667 6.333333
237 0 4.75 1.000000 4.166667
238 0 2.50 1.000000 4.000000
239 1 6.50 6.000000 6.500000
240 1 5.25 4.666667 5.666667
241 1 5.00 2.666667 6.666667
242 1 6.75 2.000000 6.833333
243 0 6.75 5.000000 7.000000
244 1 4.75 6.333333 6.166667
245 0 6.75 6.000000 6.500000
246 0 5.50 3.666667 6.166667
247 1 7.00 5.000000 7.000000
248 0 5.75 3.666667 7.000000
249 1 6.50 4.333333 6.000000
250 0 6.25 4.000000 6.666667
251 1 5.25 5.333333 7.000000
252 1 4.75 6.666667 7.000000
253 1 5.25 4.333333 6.666667
254 0 5.25 4.666667 6.166667
255 0 4.25 1.000000 4.833333
256 1 5.00 5.333333 4.833333
257 0 6.75 4.000000 5.666667
258 0 2.25 5.000000 1.833333
259 1 7.00 6.000000 6.666667
260 1 3.00 4.000000 5.000000
261 0 5.50 2.333333 5.166667
262 1 5.25 5.000000 7.000000
263 1 5.50 4.333333 6.666667
264 1 5.00 5.000000 6.000000
265 0 4.75 2.333333 5.000000
266 1 6.00 4.666667 5.333333
267 1 5.50 5.500000 3.666667
268 0 4.75 3.000000 4.833333
269 0 4.00 3.000000 4.000000
270 0 5.50 1.333333 5.000000
271 1 6.75 7.000000 6.333333
272 1 7.00 7.000000 7.000000
273 0 4.00 3.000000 3.500000
274 1 5.00 5.666667 6.333333
276 1 6.75 6.000000 6.000000
277 0 4.75 2.666667 6.000000
278 0 4.75 1.666667 3.333333
279 0 5.75 2.333333 4.666667
280 1 6.75 4.000000 7.000000
281 1 6.00 3.333333 6.000000
282 1 5.75 4.500000 6.166667
283 0 5.25 1.000000 5.666667
284 0 6.25 1.000000 5.333333
286 0 1.25 1.000000 4.000000
287 0 6.25 3.666667 5.000000
288 0 6.50 3.000000 5.500000
289 0 6.00 1.000000 7.000000
290 1 6.25 7.000000 6.500000
291 1 5.25 4.333333 5.833333
292 1 6.50 1.000000 5.666667
293 1 6.75 5.333333 6.833333
294 0 2.75 3.000000 4.000000
295 1 4.75 4.000000 5.833333
296 0 6.00 5.000000 5.666667
297 0 5.00 6.000000 7.000000
298 0 5.50 4.000000 5.500000
299 0 4.75 4.333333 4.166667
300 0 1.25 2.666667 5.500000
301 0 3.25 4.000000 6.000000
302 0 6.00 5.333333 6.666667
303 1 5.50 2.333333 6.166667
304 1 6.50 4.000000 6.500000
305 0 4.75 2.000000 5.000000
306 0 4.25 1.000000 4.500000
307 1 5.75 5.000000 6.833333
308 0 2.50 2.000000 6.500000
309 0 4.00 3.000000 5.166667
310 1 6.75 4.666667 5.000000
311 0 2.50 3.000000 6.500000
312 0 4.75 1.000000 3.833333
313 1 7.00 1.000000 6.666667
314 0 4.00 4.000000 3.333333
315 0 5.50 3.333333 4.000000
316 0 6.00 4.000000 5.333333
317 0 5.00 2.000000 4.666667
318 0 1.00 1.000000 4.166667
319 0 3.50 1.333333 3.833333
321 0 6.50 4.666667 6.166667
322 0 2.50 4.000000 3.833333
323 0 5.25 4.666667 6.333333
324 1 6.50 6.666667 6.833333
325 1 5.50 4.666667 7.000000
326 0 5.25 4.000000 5.333333
327 0 4.25 3.666667 5.500000
328 1 6.25 5.333333 6.333333
329 1 6.00 4.333333 6.500000
330 0 6.00 5.000000 5.500000
331 1 4.75 4.666667 5.833333
332 0 3.50 2.666667 5.833333
333 0 6.75 1.000000 5.333333
334 0 6.00 4.333333 6.000000
335 0 4.50 1.000000 5.833333
336 0 2.50 2.000000 4.000000
337 1 4.75 4.666667 6.500000
338 0 4.25 4.000000 6.166667
339 0 6.25 3.000000 4.000000
340 0 5.50 3.333333 6.333333
341 1 5.75 5.333333 6.333333
342 0 4.75 2.000000 4.000000
343 0 5.00 3.666667 5.666667
345 0 5.50 4.666667 5.500000
346 0 6.50 6.000000 6.833333
347 1 5.75 4.666667 5.666667
348 1 6.50 6.666667 7.000000
349 0 4.00 2.000000 4.833333
350 1 6.00 4.000000 6.166667
351 0 6.00 1.000000 6.500000
352 0 4.00 3.666667 5.833333
353 0 6.25 1.333333 5.833333
354 0 4.75 3.666667 3.666667
355 0 4.75 1.000000 1.000000
356 1 4.75 6.000000 6.833333
357 0 6.25 1.000000 4.000000
358 1 6.50 5.000000 6.166667
359 0 6.00 3.666667 6.166667
360 0 5.50 2.000000 5.333333
361 1 4.50 3.666667 5.166667
362 0 5.25 4.333333 5.333333
363 0 6.50 6.333333 6.666667
364 1 3.75 5.333333 6.666667
365 1 4.00 4.000000 6.000000
366 1 5.00 4.000000 5.000000
367 0 6.50 4.666667 6.000000
368 1 4.00 4.000000 5.000000
369 0 4.00 5.000000 5.833333
370 0 6.25 4.000000 4.500000
371 0 2.75 1.000000 4.000000
372 0 4.75 5.666667 6.166667
373 0 5.75 5.000000 6.666667
374 0 4.25 5.333333 4.500000
375 1 7.00 4.666667 6.500000
376 0 5.25 1.333333 6.166667
377 1 5.50 6.666667 6.500000
378 0 7.00 6.000000 6.166667
379 1 6.25 4.666667 4.666667
380 0 5.75 4.666667 4.833333
381 0 6.75 4.000000 6.333333
382 0 4.00 1.666667 4.833333
383 0 5.75 1.000000 6.000000
384 1 5.25 6.000000 6.833333
385 0 6.50 1.000000 5.166667
386 0 3.50 3.000000 3.666667
387 0 3.75 2.666667 4.333333
388 0 3.00 3.000000 5.166667
389 0 3.50 2.333333 6.333333
390 0 5.75 5.000000 5.666667
391 1 6.75 7.000000 6.500000
392 0 4.25 3.000000 6.000000
393 0 5.75 3.000000 4.666667
394 1 6.25 7.000000 6.500000
395 0 6.25 4.000000 6.166667
396 0 5.75 3.000000 6.000000
397 0 4.75 1.000000 5.000000
398 0 4.75 2.333333 4.333333
399 1 6.50 5.000000 5.000000
400 1 6.50 5.333333 6.166667
401 1 6.00 7.000000 7.000000
402 0 5.25 1.000000 6.000000
403 0 4.50 5.000000 4.833333
404 0 5.00 3.333333 4.833333
405 1 5.00 4.000000 5.666667
406 1 6.75 4.000000 6.666667
407 0 6.25 4.666667 6.666667
408 0 6.25 3.000000 2.500000
409 0 5.50 5.000000 4.166667
411 1 3.25 4.000000 5.833333
412 1 5.50 4.666667 6.000000
413 0 4.75 3.000000 6.166667
414 1 5.50 5.666667 6.500000
415 0 5.50 4.000000 5.833333
416 1 5.50 6.333333 5.333333
417 1 5.50 6.000000 6.166667
418 0 4.00 3.333333 4.500000
419 1 6.50 3.666667 6.333333
420 1 5.75 4.000000 6.833333
421 0 4.25 3.000000 4.333333
422 1 5.75 2.666667 6.666667
423 0 5.00 3.666667 6.166667
424 0 4.00 4.000000 4.000000
425 0 4.75 4.000000 6.000000
426 0 4.75 6.333333 6.000000
427 1 6.00 4.666667 5.500000
428 0 1.75 5.666667 4.333333
429 0 2.75 1.000000 4.000000
430 1 6.50 4.000000 5.000000
431 0 5.75 1.000000 4.000000
432 0 5.50 1.000000 5.500000
433 0 3.00 1.000000 6.500000
434 0 3.00 5.333333 5.833333
435 0 5.00 4.000000 5.500000
436 0 5.50 3.333333 4.000000
437 0 3.25 4.000000 5.166667
438 0 2.50 4.333333 5.000000
439 0 5.00 2.333333 4.500000
441 1 6.25 6.333333 6.333333
442 0 4.50 2.000000 4.000000
443 0 5.25 5.333333 6.333333
444 0 5.75 4.000000 6.833333
445 1 6.75 4.000000 6.833333
446 0 6.25 1.000000 4.000000
447 1 4.25 4.000000 4.833333
448 0 5.00 1.500000 4.000000
450 0 5.50 4.000000 5.000000
451 0 6.75 3.333333 6.833333
452 0 4.25 4.000000 6.000000
453 1 5.25 5.666667 5.833333
454 0 4.25 2.000000 3.833333
455 1 5.75 4.000000 4.666667
456 1 6.00 7.000000 7.000000
457 0 5.50 3.333333 6.333333
458 0 5.75 4.333333 6.000000
459 0 5.00 4.500000 5.333333
461 0 5.00 1.000000 3.666667
463 0 5.50 7.000000 3.500000
464 1 6.00 5.000000 6.666667
465 0 5.50 2.666667 5.333333
466 0 6.25 3.000000 4.500000
467 0 5.25 2.666667 4.500000
468 0 5.75 4.000000 6.000000
469 0 3.75 4.666667 5.500000
470 1 5.75 5.000000 5.166667
471 0 5.75 3.000000 5.166667
472 0 5.25 3.333333 6.000000
473 0 4.25 4.000000 3.833333
474 1 4.75 5.000000 6.166667
475 0 5.75 3.000000 5.333333
476 1 7.00 7.000000 7.000000
477 1 4.50 4.000000 3.833333
478 1 6.00 1.000000 6.833333
479 0 4.25 3.333333 3.333333
480 1 6.25 4.000000 6.666667
482 0 5.25 1.666667 4.666667
483 0 4.00 4.000000 4.000000
484 1 6.25 4.000000 5.500000
485 0 1.00 1.000000 4.000000
486 1 7.00 4.000000 7.000000
487 0 4.25 3.333333 5.500000
488 1 5.00 6.000000 5.833333
489 1 6.25 6.000000 6.000000
490 1 4.50 4.000000 5.500000
491 0 4.75 5.000000 5.000000
492 0 5.50 4.000000 4.000000
494 0 6.00 1.000000 3.500000
495 1 6.25 4.000000 6.333333
496 1 5.25 5.333333 6.833333
497 1 5.75 3.333333 6.666667
498 0 5.50 3.000000 5.333333
499 0 6.00 5.333333 5.500000
500 0 6.00 3.000000 4.833333
501 0 4.00 3.333333 6.000000
502 0 4.25 3.666667 4.833333
503 0 5.75 1.666667 4.666667
504 0 4.50 2.666667 3.833333
505 1 6.00 5.000000 6.666667
506 1 4.50 5.666667 6.166667
507 1 7.00 6.000000 6.833333
508 1 4.50 4.000000 4.000000
509 0 5.50 6.000000 5.500000
510 1 6.00 2.333333 6.000000
511 0 5.75 3.666667 5.666667
512 1 6.50 5.666667 7.000000
513 0 4.25 4.000000 4.000000
514 0 6.25 4.000000 5.833333
515 1 7.00 1.333333 7.000000
516 0 6.00 5.333333 6.000000
517 1 5.25 4.666667 4.166667
518 0 6.25 2.000000 7.000000
519 0 2.75 1.333333 3.333333
520 0 3.25 1.000000 4.000000
521 0 3.25 3.333333 4.666667
522 0 5.25 1.000000 3.666667
523 1 4.00 4.000000 4.500000
524 1 4.25 3.666667 5.000000
526 1 6.25 2.666667 5.000000
527 1 5.25 5.333333 6.000000
528 0 2.75 3.333333 5.000000
529 0 3.50 2.500000 5.166667
530 1 6.00 3.666667 5.000000
531 0 3.75 2.333333 2.666667
532 0 4.50 4.000000 5.166667
533 1 2.75 3.666667 4.166667
534 1 5.00 3.000000 5.166667
535 1 5.25 5.333333 5.833333
536 1 6.00 6.000000 6.000000
537 1 6.25 5.333333 7.000000
538 0 3.75 2.000000 1.666667
539 1 6.50 4.666667 5.500000
540 1 5.75 3.666667 7.000000
541 1 4.25 4.000000 4.833333
542 0 5.25 4.000000 6.500000
543 1 6.75 6.333333 6.666667
544 0 3.50 5.333333 5.333333
545 0 5.00 2.000000 6.000000
546 0 6.25 3.000000 6.500000
547 0 3.75 2.000000 3.833333
548 0 4.00 3.333333 6.166667
549 1 5.00 2.000000 2.000000
550 0 4.25 4.000000 4.333333
551 0 5.50 1.666667 4.666667
552 1 4.25 5.333333 5.500000
553 0 4.50 2.333333 4.000000
554 1 5.00 5.000000 5.833333
555 0 3.50 1.333333 4.166667
556 1 5.00 6.000000 4.833333
557 0 5.00 1.000000 7.000000
558 0 3.75 2.666667 5.833333
559 0 4.00 3.000000 3.166667
560 0 4.50 4.666667 3.833333
561 0 3.75 3.000000 3.000000
562 0 4.25 3.333333 4.666667
563 1 6.50 4.000000 5.000000
Nu kan je die data aan glm() geven in plaats van de hele
dataset barrieres.
model0 <- glm(formula=formule0, data=model1$model, family=binomial(link="logit")) # vergeet niet om het object "formule0" te maken als je dat nog niet hebt gedaan
anova() zal nu de modelvergelijking kunnen uitvoeren
zonder foutmeldingen.
anova(model0, model1, test="LRT")Analysis of Deviance Table
Model 1: intentie ~ 1
Model 2: intentie ~ pcg + psn + attitudes
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 538 726.63
2 535 559.08 3 167.54 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Je ziet een p-waarde voor de \(\chi^2\)-toets die kleiner is dan 0.05, dus
je kan de nulhypothese verwerpen op het 5% significantieniveau. Het
toevoegen van de drie predictoren leidt tot een significante verbetering
van het model. Kortom, de drie predictoren dragen wel degelijk bij tot
een beter begrip van intentie.
Area Under the Curve (AUC)
Nog een optie om het model te beoordelen is het berekenen van de “Area Under the Curve” (AUC). Dit komt niet voor in de opleidingen van de FPPW. Onderaan deze pagina vind je een gedetailleerde uitleg.
Hosmer-Lemeshow toets
Het is ook mogelijk om een formele “goodness-of-fit test” uit te voeren, met name de Hosmer-Lemeshow toets. Ook dit komt niet voor als leerstof in de opleidingen van de FPPW. Onderaan deze pagina vind je een gedetailleerde uitleg.
Antwoord op de onderzoeksvraag
Herinner je wat de onderzoeksvraag van Van Nieuwenhove & De Wever (2023) was: kan de intentie om wel of niet een opleiding te volgen in het volwassenenonderwijs worden verklaard met behulp van de drie variabelen van de theorie van gepland gedrag?
Die vraag kan je beantwoorden door de verklarende kracht van het model in zijn geheel te beoordelen, zoals we net hebben gedaan. De verschillende methoden om het model te beoordelen zijn dus ook methoden om een antwoord te formuleren op de onderzoeksvraag.
8 Parameterschattingen interpreteren
In de output van summary(model1) vind je een kolom
Estimate. Daarin staan parameterschattingen.
In logistische regressie zijn het niet de waarden van de uitkomstvariabele zelf die gemodelleerd worden in functie van predictoren, maar wel de logaritmes van de odds. In dit voorbeeld gaat het om de odds op het hebben van de intentie om een opleiding te volgen. Voor een uitgebreidere uitleg over logaritmes, odds en oddsratio’s, kan je in deze syllabus terecht.
In de kolom Estimate vind je bij elke predictor een
schatting. De interpretatie van deze waarden is niet zo eenvoudig. Deze
kolom met parameterschattingen kan je omrekenen naar odds en
oddsratio’s. De kolom maakt deel uit van het object model1
en kan je afzonderlijk selecteren met model1$coef.
Vervolgens geef je die aan de functie exp(). R berekent dan
\(e^{b_0}\), \(e^{b_1}\) enzovoort. Soms wordt dit ook
geschreven als \(exp(b_0)\), \(exp(b_1)\), enzovoort.
exp(model1$coef) (Intercept) pcg psn attitudes
0.0006343702 1.4041237061 1.7294155104 1.6615293473
De odds en oddsratio’s die je nu bekomt zijn makkelijker te interpreteren.
Onder (Intercept) vind je de geschatte odds voor iemand
met waarde 0 voor alle predictoren.
Onder de namen van de predictoren kan je telkens de geschatte
oddsratio voor die predictor aflezen. Voor pcg bijvoorbeeld
is die ongeveer gelijk aan 1.4. Voor elke toename van de predictor
pcg met 1 eenheid schatten we dus dat de odds op de
intentie om een opleiding te volgen vermenigvuldigd worden met factor
1.4 (bij gelijke waarden van de andere predictoren). Deze geschatte
factor is groter dan 1, dus de odds nemen toe naarmate de waarde voor
pcg toeneemt. Er is met andere woorden een positief
verband.
De andere parameterschattingen kan je natuurlijk analoog interpreteren. Voor meer uitleg bij de interpretatie van de parameterschattingen, zie opnieuw deze syllabus, vooral op pagina 234 en 235.
Een categorische predictor
Alle predictoren tot hiertoe waren continue variabelen. Binaire
logistische regressie kan echter even goed met categorische predictoren.
Louter om dit te kunnen illustreren bouwen we hier een nieuw model. Het
enige verschil met model1 is de extra predictor
opl.niveau. Dat is een categorische variabele met 3
niveaus: 1, 2 en 3, wat staat
voor “laagopgeleid”, “middenopgeleid” en “hoogopgeleid”.
Erg belangrijk hier is dat deze kolom in het data frame
barrieres zeker van het type factor moet zijn.
Anders bestaat het gevaar dat R de niveaus 1,
2, en 3 gaat interpreteren als echte getallen,
wat fout zou zijn.
barrieres$opl.niveau <- factor(barrieres$opl.niveau)
Nu kan je het model specifiëren en de output opvragen. Dat gebeurt op dezelfde manier als voorheen.
formule.extra <- "intentie ~ pcg + psn + attitudes + opl.niveau"
model.extra <- glm(formula=formule.extra, data=barrieres, family=binomial(link="logit"))
summary(model.extra)
Call:
glm(formula = formule.extra, family = binomial(link = "logit"),
data = barrieres)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.2471 -0.7953 -0.3663 0.8649 2.6553
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.88967 0.78838 -8.739 < 2e-16 ***
pcg 0.30604 0.11708 2.614 0.00895 **
psn 0.53003 0.08053 6.582 4.64e-11 ***
attitudes 0.40196 0.13124 3.063 0.00219 **
opl.niveau2 0.04436 0.29561 0.150 0.88072
opl.niveau3 0.81788 0.24947 3.278 0.00104 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 726.63 on 538 degrees of freedom
Residual deviance: 545.63 on 533 degrees of freedom
(24 observations deleted due to missingness)
AIC: 557.63
Number of Fisher Scoring iterations: 5
In de output zie je twee dummyvariabelen opl.niveau2 en
opl.niveau3. Deze dummies vertegenwoordigen de niveaus
2 en 3 van opl.niveau, wat
betekent dat niveau 1 het referentieniveau is. Je kan dit
ook zien met de functie contrasts().
contrasts(barrieres$opl.niveau) 2 3
1 0 0
2 1 0
3 0 1
Wat betekenen de schattingen bij de dummyvariabelen? Net als eerder is het voor de interpretatie makkelijker om de schattingen op de odds-schaal te interpreteren.
exp(model.extra$coef)(Intercept) pcg psn attitudes opl.niveau2 opl.niveau3
0.001018252 1.358032303 1.698987126 1.494744672 1.045356533 2.265682112 Hieruit lezen we af dat de schatting bij opl.niveau2
ongeveer gelijk is aan 1.05. Dit is een oddsratio. De geschatte odds op
intentie is een factor 1.05 groter voor iemand die
“middenopgeleid” is in vergelijking met iemand die “laagopgeleid” (= het
referentieniveau) is, gegeven dat de waarden voor alle andere
predictoren gelijk blijven.
In deze syllabus vind je vanaf p.242 een extra voorbeeld van de interpretatie van parameterschattingen bij een categorische predictor.
Het model.extra diende enkel ter illustratie van een
logistisch regressiemodel met een categorische predictor. In wat volgt
keren we terug naar model1, dus zonder de predictor
opl.niveau.
9 Effecten visualiseren
Verschillende functies uit het package effects laten je
toe om heel eenvoudig plots te creëren van de effecten in een logistisch
regressiemodel.
install.packages("effects") # package eenmalig installeren
library(effects) # package laden
Om meteen plots van alle effecten in het model te bouwen gebruik je
de functie predictorEffects(). Je hoeft enkel het object
model1 aan die functie te geven. Het resultaat geef je aan
de functie plot(). Het resultaat is een plot voor elke
predictor in het model.
effects.model1 <- predictorEffects(model1)
plot(effects.model1)
Op de horizontale as van elke plot verschijnen de verschillende
waarden die de predictor kan aannemen. De verticale as toont de kans
(dus niet de log(odds) of de odds!) op intentie voor elke
waarde van de predictor. Kansen aflezen op de verticale as is niet
noodzakelijk interessant. De plot toont immers de lineaire relatie
tussen de predictor en de uitkomst intentie voor
specifieke vaste waarden van de overige predictoren. Voor andere
waarden van de overige predictoren zal de rechte hoger of lager liggen.
Het interessante aan een dergelijke plot is de helling van de rechte.
Die vertelt je iets over de relatie tussen de predictor en de
uitkomst.
Om de plots correct te lezen moet je wel letten op de assen. De verschillende plots gebruiken een verschillende verticale as. Bovendien is de verticale as herschaald. De streepjes bij 0.1 en 0.2 liggen niet even ver uit elkaar als de streepjes tussen 0.2 en 0.3, enzovoort. Dat komt omdat in een logistisch regressiemodel de kans geen lineaire functie is van de predictoren.
Je kan ook plots maken voor afzonderlijke predictoren. Daarvoor
gebruik je de functie Effect, waarbij je de naam van de
predictor meegeeft samen met het object model1.
effect.pcg <- Effect("pcg", model1)
plot(effect.pcg)
Bij deze laatste plot gelden natuurlijk dezelfde opmerkingen over de verticale as die we zonet al hebben uiteengezet.
10 Toetsen
Toets voor een parameter
Om het effect van een predictor in een logistische regressiemodel te toetsen zijn er twee veelgebruikte methoden: de Wald-toets en de Likelihood Ratio Toets (LRT). Beide toetsen zijn niet equivalent. De voorkeur gaat uit naar de LRT, maar voor grote steekproeven is er weinig verschil.
Wald toets
We kijken opnieuw naar de output van
summary(model1).
summary(model1)
Call:
glm(formula = formule1, family = binomial(link = "logit"), data = barrieres)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.1605 -0.8501 -0.3922 0.8907 2.6769
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.36288 0.76557 -9.617 < 2e-16 ***
pcg 0.33941 0.11529 2.944 0.00324 **
psn 0.54778 0.07876 6.955 3.53e-12 ***
attitudes 0.50774 0.12637 4.018 5.87e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 726.63 on 538 degrees of freedom
Residual deviance: 559.08 on 535 degrees of freedom
(24 observations deleted due to missingness)
AIC: 567.08
Number of Fisher Scoring iterations: 5
Je ziet in de output van summary(model1) een Wald-toets
op elke rij. De toets bij het (Intercept) is inhoudelijk
niet echt interessant. De toetsen bij de predictoren daarentegen zijn
dat wel. De nulhypothese bij die toetsen houdt telkens in dat de
parameter gelijk is aan 0 in de populatie, \(\beta_1 = 0\). Als dat klopt, dan zal de
log(odds) op intentie gemiddeld gezien gelijk blijven bij
elke toename met 1 eenheid van de predictor.
Omdat \(\beta_1\) niet zo makkelijk te interpreteren is, veranderen we de formulering van deze toets best in termen van \(exp(\beta_1)\). De nulhypothese \(\beta_1 = 0\) wordt dan \(exp(\beta_1)=1\). Beide zijn perfect equivalent want \(exp(0) = 1\). Het voordeel van de nieuwe formulering is dat \(exp(\beta_1)\) een oddsratio is, wat makkelijker is voor de interpretatie.
Als de nulhypothese \(exp(\beta_1)=1\) waar is, dan blijven de
odds op intentie gelijk bij elke toename met 1 eenheid van
de predictor pcg, gegeven dat de waarden voor de andere
predictoren gelijk blijven.
Als daarentegen \(exp(\beta_1)\neq1\), dan veranderen de odds
op intentie wel bij elke toename met 1 eenheid van
pcg.
In dit geval zien we in de output bij pcg een p-waarde
van 0.0032409. Dat is kleiner dan 0.05, dus we verwerpen de nulhypothese
op het 5%-significantieniveau.
Likelihood Ratio Toets (LRT)
De LRT voer je uit door twee geneste modellen te construeren. Het ene
model (hier model1) bevat alle predictoren, het andere
bevat dezelfde predictoren behalve de ene die je wil toetsen. Vervolgens
voer je een modelvergelijking uit door beide modellen aan de functie
anova() te geven.
Stel dat je het effect van pcg wil toetsen, dan bouw je
eerst een model zonder die predictor. Om modellen te kunnen vergelijken
met anova() is het belangrijk dat ze gefit zijn op basis
van dezelfde data. Daar kan je voor zorgen door de data van het
complexere model te hergebruiken met model1$model, net
zoals toen je model1 wilde
vergelijken met het nulmodel.
formule.min.pcg <- "intentie ~ psn + attitudes"
model.min.pcg <- glm(formula=formule.min.pcg, data=model1$model, family=binomial(link="logit"))
De nulhypothese van de toets stelt dat beide modellen niet van elkaar verschillen. (Technisch gesproken: “de deviance is hetzelfde”.) De alternatieve hypothese houdt in dat het meer uitgebreide model beter is dan het minder uitgebreide model. (“De deviance is lager”.)
anova(model.min.pcg, model1, test="LRT") # let op de volgorde waarin de modellen aan de functie worden gegevenAnalysis of Deviance Table
Model 1: intentie ~ psn + attitudes
Model 2: intentie ~ pcg + psn + attitudes
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 536 568.13
2 535 559.08 1 9.0508 0.002626 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
De output geeft je een heel kleine p-waarde. Je kan de nulhypothese
verwerpen op het 5% significantieniveau. Je concludeert dat het model
met pcg beter is dan het model zonder pcg.
Merk trouwens op dat de p-waarde licht verschilt van die bij de Wald-toets. Beide toetsen zijn dus inderdaad niet equivalent (maar het maakt in dit geval niet uit voor de conclusie over de nulhypothese).
Toets voor meerdere predictoren
Soms zijn onderzoekers geïnteresseerd in het effect van een verzameling predictoren. Dit kan op een analoge manier getoetst worden, via een LRT waarbij je twee geneste modellen vergelijkt. Het eerste model bevat alle predictoren, het andere model bevat dezelfde predictoren behalve degene die je wil toetsen.
Stel dat je het effect van pcg en attitudes
wil toetsen, dan bouw je eerst een model zonder die predictoren. Opnieuw
geldt dat de datasets van beide modellen die je wil vergelijken even
groot moeten zijn. Dit kan door de “dataset” model1$model
te gebruiken, net als bij de toets voor het volledige model en bij de
toets voor één predictor.
formule.min.pcg.attitudes <- "intentie ~ psn"
model.min.pcg.attitudes <- glm(formula=formule.min.pcg.attitudes, data=model1$model, family=binomial(link="logit"))
De nulhypothese van de toets stelt dat beide modellen niet van elkaar verschillen. (Technisch gesproken: de deviance is hetzelfde.) De alternatieve hypothese houdt in dat het meer uitgebreide model beter is dan het minder uitgebreide model. (De deviance is lager.)
anova(model.min.pcg.attitudes, model1, test="LRT") # let op de volgorde waarin de modellen aan de functie worden gegevenAnalysis of Deviance Table
Model 1: intentie ~ psn
Model 2: intentie ~ pcg + psn + attitudes
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 537 608.21
2 535 559.08 2 49.131 2.144e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
De output geeft je een heel kleine p-waarde. Je kan de nulhypothese
verwerpen op het 5% significantieniveau. Je concludeert dat het model
met pcg en attitudes beter is dan het model
zonder die twee predictoren. De verklarende kracht van het model is
hoger ten gevolge van het toevoegen van de predictoren pcg
en attitudes.
Deze manier van toetsen is trouwens ook bruikbaar om het effect van
een categorische predictor met meer dan twee niveaus te toetsen. Een
dergelijke predictor zal gecodeerd zijn door twee of meer
dummyvariabelen. Je kan een modelvergelijking zoals degene hierboven
gebruiken om een model met die dummyvariabelen te vergelijken met een
model zonder die dummyvariabelen. We illustreren dit voor het extra
model dat eerder al werd opgesteld. We voeren
een toets uit voor de categorische predictor
opl.niveau.
anova(model1, model.extra, test="LRT")Analysis of Deviance Table
Model 1: intentie ~ pcg + psn + attitudes
Model 2: intentie ~ pcg + psn + attitudes + opl.niveau
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 535 559.08
2 533 545.63 2 13.457 0.001196 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1De output geeft je een heel kleine p-waarde. Je kan de nulhypothese verwerpen op het 5% significantieniveau.
Toets voor het hele model
Ten slotte kan je ook het hele model vergelijken met het nulmodel via een LRT. Dat zijn we eerder al tegengekomen, bij de beoordeling van het hele model.
11 Predicties
Dankzij de schattingen die we hebben bekomen in de output van
summary(model1) is het nu mogelijk om predicties te
berekenen. Voor een persoon met bepaalde waarden voor de predictoren kan
je berekenen wat de geschatte waarde is voor de uitkomstvariabele.
Log-schaal
In logistische regressie modelleren we niet intentie (de
uitkomstvariabele) zelf in functie van de predictoren, maar wel de
logaritme van de odds op intentie. Wanneer je waarden voor
de predictoren in de regressievergelijking stopt, dan kom je natuurlijk
waarden voor de logaritme van de odds uit, niet voor de
uitkomstvariabele zelf. Voor meer uitleg kan je opnieuw terecht
in
deze syllabus.
Je kan de logaritme van de odds zelf berekenen door waarden voor de predictoren in te vullen in de regressievergelijking.
\[log\left(\frac{\pi}{1-\pi}\right) = \beta_0 + \beta_1 \: pcg+ \beta_2 \: psn+ \beta_3 \: attitudes\]
Je kan dit makkelijk opvragen in R voor alle observaties in je
dataset ineens. Gebruik daarvoor de functie predict().
predict(model1) 1 2 3 4 5 6 7 9 10
0.654663502 1.006043610 0.810331044 -1.059501175 -0.766508012 -0.004670030 1.123508321 -1.345521221 0.146003933
11 12 13 14 15 16 17 18 19
0.289704779 -2.324544599 0.329749803 2.036480814 0.204851426 -2.318032030 1.078930266 1.599326257 1.136396114
20 21 22 23 24 25 26 28 29
0.935689969 -1.425380995 -2.637952150 1.052003309 0.503759265 -0.003979207 -1.848726663 0.881836055 -1.386026794
30 31 32 33 34 35 36 37 38
-1.693519670 2.036480814 -0.374852606 0.668011843 -0.792974421 1.038654968 -0.258769540 -2.072057009 0.544034563
39 40 41 42 43 44 45 46 47
0.614618478 1.710415744 0.008448037 -1.986973381 1.052233583 -0.088372012 -2.358077054 0.793140495 -2.990253355
48 49 50 51 52 53 54 55 56
-0.076405316 0.967149956 1.697527951 -0.512869050 -3.681046874 0.393130326 -4.585664604 -1.302094538 -1.783136327
57 58 59 60 61 62 63 64 65
0.093531664 1.671291817 -0.174146461 -0.264684215 -1.573384597 -1.275397855 0.388367022 1.515163727 -1.275167581
66 67 68 69 70 71 72 73 74
1.150205004 0.218199767 -0.148140601 0.752865196 -0.343392619 1.782611578 0.485187071 -2.122549738 -1.261819239
75 76 77 78 79 80 81 82 83
0.536968517 1.038654968 -0.682345482 -0.205606448 -0.199691773 0.908302464 -1.934270838 -2.357616506 0.119767799
84 85 86 87 88 89 90 91 92
-0.089753657 0.289935053 0.517337880 0.628657642 -3.425795993 -2.201488416 -0.903833634 -0.662621917 -1.373599549
93 94 95 96 97 98 99 100 101
-1.679019957 -1.139591224 -2.226573179 0.556461808 -2.107359203 -0.583913513 -0.037051114 -0.962681126 -1.732643597
102 103 106 107 108 109 110 111 112
-5.460203993 -0.812697986 -0.793434969 0.627506271 -0.935293621 0.771207116 1.514703178 -0.766968560 -1.985822010
114 115 116 117 118 119 120 121 122
-2.599058497 -0.102641450 -1.411572105 1.417422581 0.849454971 -0.662852191 0.177233646 -2.312808177 0.360058419
123 124 125 126 127 128 129 130 131
-2.030169790 0.119767799 -1.458521575 -4.951544425 -3.426486815 -0.108556125 -1.698513248 -2.206712268 -0.682806030
132 133 134 135 136 137 138 139 140
-0.721469409 -2.499475157 -1.979079167 -2.663497462 1.978093870 0.374788406 -0.713575195 1.020082774 0.516877332
142 143 144 146 147 148 149 150 151
-0.910208857 0.725938239 -2.148095050 2.219075313 -0.851821913 0.758319323 0.699241556 1.306102819 0.204160603
152 153 154 155 156 157 158 159 160
-2.916906150 2.231963106 -0.239506523 -2.904018357 -0.976720290 -2.666260753 -1.621323836 1.386192868 1.163092797
161 162 163 165 166 168 169 170 171
1.514472904 -2.103516996 -1.484228488 -1.262049513 0.614618478 -2.172121371 1.795038823 -1.783136327 -3.341863737
172 173 174 175 176 177 178 179 180
1.234137259 0.713050446 -1.483767940 -0.303347594 -0.409213229 -0.486402642 -0.302887046 -0.428476246 -0.812237437
181 182 183 184 185 187 188 189 190
-0.148831424 0.093761939 1.515394001 -0.851821913 0.615539575 2.049598881 -1.092112533 1.541169587 -0.896630242
191 192 193 194 195 196 197 198 199
-0.276650911 -0.009893883 0.151688334 0.992464994 0.797673525 0.713510995 -0.355589589 -0.720778586 1.006964706
200 201 202 203 205 206 207 208 209
-0.290459801 1.625792666 0.472759826 -0.682115207 -1.783136327 0.250120303 -0.258078717 -0.989147535 0.035374994
210 211 212 213 214 215 216 217 218
1.892779969 -0.395634613 1.684409884 -0.076635590 -1.920692223 -0.807243859 0.797443251 0.080643871 1.064430554
219 220 221 222 223 224 225 226 228
0.856290743 -0.467369899 1.403843965 -0.539335459 0.276356437 1.175750315 1.501815385 1.260373394 0.445372320
229 230 231 232 233 234 235 236 237
-0.850900816 1.274182284 0.081104420 -1.255213741 0.094683036 -0.186573706 2.316816459 0.237232510 -3.087303678
238 239 240 241 242 243 244 245 246
-3.935606932 1.430310374 -0.147449779 -0.820131652 -0.506724101 1.221249466 0.849685245 1.515163727 -0.356510686
247 248 249 250 251 252 253 254 255
1.306102819 0.151458060 0.263468644 0.334513107 0.894723848 1.455395137 0.177694195 0.106419457 -2.918518069
256 257 258 259 260 261 262 263 264
-0.290229527 -0.003518659 -2.929426323 1.684640158 -1.614811266 -1.594627153 0.712129349 0.262547547 0.119537525
265 266 267 268 269 270 271 272 273
-1.933810290 -0.062135877 -0.621587122 -1.653244371 -2.330919823 -2.227033727 1.978324144 2.401669812 -2.584789059
274 276 277 278 279 280 281 282 283
0.653972679 1.261294491 -1.243477319 -3.145230074 -1.763643036 0.673465970 -0.454021558 0.184828914 -2.155989265
284 286 287 288 289 290 291 292 293
-1.985822010 -4.359873697 -0.694312178 -0.720778586 -1.224444577 1.893240517 -0.245421199 -1.731722500 1.319220887
294 295 296 297 298 299 300 301 302
-2.755186587 -0.597722403 0.289704779 1.175059493 -0.512408502 -1.261358691 -2.685293495 -1.022219442 0.980037749
303 304 305 306 307 308 309 310 311
-1.086888681 0.334743382 -2.116404789 -3.087764226 0.797212976 -2.118477257 -1.738558272 0.023178024 -1.570693761
312 313 314 315 316 317 318 319 321
-3.256549835 -1.054277323 -2.121628641 -1.639205207 -0.427324875 -2.200797593 -4.360103971 -3.498222101 0.530686222
322 323 324 325 326 327 328 329 330
-2.376879523 0.191042536 1.964745529 0.614388203 -0.681884933 -1.119269765 0.895644945 0.347631175 0.205081700
331 332 333 334 335 336 337 338 339
-0.232533406 -1.752367162 -1.816115305 0.093761939 -2.325926244 -3.387823436 0.105958909 -0.598182952 -1.567239647
340 341 342 343 345 346 347 348 349
-0.454482106 0.725938239 -2.624143260 -0.780086628 -0.147219504 1.599556531 0.022256927 2.049368607 -2.455587926
350 351 352 353 354 355 356 357 358
-0.004209481 -1.478313813 -1.034876960 -1.549358276 -1.880416924 -4.695142172 1.005583061 -2.662806639 0.713280720
359 360 361 362 363 364 365 366 367
-0.186803980 -1.692598573 -1.203662569 -0.499290435 1.697527951 0.216357574 -0.767659383 -0.935984443 0.446063143
368 369 370 371 372 373 374 375 376
-1.275397855 -0.304498966 -0.765586915 -3.850753580 0.484496248 0.712589898 -0.714035743 0.869639085 -1.719525530
377 378 379 380 381 382 383 384 385
1.456085960 1.430770922 -0.315774839 -0.400858466 0.334973656 -2.638182425 -1.817036402 1.175289767 -1.985591736
386 387 388 389 390 391 392 393 394
-2.669872686 -2.429121517 -2.077971684 -1.681092425 0.204851426 2.062947223 -1.230589526 -1.398454038 1.893240517
395 396 397 398 399 400 401 402 403
0.080643871 -0.721469409 -2.664188285 -2.272302604 0.120919170 0.895875219 2.062256400 -1.986743107 -0.642530732
404 405 406 407 408 409 411 412 413
-1.385796519 -0.597492129 0.504219813 0.699702105 -2.328847355 -0.641609635 -1.106842520 0.106649732 -0.976259742
414 415 416 417 418 419 420 421 422
0.908302464 -0.343162344 0.681129911 0.921650805 -1.894456088 0.067525804 0.249429480 -2.076820313 -0.565571593
423 424 425 426 427 428 429 430 431
-0.526217392 -1.783136327 -0.513099324 0.765062167 0.022487201 -1.464597852 -3.850753580 -0.426864326 -2.832513345
432 433 434 435 436 437 438 439 441
-2.155758990 -2.496554047 -0.461317878 -0.682115207 -1.639205207 -1.445334835 -1.601923473 -2.102826173 1.443428441
442 443 444 445 446 447 448 450 451
-2.708996613 0.556231533 0.249429480 0.588842892 -2.662806639 -1.275167581 -2.813181656 -0.766277738 0.223653894
452 453 454 455 456 457 458 459 461
-0.682806030 0.484956796 -2.878473045 -0.850670542 2.062256400 -0.454482106 0.008908586 -0.492846538 -3.256319561
463 464 465 466 467 468 469 470 471
0.115465043 0.797443251 -1.327409575 -1.313370411 -1.835378321 -0.173685913 -0.741192974 -0.049017810 -1.144584802
472 473 474 475 476 477 478 479 480
-0.708581616 -1.782906052 0.119307250 -1.059961724 2.401669812 -1.698052700 -1.309067655 -2.401964286 0.334513107
482 483 484 485 486 487 488 489 490
-2.298538739 -1.783136327 -0.257848443 -4.444727050 0.758319323 -1.301864264 0.582697942 1.091587785 -0.851821913
491 492 494 495 496 497 498 499 500
-0.473054300 -1.274016209 -3.001529228 0.165266950 0.810100769 -0.200382596 -1.144815077 0.387676199 -1.228977607
501 502 503 504 505 506 507 508 509
-1.132848381 -1.457762079 -2.128832033 -2.428430694 0.797443251 0.399642895 1.769263237 -1.613429621 0.583158491
510 511 512 513 514 515 516 517 518
-1.001805054 -0.525526569 1.501585111 -1.698282974 -0.088602286 -0.702436667 0.641545435 -0.909057486 -0.591807728
519 520 521 522 523 524 526 527 528
-4.006651395 -3.681046874 -2.064393068 -3.171466208 -1.529267091 -1.373139001 -1.242095674 0.386985376 -2.064853617
529 530 531 532 533 534 535 536 537
-2.182156726 -0.779165531 -3.457946802 -1.021068070 -2.305374511 -1.399144861 0.302362298 1.006734432 1.234137259
538 539 540 541 542 543 544 545 546
-4.148279773 0.192193907 0.151458060 -1.275167581 -0.089523383 1.782381304 -0.545480408 -1.523812964 -0.297893468
547 548 549 550 551 552 553 554 555
-3.048179751 -1.048225302 -3.554766851 -1.529036817 -2.213685386 -0.206297271 -2.526402115 0.034914446 -3.328975944
556 557 558 559 560 561 562 563
0.074959470 -1.563857988 -1.667513809 -2.754035216 -1.332863702 -2.923511648 -1.724979657 -0.426864326
Bijvoorbeeld, voor de 54e observatie in de dataset kan je aflezen dat \(log\left(\frac{\pi}{1-\pi}\right) = -4.5856646\). Als je dat wil kan je altijd narekenen dat dit inderdaad de juiste geschatte log(odds) is voor een persoon met volgende waarden voor de predictoren:
attitudes pcg psn
1.833333 2.75 1.666667
Inderdaad,
\[ \begin{aligned} b_0 + b_1 \: pcg_{54}+ b_2 \: psn_{54} + b_3 \: attitudes_{54} &= -7.36288 + 0.33941 \times 2.75 + 0.54778 \times 1.6666667 + 0.50774 \times 1.8333333 \\ &\approx -4.5856646 \end{aligned} \]
De log(odds) is niet zo interessant om te rapporteren. Je rekent predicties op de log(odds)-schaal best om naar odds of naar kansen. Hieronder tonen we hoe je dat kan doen.
Odds-schaal
Predicties op de odds-schaal kan je bekomen door het getal \(e\) te verheffen tot de macht \(\beta_0 + \beta_1 \: pcg+ \beta_2 \: psn+ \beta_3 \: attitudes\).
Opnieuw kan je dit laten berekenen voor alle observaties in je
dataset. Het commando predict(model1) leverde log(odds) op.
Om hieruit de odds te berekenen in R gebruik je de functie
exp().
log.odds <- predict(model1)
exp(log.odds) 1 2 3 4 5 6 7 9 10
1.924494820 2.734759806 2.248652265 0.346628674 0.464632731 0.995340858 3.075625578 0.260403945 1.157200740
11 12 13 14 15 16 17 18 19
1.336033005 0.097827985 1.390620156 7.663592084 1.227342700 0.098467176 2.941531213 4.949696477 3.115520128
20 21 22 23 24 25 26 28 29
2.548971563 0.240416848 0.071507556 2.863381615 1.654930915 0.996028699 0.157437510 2.415330317 0.250066901
30 31 32 33 34 35 36 37 38
0.183871218 7.663592084 0.687390589 1.950355851 0.452496876 2.825414184 0.772000918 0.125926483 1.722944185
39 40 41 42 43 44 45 46 47
1.848951049 5.531260589 1.008483823 0.137109777 2.864041054 0.915420266 0.094601963 2.210327059 0.050274698
48 49 50 51 52 53 54 55 56
0.926440630 2.630436906 5.460432239 0.598775196 0.025196583 1.481611469 0.010196971 0.271961562 0.168110071
57 58 59 60 61 62 63 64 65
1.098045372 5.319034579 0.840173836 0.767448260 0.207342224 0.279319817 1.474570884 4.550166050 0.279384145
66 67 68 69 70 71 72 73 74
3.158840417 1.243835521 0.862309865 2.123074335 0.709359649 5.945362943 1.624478873 0.119725969 0.283138461
75 76 77 78 79 80 81 82 83
1.710812693 2.825414184 0.505430124 0.814153425 0.818983147 2.480108883 0.144529616 0.094645542 1.127235076
84 85 86 87 88 89 90 91 92
0.914156354 1.336340694 1.677555846 1.875091844 0.032523382 0.110638360 0.405014004 0.515497968 0.253193933
93 94 95 96 97 98 99 100 101
0.186556720 0.319949783 0.107897543 1.744489230 0.121558554 0.557711479 0.963626879 0.381867677 0.176816360
102 103 106 107 108 109 110 111 112
0.004252688 0.443659463 0.452288527 1.872934160 0.392470611 2.162374917 4.548070960 0.464418794 0.137267732
114 115 116 117 118 119 120 121 122
0.074343540 0.902450489 0.243759766 4.126471077 2.338372024 0.515379276 1.193910013 0.098982899 1.433413151
123 124 125 126 127 128 129 130 131
0.131313224 1.127235076 0.232579873 0.007072478 0.032500922 0.897128542 0.182955331 0.110061908 0.505197402
132 133 134 135 136 137 138 139 140
0.486037543 0.082128092 0.138196435 0.069704008 7.228950525 1.454683580 0.489889612 2.773424321 1.676783428
142 143 144 146 147 148 149 150 151
0.402440163 2.066669221 0.116706266 9.198820902 0.426636928 2.134685487 2.012225968 3.691758192 1.226495117
152 153 154 155 156 157 158 159 160
0.054100808 9.318140634 0.787016139 0.054802561 0.376544031 0.069511661 0.197636887 3.999594048 3.199814363
161 162 163 165 166 168 169 170 171
4.547023777 0.122026506 0.226677157 0.283073269 1.848951049 0.113935661 6.019708420 0.168110071 0.035370974
172 173 174 175 176 177 178 179 180
3.435413372 2.040205313 0.226781577 0.738342408 0.664172596 0.614834200 0.738682529 0.651501066 0.443863837
181 182 183 184 185 187 188 189 190
0.861714367 1.098298252 4.551213956 0.426636928 1.850654897 7.764785882 0.335506975 4.670049107 0.407942012
191 192 193 194 195 196 197 198 199
0.758319171 0.990154901 1.163797464 2.697876532 2.220369281 2.041145142 0.700760161 0.486373425 2.737279945
200 201 202 203 205 206 207 208 209
0.747919594 5.082446121 1.604415998 0.505546525 0.168110071 1.284179898 0.772534418 0.371893582 1.036008133
210 211 212 213 214 215 216 217 218
6.637795921 0.673252649 5.389269700 0.926227319 0.146505513 0.446085849 2.219858046 1.083984790 2.899187582
219 220 221 222 223 224 225 226 228
2.354411361 0.626648251 4.070818013 0.583135641 1.318317678 3.240573484 4.489832452 3.526738104 1.561071306
229 230 231 232 233 234 235 236 237
0.427030083 3.575776244 1.084484132 0.285014922 1.099310358 0.829797399 10.143331138 1.267735845 0.045624808
238 239 240 241 242 243 244 245 246
0.019533840 4.179996352 0.862905774 0.440373675 0.602465968 3.391422556 2.338910553 4.550166050 0.700114990
247 248 249 250 251 252 253 254 255
3.691758192 1.163529502 1.301436487 1.397259904 2.446660046 4.286176759 1.194459993 1.112288336 0.054013673
256 257 258 259 260 261 262 263 264
0.748091841 0.996487524 0.053427680 5.390510853 0.198928211 0.202984197 2.038326951 1.300238290 1.126975533
265 266 267 268 269 270 271 272 273
0.144596194 0.939755187 0.537091331 0.191427837 0.097206293 0.107847862 7.230615358 11.041598386 0.075411985
274 276 277 278 279 280 281 282 283
1.923165794 3.529988068 0.288379684 0.043057017 0.171419239 1.961022401 0.635069043 1.203012604 0.115788588
284 286 287 288 289 290 291 292 293
0.137267732 0.012780002 0.499417841 0.486373425 0.293920905 6.640853652 0.782374934 0.176979300 3.740505963
294 295 296 297 298 299 300 301 302
0.063597152 0.550063032 1.336033005 3.238335595 0.599051024 0.283268890 0.068201174 0.359795508 2.664556825
303 304 305 306 307 308 309 310 311
0.337264200 1.397581694 0.120463944 0.045603800 2.219346929 0.120214545 0.175773636 1.023448722 0.207900899
312 313 314 315 316 317 318 319 321
0.038521073 0.348444149 0.119836299 0.194134278 0.652251618 0.110714818 0.012777059 0.030251119 1.700098553
322 323 324 325 326 327 328 329 330
0.092839831 1.210510941 7.133097190 1.848525332 0.505662953 0.326518143 2.448914695 1.415710004 1.227625358
331 332 333 334 335 336 337 338 339
0.792523274 0.173363079 0.162656396 1.098298252 0.097692915 0.033782126 1.111776191 0.549809760 0.208620254
340 341 342 343 345 346 347 348 349
0.634776630 2.066669221 0.072501846 0.458366302 0.863104501 4.950836396 1.022506460 7.762998058 0.085812729
350 351 352 353 354 355 356 357 358
0.995799366 0.228021852 0.355270091 0.212384222 0.152526501 0.009139568 2.733500606 0.069752178 2.040675173
359 360 361 362 363 364 365 366 367
0.829606340 0.184040659 0.300093085 0.606961186 5.460432239 1.241546244 0.464098074 0.392199577 1.562150102
368 369 370 371 372 373 374 375 376
0.279319817 0.737492791 0.465060900 0.021263707 1.623357033 2.039265915 0.489664046 2.386049536 0.179151130
377 378 379 380 381 382 383 384 385
4.289138770 4.181921887 0.729223624 0.669744846 1.397903558 0.071491092 0.162506642 3.239081386 0.137299345
386 387 388 389 390 391 392 393 394
0.069261043 0.088114205 0.125183867 0.186170488 1.227342700 7.869127742 0.292120314 0.246978488 6.640853652
395 396 397 398 399 400 401 402 403
1.083984790 0.486037543 0.069655871 0.103074567 1.128533689 2.449478682 7.863693447 0.137141354 0.525959676
404 405 406 407 408 409 411 412 413
0.250124491 0.550189712 1.655693266 2.013152909 0.097407959 0.526444359 0.330601181 1.112544497 0.376717488
414 415 416 417 418 419 420 421 422
2.480108883 0.709523015 1.976109298 2.513436160 0.150400117 1.069857866 1.283293064 0.125328083 0.568035369
423 424 425 426 427 428 429 430 431
0.590835648 0.168110071 0.598637329 2.149127975 1.022741944 0.231170938 0.021263707 0.652552080 0.058864720
432 433 434 435 436 437 438 439 441
0.115815255 0.082368348 0.630452238 0.505546525 0.194134278 0.235667154 0.201508549 0.122110834 4.235191058
442 443 444 445 446 447 448 450 451
0.066603602 1.744087565 1.283293064 1.801902213 0.069752178 0.279384145 0.060013745 0.464739736 1.250638091
452 453 454 455 456 457 458 459 461
0.505197402 1.624104840 0.056220544 0.427128428 7.863693447 0.634776630 1.008948385 0.610885009 0.038529944
463 464 465 466 467 468 469 470 471
1.122395278 2.219858046 0.265163258 0.268912183 0.159553128 0.840560865 0.476545070 0.952164172 0.318356071
472 473 474 475 476 477 478 479 480
0.492342032 0.168148787 1.126716049 0.346469072 11.041598386 0.183039610 0.270071739 0.090539932 1.397259904
482 483 484 485 486 487 488 489 490
0.100405455 0.168110071 0.772712333 0.011740310 2.134685487 0.272024195 1.790863565 2.979000332 0.426636928
491 492 494 495 496 497 498 499 500
0.623096236 0.279706005 0.049710991 1.179708000 2.248134518 0.818417570 0.318282770 1.473552569 0.292591568
501 502 503 504 505 506 507 508 509
0.322114444 0.232756583 0.118976173 0.088175098 2.219858046 1.491292055 5.866529526 0.199203250 1.791688535
510 511 512 513 514 515 516 517 518
0.367215998 0.591243952 4.488798679 0.182997466 0.915209493 0.495376764 1.899414033 0.402903787 0.553326117
519 520 521 522 523 524 526 527 528
0.018194219 0.025196583 0.126895284 0.041942057 0.216694426 0.253310568 0.288778398 1.472534957 0.126836856
529 530 531 532 533 534 535 536 537
0.112797994 0.458788697 0.031494360 0.360210005 0.099721447 0.246807929 1.353051345 2.736649692 3.435413372
538 539 540 541 542 543 544 545 546
0.015791558 1.211905491 1.163529502 0.279384145 0.914366884 5.943994036 0.579563289 0.217879534 0.742380423
547 548 549 550 551 552 553 554 555
0.047445208 0.350559334 0.028588039 0.216744331 0.109297103 0.813591183 0.079946141 1.035531111 0.035829778
556 557 558 559 560 561 562 563
1.077840465 0.209326931 0.188715666 0.063670418 0.263720961 0.053744623 0.178176676 0.652552080
Als je dat wil kan je dit opnieuw zelf narekenen voor een bepaalde observatie, bijvoorbeeld nummer 54.
\[\frac{\pi}{1-\pi} = e^{-4.5856646} \approx 0.010197\]
Kansen
Als je eenmaal de odds hebt, dan is het ook mogelijk om de voorspelde kansen \(\pi\) te berekenen.
\[ \begin{aligned} odds &= \frac{\pi}{1-\pi} \\ & \Updownarrow \\ \pi &= \frac{odds}{1+odds} \end{aligned} \]
In R kan je op basis van die formule een object (hier
kansen) maken dat voor elke observatie de voorspelde kans
bevat.
kansen <- odds/(1+odds)
Je kan de kansen ook veel sneller bekomen via de functie
predict(), door een extra argument mee te geven, of nog
korter via de functie fitted().
kansen <- predict(model1, type="response")
# of
kansen <- fitted(model1)
De 54e observatie scoorde vrij laag op de drie TGG-predictoren. Nu
zie je in het object kans dat voor deze persoon de
verwachte kans op intentie erg laag is. Eigenlijk kon je
dit ook al aflezen aan de voorspelde odds en zelfs aan de log(odds). De
lage voorspelde kans is niet zo verwonderlijk. Uit de
parameterschattingen van het model wist je immers al dat er een positief
verband is tussen elke predictor en intentie. Bij lage
waarden van de predictoren zal de voorspelde kans dan ook eerder klein
zijn.
12 Extra informatie
In deze sectie vind je meer gedetailleerde uitleg bij twee technieken die nuttig zijn in de context van logistische regressie, maar die geen deel uitmaken van de leerstof in opleidingen van de FPPW. We verwezen hier al naar in deze subsectie.
Area Under the Curve (AUC)
De AUC is een maat die de voorspellende kracht van een logistisch
regressiemodel uitdrukt. Om deze maat te berekenen zal je het package
pROC nodig hebben.
install.packages("pROC") # eenmalig
library(pROC) # in elke R-sessie
Logistische regressie probeert op basis van waarden voor predictoren
om gevallen zo goed mogelijk te klasseren als “succes” of “mislukking”3. In deze
demonstratie probeerde je te voorspellen of iemand de intentie heeft om
een opleiding te volgen door informatie over pcg,
psn en attitudes van die persoon te gebruiken.
Als je model extreem goed is, dan is elke voorspelling van een succes
ook in werkelijkheid een succes, en is elke voorspelling van een
mislukking daadwerkelijk een mislukking. Het zou nuttig zijn om in een
getal uit te drukken hoe goed de voorspellingen op basis van je model
aansluiten bij de geobserveerde intentie. Alleen blijkt dat
niet zo eenvoudig.
Je model geeft immers niet vanzelf voorspellingen in de vorm van
“succes” en “mislukking” waarmee je de geobserveerde
intentie kan vergelijken. In plaats daarvan krijg je voor
een bepaalde observatie een geschat logaritme van de odds, waaruit je
vervolgens een geschatte odds en een geschatte kans kan berekenen. Een
meer uitgebreide uitleg vind je wat hoger op deze pagina, onder
Predicties. In R kan je de voorspelde kansen
makkelijk laten berekenen met de functie fitted().
kansen <- fitted(model1)
Eenmaal je de kansen hebt kan je overgaan naar een uitdrukkelijke voorspelling van “succes” of “mislukking” door een cutoff te kiezen, bv. 0.50. Als de kans voor een bepaalde persoon in de dataset kleiner of gelijk is aan 0.50, dan is je voorspelling dat die persoon niet de intentie heeft om een opleiding te volgen (“mislukking”). Als de kans groter is dan 0.50, dan voorspel je dat die persoon wel de intentie heeft (“succes”).
In R kan je die uitdrukkelijke voorspelling van “succes” of
“mislukking” maken met onderstaande lijn code. In plaats van de woorden
succes en mislukking te gebruiken, coderen we
dit als een 1 en een 0.
# We creëren een vector genaamd "intentie.voorspeld". Daarin komt een waarde 1 voor bij elke observatie die voldoet aan de "test" en een waarde 0 bij elke observatie die niet voldoet aan de "test".
intentie.voorspeld <- ifelse(test=kansen > 0.50, yes=1, no=0)
Om de voorspelde kansen en intenties naast elkaar te zien kan je ze samen in een data frame stoppen. En inderdaad, in dat data frame zie je dat de voorspelde intentie gelijk is aan 1, overal waar de kans groter is dan 0.50.
data.frame(kansen, intentie.voorspeld) kansen intentie.voorspeld
1 0.658060601 1
2 0.732245164 1
3 0.692180043 1
4 0.257404792 0
5 0.317234977 0
6 0.498832495 0
7 0.754638894 1
9 0.206603562 0
10 0.536436280 1
11 0.571923856 1
12 0.089110486 0
13 0.581698499 1
14 0.884574436 1
15 0.551034513 1
16 0.089640526 0
17 0.746291493 1
18 0.831924199 1
19 0.757017347 1
20 0.718228230 1
21 0.193819399 0
22 0.066735466 0
23 0.741159404 1
24 0.623342365 1
25 0.499005200 0
26 0.136022471 0
28 0.707202552 1
29 0.200042814 0
30 0.155313530 0
31 0.884574436 1
32 0.407368983 0
33 0.661057835 1
34 0.311530361 0
35 0.738590398 1
36 0.435666207 0
37 0.111842545 0
38 0.632750460 1
39 0.648993618 1
40 0.846890200 1
41 0.502111997 1
42 0.120577432 0
43 0.741203578 1
44 0.477921364 0
45 0.086425903 0
46 0.688505258 1
47 0.047868141 0
48 0.480907958 0
49 0.724551059 1
50 0.845211595 1
51 0.374521194 0
52 0.024577319 0
53 0.597036034 1
54 0.010094042 0
55 0.213812721 0
56 0.143916293 0
57 0.523365885 1
58 0.841747978 1
59 0.456573080 0
60 0.434212575 0
61 0.171734426 0
62 0.218334629 0
63 0.595889531 1
64 0.819825210 1
65 0.218373931 0
66 0.759548360 1
67 0.554334535 1
68 0.463032431 0
69 0.679802690 1
70 0.414985606 0
71 0.856019043 1
72 0.618971976 1
73 0.106924348 0
74 0.220660879 0
75 0.631106936 1
76 0.738590398 1
77 0.335738017 0
78 0.448778705 0
79 0.450242295 0
80 0.712652669 1
81 0.126278616 0
82 0.086462273 0
83 0.529906210 1
84 0.477576637 0
85 0.571980233 1
86 0.626525063 1
87 0.652185024 1
88 0.031498931 0
89 0.099616909 0
90 0.288263322 0
91 0.340150881 0
92 0.202038908 0
93 0.157225286 0
94 0.242395421 0
95 0.097389460 0
96 0.635633476 1
97 0.108383601 0
98 0.358032592 0
99 0.490738281 0
100 0.276341710 0
101 0.150249747 0
102 0.004234679 0
103 0.307315869 0
106 0.311431591 0
107 0.651923802 1
108 0.281851989 0
109 0.683781959 1
110 0.819757172 1
111 0.317135232 0
112 0.120699575 0
114 0.069199038 0
115 0.474362142 0
116 0.195986212 0
117 0.804934040 1
118 0.700452798 1
119 0.340099198 0
120 0.544192791 1
121 0.090067734 0
122 0.589054576 1
123 0.116071501 0
124 0.529906210 1
125 0.188693551 0
126 0.007022809 0
127 0.031477863 0
128 0.472887589 0
129 0.154659543 0
130 0.099149343 0
131 0.335635314 0
132 0.327069491 0
133 0.075894982 0
134 0.121417033 0
135 0.065161958 0
136 0.878477821 1
137 0.592615517 1
138 0.328809335 0
139 0.734988722 1
140 0.626417293 1
142 0.286957100 0
143 0.673913315 1
144 0.104509368 0
146 0.901949450 1
147 0.299050809 0
148 0.680988729 1
149 0.668019594 1
150 0.786860286 1
151 0.550863601 1
152 0.051324131 0
153 0.903083314 1
154 0.440407964 0
155 0.051955278 0
156 0.273543034 0
157 0.064993832 0
158 0.165022378 0
159 0.799983761 1
160 0.761894238 1
161 0.819723145 1
162 0.108755458 0
163 0.184789581 0
165 0.220621282 0
166 0.648993618 1
168 0.102282084 0
169 0.857543941 1
170 0.143916293 0
171 0.034162609 0
172 0.774541871 1
173 0.671074846 1
174 0.184858969 0
175 0.424739341 0
176 0.399100789 0
177 0.380741379 0
178 0.424851873 0
179 0.394490249 0
180 0.307413916 0
181 0.462860674 0
182 0.523423327 1
183 0.819859222 1
184 0.299050809 0
185 0.649203416 1
187 0.885907082 1
188 0.251220683 0
189 0.823634684 1
190 0.289743476 0
191 0.431275040 0
192 0.497526550 0
193 0.537849537 1
194 0.729574530 1
195 0.689476606 1
196 0.671176497 1
197 0.412027620 0
198 0.327221556 0
199 0.732425717 1
200 0.427891304 0
201 0.835592461 1
202 0.616036762 1
203 0.335789374 0
205 0.143916293 0
206 0.562206111 1
207 0.435836061 0
208 0.271080488 0
209 0.508842826 1
210 0.869072176 1
211 0.402361621 0
212 0.843487590 1
213 0.480850474 0
214 0.127784394 0
215 0.308478123 0
216 0.689427302 1
217 0.520150049 1
218 0.743536319 1
219 0.701885102 1
220 0.385238942 0
221 0.802793159 1
222 0.368342185 0
223 0.568652731 1
224 0.764182839 1
225 0.817845078 1
226 0.779090379 1
228 0.609538400 1
229 0.299243925 0
230 0.781457845 1
231 0.520264998 1
232 0.221798920 0
233 0.523653091 1
234 0.453491408 0
235 0.910260228 1
236 0.559031533 1
237 0.043634014 0
238 0.019159580 0
239 0.806949671 1
240 0.463204197 0
241 0.305735715 0
242 0.375961786 0
243 0.772283358 1
244 0.700501111 1
245 0.819825210 1
246 0.411804492 0
247 0.786860286 1
248 0.537792298 1
249 0.565488769 1
250 0.582857078 1
251 0.709864046 1
252 0.810827362 1
253 0.544307026 1
254 0.526579784 1
255 0.051245704 0
256 0.427947676 0
257 0.499120336 0
258 0.050717938 0
259 0.843517987 1
260 0.165921704 0
261 0.168733885 0
262 0.670871497 1
263 0.565262432 1
264 0.529848847 1
265 0.126329438 0
266 0.484471027 0
267 0.349420571 0
268 0.160670946 0
269 0.088594364 0
270 0.097348983 0
271 0.878502402 1
272 0.916954546 1
273 0.070123810 0
274 0.657905138 1
276 0.779248867 1
277 0.223831288 0
278 0.041279639 0
279 0.146334662 0
280 0.662278813 1
281 0.388405031 0
282 0.546076133 1
283 0.103772874 0
284 0.120699575 0
286 0.012618734 0
287 0.333074495 0
288 0.327221556 0
289 0.227155234 0
290 0.869124571 1
291 0.438950817 0
292 0.150367386 0
293 0.789052053 1
294 0.059794399 0
295 0.354864945 0
296 0.571923856 1
297 0.764058325 1
298 0.374629086 0
299 0.220740090 0
300 0.063846751 0
301 0.264595306 0
302 0.727115707 1
303 0.252204613 0
304 0.582913065 1
305 0.107512558 0
306 0.043614800 0
307 0.689377994 1
308 0.107313858 0
309 0.149496153 0
310 0.505794247 1
311 0.172117513 0
312 0.037092240 0
313 0.258404584 0
314 0.107012337 0
315 0.162573239 0
316 0.394765307 0
317 0.099678888 0
318 0.012615865 0
319 0.029362860 0
321 0.629643148 1
322 0.084952825 0
323 0.547615901 1
324 0.877045610 1
325 0.648941160 1
326 0.335840735 0
327 0.246146760 0
328 0.710053716 1
329 0.586043027 1
330 0.551091481 1
331 0.442127188 0
332 0.147748878 0
333 0.139900659 0
334 0.523423327 1
335 0.088998401 0
336 0.032678187 0
337 0.526464971 1
338 0.354759516 0
339 0.172610258 0
340 0.388295635 0
341 0.673913315 1
342 0.067600672 0
343 0.314301216 0
345 0.463261455 0
346 0.831956395 1
347 0.505564002 1
348 0.885883804 1
349 0.079030874 0
350 0.498947631 0
351 0.185682243 0
352 0.262139697 0
353 0.175178972 0
354 0.132340992 0
355 0.009056793 0
356 0.732154858 1
357 0.065204053 0
358 0.671125673 1
359 0.453434338 0
360 0.155434408 0
361 0.230824307 0
362 0.377707434 0
363 0.845211595 1
364 0.553879380 1
365 0.316985646 0
366 0.281712179 0
367 0.609702804 1
368 0.218334629 0
369 0.424458044 0
370 0.317434517 0
371 0.020820975 0
372 0.618809035 1
373 0.670973180 1
374 0.328707703 0
375 0.704670593 1
376 0.151932288 0
377 0.810933302 1
378 0.807021406 1
379 0.421705796 0
380 0.401106102 0
381 0.582969050 1
382 0.066721126 0
383 0.139789861 0
384 0.764099835 1
385 0.120724017 0
386 0.064774681 0
387 0.080978821 0
388 0.111256365 0
389 0.156950868 0
390 0.551034513 1
391 0.887249341 1
392 0.226078261 0
393 0.198061547 0
394 0.869124571 1
395 0.520150049 1
396 0.327069491 0
397 0.065119889 0
398 0.093442973 0
399 0.530193013 1
400 0.710101122 1
401 0.887180214 1
402 0.120601852 0
403 0.344674688 0
404 0.200079667 0
405 0.354917665 0
406 0.623450489 1
407 0.668121722 1
408 0.088761849 0
409 0.344882770 0
411 0.248460009 0
412 0.526637190 1
413 0.273634563 0
414 0.712652669 1
415 0.415041511 0
416 0.663990835 1
417 0.715378349 1
418 0.130737223 0
419 0.516875039 1
420 0.562036072 1
421 0.111370262 0
422 0.362259283 0
423 0.371399553 0
424 0.143916293 0
425 0.374467253 0
426 0.682451775 1
427 0.505621563 1
428 0.187765103 0
429 0.020820975 0
430 0.394875350 0
431 0.055592295 0
432 0.103794292 0
433 0.076100107 0
434 0.386673233 0
435 0.335789374 0
436 0.162573239 0
437 0.190720578 0
438 0.167712955 0
439 0.108822435 0
441 0.808985004 1
442 0.062444569 0
443 0.635580142 1
444 0.562036072 1
445 0.643099607 1
446 0.065204053 0
447 0.218373931 0
448 0.056616007 0
450 0.317284856 0
451 0.555681563 1
452 0.335635314 0
453 0.618917665 1
454 0.053228035 0
455 0.299292215 0
456 0.887180214 1
457 0.388295635 0
458 0.502227132 1
459 0.379223226 0
461 0.037100466 0
463 0.528834233 1
464 0.689427302 1
465 0.209588175 0
466 0.211923399 0
467 0.137598808 0
468 0.456687351 0
469 0.322743328 0
470 0.487748001 0
471 0.241479581 0
472 0.329912327 0
473 0.143944666 0
474 0.529791483 1
475 0.257316769 0
476 0.916954546 1
477 0.154719765 0
478 0.212642901 0
479 0.083023032 0
480 0.582857078 1
482 0.091244054 0
483 0.143916293 0
484 0.435892682 0
485 0.011604075 0
486 0.680988729 1
487 0.213851432 0
488 0.641687966 1
489 0.748680594 1
490 0.299050809 0
491 0.383893587 0
492 0.218570518 0
494 0.047356836 0
495 0.541222953 1
496 0.692130977 1
497 0.450071306 0
498 0.241437404 0
499 0.595723167 1
500 0.226360419 0
501 0.243635826 0
502 0.188809848 0
503 0.106325922 0
504 0.081030248 0
505 0.689427302 1
506 0.598601859 1
507 0.854366023 1
508 0.166113000 0
509 0.641793851 1
510 0.268586674 0
511 0.371560848 0
512 0.817810771 1
513 0.154689652 0
514 0.477863908 0
515 0.331272209 0
516 0.655102725 1
517 0.287192743 0
518 0.356220185 0
519 0.017869104 0
520 0.024577319 0
521 0.112606101 0
522 0.040253733 0
523 0.178100944 0
524 0.202113167 0
526 0.224071414 0
527 0.595556780 1
528 0.112560088 0
529 0.101364304 0
530 0.314499761 0
531 0.030532750 0
532 0.264819406 0
533 0.090678823 0
534 0.197951844 0
535 0.575019898 1
536 0.732380586 1
537 0.774541871 1
538 0.015546062 0
539 0.547901118 1
540 0.537792298 1
541 0.218373931 0
542 0.477634090 0
543 0.855990660 1
544 0.366913623 0
545 0.178900727 0
546 0.426072523 0
547 0.045296124 0
548 0.259566037 0
549 0.027793478 0
550 0.178134655 0
551 0.098528251 0
552 0.448607818 0
553 0.074027896 0
554 0.508727725 1
555 0.034590411 0
556 0.518731098 1
557 0.173093748 0
558 0.158755934 0
559 0.059859160 0
560 0.208686070 0
561 0.051003462 0
562 0.151230864 0
563 0.394875350 0
Nog interessanter is om meteen ook de geobserveerde intentie
(intentie) hierbij te zien, zodat we geval per geval kunnen
nagaan of onze voorspelling van de intentie ook juist is. Daartoe voegen
we kans en intentie.voorspeld toe aan het
bestaande data frame barrieres. Let op! Dit zal enkel
lukken als het aantal rijen van barrieres,
kansen en intentie.voorspeld gelijk is. Het
oorspronkelijke data frame barrieres zal je dus moeten
inkorten. Er waren immers ontbrekende data, waardoor sommige observaties
uit barrieres geen waarde hebben bij kansen en
bij intentie.voorspeld. Het inkorten van het data frame kan
je als volgt doen.
barrieres <- barrieres[,c("pcg", "psn", "attitudes", "intentie")] # let op: we overschrijven hier de oorspronkelijke dataset. Wees daar steeds voorzichtig mee
barrieres <- na.omit(barrieres)
Nu kan je barrieres, kans en
intentie.voorspeld samenvoegen tot één data frame zonder
foutmeldingen te krijgen.
barrieres <- data.frame(barrieres, kansen, intentie.voorspeld)
head(barrieres) # toont de eerste zes rijen pcg psn attitudes intentie kansen intentie.voorspeld
1 5.75 5.666667 5.833333 1 0.6580606 1
2 5.25 6.000000 6.500000 1 0.7322452 1
3 5.50 5.333333 6.666667 0 0.6921800 1
4 6.25 3.000000 5.000000 0 0.2574048 0
5 5.25 4.000000 5.166667 0 0.3172350 0
6 5.50 4.000000 6.500000 1 0.4988325 0
Deze eerste zes rijen bevatten meteen een voorbeeld van alle mogelijke situaties met betrekking tot voorspelde en geobserveerde intentie.
-
Wanneer
intentieenintentie.voorspeldallebei gelijk zijn aan 1, dan is de voorspelling correct positief (“true positive”). -
Wanneer
intentieenintentie.voorspeldallebei gelijk zijn aan 0, dan is de voorspelling correct negatief (“true negative”). -
Wanneer
intentie = 0enintentie.voorspeld = 1, dan hebben we een vals positieve voorspelling gemaakt (“false positive”). -
Wanneer
intentie = 1enintentie.voorspeld = 0, dan hebben we een vals negatieve voorspelling gemaakt (“false negative”).
Nu zou je denken dat je een simpele berekening kan uitvoeren met het aantal correcte en foute voorspellingen die je model maakt. Daarmee zou je dan het model kunnen beoordelen. Jammer genoeg ligt het niet zo eenvoudig. Het probleem is dat de cutoff van 0.50 eigenlijk een willekeurige keuze is. Er kunnen goede redenen zijn om voor bijvoorbeeld 0.20 of 0.80 te kiezen in plaats van 0.50.
- Een lagere cutoff van 0.20 heeft als voordeel dat je minder risico loopt om een observatie te klasseren als “mislukking” terwijl het in werkelijkheid een “succes” is. Het nadeel van deze lage cutoff is het omgekeerde: het risico is groter dat je een observatie klasseert als “succes” terwijl het in werkelijkheid een “mislukking” is.
- Een hogere cutoff van 0.80 heeft als voordeel dat je minder risico loopt om een observatie te klasseren als “succes” terwijl het in werkelijkheid een “mislukking” is. Het nadeel van deze hoge cutoff is het omgekeerde: het risico is groter dat je een observatie klasseert als “mislukking” terwijl het in werkelijkheid een “succes” is.
Afhankelijk van welk risico je beter kan tolereren kan je kiezen voor een ander cutoff percentage. Er is met andere woorden geen absoluut juiste of foute cutoff.4
Een gevolg daarvan is dat de vector intentie.voorspeld
er helemaal anders kan uitzien als je kiest voor een andere cutoff. Het
is dus niet mogelijk om intentie.voorspeld te gebruiken om
een algemene beoordeling te maken van de voorspellende kracht van je
model.
Je bent voorlopig nog niet dichter gekomen bij een manier om de voorspellende kracht van het model te kwantificeren. Daarvoor heb je een perspectief nodig dat niet gebonden is aan één specifieke keuze voor een cutoff waarde. Dat kan je bekomen door de cutoff in hele kleine sprongen te laten variëren van 1 naar 0.
Om de uitleg goed te kunnen volgen is het een goed idee om nu eerst
de data te rangschikken volgens oplopende voorspelde kans. De functie
order() geeft je eerst het rijnummer van de laagste kans,
dan het rijnummer van de tweede laagste kans, dan het rijnummer van de
derde laagste kans, enzovoort. Die rijnummers stop je in de vector
volgorde. Vervolgens herordenen we het data frame in die
volgorde. Je ziet dat het data frame inderdaad gerangschikt
is volgens de kolom kansen. (De kolom
intentie.voorspeld die gebaseerd was op een cutoff van
0.5 laten we achterwege.)
volgorde <- order(barrieres$kansen)
barrieres <- barrieres[volgorde,]
barrieres <- subset(barrieres, select= -c(intentie.voorspeld)) # kolom intentie.voorspeld weglaten
head(barrieres) pcg psn attitudes intentie kansen
102 1.00 1.000000 2.000000 0 0.004234679
126 2.00 1.000000 2.333333 0 0.007022809
355 4.75 1.000000 1.000000 0 0.009056793
54 2.75 1.666667 1.833333 0 0.010094042
485 1.00 1.000000 4.000000 0 0.011604075
318 1.00 1.000000 4.166667 0 0.012615865
Nu ben je klaar om de cutoff in kleine sprongen te laten variëren van 1 naar 0.
Wat gebeurt er bij een cutoff gelijk aan 1? Dan zal je pas
voorspellen dat de intentie er is (“succes”) als de voorspelde kans
precies gelijk is aan 1. In alle andere gevallen is de voorspelling dat
de intentie er niet is (“mislukking”). Het is niet zo verrassend dat je
met zo’n hoge cutoff nooit een voorspelling “succes” zal maken. De
voorspelde kans is namelijk nooit exact gelijk aan 1. Dat zie je aan de
onderkant van het geordende data frame barrieres, waar de
hoogste voorspelde kansen te vinden zijn.
tail(barrieres) # toont de laatste zes rijen pcg psn attitudes intentie kansen
391 6.75 7.000000 6.5 1 0.8872493
146 7.00 6.666667 7.0 1 0.9019494
153 6.50 7.000000 7.0 0 0.9030833
235 6.75 7.000000 7.0 1 0.9102602
272 7.00 7.000000 7.0 1 0.9169545
476 7.00 7.000000 7.0 1 0.9169545Die extreem hoge cutoff heeft als voordeel dat je nooit een vals positieve voorspelling maakt. Het nadeel is dat je ook nooit een correct positieve voorspelling maakt.
De proportie vals positieve en correct positieve voorspellingen kan je in een assenstelsel plaatsen. De horizontale as stelt de proportie vals positieve voorspellingen voor5. Op de verticale as lees je de proportie correct positieve voorspellingen af6. Bij een cutoff gelijk aan 1 kom je in het punt helemaal links onderaan terecht.
Als je nu de cutoff geleidelijk laat zakken, dan zal je van steeds
meer gevallen (observaties) voorspellen dat het een “succes” zal zijn.
Voor elk individueel geval kan dat ofwel terecht zijn (correct positief)
ofwel verkeerd (vals positief). Op het moment dat je cutoff zakt onder
0.9169545 dan maak je een eerste positieve voorspelling. Dat kan je zien
in de output van tail(barrieres). Dat is een correct
positieve voorspelling, want de geobserveerde intentie is
gelijk aan 1.
De proportie correct positieve voorspellingen is dus toegenomen, terwijl de proportie vals positieve voorspellingen nog steeds gelijk is aan 0. Op de grafiek hieronder (waarin we inzoomen op de linkeronderkant van de grafiek hierboven) zie je dat er daardoor een sprong naar boven wordt gemaakt. Er is zelfs meteen een dubbele sprong, want er zijn twee observaties met identiek dezelfde voorspelde kans. Van zodra de cutoff onder die kans duikt, wordt natuurlijk voor beide gevallen een positieve voorspelling gemaakt. Het blijkt in beide gevallen om een correct positieve voorspelling te gaan, dus we maken een dubbele sprong omhoog.
De curve die je nu stapsgewijs aan het maken bent heet een ROC-curve.
Wanneer je de cutoff verder laat zakken tot 0.9102602 komt er nog een correct positieve voorspelling bij. Dat vertaalt zich in nog een sprong naar boven op de grafiek.
Bij een volgende verlaging van de cutoff tot 0.9030833 maak je een
bijkomende positieve voorspelling, maar deze keer is dat een vals
positieve voorspelling (want de geobserveerde intentie is
gelijk aan 0). De proportie vals positieve voorspellingen
is nu niet meer gelijk aan 0, maar iets hoger. Visueel vertaalt dat zich
in een sprong naar rechts.
Dit proces van dalende cutoffs blijf je herhalen tot de cutoff gezakt is naar 0. Dan worden alleen nog positieve voorspellingen gemaakt. De proportie correct positieve voorspellingen is gelijk aan 1 en hetzelfde geldt voor de proportie vals positieve voorspellingen. Je bent terechtgekomen in het punt helemaal rechtsboven op de grafiek. De resulterende curve zie je hieronder.
Je hoeft die ROC-curve natuurlijk niet manueel en stapsgewijs op te
bouwen. Met behulp van het package pROC kan je die
makkelijk plotten. Dat gebeurt in twee stappen. Eerst gebruik je de
functie roc(). Die heeft twee argumenten nodig: enerzijds
de namen van de kolommen met de geobserveerde intentie en
de berekende kansen, anderzijds het data frame waarin die
kolommen te vinden zijn (barrieres).
mijn.roc <- roc(formula = intentie ~ kansen, data = barrieres)
In een tweede stap geef je het object mijn.roc aan de
(basis)functie plot().
plot(mijn.roc)
Dat is misschien allemaal heel boeiend, maar wat heeft dit nu te maken met de voorspellende kracht van het model? Gelijk welk model zal evolueren op de grafiek van linksonder naar rechtsboven wanneer de cutoff afneemt van 1 naar 0. Waarin onderscheiden goede modellen zich dan van slechte? De sleutel zit in het traject dat een model aflegt van linksonder naar rechtsboven.
Een goed model is een model dat relatief hoge kansen voorspelt voor
gevallen waar intentie gelijk is aan 1. Een
dergelijk model zal bij afnemende cutoffwaarden van 1 naar 0 eerst
vooral correct positieve voorspellingen maken. Visueel betekent
dit dat de curve eerst steil omhoog zal gaan. Tevens zal een goed model
lagere kansen voorspellen voor gevallen waar intentie
gelijk is aan 0. De vals positieve voorspellingen
(en de bijhorende sprongen naar rechts) zullen dus pas komen nadat de
meeste correct positieve voorspellingen (en bijhorende sprongen omhoog)
al gepasseerd zijn.
Het traject van een goed model zal dus dicht tegen de linkerzijde en bovenzijde van de grafiek lopen.
Een slecht model doet natuurlijk het omgekeerde. Dan worden juist
relatief hoge kansen voorspeld bij gevallen waar de geobserveerde
intentie gelijk is aan 0. Bij afnemende cutoff
zullen dan eerst vooral vals positieve voorspellingen worden
gemaakt. Dat betekent dat de curve eerst horizontale sprongen maakt.
Daarna volgen pas de correct positieve voorspellingen en de bijhorende
sprongen naar boven.
Het traject van een slecht model loopt bijgevolg dicht tegen de onderzijde en rechterzijde van de grafiek.
Waar bevindt jouw model zich precies tussen die twee uitersten? Dat kan je beoordelen aan de hand van de oppervlakte onder de curve (vandaar de naam “Area Under the Curve”).
Om de uitleg te verduidelijken zie je hieronder dezelfde curve, waarbij de oppervlakte onder de curve groen is gekleurd en de oppervlakte erboven rood. Je kan nu berekenen hoe groot de groene oppervlakte onder de curve is ten opzichte van het totaal van de groene plus de rode oppervlakte. In een ideaal model neemt de groene oppervlakte de hele grafiek in. Dan is de AUC gelijk aan 1. Bij het slechtst denkbare model is de groene oppervlakte onbestaande. De AUC is dan gelijk aan 0. De verhouding van de groene oppervlakte ten opzichte van de totale oppervlakte in jouw model ligt ergens tussen 0 en 1. Hoe dichter bij 1, hoe beter.
Met de functie auc() uit het package pROC
kan je de AUC berekenen. Daarvoor heb je enkel het object
mijn.roc nodig dat je al eerder maakte.
auc(mijn.roc)Area under the curve: 0.8108
De AUC interpreteren
De vraag is nu of deze waarde voor de AUC hoog of laag is. Met andere woorden, is jouw model nu goed of niet? Jammer genoeg is er geen onderbouwde grenswaarde die goede modellen duidelijk scheidt van slechte modellen. Toch zijn er enkele zaken die het vermelden waard zijn en die als richtlijnen kunnen dienen bij de interpretatie van de AUC.
-
Als je gewoon zou gokken wat de geobserveerde
intentieis, dan mag je verwachten om een AUC te bekomen van 0.50. De berekende AUC van 0.8108 moet je dus interpreteren als een waarde op een spectrum van 0.50 tot 1 en niet op een spectrum van 0 tot 1. - Sommige auteurs hebben wel vuistregels voorgesteld voor een grenswaarde tussen goede en slechte modellen. Hosmer, Lemeshow & Sturdivant (2013) beschouwen een AUC vanaf 0.70 als aanvaardbaar. De auteurs waarschuwen echter zelf dat dit niet als een magische grens mag worden beschouwd en dat de context relevant is. Als je doel bijvoorbeeld is om een test te ontwikkelen die een levensbedreigende ziekte kan detecteren, dan is 0.70 onaanvaardbaar laag.
- De AUC is waarschijnlijk nog het meest nuttig wanneer je modellen gaat vergelijken. Dan focus je vooral op de vraag in welke mate de AUC toeneemt wanneer je bijvoorbeeld een predictor toevoegt aan je model.
Hosmer-Lemeshow
Hoe kan je een “goodness-of-fit” toets uitvoeren voor een logistisch regressiemodel? Hosmer-Lemeshow is hiervoor de meest gebruikte toets. Die steunt op een vergelijking van geobserveerde en verwachte frequenties.
Je kan daarvoor het package generalhoslem gebruiken.
install.packages("generalhoslem") # eenmalig installeren
library(generalhoslem) # bij de start van elke sessie laden
De functies uit het package maken het uitvoeren van de toets makkelijk, maar omdat het ook belangrijk is om te begrijpen wat de toets inhoudt vind je hieronder eerst wat theoretische achtergrond. We gaan daar dieper op in onder meer omdat Hosmer-Lemeshow geen leerstof is in opleidingen aan de FPPW.
In een eerste stap worden de gegevens opgedeeld in groepen, bijvoorbeeld tien. Je kan dan een tabel opstellen die er als volgt uit ziet.
| Groep 1 | Groep 2 | Groep 3 | Groep 4 | Groep 5 | Groep 6 | Groep 7 | Groep 8 | Groep 9 | Groep 10 | |
|---|---|---|---|---|---|---|---|---|---|---|
| Intentie | ||||||||||
| Geen intentie |
Het is een \(2 \times 10\) tabel. De twee rijen staan voor de twee mogelijke uitkomsten: intentie om een opleiding te volgen of niet (ook wel “succes” en “mislukking” genoemd). Er zijn tien kolommen omdat de data in dit voorbeeld zijn opgedeeld in tien groepen. Dat aantal is in principe vrij te kiezen, maar doorgaans werkt men met tien groepen.
De opdeling in groepen gebeurt door de gegevens te rangschikken volgens de voorspelde kans en vervolgens tien groepen van gelijke grootte te maken.
Nu kan je de twintig cellen van de bovenstaande tabel op twee manieren opvullen.
Je kan in elke cel van de tabel de geobserveerde frequentie voor die cel \(f_{obs, ij}\) noteren.
Groep 1 Groep 2 Groep 3 Groep 4 Groep 5 Groep 6 Groep 7 Groep 8 Groep 9 Groep 10
Intentie 2 3 9 17 13 25 34 29 41 44
Geen intentie 52 51 45 38 40 28 20 25 14 9
Anderzijds kan je de twintig cellen in de tabel ook opvullen met voorspelde frequenties voor elke cel \(f_{exp, ij}\). Die zijn gebaseerd op de voorspelde kansen op basis van het model.
Groep 1 Groep 2 Groep 3 Groep 4 Groep 5 Groep 6 Groep 7 Groep 8 Groep 9 Groep 10
Intentie 2.244678 5.498507 8.860662 13.40113 17.9386 22.92589 28.47789 33.4737 39.48311 44.695836
Geen intentie 51.755322 48.501493 45.139338 41.59887 35.0614 30.07411 25.52211 20.5263 15.51689 8.304164
Beide tabellen zullen doorgaans niet identiek zijn. De verschillen per cel kan je gebruiken om het model te beoordelen. Bij een goed model sluiten de voorspelde frequenties dicht aan bij de geobserveerde frequenties. De verschillen zijn met andere woorden klein. Bij een slecht model geldt het omgekeerde.
Om de “goodness-of-fit” van het model te toetsen bereken je een teststatistiek aan de hand van volgende formule:
\[\sum^2_{i=1} \sum^N_{j=1} \frac{
(f_{obs, ij} - f_{exp, ij})^2 }{ f_{exp, ij} }\]
In de formule zie je dat de verschillen tussen de geobserveerde en de voorspelde frequenties inderdaad centraal staan (\(f_{obs, ij} - f_{exp, ij}\)). Als die verschillen over het algemeen groot zijn, dan neemt de teststatistiek een hoge waarde aan. Als die verschillen over het algemeen klein zijn, dan krijgt de teststatistiek een lage waarde.
De nulhypothese stelt dat er geen verschillen zijn tussen de geobserveerde en voorspelde frequenties. De alternatieve hypothese stelt dat er minstens één verschil is (in de twintig cellen).
De teststatistiek benadert een \(\chi^2\)-verdeling met \(N-2\) vrijheidsgraden, waarbij \(N\) het aantal groepen voorstelt (hier 10). Aan de hand van de waarde van de teststatistiek in jouw steekproef kan je afleiden of er bewijs is tegen de nulhypothese.
In R kan je de toets uitvoeren met de functie logitgof()
uit het package generalhoslem. Deze functie vereist twee
stukken informatie: de geobserveerde intentie en de
voorspelde kansen.
In R kan je de voorspelde kansen laten berekenen met de functie
fitted(). Een meer uitgebreide uitleg vind je wat hoger op
deze pagina, onder Predicties. Hieronder voegen
we die kansen meteen toe aan het data frame barrieres.7
barrieres$kansen <- fitted(model1)
Nu heb je alles om de functie logitgof() te
gebruiken.
logitgof(obs=barrieres$intentie, exp=barrieres$kansen)
Hosmer and Lemeshow test (binary model)
data: barrieres$intentie, barrieres$kansen
X-squared = 9.0725, df = 8, p-value = 0.3362
Je stelt een p-waarde vast die groter is dan 0.05. Je besluit dat je de nulhypothese niet kan verwerpen op het 5% significantieniveau.
Enkele bedenkingen bij de Hosmer-Lemeshow toets
De Hosmer-Lemeshow toets werkt niet in alle omstandigheden perfect.
Zo is er in bepaalde gevallen een lage power. Dat betekent dat vaak verkeerdelijk besloten wordt dat het model goed fit, terwijl dat niet het geval is.
Daarnaast is de indeling in groepen niet altijd objectief te verantwoorden. Zowel het aantal groepen als de manier waarop de data worden verdeeld kunnen verschillen ten gevolge van arbitraire keuzes. Dat kan een invloed hebben op de waarde van de teststatistiek en dus op de conclusie van de toets.
13 Referenties
Van Nieuwenhove L. & De Wever B. (2023). Psychosocial barriers to adult learning and the role of prior learning experiences: A comparison based on educational level. Adult Education Quarterly 74 (1), pp.62-87. DOI: 10.1177/07417136231147491
14 Voetnoten
Elk van deze drie variabelen is berekend door meerdere items uit een enquête te combineren.↩︎
In sommige vakgebieden, zoals de psychometrie, geven onderzoekers vaak de voorkeur aan
link="probit". Dat heeft voordelen in bepaalde gevallen maar hier gaan we daar niet dieper op in.↩︎Deze termen hebben natuurlijk geen enkele normatieve betekenis. “Succes” betekent enkel dat een bepaalde gebeurtenis zich voordoet en “mislukking” dat die gebeurtenis zich niet voordoet. In deze demonstratie slaat die gebeurtenis op het hebben van de intentie om een opleiding te volgen.↩︎
Wel is het zo dat cutoffs van 0 of van 1 weinig zinvol zijn.↩︎
\(\frac{Aantal \; vals \; positieve \; voorspellingen}{Aantal \; mislukkingen}\)↩︎
\(\frac{Aantal \; correct \; positieve \; voorspellingen}{Aantal \; successen}\)↩︎
Dat hadden we ook al gedaan bij de uitleg over de AUC. Hier herhalen we het voor de duidelijkheid. Merk op dat we bij de uitleg over de AUC de volgorde van het data frame
barriereshebben veranderd. Als je dat ook hebt gedaan, dan zal onderstaande code niet tot het correcte resultaat leiden.↩︎