Logistische regressie in R – tutorial

Logistische regressie is een eenvoudig maar krachtig model om binaire uitkomsten te voorspellen. Met andere woorden: of iets wel of niet gebeurt. Het is een type classificatiemodel voor supervised machine learning.

Logistische regressie wordt in bijna elke sector gebruikt—marketing, gezondheidszorg, sociale wetenschappen en meer—en is een essentieel onderdeel van de gereedschapskist van elke data scientist.

Om het meeste uit deze tutorial te halen, heb je basiskennis van R nodig. Het helpt ook als je iets weet over een verwant modeltype: lineaire regressie. Lees de Linear Regression in R Tutorial om daar meer over te weten te komen.

Een overzicht van logistische regressie

Stel dat je wilt voorspellen of het vandaag een zonnige dag wordt of niet. Er zijn twee mogelijke uitkomsten: "zonnig" of "niet zonnig". De uitkomstvariabele staat ook wel bekend als een "doelvariabele" of een "afhankelijke variabele".

Er zijn veel variabelen die de uitkomst kunnen beïnvloeden, zoals ‘temperatuur van de dag ervoor’, ‘luchtdruk’ enz. De beïnvloedende variabelen staan bekend als features, onafhankelijke variabelen of voorspellers—al deze termen betekenen hetzelfde. 

Andere voorbeelden zijn of een klant je product wel of niet zal kopen, of een e-mail spam is of niet, of een transactie frauduleus is of niet, en of een geneesmiddel een patiënt wel of niet geneest.  

Hoe werkt logistische regressie?

Logistische regressie zoekt de best mogelijke fit tussen de voorspellers en de doelvariabele om de kans te voorspellen dat de doelvariabele tot een gelabelde klasse/categorie behoort. 

Lineaire regressie probeert de beste rechte lijn te vinden die de uitkomst op basis van de features voorspelt. Die vormt een vergelijking zoals

y_predictions = intercept + slope * features

en gebruikt optimalisatie om de best mogelijke waarden voor intercept en helling te vinden.

Logistische regressie werkt vergelijkbaar, behalve dat het regressie uitvoert op de kansen dat de uitkomst tot een categorie behoort. Het gebruikt een sigmoidfunctie (de cumulatieve verdelingsfunctie van de logistische verdeling) om de rechterkant van die vergelijking te transformeren.

y_predictions = logistic_cdf(intercept + slope * features)

Ook hier gebruikt het model optimalisatie om de best mogelijke waarden voor intercept en helling te vinden.

Omdat het algoritme voor logistische regressie sterk lijkt op de vergelijking voor lineaire regressie, valt het onder de familie van modellen die "generalized linear models" worden genoemd. Daarom staat er "regressie" in de naam, ook al is het een classificatiemodel.

De sigmoidfunctie lijkt op een S-vormige curve in de onderstaande afbeelding. Ze neemt reëelwaardige invoer en zet die om naar waarden tussen 0 en 1 (door aan beide kanten te krimpen, d.w.z. negatieve waarden naar 0 en zeer hoge positieve naar 1). Bovendien is de drempelwaarde de doorslaggevende factor die de uitvoer indeelt in categorieën of klassen wanneer die op deze kansen wordt toegepast.

Complexe concepten begrijp je het best met voorbeelden, dus laten we een analogie gebruiken om het LR-algoritme uit te leggen. Stel dat het LR-model de taak krijgt om een frauduleuze transactie te identificeren op basis van meerdere fraude-indicatoren, zoals de locatie van de gebruiker, het aankoopbedrag, IP-adres, enz. Het doel is om de waarschijnlijkheid te bepalen dat een bepaalde transactie legitiem of frauduleus is – dat is de doelvariabele.

Het model kent gewichten toe aan de voorspellers op basis van hun impact op de doelvariabele en combineert ze om de genormaliseerde score of kans op fraude te berekenen. 

De dataset

We gebruiken een dataset van een direct marketingcampagne van een Portugese bank die via telefoontjes klanten benadert. De campagne heeft als doel om abonnementen op een bankdeposito te verkopen, weergegeven door de variabele y (wel of geen abonnement). Het doel van het logistische regressiemodel is om te voorspellen of een klant een abonnement zal nemen op basis van de voorspellers, oftewel kenmerken van de klant, zoals demografische informatie.

De datadictionary voor deze dataset en veel andere nuttige datasets vind je op de website van Datacamp.

Workflow voor logistische regressie

De volledige machine learning-workflow wordt behandeld in de infographic A Beginner's Guide to The Machine Learning Workflow. Hier richten we ons op de stappen voor datapreparatie en modellering. We behandelen specifiek:

• De data splitsen in trainings- en testsets
• Het model fitten (trainen)
• Voorspellingen maken
• De modelprestatie evalueren
• Hyperparametertuning

Packages voor logistische regressie

In R zijn er twee populaire workflows voor het modelleren van logistische regressie: base-R en tidymodels.

De base-R-workflow is eenvoudiger en bevat functies zoals glm() en summary() om het model te fitten en een modeloverzicht te genereren.

