17  Estimations

Objectives

This chapter gives the code to perform the classification of individuals’ mental health given their characteristics and consumption patterns. This is a robustness check that replicates the analysis performed in Chapter 5, using the self-assessed health (SAH) instead of the declarative question on depression to define the imaginary healthy. The MHI-5 score is used here to define the imaginary healthy.

Let us load the data obtained after cleaning (see the chapter where the data is cleaned, Chapter 15).

Let us first load {tidyverse}:

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

The cleaned data:

load("../data/df_clean_sah.rda")
dim(df_clean)

[1] 5380   94

17.1 Samples

We only focus on people who state they did not experience a depression episode during the previous year. We consider two types of individuals: those whose MHI-5 score is considered low (less than or equal to 60, i.e., the first quartile of the distribution of MHI-5 scores among all respondents), and those whose score is considered high (greater than 60).

Let us create our target variable based on this. First, we need to load some packages.

library(caret)

Loading required package: lattice

Attaching package: 'caret'

The following object is masked from 'package:purrr':

    lift

library(MLmetrics)

Attaching package: 'MLmetrics'

The following objects are masked from 'package:caret':

    MAE, RMSE

The following object is masked from 'package:base':

    Recall

library(parallel) 
library(doParallel)

Loading required package: foreach

Attaching package: 'foreach'

The following objects are masked from 'package:purrr':

    accumulate, when

Loading required package: iterators

We split the dataset into two subsets. One for training the model (df_train), with 80% observations and a second for validating it (df_test) with the remaining 20% observations.

set.seed(123)
ind_build <- createDataPartition(df_clean$status, p = .8, list = FALSE)

df_train <- df_clean[ind_build,] |> as.data.frame()
df_test <- df_clean[-ind_build,] |> as.data.frame()

Here are the two levels of the target variable:

• Level 1: MHI5 > q1 (healthy) (Not_D_and_sup_Q1).
• Level 2 (positive class): MHI5 < q1 (imaginary healthy) (Not_D_and_inf_Q1)

table(df_train$status) |> prop.table()

Not_D_and_inf_Q1 Not_D_and_sup_Q1 
        0.295935         0.704065 

table(df_test$status) |> prop.table()

Not_D_and_inf_Q1 Not_D_and_sup_Q1 
        0.295814         0.704186 

Let us save the train and test samples.

save(df_train, file ="../data/out/df_train_sah.rda")
save(df_test, file ="../data/out/df_test_sah.rda")

17.2 Grid Search

We will train different classifiers. For models that require hyperparameters, we will use the {caret} package to perform a grid search to find the optimal parameters. We will perform repeated 5-fold cross validation (10 repetitions) to find the hyperparameters. We will optimize the sensitivity (true positive rate, returned by the twoClassSummary() function from {caret}).

control <- trainControl(
  method = "repeatedcv",
  repeats = 10,
  number = 5, #number of folds
  summaryFunction = twoClassSummary,
  search = "grid",
  classProbs = TRUE,
  verboseIter = TRUE,
  allowParallel = TRUE,
  savePredictions = TRUE
)

Some columns need to be removed from the estimation:

ind_remove <- which(
  colnames(df_train) %in% c("status", "id", "PERSONNE_statut")
)

The formula for the models:

formula <- as.formula(
  paste("status ~", paste(colnames(df_train[-ind_remove]), collapse = "+")))
formula

