Peut-on prédire le niveau des dépassements d’honoraires ?

Résumé

En mobilisant les données de l’Assurance maladie, nous cherchons à modéliser les dépassements d’honoraires des professionnels de santé en France en 2018. Le meilleur algorithme parvient à prédire le niveau de ces dépassements avec une erreur moyenne d’environ 1150 euros.

Sommaire

• Objectif de l’étude
• Données
• Transformation du jeu de données
• Modélisation
• Evaluation de la modélisation
• Conclusion

Objectif de l’étude

Cette étude vise à modéliser le niveau des dépassements d’honoraires des professionnels de santé en France. Ce modèle peut être utile à l’Assurance maladie pour anticiper les dépassements d’honoraires, pour détecter des tendances ou contribuer à repérer des anomalies dans l’évolution des dépassements.

Données

Nous recourons aux données mises à disposition par l’Assurance maladie sur le site d’ouverture des données publiques https://www.data.gouv.fr

Au moment de la réalisation de cette étude, les données les plus récentes sont celles du mois d’août 2018. Toutes les données sont disponibles à l’adresse suivante, ainsi que le dictionnaire des variables : https://www.data.gouv.fr/fr/datasets/depenses-d-assurance-maladie-hors-prestations-hospitalieres-donnees-nationales/

Le code qui suit télécharge les extensions requises ainsi que les données, puis prépare ces dernières en modifiant la nature de certaines variables et en harmonisant la syntaxe de séparation des décimales, la rendant compatible avec R. De plus, certaines colonnes sont supprimées :

• certaines constituent de fait, sous forme de codage, des doublons avec d’autres colonnes : statistique mensuelle, etc.
• certaines contiennent de nombreux niveaux, nous préférons conserver les colonnes qui agrègent ces niveaux : spécialité prescripteur, etc.
• la colonne remboursement final est abandonnée, nous préférons conserver celle qui correspond à la base du remboursement

Dès lors qu’un niveau plus agrégé a été retenu pour certaines variables, il apparait qu’environ 10 % des observations sont en double ; elles sont donc retirées du jeu de données, qui contient au final environ 185 000 observations et 15 colonnes, dont la variable à prédire (dépassements d’honoraires, “dep_mon”).

La visualisation des données est rendue difficile par la grande dispersion de la distribution de certaines d’entre elles, en particulier des dépassements d’honoraires, dont la médiane est 0 tandis que la moyenne est de 1789 euros et l’écart type de 108 262 euros ; le maximum s’élève à 31 524 586 euros, le minimum à -3900 euros, soit un écart extrêmement élevé qui rendra difficile la prédiction. Pour 83 % des observations, aucun dépassement d’honoraire n’est constaté sur ce mois d’août 2018. A noter que la présence de valeurs négatives dans le jeu de données pour plusieurs colonnes s’explique par des mesures de régularisation pour des soins anciens (source : discussion relative à ce jeu de données sur le site data.gouv.fr) ; il convient de conserver ces valeurs négatives puisque de telles régularisations sont opérées chaque mois.

# required packages
library(readr)
library(caret)
library(h2o)

August2018 <- read_delim("n201808.csv", ";", escape_double = FALSE, locale = locale(encoding = "ISO-8859-1"), trim_ws = TRUE)

August2018 <- data.frame(August2018)

# change the nature of some columns
August2018$l_serie <- factor(August2018$l_serie)
August2018$l_prs_nat <- factor(August2018$l_prs_nat)
August2018$sns_date <- factor(August2018$sns_date)
August2018$l_asu_nat <- factor(August2018$l_asu_nat)
August2018$l_cpl_cod <- factor(August2018$l_cpl_cod)
August2018$l_pre_spe1 <- factor(August2018$l_pre_spe1)
August2018$l_exe_spe1 <- factor(August2018$l_exe_spe1)
August2018$l_pre_stj1 <- factor(August2018$l_pre_stj1)
August2018$l_exe_stj1 <- factor(August2018$l_exe_stj1)
August2018$rem_date <- factor(August2018$rem_date)

# syntax
August2018$dep_mon <- gsub('.', "", fixed = TRUE, x = August2018$dep_mon)
August2018$dep_mon <- gsub(',', ".", fixed = TRUE, x = August2018$dep_mon)
August2018$dep_mon <- as.numeric(August2018$dep_mon)