De tidymodels-workflow maakt het makkelijker om meerdere modellen te beheren en biedt een consistente interface voor verschillende modeltypen.

In deze tutorial gebruiken we de tidymodels-workflow.

Data visualiseren

Importeer het tidymodels-package met de functie library(). 

De dataset staat in een CSV-bestand met Europese opmaak (komma’s als decimaalteken en puntkomma’s als scheidingsteken). We lezen het in met read_csv2() uit het readr-package. 

Zet de doelvariabele y om naar een factorvariabele voor het modelleren. 

Gebruik de functie ggplot() om het aantal per beroepsgroep te plotten ten opzichte van y.

library(readr)
library(tidymodels)

# Read the dataset and convert the target variable to a factor
bank_df <- read_csv2("bank-full.csv")
bank_df$y = as.factor(bank_df$y)

# Plot job occupation against the target variable
ggplot(bank_df, aes(job, fill = y)) +
    geom_bar() +
    coord_flip()

Data splitsen

Laten we de dataset splitsen in een trainingsset om het model te fitten en een testset om het model te evalueren, zodat we zeker weten dat het getrainde model ook werkt op onzichtbare data.

Data splitsen in een train- en testset kan met de functie initial_split() en het attribuut prop dat het aandeel voor de traindata bepaalt.

# Split data into train and test
set.seed(421)
split <- initial_split(bank_df, prop = 0.8, strata = y)
train <- split %>% 
         training()
test <- split %>% 
        testing()

Het model fitten en evalueren

Om het model te maken, declareer je een logistic_reg()-model. Dit heeft de argumenten mixture en penalty nodig die de hoeveelheid regularisatie bepalen. Een mixture-waarde van 1 staat voor een lasso-model en 0 voor ridge-regressie. Waarden daartussen zijn ook toegestaan. Het argument penalty geeft de sterkte van de regularisatie aan.

Let op dat je "double" floatingpointgetallen moet doorgeven aan mixture en penalty

Stel de "engine" (de backendsoftware die de berekeningen uitvoert) in met set_engine(). Er zijn meerdere keuzes: de standaardengine is "glm", die een klassieke logistische regressie uitvoert. Dit heeft vaak de voorkeur van statistici omdat je p-waarden voor elke coëfficiënt krijgt, wat het makkelijker maakt om het belang van elke coëfficiënt te begrijpen.

Hier gebruiken we de "glmnet"-engine. Deze heeft de voorkeur van machinelearningspecialisten omdat je hiermee kunt regulariseren, wat de voorspellingen kan verbeteren, vooral als je veel features hebt. (In Python neemt het scikit-learn-package standaard enige regularisatie op in logistische regressie.)

Roep de methode fit() aan om het model te trainen op de trainingsdata uit de vorige stap. Deze neemt een formule als eerste argument. Aan de linkerkant van de formule gebruik je de doelvariabele (hier y). Aan de rechterkant kun je elke gewenste feature opnemen. Een punt betekent: "gebruik alle variabelen die niet aan de linkerkant van de formule staan". Voor meer informatie over het schrijven van formules, lees de R Formula Tutorial.

# Train a logistic regression model
model <- logistic_reg(mixture = double(1), penalty = double(1)) %>%
  set_engine("glmnet") %>%
  set_mode("classification") %>%
  fit(y ~ ., data = train)

# Model summary
tidy(model)

De output staat hieronder, waarbij de kolom estimate de coëfficiënten van de voorspellers weergeeft.

# A tibble: 43 × 3
   term              estimate penalty
   <chr>                <dbl>   <dbl>
 1 (Intercept)      -2.59           0
 2 age              -0.000477       0
 3 jobblue-collar   -0.183          0
 4 jobentrepreneur  -0.206          0
 5 jobhousemaid     -0.270          0
 6 jobmanagement    -0.0190         0
 7 jobretired        0.360          0
 8 jobself-employed -0.101          0
 9 jobservices      -0.105          0
10 jobstudent        0.415          0
# ... with 33 more rows
# ℹ Use `print(n = ...)` to see more rows

Voorspellingen maken

Maak voorspellingen op de testdata met de functie predict(). Je kunt kiezen welk type voorspelling je wilt.

• type = "class" geeft voor elke observatie de meest waarschijnlijke doelwaarde terug. Hier zal dat "yes" of "no" zijn, afhankelijk van of het model denkt dat de klant een termijndeposito zal nemen of niet.
• type = "prob" geeft voor elke observatie de kans op elke doelwaarde terug. Hier krijg je de kans op "yes" en de kans op "no" (die per observatie optellen tot één).

# Class Predictions
pred_class <- predict(model,
                      new_data = test,
                      type = "class")

# Class Probabilities
pred_proba <- predict(model,
                      new_data = test,
                      type = "prob")

Evalueer het model met de functie accuracy() met als truth-argument y en als estimate-argument de voorspellingen uit de vorige stap.

results <- test %>%
           select(y) %>%
           bind_cols(pred_class, pred_proba)

accuracy(results, truth = y, estimate = .pred_class)

Hyperparametertuning

In plaats van specifieke waarden door te geven aan de argumenten mixture en penalty (de "hyperparameters"), kun je de voorspellende kracht van het model optimaliseren door te tunen. 

