Régressions linéaires et bootstrap sous R

Consolidation de modèles

Comment obtenir un modèle de régression solide, validé, dont la fiabilité serait connue avec un risque minimalisé de sur-interprétation

L'essentiel de cette page !

Lorsqu'on fait un modèle de régression, il est pertinent d'établir la confiance que l'on peut avoir dans son modèle (pour éviter surinterprétation et sous-interprétation) en estimant dans quelles mesures mon échantillon de donner va le devier. Pour cela on fait de l'échantillonnage aléatoire dans ses propres valeurs.

Dans un deuxième temps, on va pouvoir ainsi déterminer sur un échantillon test, l'erreur moyenne de mon modèle mais aussi le modèle optimisé pour éviter sa déformation par effet de levier.

Pour ne pas avoir à lire toute la page dessous, une fonction permet de valider automatiquement un modèle de régression par bootstrap : bootreg() du package {KefiR}.

  • Installer valreg() (package {KefiR}

install.packages("devtools") ; require(devtools) # Risque d'erreur si RTools non installé.

devtools::install_github("Antoine-Masse/KefiR")

  • Exécuter le package {KefiR}

library("KefiR")

  • Exemple d'utilsation de valreg()

bootreg(reg, plot=T, analyse=T)

  • plot : permet d'activer le bilan graphique.

  • analyse : permet d'obtenir en sortie le tableau de la distribution des p-values et des coefficients

1- Ouvrons un jeu de données (c'est toujours la première étape) - cliquer sur le titre pour afficher le code

Ouvrons le jeu de données indiqué ici - suivre le lien.

Ce jeu de données illustre les analyses sanguines et d'un patient et l'évaluation de son état de santé.

L'idéal serait de pouvoir prédire l'état de Santé d'une personne en fonction de ses futurs analyses.

library("openxlsx")

data <- read.xlsx(file.choose())

data <- data[c(-1,-2),]

colnames(data)

View(data)

plot(data[,c(2,3,6,7,8,15)])

On va essayer de construire le modèle "parfait" pour prédire la santé du patient.

Attention : l'objectif est d'éviter de modéliser le "bruit". Notre modèle ne doit pas être capable de prédire l'ensemble des résultats actuels car certains sont peut-être faux... Il doit avant tout être juste...

2- Création d'un modèle de régression global

# Faire une régression

global_reg <- lm(Santé~VGM+TCMH+Plaquettes+Leucocytes+Neutrophiles,data= data)

# Graphique

plot.new()

par(mar = c(4, 8, 2, 2)) #par(mar=c(marge basse, marge gauche, marge haute, marge droite))

barplot(global_reg$coefficients[-1],horiz=T,col="#A8372A",las=2)


summary(global_reg)


Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -64.22858 21.36586 -3.006 0.0169 *

VGM 0.33517 0.21641 1.549 0.1600

TCMH 0.96400 1.02844 0.937 0.3760

Plaquettes 0.02173 0.01761 1.234 0.2523

Leucocytes -1.17781 0.66893 -1.761 0.1163

Neutrophiles 2.05687 0.74348 2.767 0.0244 *

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

On pourrait décrire ce graphique en disant que les paramètres qui influencent le plus la santé sont les neutrophiles et les leucocytes.

Il n'en est rien car chaque jeu de valeurs ne se mesure pas avec les mêmes unités.

3- Refaire la régression en centrant et réduisant les données d'abord

Oui mais moi, je veux savoir ce qui influence le plus la santé du patient ! Alors je fais quoi ?

Bien, déjà je peux centrer réduire les valeurs...

centration <- function(x) { # Voici une fonction pour centrer-réduire

temp <- x

temp <- (temp-mean(temp,na.rm=T))/sd(temp,na.rm=T)

return(temp)}

data <- data.frame(apply(data,2,centration))

Et on recommence la partie 1...

Bon, ça semble se confirmer que la santé du patient dépend des neutrophiles et leucocytes, l'un ayant un impact positif et l'autre un impact négatif.

Plaquettes, TCMG et VGM auraient un effet moindre.

Et bien pas nécessairement !

Il faut d'abord vérifier si certains paramètres ne mériterait pas d'être éliminés car leur impact bruité contribuerait à faire de l'overfitting, du surapprentissage/surinterprétation.

Exemple : si j'apprends à prédire le sexe des garçons et filles à partir d'un jeu de photos d'apprentissage, si j'ai des photos de blondes et pas de blonds, mon modèle risque de croire que tout ce qui est blond est une fille...

On peut identifier les paramètres pertinents en regardant valeur Pr. Si cette valeur est supérieure à 0.05, c'est que le paramètre n'a a priori pas d'effet significatif sur la Santé.

Cela dépend toutefois de mes choix de paramètres et des valeurs à effet levier qui peuvent fausser totalement mon modèle.

On peut éviter ces problèmes en faisant du bootstrap.

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -3.762e-15 1.310e-01 0.000 1.0000

VGM 3.698e-01 2.388e-01 1.549 0.1600

TCMH 3.031e-01 3.234e-01 0.937 0.3760

Plaquettes 2.264e-01 1.835e-01 1.234 0.2523

Leucocytes -1.103e+00 6.263e-01 -1.761 0.1163

Neutrophiles 1.716e+00 6.201e-01 2.767 0.0244 *

4- Trouver le meilleur modèle par bootstrap

Le bootstrap constite à fabriquer deux paquets :

  • Échantillon d'apprentissage. Je construits mon modèle à partir d'un lot d'individus. Je construits ce lot en prélevant dans mes individus autant d'individus que j'ai déjà. On fait du tirage avec remise. Certains individus peuvent être tirés plusieurs fois.

  • Echantillon test. Les individus non tirés au sort (s'il y en a) serviront à vérifier la fiabilité du modèle.

Ce que l'on fait par la suite :

  • On fabrique un grand nombre de modèles à partir de couples d'échantillons apprentissage/test différents.

  • On va ainsi avoir plusieurs valeur Pr pour chaque paramètre (pour chaque variable) et ainsi estime si son effet sur la Santé semble solide ou fragile.

  • On va avoir plusieurs valeurs de chaque coefficient directeur et ainsi récupérer la valeur la plus représentative.

4.1- Mon bootstrap.

indices <- c(1:nrow(data)) # numérotation des individus

enregistrement <- 0

erreur <- c() ; err_prediction <- c() ; err_verity <- c()

for (i in c(1:1000)) {

indices_training <- sample(indices,size=nrow(data),replace=T)

indices_test <- setdiff(indices ,indices_training)

if (length(indices_test) > 0) {

enregistrement <- enregistrement + 1

training <- data[indices_training,]

test <- data[indices_test,]

reg <- lm(Santé~VGM+TCMH+Plaquettes+Leucocytes+Neutrophiles,data= training) # Valeur à adapter manuellement pour proposer une structure à mon modèle : quelles variables j'utilise !

if (enregistrement == 1) {

p_values <- summary(reg)[[4]][,4] # p_values correspond aux valeurs Pr récupérés pour chaque variable et chaque cycle du bootstrap (ici 1000)

# Pour extraction des p-values sur glm() : mettre cela : coef(summary(model2))[,4]

coeff <- reg$coefficients

} else {

p_values <- rbind(p_values,summary(reg)[[4]][,4])

coeff <- rbind(coeff,reg$coefficients) # Coeff correspond aux coefficients récupérés pour chaque variable et chaque cycle du bootstrap (ici 1000)

}

# Test

err_prediction <- c(err_prediction,predict(reg,test) )

err_verity <- c(err_verity,test$Santé)

erreur <- c(erreur,(predict(reg,test)-test$Santé)/test$Santé)

}}

# erreur correspond à l'erreur faite pour chaque prédiction sur l'échantillon test

4.2- Trouver les modèles les plus fiables.

Une petite fonction utile boxplot_Pr() :

boxplot_Pr <- function(x) {

my_min <- min(c(apply(x,2,quantile)[2,]),0.0009)

boxplot(x,log="y",ylim=c(my_min,1))

abline(h=0.05,col="red",lwd=2)

abline(h=0.01,col="orange",lwd=2)

abline(h=0.001,col="green",lwd=2)

}

boxplot_Pr(p_values)


Ce graphique indiquer que les Neutrophiles ont bien un effet significatif.

TCMH et Plaquettes n'ont pas f'effet pertinent.

En revanche, VGM et Leucocytes sont sur la marge, ne les éliminons pas tout de suite, éliminons d'abord TCMH et Plaquettes.

Le trait rouge indique le seuil de significativité à 0.05, orange à 0.01 et vert à 0.001

Recommençons le bootstrap en mettant cette fois-ci le modèle suivant :

1) # Ligne à mettre à la place de celle en gras dans le code partie 4.1

reg <- lm(Santé~VGM+Leucocytes+Neutrophiles,data= training)

boxplot_Pr(p_values)

Cette fois-ci les leucocytes semblent à éliminer ainsi que la valeur à l'intersection (qui s'élimine du fait que l'on a centré-réduit nos valeurs)

2) # Ligne à mettre à la place de celle en gras dans le code partie 4.1

reg <- lm(Santé~VGM+Neutrophiles+0,data= training)

boxplot_Pr(p_values)

Ca y est, on a trouvé 2 paramètres très significatifs !

VGM et Neutrophiles agissent bien sur la Santé de mon patient. J'ai donc un modèle simple qui m'évite la surinterprétation.

Comment trancher si une de mes variables à ses p-values (valeur Pr) qui dépassent le seuil de significativité (ce qui est souvent le cas pour quelques valeurs) ?

Il suffit de calculer un intervalle de confiance ! Oui... Mais les p-values ne sont pas réparties selon la loi normale...

Alors il y a plus simple, je prends la valeur dans la tranche 99% en partant de la plus basse... et j'aurais mon intervalle à 99%.

La fonction qui le fera :

confiance <- function(x,seuil=0.99) { # seuil

temp = sort(x) ; valeur_seuil = round(length(x)*seuil)

temp <- temp[valeur_seuil]

return(temp)

}

confiance(p_values[,1],0.99) # Valeur max du coefficient directeur de VGM (confiance à 99%)

confiance(p_values[,2],0.99) # Valeur max du coefficient directeur de Neutrophiles (confiance à 99%)

Ainsi, Neutrophiles à un effet significatif avec une valeur Pr maximale à 99% de 0.005.

En revanche, le VGM a une valeur à 99% qui va à 0.08... Il faudra donc réfléchir à deux fois avant de le conserver, mais un confiance à 95% (Pr de 0.02) pourrait peut-être suffire ? Faire un choix s'impose...

Je dois maintenant voir quel taux d'erreur, implique mon modèle !

4.3- Analyser le taux d'erreur, la fiabilité de mon modèle

erreur contient ici l'erreur commise à chaque prédiction faite sur mon échantillon test durant les 1000 iterations.

Que voit-on ?

Pour bien visualiser cela, je vais utiliser quelques petites fonctions maison :

mode <- function(x) {

densite <- density(x)

mode <- densite$x[which(densite$y==max(densite$y))]

return(mode)}

  • Pour visualiser un histogramme des valeurs en superposant, densité, moyenne, médiane et mode.

hist_decrypt <- function(x,breaks=20) {

densite <- density(x) # créer une liste des densités

hist(x,freq=F,col="#AAFFAA",ylim=c(0,max(densite$y)),breaks=breaks) # il faut que freq=F

lines(densite, col = "red",lwd=3) # Superposer la ligne

abline(v=mean(x),col="red",lwd=2)

abline(v=median(x),col="orange",lwd=2)

abline(v=mode(x),col="green",lwd=2)}

Maintenant, je peux visualiser mes erreurs :

erreur <- (err_prediction-err_verity)/err_verity

hist_decrypt(erreur)

# Erreur moyenne

mean(erreur)

int.ech(erreur)

Faire une estimation de la Santé en fonction de VGM et du taux de Neutrophiles semble ainsi apporter 10% d'erreur +/- 3% avec un intervalle de confiance à 95%.

4.4- Récupérer le meilleur lot de coefficients directeurs.

Les coefficients directeurs proposés par lm dépendent énormément de certains valeurs potentiellement fausses ou extrêmes susceptibles de déséquilibrer le modèle.

Autant donc récupérer les valeurs de coefficients les plus pertinents pamis tout ceux testés durant le boostrap.

La commande hist_decrypt permet de voir cela :

  • Essayons ici de voir la répartition des valeurs du coefficient Neutrophile (coefficient 2) :

hist_decrypt(coeff[,2])

La moyenne apparaît en rouge, la médiane en orange et le mode en vert. Il apparaît ici que le mode semble plus pertinent.

On va donc a priori récupérer le mode de chaque coefficient.

vgm <- mode(coeff[,1]) # Coefficient directeur de VGM

neutro <- mode(coeff[,2]) # Coefficient directeur de Neutrophile

  • Attention : ce modèle fera plus d'"erreur" sur notre jeu de données de base car il ne modélise que peu le bruit.

  • Il sera en revanche beaucoup plus fiable pour faire des prédictions sur de nouvelles valeurs !

Voyons son taux d'erreur sur nos données actuelles :

prediction <- data$VGM * vgm + data$Neutrophile * neutro

plot(data$Santé, prediction , xlab = "Données de Santé réelles", ylab = "Données prédites",pch=16,col="orange",cex=2)

Ma prédiction est donc plutôt bonne, à part le point en bas droite, j'ai plutôt bien prédit la santé de mon patient sans faire l'erreur d'y impliquer tous les paramètres sanguins non-pertinents !