Embedded Files
L'essentiel de cette page !
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
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
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
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
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.
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.
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
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 :
Pour calculer un intervalle de confiance
Je dois aller copier-coller le code de la fonction int.ech()
Pour calculer le mode
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.
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 !
Google Sites
Report abuse