August2018$rec_mon <- gsub('.', "", fixed = TRUE, x = August2018$rec_mon)
August2018$rec_mon <- gsub(',', ".", fixed = TRUE, x = August2018$rec_mon)
August2018$rec_mon <- as.numeric(August2018$rec_mon)

August2018$act_coe <- gsub('.', "", fixed = TRUE, x = August2018$act_coe)
August2018$act_coe <- gsub(',', ".", fixed = TRUE, x = August2018$act_coe)
August2018$act_coe <- as.numeric(August2018$act_coe)

# remove duplicated or useless columns
August2018$asu_nat <- NULL
August2018$prs_nat <- NULL
August2018$SERIE <- NULL
August2018$cpl_cod <- NULL
August2018$pre_stj1 <- NULL
August2018$exe_stj1 <- NULL
August2018$l_pre_spe <- NULL # many levels, use the grouping column instead (l_pre_spe1)
August2018$pre_spe <- NULL  # many levels, use the grouping column instead (l_pre_spe1)
August2018$pre_spe1 <- NULL
August2018$l_exe_spe <- NULL # many levels, use the grouping column instead (l_exe_spe1)
August2018$exe_spe <- NULL # many levels, use the grouping column instead (l_exe_spe1)
August2018$exe_spe1 <- NULL 
August2018$rem_mon <- NULL # keep instead the reimbursement base and the rate

# duplicated rows
cat("The number of duplicated rows is", nrow(August2018) - nrow(unique(August2018))) 
# There are now 20725 duplicated rows (approximately 10 % of the dataset) because we removed some columns to keep a more aggregated level (l_exe_spe1, l_pre_spe1)

August2018 <- unique(August2018) # remove the duplicated rows

# target variable
summary(August2018$dep_mon)
sd(August2018$dep_mon)
sum(August2018$dep_mon == 0) / nrow(August2018)

Transformation du jeu de données

Avant d’opérer la modélisation des dépassements d’honoraires, le jeu de données est coupée aléatoirement en deux, l’un pour entraîner les modèles (avec 75 % des observations), l’autre pour affiner leurs paramètres et réduire le sur-apprentissage (overfitting).

Les variables discrètes sont encodées de manière binaire (one hot encoding). Les variables prédictives sont transformées en les centrant et les réduisant et en appliquant la méthode de Yeo-Johnson. En outre, les nombreuses variables binaires dont la variance est très faible sont retirées des jeux de données. Il demeure 33 variables, dont la variable prédictive qui n’a subi aucune transformation.

Enfin, les noms des colonnes sont modifiées afin d’éviter les caractères spéciaux, qui rendent leur lecture difficile voire impossible dans l’interface d’h2o à laquelle nous recourerons pour modéliser.

# Split the training set into a train and a test set

        # Step 1: Get row numbers for the training data
set.seed(1)
trainRowNumbers <- createDataPartition(August2018$dep_mon, p = 0.75, list = FALSE)

        # Step 2: Create the training dataset
trainData <- August2018[trainRowNumbers,]

        # Step 3: Create the test dataset
testData <- August2018[-trainRowNumbers,]

        # Store Y for later use.
y <- trainData$dep_mon
y_test <- testData$dep_mon

# pre process the training data

        # One hot encoding
dummies_model <- dummyVars(dep_mon ~ ., data = trainData)
trainData <- predict(dummies_model, newdata = trainData) 
trainData <- data.frame(trainData)

        # Transform the training data
preProcess_yeojohnson_model <- preProcess(trainData, method = c("center", "scale", "YeoJohnson", "nzv"))
trainData <- predict(preProcess_yeojohnson_model, newdata = trainData)

        # Append the training response variable, that remains untransformed
trainData$dep_mon <- y


# pre process the test data
testData <- predict(dummies_model, newdata = testData) # One hot encoding
testData <- data.frame(testData)
testData <- predict(preProcess_yeojohnson_model, newdata = testData) # Transform
testData$dep_mon <- y_test # Append the response variable


# change the columns'names to avoid special characters
        # training set