status ~ PERSONNE_pb_asthm + PERSONNE_pb_bronchit + PERSONNE_pb_infarctus + 
    PERSONNE_pb_coronair + PERSONNE_pb_hypertens + PERSONNE_pb_avc + 
    PERSONNE_pb_arthros + PERSONNE_pb_lombalgi + PERSONNE_pb_cervical + 
    PERSONNE_pb_diabet + PERSONNE_pb_allergi + PERSONNE_pb_cirrhos + 
    PERSONNE_pb_urinair + PERSONNE_age + PERSONNE_sexe + PERSONNE_couple + 
    PERSONNE_ss + PERSONNE_regime + PERSONNE_rap_pcs8 + PERSONNE_ald + 
    SOINS_ald_am + MENAGE_revucinsee + MENAGE_tu + MENAGE_nbpers + 
    ensol_2011 + MUTUELLE_assu + MUTUELLE_typcc + SOINS_remomn + 
    SOINS_remspe + SOINS_rempha + SOINS_remkin + SOINS_reminf + 
    SOINS_remden + SOINS_remmat + SOINS_remtra + SOINS_remopt + 
    SOINS_rempro + SOINS_remurg + SOINS_tmomn + SOINS_tmspe + 
    SOINS_tmpha + SOINS_tmkin + SOINS_tminf + SOINS_tmden + SOINS_tmmat + 
    SOINS_tmtra + SOINS_tmopt + SOINS_tmpro + SOINS_tmurg + SOINS_dpaomn + 
    SOINS_dpaspe + SOINS_dpapha + SOINS_dpakin + SOINS_dpainf + 
    SOINS_dpaden + SOINS_dpamat + SOINS_dpatra + SOINS_dpaopt + 
    SOINS_dpapro + SOINS_dpaurg + SOINS_pf_fromn + SOINS_pf_frspe + 
    SOINS_pf_frpha + SOINS_pf_frkin + SOINS_pf_frinf + SOINS_pf_frden + 
    SOINS_pf_frtra + SOINS_pf_frurg + SOINS_seac_omn + SOINS_seac_spe + 
    OPINION1_renonc_cons + OPINION1_renonc_dent + OPINION1_renonc_fin + 
    OPINION1_renonc_loin + OPINION1_renonc_long + QST_ct_depech + 
    QST_ct_liberte + QST_ct_apprend + QST_ct_aidecol + QST_ct_travnuit + 
    QST_ct_repet + QST_ct_lourd + QST_ct_posture + QST_ct_produit + 
    QES_association + QES_tpsami + QES_tpsasso + QES_tpscolleg + 
    QES_tpsfamil + QES_mere_etude + QES_pere_etude

17.3 Training

17.3.1 Random Forest

This section presents the code used to classify the respondents accordingly with their label, using a random forest (Breiman (2001)).

We will need the following libraries:

library(randomForest)

Warning: package 'randomForest' was built under R version 4.4.1

randomForest 4.7-1.2

Type rfNews() to see new features/changes/bug fixes.

Attaching package: 'randomForest'

The following object is masked from 'package:dplyr':

    combine

The following object is masked from 'package:ggplot2':

    margin

library(ranger)

Attaching package: 'ranger'

The following object is masked from 'package:randomForest':

    importance

We are relying on the {ranger} package to classify individuals as imaginary healthy (Not_D_and_inf_Q1) or not (Not_D_and_sup_Q1). We use the ranger() function, that performs the classification using a random forest algorithm. Some parameters are needed to estimate the model:

• the number of trees to grow (num.trees)
• the number of variables to possibly split at each node of the trees (mtry)
• the splitting rule (splitrule)
• the minimal node size (min.node.size).

The choice of these hyper-parameters can greatly affect the quality of fit. Hence, we loop over a grid of possible combinations of these hyperparameters. More specifically:

• we set the number of trees to 500 (default value in ranger())
• we let the number of variables to possibly split at each node vary between 3 to the square root of the number of features (3 to 8)
• we select the Gini index as the splitting rule
• we let the minimum number of observations in terminal nodes vary in the following set: {50, 75, 100, 150}.

The grid is created as follows:

tune_grid_rf <- expand.grid(
  mtry = seq(3, sqrt(ncol(df_train[, -ind_remove]))),
  splitrule = "gini",
  min.node.size = c(50, 75, 100, 150)
)

We iterate over this grid. At each iteration, we estimate the model and compute the sensitivity.

The following code is not run, we load the results obtained when running this code.

no_cores <- detectCores() - 1
cl <- makePSOCKcluster(no_cores)
registerDoParallel(cl)

grid_search_rf <- train(
  form = formula,
  data = df_train,
  method = "ranger",
  trControl = control,
  tuneGrid = tune_grid_rf,
  metric = "Sens",
  # Arguments passed to ranger()
  num.trees = 500,
  importance = "permutation",
  respect.unordered.factors = "partition",
  classification = T,
  # slows the process greatly
  # keep.inbag = T,
  class.weights = 1 / prop.table(table(df_train$status))
)
stopCluster(cl)