Het idee is dat je het model vaak draait met verschillende waarden voor de hyperparameters en bekijkt welke de beste voorspellingen geeft.

• Zet in de functie logistic_reg() de argumenten mixture en penalty op een aanroep van tune().
• Gebruik de functie grid_regular() om een raster met mogelijke waarden voor mixture en penalty te definiëren.
• De functie workflow() maakt een object om de modeldetails op te slaan, wat nodig is wanneer je het vaak draait. 
• Kies het beste model met de functie select_best(). Je kunt kiezen uit verschillende metrics die bepalen wat "het beste" is. Hier gebruiken we de metric "receiver operating characteristic area under the curve" (ROC AUC).

# Define the logistic regression model with penalty and mixture hyperparameters
log_reg <- logistic_reg(mixture = tune(), penalty = tune(), engine = "glmnet")

# Define the grid search for the hyperparameters
grid <- grid_regular(mixture(), penalty(), levels = c(mixture = 4, penalty = 3))

# Define the workflow for the model
log_reg_wf <- workflow() %>%
  add_model(log_reg) %>%
  add_formula(y ~ .)

# Define the resampling method for the grid search
folds <- vfold_cv(train, v = 5)

# Tune the hyperparameters using the grid search
log_reg_tuned <- tune_grid(
  log_reg_wf,
  resamples = folds,
  grid = grid,
  control = control_grid(save_pred = TRUE)
)

select_best(log_reg_tuned, metric = "roc_auc")

# A tibble: 1 × 3
       penalty mixture .config              
         <dbl>   <dbl> <chr>                
1 0.0000000001       0 Preprocessor1_Model01

Meer evaluatiemetrics

Met de beste hyperparameters:

1. Trainen we een model voor logistische regressie
2. Gebruiken we het om voorspellingen op de testset te genereren
3. Maken we een confusion matrix met de echte waarden en de schattingen

# Fit the model using the optimal hyperparameters
log_reg_final <- logistic_reg(penalty = 0.0000000001, mixture = 0) %>%
                 set_engine("glmnet") %>%
                 set_mode("classification") %>%
                 fit(y~., data = train)

# Evaluate the model performance on the testing set
pred_class <- predict(log_reg_final,
                      new_data = test,
                      type = "class")
results <- test %>%
  select(y) %>%
  bind_cols(pred_class, pred_proba)

# Create confusion matrix
conf_mat(results, truth = y,
         estimate = .pred_class)

Truth
Prediction   no  yes
       no  7838  738
       yes  147  320

Je kunt de precisie (positive predictive value, het aantal true positives gedeeld door het aantal voorspelde positives) berekenen met de functie precision().

precision(results, truth = y,
          estimate = .pred_class)

# A tibble: 1 × 3
  .metric   .estimator .estimate
  <chr>     <chr>          <dbl>
1 precision binary         0.914

Op dezelfde manier kun je de recall (sensitiviteit, het aantal true positives gedeeld door het aantal werkelijke positives) berekenen met de functie recall().

recall(results, truth = y,
          estimate = .pred_class)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 recall  binary         0.982

Laten we de variabelen begrijpen die de beslissing om te abonneren beïnvloeden. In een scenario met logistische regressie bepalen de coëfficiënten hoe gevoelig de doelvariabele is voor individuele voorspellers. Hoe hoger de waarde van de coëfficiënten, hoe groter hun belang. Sorteer de variabelen in aflopende volgorde van de absolute waarde van hun coëfficiënt en toon alleen de coëfficiënten met een absolute waarde groter dan 0,5.

coeff <- tidy(log_reg_final) %>% 
  arrange(desc(abs(estimate))) %>% 
  filter(abs(estimate) > 0.5)

# A tibble: 10 × 3
   term            estimate      penalty
   <chr>              <dbl>        <dbl>
 1 (Intercept)       -2.59  0.0000000001
 2 poutcomesuccess    2.08  0.0000000001
 3 monthmar           1.62  0.0000000001
 4 monthoct           1.08  0.0000000001
 5 monthsep           1.03  0.0000000001
 6 contactunknown    -1.01  0.0000000001
 7 monthdec           0.861 0.0000000001
 8 monthjan          -0.820 0.0000000001
 9 housingyes        -0.550 0.0000000001
10 monthnov          -0.517 0.0000000001

Plot de feature importance met de functie ggplot().

ggplot(coeff, aes(x = term, y = estimate, fill = term)) + geom_col() + coord_flip()

Ga een stap verder

Hiermee zijn we aan het einde gekomen van deze tutorial, waarin we hebben laten zien hoe je met het tidymodels-package een model voor logistische regressie traint en evalueert. We hebben ook behandeld hoe je de modelresultaten interpreteert en ze visualiseert, bijvoorbeeld als feature importance. 

Ben je benieuwd om meer te leren over het modelleren van AI-modellen met tidymodels, bekijk dan de cursus Modeling with tidymodels in R. Wil je modelleren met de base-R-aanpak, bekijk dan de cursussen Introduction to Regression in R, Intermediate Regression in R en Generalized Linear Models in R.