colnames(trainData)[1:21] <- c("lserieAMIsoinsinfirmiers", "lserieCSpecialistes",
                               "lserieTotalfraisdedeplacementdesauxiliairesmedicaux", "remdate201807", 
                               "remdate201808", "snsdate201804", "snsdate201805", "snsdate201806", 
                               "snsdate201807", "snsdate201808", "lasunatAccidentduTravailMaladieProfessionnelle",
                               "lasunatMaladie", "lasunatMaternite", "lcplcodMajorationdenuit", "lcplcodSansmajoration",
                               "REMTAU", "lprespe1Autres", "lprespe1MedecinsOmnipraticiensliberaux",
                               "lprespe1MedecinsSpecialistesliberaux", "lprestj1Liberal", "lprestj1Salarie"  )

colnames(trainData)[22:33] <- c("lexespe1Autres", "lexespe1Infirmiers", "lexespe1Laboratoires",
                                "lexespe1MasseursKinesitherapeutes", "lexespe1MedecinsOmnipraticiens",
                                "lexespe1MedecinsSpecialistes", "lexestj1Liberal", "lexestj1Salarie", "recmon",
                                "actdnb", "actcoe", "depmon")

        # test set
colnames(testData)[1:21] <- c("lserieAMIsoinsinfirmiers", "lserieCSpecialistes",
                               "lserieTotalfraisdedeplacementdesauxiliairesmedicaux", "remdate201807", 
                               "remdate201808", "snsdate201804", "snsdate201805", "snsdate201806", 
                               "snsdate201807", "snsdate201808", "lasunatAccidentduTravailMaladieProfessionnelle",
                               "lasunatMaladie", "lasunatMaternite", "lcplcodMajorationdenuit", "lcplcodSansmajoration",
                               "REMTAU", "lprespe1Autres", "lprespe1MedecinsOmnipraticiensliberaux",
                               "lprespe1MedecinsSpecialistesliberaux", "lprestj1Liberal", "lprestj1Salarie")

colnames(testData)[22:33] <- c("lexespe1Autres", "lexespe1Infirmiers", "lexespe1Laboratoires",
                                "lexespe1MasseursKinesitherapeutes", "lexespe1MedecinsOmnipraticiens",
                                "lexespe1MedecinsSpecialistes", "lexestj1Liberal", "lexestj1Salarie", "recmon",
                                "actdnb", "actcoe", "depmon")

Modélisation

Nous choisissons de traiter la modélisation des dépassements d’honoraires en problème de régression car il existe une très grande diversité dans le niveau des dépassements. Autrement dit, la question est moins de savoir si un professionnel de santé recourt aux dépassements d’honoraires, que de savoir dans quelle mesure le prix qu’il pratique est supérieur au tarif de convention fixé pour chaque acte médical par la sécurité sociale.

Les algorithmes suivants sont utilisés :

• Modèle linéaire généralisé
• Forêts aléatoires
• Réseaux de neurones
• Gradient boosting

Comme métrique de performance, nous retenons l’erreur moyenne absolue (mean absolute error, MAE), dont l’interprétation est simple.

Modèle linéaire généralisé

Sans surprise, les résultats obtenus du modèle linéaire généralisé sont médiocres: la MAE s’élève à 3639 euros, certes sans fort sur-apprentissage. Des modèles plus complexes sont donc nécessaires.

h2o.init(nthreads = -1, max_mem_size = '4g')

training_data_yeojohnson_h2o <- as.h2o(trainData)
testing_data_yeojohnson_h2o <- as.h2o(testData)

Y <- "depmon"
X <- setdiff(names(training_data_yeojohnson_h2o), Y)

# train the model
glm1 <- h2o.glm(y = Y,
                x  =  X,
                training_frame = training_data_yeojohnson_h2o,
                family = "gaussian",
                nfolds = 5, 
                alpha = 0.5,
                seed = 1)
glm1
# 3229 on both training and cross validation data

# train the model
glm2 <- h2o.glm(y = Y,
                x  =  X,
                training_frame = training_data_yeojohnson_h2o,
                family = "gaussian",
                nfolds = 5,
                seed = 1)
glm2
# same