From the code above, it can be noted that:

• the variable importance mode we select is "permutation"
• we set the respect.unordered.factors argument to "partition", so that all possible bi partitions are considered for splitting
• the argument class.weights is set to the inverse proportion of our target variable in the training sample, to account for the imbalance of the target variable in the dataset.

We save the results:

dir.create("../data/out/estim/", recursive = TRUE)
save(
  grid_search_rf, tune_grid_rf, 
  file = "../data/out/estim/v3/grid_search_rf_sah.rda"
)

17.3.2 Extreme Gradient Boosting (XGBoost)

Let us now build a classifier on the same data using extreme gradient boosting (J. H. Friedman (2001)). We need to load the {xgboost} package.

library(xgboost)

Attaching package: 'xgboost'

The following object is masked from 'package:dplyr':

    slice

The grid to perform the grid search:

tune_grid_xgb <- expand.grid(
  nrounds = 500, # no. boosting iterations
  max_depth = 3:6, # max tree depth
  colsample_bytree = (1:10)/10, # pct columns used
  eta = 0.01, # learning rate
  gamma = c(0, 5, 10), # complexity parameter
  min_child_weight = c(50, 100, 150), # Minimum sum of instance weight in a child
  subsample = c(.7, .8, .9, 1) # Subsample ratio of the training instances
)

And the estimation (not run in this notebook):

cl <- makePSOCKcluster(no_cores)
registerDoParallel(cl)

set.seed(123)
grid_search_xgb <- train(
  formula,
  data = df_train,
  trControl = control,
  tuneGrid = tune_grid_xgb,
  method = "xgbTree",
  scale_pos_weight = sum(df_train$status != "Not_D_and_inf_Q1") /
    sum(df_train$status == "Not_D_and_inf_Q1"),
  metric = "Sens"
)
stopCluster(cl)
registerDoSEQ()

We save the results:

save(
  grid_search_xgb, tune_grid_xgb, 
  file = "../data/out/estim/v3/grid_search_xgb_sah.rda"
)

17.3.3 Support Vector Machine model (SVM)

Let us now build an SVM (Ben-Hur et al. (2002)) with the {e1071} package.

library(e1071)

cl <- makePSOCKcluster(no_cores)
registerDoParallel(cl)
svm_model <- train(
  formula,              # Formula
  data = df_train,      # Data
  method = "svmRadial", # SVM with Radial kernel
  trControl = control,  # Training control
  preProcess = c("center", "scale"), # Preprocessing
  tuneLength = 100      # Number of hyperparameter values to try  
)
stopCluster(cl)
registerDoSEQ()

Let us save the results:

save(svm_model, file = "../data/out/estim/v3/svm_model_sah.rda")

17.3.4 Penalized Logistic Regression Model (glmnet)

Let us now estimate a penalized logistic regression model (J. Friedman, Hastie, and Tibshirani (2010)). Again, the following code is not run in this chapter.

cl <- makePSOCKcluster(no_cores)
registerDoParallel(cl)

glmnet_model <- train(
  formula,             # Formula
  data = df_train,     # Data
  method = "glmnet",   # Penalized logistic regression
  trControl = control, # Training control
  preProcess = c("center", "scale"), # Preprocessing
  tuneLength = 100     # Number of hyperparameter values to try
)

stopCluster(cl)
registerDoSEQ()

Let us save the results:

save(glmnet_model, file = "../data/out/estim/v3/glmnet_model_sah.rda")

Ben-Hur, Asa, David Horn, Hava T. Siegelmann, and Vladimir Vapnik. 2002. “Support Vector Clustering.” J. Mach. Learn. Res. 2 (March): 125–37.

Breiman, Leo. 2001. “Random Forests.” Machine Learning 45 (1): 5–32. https://doi.org/10.1023/A:1010933404324.

Friedman, Jerome H. 2001. “Greedy function approximation: A gradient boosting machine.” The Annals of Statistics 29 (5): 1189–1232. https://doi.org/10.1214/aos/1013203451.

Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. 2010. “Regularization Paths for Generalized Linear Models via Coordinate Descent.” Journal of Statistical Software 33 (1). https://doi.org/10.18637/jss.v033.i01.