Chapter 12 Logistische regressie
Logistische regressie kan worden gebruikt als de afhankelijke variabele binair is, dat wil zeggen, alleen de waarden 0 en 1 heeft. Er bestaan ook logistische regressiemodellen voor categorische afhankelijke variabelen (meer dan twee categorieën); hierover later meer. Wanneer de afhankelijke variabele binair is, is een lineair regressiemodel niet geschikt om de scores van de afhankelijke variabele te voorspellen, en wordt een logistisch model gebruikt. De onafhankelijke variabelen moeten minimaal van intervalniveau zijn of dummy-variabelen. De procedure loopt via glm(). en summary().
De functie die wordt gebruikt om het model te schatten isglm(); dit staat voor Generalized Linear Model. De structuur is hetzelfde als bij de ‘gewone’ lineaire regressie functie lm(), waarin je de regressievergelijking moet opgeven met een ~ tussen de afhankelijke en onafhankelijke variabelen (bijv. y ~ x1 + x2). Voor de glm() functie moet een extra argument worden gebruikt waarmee je aangeeft welk type model (familie van modellen) er moet worden geschat. Dit argument is family = [naam verdeling] (link = “[naam link]”); de namen in tussen vierkante haken [ ] moeten worden ingevuld. Voor logistische regressie is de code: family = binomial(link = “logit”). Als voorbeeld zal de variabele cut uit de diamonds dataframe worden gebruikt. Deze moet eerst binair worden gemaakt om te worden gebruikt in logistische regressie.
diamonds <- diamonds %>% mutate(
cut_bin = case_when(
cut == "Fair" ~ 0,
cut == "Good" ~ 0,
cut == "Very Good" ~ 1,
cut == "Premium" ~ 1,
cut == "Ideal" ~ 1))De nieuwe variabele cut_bin geeft nu aan of de kwaliteit goed is (cut_bin = 1) of niet (cut_bin = 0). Met een logistische regressie kunnen we nu schatten wat de kans is dat de cut van goede kwaliteit is, gegeven het gewicht van de diamant (carat).
LR_model1 <- glm(cut_bin ~ carat,
family = binomial (link = "logit"),
data = diamonds)
summary(LR_model1)##
## Call:
## glm(formula = cut_bin ~ carat, family = binomial(link = "logit"),
## data = diamonds)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.37456 0.02628 90.35 <2e-16 ***
## carat -0.46461 0.02578 -18.02 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 39756 on 53939 degrees of freedom
## Residual deviance: 39445 on 53938 degrees of freedom
## AIC: 39449
##
## Number of Fisher Scoring iterations: 4We hebben de analyse nu uitgevoerd, de output bevat de volgende zaken:
- Call: hier wordt de code van de uitgevoerde analyse weergeven.
- Deviance Residuals: hier wordt de verdeling van de residuen weergeven.
- Coefficients: voor de onafhankelijke variabelen wordt in de output de regressiecoëfficiënt, de standaardfout, de z-waarde en tweezijdige overschrijdingskans voor het toetsen of de coëfficiënt van nul verschilt weergeven (ook wel de Waldtoets).
Vanuit de output kan herleid worden dat het model om de prijs te van een diamant te bepalen vanuit het gewicht van de diamand er als volgt uitziet: \(log [ \frac{\hat{p}(y=1)}{1- \hat{p}(y=1)}] = 2.375 - 0.465* carat\) met y = kwaliteit van de diamand. De negatieve waarde van de helling laat zien dat naarmate het gewicht van een diamant toeneemt, de logodds (\(log [ \frac{\hat{p}(y=1)}{1- \hat{p}(y=1)}]\)) kleiner wordt, daarmee worden de odds (\(\hat{p}(y=1)\)) kleiner en dus wordt de kans op hebben van een kwalitatief goede diamand kleiner. De coëfficiënten van het regressiemodel kunnen ook opgevraagd worden met de code coef()
## (Intercept) carat
## 2.3745562 -0.464612112.1 Betrouwbaarheidsinterval voor regressiecoëfficiënten
Het berekenen van de betrouwbaarheidsintervallen van regressiecoëfficieënten kan met de code confint.default().
## 2.5 % 97.5 %
## (Intercept) 2.3230432 2.4260692
## carat -0.5151443 -0.414079812.2 Odds ratio
De odds ratio geeft de factor waarmee de odds vermenigvuldigd wordt als een onafhankelijke variabele 1 waarde stijgt. De odds-ratio wordt berekend door \(e^{helling}\) of in R kan je dit berekenen met de code exp().
De odds-ratio van de variabele carat is 0.628 (wat overeenkomt met \(e^{-0.4646121}\)), wat betekent dat als het gewicht van de diamant met één eenheid stijgt, de odds op het hebben van een kwalitatief goede diamant 0.628 keer zo klein wordt.
## (Intercept) carat
## 10.7462431 0.6283788Het betrouwbaarheidsinterval van de odds-ratio kan ook berekend worden emt de functie exp()
## 2.5 % 97.5 %
## (Intercept) 10.2066883 11.3143204
## carat 0.5974144 0.660948212.3 Geschatte kansen
In logistische regressie kunnen de geschatte kansen, ofwel de kans dat de waarde van de afhankelijke variabele gelijk is aan 1, berekend worden met de functie predict(). In deze functie gebruik je het object met het geschatte model en het extra argument type = “response” (anders worden voorspelde log-odds scores gegeven i.p.v. geschatte kansen): p_hat = predict(LR_model1, type = “response”). De geschatte kansen (p_hat) kunnen in de dataset worden bewaard (als extra variabele) met behulp van de functie mutate().
12.4 Classificatietabellen
Classificatietabellen weergeven in hoeverre de score (0 of 1) op de afhankelijke variabele goed voorspeld kan worden op basis van de geschatte kansen. Hierbij wordt er gebruik gemaakt van een cut-value, wat de grens tussen score 0 en 1 aangeeft. Bij een cut-value van 0.5 worden kansen onder de 0.5 geclassificeerd als score 0 en boven de 0.5 als score 1. Voor het maken van classificatietabellen is devolgde functie geschreven.
# Functie voor classificatietabellen
Class.table <- function(LOGMOD = NULL){
library(tidyverse)
DATSET <- LOGMOD$data
DATSET <- mutate(DATSET,
p_hat = predict(LOGMOD, type = "response"),
y_hat = as.factor(ifelse(p_hat >= 0.5, 1, 0)))
DV <- LOGMOD$formula[[2]]
Class_tmp <- table(DATSET[[DV]], DATSET$y_hat)
Class_tab <- matrix(data = NA, nrow = 3, ncol = 3)
# Controle of alle cases in 1 groep geclassificeerd worden
C1 <- ifelse(Class_tmp[1]+Class_tmp[2] == length(DATSET[[DV]]), T, F)
# Controle of dit de 0 groep is (anders 1)
C2 <- ifelse(dimnames(Class_tmp)[[2]] == "0", T, F)
if(C1){
Class_tab[1,1] <- ifelse(C2, Class_tmp[1], 0)
Class_tab[1,2] <- ifelse(C2, 0, Class_tmp[1])
Class_tab[2,1] <- ifelse(C2, Class_tmp[2], 0)
Class_tab[2,2] <- ifelse(C2, 0, Class_tmp[2])
} else {
Class_tab[1:2,1:2] <- Class_tmp
}
Class_tab[1,3] <- round(Class_tab[1,1] / (Class_tab[1,1] + Class_tab[1,2]), 4)
Class_tab[2,3] <- round(Class_tab[2,2] / (Class_tab[2,1] + Class_tab[2,2]), 4)
Class_tab[3,3] <- round((Class_tab[1,1] + Class_tab[2,2]) / length(DATSET[[DV]]), 4)
rownames(Class_tab) <- c("Obs0", "Obs1", "Tot")
colnames(Class_tab) <- c("Exp0", "Exp1", "Tot")
return(Class_tab)
}De functie heeft één argumenten (input), namelijk het model. In het onderstaande voorbeeld is een classificatietabel van model 1 (kwaliteit van diamant schatten vanuit het gewicht) weergeven. In dit geval kan 87.92% van de cases goed voorspeld worden. Door van ieder model (incl. het lege model) een classificatiemodel te maken, kan de invloed van toegevoegde variabelen op het goed voorspellen van de scores op de afhankelijke variabelen onderzocht worden. Daarnaast geven classificatietabellen inzicht in hoe goed het model in staat is om de waarden op de afhankelijke variabelen te voorspellen.
## Exp0 Exp1 Tot
## Obs0 0 6516 0.0000
## Obs1 0 47424 1.0000
## Tot NA NA 0.8792Voor het toetsen van de fit van het model is het noodzakelijk om het zogenoemde lege model (met alleen een constante en geen predictoren) te schatten. Dit model wordt geschat door na de ~ in de code voor een logistische regressie een 1 te plaatsen, wordt het lege model geschat.
LR_model0 <- glm(cut_bin ~ 1,
family = binomial (link = "logit"),
data = diamonds)
summary(LR_model0)##
## Call:
## glm(formula = cut_bin ~ 1, family = binomial(link = "logit"),
## data = diamonds)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.98487 0.01321 150.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 39756 on 53939 degrees of freedom
## Residual deviance: 39756 on 53939 degrees of freedom
## AIC: 39758
##
## Number of Fisher Scoring iterations: 412.5 Likelihood ratio toets
De likelihood ratio toets wordt gebruikt om de -2LogLikelihood of Deviance van verschillende modellen tegen elkaar te toetsen. Met deze toets kan een directe vergelijking tussen twee modellen gemaakt worden, waarbij getoetst wordt of het complete model beter in staat is om de score op de afhankelijke variabele te schatten dat het gereduceerde model. De Likelihood ratio toets wordt uitgevoerd met de code anova(). De functie anova() is algemener dan alleen de variantieanalyse, zoals gebruikt in het lineaire model. Het wordt ook gebruikt in sequentiële analyses, om te toetsen of geneste modellen significant van elkaar verschillen; niet alleen lineaire modellen, maar alle generalized linear models.
## Analysis of Deviance Table
##
## Model 1: cut_bin ~ 1
## Model 2: cut_bin ~ carat
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 53939 39756
## 2 53938 39445 1 310.35 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In de output is onder Resid. Dev de deviance score van beide modellen weergeven. Onder Deviance is het verschil in deviance tussen de twee modellen weergeven. Tot slot wordt onder Pr(>Chi) de toets voor het verschil in Deviance weergeven.
12.6 Hosmer Lemeshow toets
De Hosmer Lemeshow toets toetst met een chi-kwadraat toets of de geobserveerde en geschatte kansen significant van elkaar verschillen. De Hosmer Lemeshow toets wordt uitgevoerd met de code hltest() uit de glmtoolbox-package.
##
## The Hosmer-Lemeshow goodness-of-fit test
##
## Group Size Observed Expected
## 1 5442 4692 4471.6470
## 2 5491 4908 4699.5379
## 3 4327 3959 3755.4359
## 4 5078 3769 4426.6118
## 5 5146 4323 4532.9022
## 6 5258 4463 4666.3486
## 7 5524 4952 4938.0205
## 8 5650 5133 5077.8539
## 9 5572 5305 5025.1963
## 10 5414 4965 4890.3926
## 11 1038 955 940.0533
##
## Statistic = 1309.84
## degrees of freedom = 9
## p-value = < 2.22e-16Uit de toetsresultaat van de Hosmerlemeshow toets blijkt dat de p-waarde lager is dan het significantieniveau 0.05, de nulhypothese wordt dus verworpen. De geobserveerde en geschatte kansen zijn dus niet gelijk, wat aangeeft dat de fit van het model niet goed is.