glm2.predicted <- h2o.predict(glm2, newdata = testing_data_yeojohnson_h2o)
glm2.err <- mean(abs(testing_data_yeojohnson_h2o$depmon- glm2.predicted))
glm2.err # MAE on test set : 3639.823

Forêts aléatoires

Les différentes forêts aléatoires apportent de meilleurs résultats, en divisant par deux l’erreur moyenne du modèle linéraire généralisé, qui tombe à environ 1800 euros, sans sur-apprentissage d’après la validation croisée (cross validation). La manipulation des hyperparamètres (nombre d’arbres, etc.) n’a pas permis d’améliorer l’erreur moyenne.

h2o.init(nthreads = -1, max_mem_size = '4g')

# prepare the datasets
training_data_yeojohnson_h2o <- as.h2o(trainData)
testing_data_yeojohnson_h2o <- as.h2o(testData)

Y <- "depmon"
X <- setdiff(names(training_data_yeojohnson_h2o), Y)

# train the model
set.seed(1)
drf1 <- h2o.randomForest(training_frame = training_data_yeojohnson_h2o,
                nfolds = 3,
                y = Y,
                x = X,
                ntrees = 50,
                max_depth = 15,
                min_rows = 2,
                mtries = -1,
                seed = 1)
drf1 
# MAE : 2019.645 on cross validation data

# change the parameters
drf2 <- h2o.randomForest(training_frame = training_data_yeojohnson_h2o,
                nfolds = 5,
                y = Y,
                x = X,
                ntrees = 100,
                max_depth = 15,
                min_rows = 2,
                mtries = -1,
                seed = 2)
drf2
# MAE : 1863.483 on cross validation data, 1818.103 on training

Réseaux de neurones

Plusieurs réseaux de neurones relativement simples ont été essayés à travers trois fonctions d’activations différentes (Tanh, Maxout, ReLU), des architectures plus ou moins profondes, etc. Les meilleurs réseaux sont parvenus à diminuer l’erreur moyenne à moins de 1700 euros. Des tentatives (non apparentes ici) de régularisation n’ont pas permis de réduire le sur-apprentissage, patent au regard de la validation croisée.

h2o.init(nthreads = -1, max_mem_size = '4g')

training_data_yeojohnson_h2o <- as.h2o(trainData)
testing_data_yeojohnson_h2o <- as.h2o(testData)

Y <- "depmon"
X <- setdiff(names(training_data_yeojohnson_h2o), Y)

# train the model
DL.h2o.1 <- h2o.deeplearning(x = X, 
                             y = Y,
                             training_frame = training_data_yeojohnson_h2o,
                             nfolds = 5, 
                             distribution = "gaussian",
                             activation = "Maxout",
                             hidden = c(32,32,32),
                             sparse = TRUE, 
                             l1 = 0,
                             epochs = 100, 
                             variable_importances = TRUE,
                             standardize = FALSE,
                             seed = 1)
DL.h2o.1
# MAE : 1402 on training and 2005 on cross validation data

DL.h2o.2 <- h2o.deeplearning(x = X, 
                             y = Y,
                             training_frame = training_data_yeojohnson_h2o,
                             nfolds = 3, 
                             distribution = "gaussian",
                             activation = "Tanh",
                             hidden = c(32,32,32),
                             sparse = TRUE, 
                             l1 = 0,
                             epochs = 100, 
                             variable_importances = TRUE,
                             standardize = FALSE,
                             seed = 2)
DL.h2o.2
# MAE : 1060.002 on training and 1690.504 on cross validation data

DL.h2o.3 <- h2o.deeplearning(x = X, 
                             y = Y,
                             training_frame = training_data_yeojohnson_h2o,
                             nfolds = 3, 
                             distribution = "gaussian",
                             activation = "Rectifier",
                             hidden = c(32,32,32),
                             sparse = TRUE, 
                             l1 = 0,
                             epochs = 100, 
                             variable_importances = TRUE,
                             standardize = FALSE,
                             seed = 2)
DL.h2o.3
# MAE : 1360.989 on training and 1794.588 on cross validation data

DL.h2o.4 <- h2o.deeplearning(x = X, 
                             y = Y,
                             training_frame = training_data_yeojohnson_h2o,
                             nfolds = 3, 
                             distribution = "gaussian",
                             activation = "Tanh",
                             hidden = c(32,32,32,32,32),
                             sparse = TRUE, 
                             l1 = 0,
                             epochs = 150, 
                             variable_importances = TRUE,
                             standardize = FALSE,
                             seed = 4)
DL.h2o.4
# MAE : 1691.333 on cross validation data

DL.h2o.5 <- h2o.deeplearning(x = X, 
                             y = Y,
                             training_frame = training_data_yeojohnson_h2o,
                             nfolds = 3, 
                             distribution = "gaussian",
                             activation = "Tanh",
                             hidden = c(100,100),
                             sparse = TRUE, 
                             l1 = 0,
                             epochs = 100, 
                             variable_importances = TRUE,
                             standardize = FALSE,
                             seed = 5)
DL.h2o.5
# MAE:  474.2974 on training and 1684.774 on cross validation data

Gradient boosting

Enfin, la méthode ensembliste du gradient boosting est mobilisée et donne des résultats analogues aux meilleurs réseaux de neurones, avec une MAE d’environ 1700 euros. Le sur-apprentissage est de nouveau élevé, malgré des tentatives de régularisation.

h2o.init(nthreads = -1, max_mem_size = '4g')

# prepare the datasets
training_data_yeojohnson_h2o <- as.h2o(trainData)
testing_data_yeojohnson_h2o <- as.h2o(testData)

Y <- "depmon"
X <- setdiff(names(training_data_yeojohnson_h2o), Y)

# train the model
set.seed(5)
gbm.h2o.1 <- h2o.gbm(training_frame = training_data_yeojohnson_h2o,
                    x = X,
                    y = Y,
                    distribution = "gaussian",
                    ntrees = 50,
                    max_depth = 10,
                    learn_rate = 0.2,
                    min_rows = 2,
                    nfolds = 5)
gbm.h2o.1
# MAE : 483 on training ; 2023 on cross validation data

gbm.h2o.2 <- h2o.gbm(training_frame = training_data_yeojohnson_h2o,
                    x = X,
                    y = Y,
                    distribution = "gaussian",
                    ntrees = 50,
                    max_depth = 10,
                    learn_rate = 0.1,
                    min_rows = 2,
                    nfolds = 5)
gbm.h2o.2
# 593 on training ; 1886.904 on cross validation data

gbm.h2o.3 <- h2o.gbm(training_frame = training_data_yeojohnson_h2o,
                    x = X,
                    y = Y,
                    distribution = "gaussian",
                    ntrees = 50,
                    max_depth = 15,
                    learn_rate = 0.1,
                    min_rows = 2,
                    nfolds = 5)
gbm.h2o.3
# 288.4684 on training ; 1855.77 on cross validation data

gbm.h2o.4 <- h2o.gbm(training_frame = training_data_yeojohnson_h2o,
                    x = X,
                    y = Y,
                    distribution = "gaussian",
                    ntrees = 30,
                    max_depth = 15,
                    learn_rate = 0.1,
                    min_rows = 2,
                    nfolds = 5)
gbm.h2o.4
# 463.5689 on training ; 1756.26 on cross validation data

gbm.h2o.5 <- h2o.gbm(training_frame = training_data_yeojohnson_h2o,
                    x = X,
                    y = Y,
                    distribution = "gaussian",
                    ntrees = 30,
                    max_depth = 15,
                    learn_rate = 0.1,
                    min_rows = 2,
                    col_sample_rate = 0.8,
                    nfolds = 5)
gbm.h2o.5
# MAE : 1706.235 on cross validation data

gbm5.predicted <- h2o.predict(gbm.h2o.5, newdata = testing_data_yeojohnson_h2o)
gbm5.err <- mean(abs(testing_data_yeojohnson_h2o$depmon - gbm5.predicted))
gbm5.err # MAE : 2141.069 on test set

Entrainement final du meilleur modèle

Nous choisissons le dernier modèle ensembliste (gbm.h2o.5) comme meilleur modèle dans la mesure où il parvient à des résultats très comparables aux meilleurs réseaux de neurones mais que son temps de calcul est bien plus rapide.

Aussi, nous pouvons désormais rassembler en un jeu de données unique les deux jeux qui nous ont permis de choisir ce meilleur modèle, et entraîner le modèle sur ce grand jeu de données du mois d’août 2018 afin d’améliorer sa puissance prédictive.

h2o.init(nthreads = -1, max_mem_size = '4g')

# merge the datasets
TrainingAugust2018 <- rbind.data.frame(trainData, testData)

# prepare the full dataset
TrainingAugust2018_h2o <- as.h2o(TrainingAugust2018)
Y <- "depmon"
X <- setdiff(names(TrainingAugust2018_h2o), Y)

# train the final model
gbm.h2o.final <- h2o.gbm(training_frame = TrainingAugust2018_h2o,
                    x = X,
                    y = Y,
                    distribution = "gaussian",
                    ntrees = 30,
                    max_depth = 15,
                    learn_rate = 0.1,
                    min_rows = 2,
                    col_sample_rate = 0.8)

h2o.varimp_plot(gbm.h2o.final, 32) # most important variables to predict the response

Evaluation de la modélisation

Enfin, nous utilisons notre meilleur modèle, désormais entraîné sur un grand jeu de données, et l’appliquons sur un nouveau jeu de données, celui du mois d’août 2017, disponible sur le site data.gouv.fr, afin d’évaluer la capacité de notre modèle à se généraliser. Nous retenons le même mois car il est possible qu’existent certaines régularités mensuelles chaque année.

Le code suivant prépare les données du mois d’août 2017 pour qu’elles soient dans le même format que celles du mois d’août 2018, sans quoi le modèle ne trouverait pas à s’appliquer puisque les variables seraient différentes.

Etant données les observations des variables prédictives, notre modèle parvient à prédire le niveau des dépassements d’honoraires avec une erreur moyenne de 1156 euros par rapport aux dépassements réellement constatés.

# load the new test dataset 
August2017 <- read_delim("N201708.csv", ";", escape_double = FALSE, locale = locale(encoding = "ISO-8859-1"), trim_ws = TRUE)
August2017 <- data.frame(August2017)

# remove duplicated or useless columns
August2017$asu_nat <- NULL
August2017$prs_nat <- NULL
August2017$SERIE <- NULL
August2017$cpl_cod <- NULL
August2017$pre_stj1 <- NULL
August2017$exe_stj1 <- NULL
August2017$l_pre_spe <- NULL 
August2017$pre_spe <- NULL 
August2017$pre_spe1 <- NULL
August2017$l_exe_spe <- NULL 
August2017$exe_spe <- NULL 
August2017$exe_spe1 <- NULL 
August2017$rem_mon <- NULL 

# change the nature of columns
August2017$l_serie <- factor(August2017$l_serie)
August2017$l_prs_nat <- factor(August2017$l_prs_nat)
August2017$sns_date <- factor(August2017$sns_date)
August2017$l_asu_nat <- factor(August2017$l_asu_nat)
August2017$l_cpl_cod <- factor(August2017$l_cpl_cod)
August2017$l_pre_spe1 <- factor(August2017$l_pre_spe1)
August2017$l_exe_spe1 <- factor(August2017$l_exe_spe1)
August2017$l_pre_stj1 <- factor(August2017$l_pre_stj1)
August2017$l_exe_stj1 <- factor(August2017$l_exe_stj1)
August2017$rem_date <- factor(August2017$rem_date)

# syntax
August2017$dep_mon <- gsub('.', "", fixed = TRUE, x = August2017$dep_mon)
August2017$dep_mon <- gsub(',', ".", fixed = TRUE, x = August2017$dep_mon)
August2017$dep_mon <- as.numeric(August2017$dep_mon)

August2017$rec_mon <- gsub('.', "", fixed = TRUE, x = August2017$rec_mon)
August2017$rec_mon <- gsub(',', ".", fixed = TRUE, x = August2017$rec_mon)
August2017$rec_mon <- as.numeric(August2017$rec_mon)

August2017$act_coe <- gsub('.', "", fixed = TRUE, x = August2017$act_coe)
August2017$act_coe <- gsub(',', ".", fixed = TRUE, x = August2017$act_coe)
August2017$act_coe <- as.numeric(August2017$act_coe)

# store the target variable for later use
y_test <- August2017$dep_mon

# pre process the new test data
        # dummy variables
dummies_model <- dummyVars(dep_mon ~ ., data = August2017)
August2017 <- predict(dummies_model, newdata = August2017)
August2017 <- data.frame(August2017)

        # transformation
preProcess_yeojohnson_model <- preProcess(August2017, method = c("center", "scale", "YeoJohnson", "nzv"))
August2017 <- predict(preProcess_yeojohnson_model, newdata = August2017)

        # append the target, that remains unchanged
August2017$dep_mon <- y_test

str(August2017) 
# the variables after this pre-processing are exactly the same as the ones used during training
# this is a good sign that the model is going to generalize well

        # change the columns' names
                # note : we keep the same names as the ones used during training, even for the dates
                # this is not a problem since the variables are the same
colnames(August2017)[1:21] <- c("lserieAMIsoinsinfirmiers", "lserieCSpecialistes",
                               "lserieTotalfraisdedeplacementdesauxiliairesmedicaux", "remdate201807", 
                               "remdate201808", "snsdate201804", "snsdate201805", "snsdate201806", 
                               "snsdate201807", "snsdate201808", "lasunatAccidentduTravailMaladieProfessionnelle",
                               "lasunatMaladie", "lasunatMaternite", "lcplcodMajorationdenuit", "lcplcodSansmajoration",
                               "REMTAU", "lprespe1Autres", "lprespe1MedecinsOmnipraticiensliberaux",
                               "lprespe1MedecinsSpecialistesliberaux", "lprestj1Liberal", "lprestj1Salarie"  )

colnames(August2017)[22:33] <- c("lexespe1Autres", "lexespe1Infirmiers", "lexespe1Laboratoires",
                                "lexespe1MasseursKinesitherapeutes", "lexespe1MedecinsOmnipraticiens",
                                "lexespe1MedecinsSpecialistes", "lexestj1Liberal", "lexestj1Salarie", "recmon",
                                "actdnb", "actcoe", "depmon")

# h2o format
August2017_h2o <- as.h2o(August2017)

# use the gbm model that was trained on the entire August2018
gbm.final.predictions <- h2o.predict(gbm.h2o.final, newdata = August2017_h2o)
gbm.err <- mean(abs(August2017_h2o$depmon - gbm.final.predictions))
gbm.err # new test MAE : 1156

gbm.final.predictions_df <- as.data.frame(gbm.final.predictions)

# summary of the predictions
summary(gbm.final.predictions_df) 
# 1st Quartile : 75 ; median : 77 ; 3rd Quartile : 103
# mean : 1841 
# minimum : -576188
# maximum : 19092275

Conclusion

Notre meilleur algorithme a modélisé les dépassements d’honoraires des professionnels de santé et son évaluation sur les données d’août 2017 fait apparaître une erreur moyenne de 1156 euros. Autrement dit, pour chacune des observations du jeu de données d’août 2017, la valeur prédite par l’algorithme se trompe, en moyenne, de 1156 euros. Comment interpréter ce résultat ?

D’un côté, cette erreur moyenne peut sembler élevée dans la mesure où il n’existe en fait aucun dépassement d’honoraires pour de très nombreux actes médicaux.

D’un autre côté, cette erreur moyenne de 1156 euros est faible au regard de la très grande dispersion des observations, comprises entre -3900 euros et plus de 31 000 000 euros, ou de l’écart type, d’environ 108 000 euros.

Par ailleurs, l’algorithme a mis en évidence les variables les plus importantes pour prédire le niveau, qu’il soit négatif, nul, ou élevé, des dépassements d’honoraires. Il s’agit principalement du montant de base du remboursement, du nombre d’actes médicaux pratiqués, du taux de remboursement, du coefficient global. De plus, sans surprise, la spécialité du prescripteur en “médecin spécialiste libéral” ou en “médecin omnipraticien libéral”, le mode d’exercice des exécutants en “médecin spécialiste” et le statut “libéral” de l’exécutant sont également utiles pour prédire le niveau des dépassements. Il ne s’agit toutefois pas là d’un lien de causalité, qui ne pourrait être mis en évidence que par des méthodes économétriques, tandis que les algorithmes d’apprentissage des données recherchent plutôt des corrélations.