Validation de modèle de régression linéaire avec R

L'essentiel de cette page

De nombreux tests permettent de valider un modèle de régression linéaire simple ou multiple et sont expliqués dans la page supérieure à celle-ci. On peut aussi les appliquer avec ma fonction valreg() de mon package {KefiR}.

Pourtant ces tests ne sont pas toujours faciles à comprendre. Voici des exemples qui montre pour un modèle de régression linéaire peut avoir l'air bon et se révéler complètement foireux.

Très difficile de modéliser un modèle qui ne planterait que sur le test de rainbow. On retiendra que ce test permet de détecter une mauvaise linéarité.

Partons d'un exemple :

Prenons un  exemple simple, j'essaye de prédire la surface d'une pizza en fonction de son rayon. Le modèle ne sera pas linéaire puisque la surface varie en fonction du carré du rayon.

X <- rnorm(20,3,2)

Y <- X*X+3*X+rnorm(20,1,3)+ifelse(X>2,0,1)*3

plot(X,Y,cex=2,pch=16)

reg <- lm(Y~X)

summary(reg)

valreg(reg)

Cela peut avoir des conséquences minimes ou très importantes. Un modèle en théorie fiable doit en général taper juste et se tromper parfois en excès parfois en moins. Les résidus (différence entre Y réel et Y prédits) doivent donc suivre la loi normale et se répartir selon une cloche. En revanche, dans le cas d'une distribution non normale des résidus, cela implique qu'un modèle en moyenne semble parfaitement juste mais qu'il se trompe systématiquement : ou radicalement vers le haut, ou radicalement vers le bas.

Partons d'un exemple :

J'ai des données qui ne sont pas homogènes. Par exemple, j'ai mesuré la puissance fourni par des moteurs en fonction de la quantité d'énergie consommée alors qu'il y a 2 types de moteurs. L'un ayant une performance de base meilleure et l'autre plus basse.

X1 <- rnorm(10,3,4)

X2 <- rnorm(10,3,4)

Y1 <- 2*X1+rnorm(10,0.5,0.5)+4

Y2 <- 2*X2+rnorm(10,0.5,0.5)-4

X <- c(X1,X2) ; Y <- c(Y1,Y2)

plot(X,Y,cex=2,col="green",pch=16)

reg <- lm(Y~X)

abline(reg,lwd=3,col="red")

summary(reg)

plot(Y,(Y-predict(lm(Y~X),data.frame(X))))

Un test de Shapiro l'aurait vu ou valreg()

library(KefiR)

valreg(reg)

R² élevé, p-value au TOP. On a un modèle en moyenne qui tient la route, mais qui ne prédit jamais juste, ou trop haut ou trop bas. Il aurait fallu séparé les 2 types de données dans cet exemple/

On peut avoir des données où Y semble bien dépendre de X pour les valeurs basses, mais pas sur les valeurs élevées.

On aura un modèle incapable de prédire sur les valeurs hautes dans cet exemple.

Exemple de code

X <- sort(rnorm(20,10,3))

SD <- round(X/3,0)-1.5

Y <- c()

for (i in 1:length(X)) {Y[i] <- 2*X[i] + rnorm(1,0,SD[i])}

plot(X,Y,pch=16,cex=2)

reg <- lm(Y~X)

abline(reg,lwd=3,col="red")

Cela, le test de Breush-Pagan l'aurait détecté ou la fonction valreg().

library(KefiR)

valreg(reg)

On peut avoir un modèle R qui semble excellent, mais celui-ci se révélera complètement faux à cause d'un effet de levier, c'est à dire un ou plusieurs points extrêmes qui donneront l'illusion d'avoir un modèle linéaire cohérent.

Partons d'un exemple

Je fais un modèle qui prédit le poids des consommateurs en fonction du nombre de cacahuètes qu'ils ont mangé la semaine d'avant, mais sans me rendre compte que parmi les consommateurs, il y a un éléphant qui s'est glissé (le cas particulier aux données extrêmes). Je vais faire un modèle qui a l'air hyper costaud : genre, si je mange, 10 kg de cacahuète, je vais être gros (l'éléphant), si je mange entre 0 et 200 g, je vais être léger (un humain en somme).

En réalité, si je dégage l'éléphant, je ne trouverai aucune relation entre la consommation de cacahuètes et le poids des gens : aucun modèle ne pouvant prédire efficacement le poids de qqn en fonction de ce critère. Ainsi, en cas d'effet de levier, où j'écarte la donnée extrêmes (grâce à la distance de Cook ou encore au test de Dixon pour les valeurs aberrantes) si je peux me le permettre , ou bien je trouve un autre modèle. 

Cacahuètes=c(abs(rnorm(10,0),10)) ; Poids=c(rnorm(10,60,3),4000)

plot(Cacahuètes,Poids,cex=2,pch=16) 

reg <- lm(Poids~Cacahuètes)

abline(reg,lwd=3,col="red")

summary(reg)

R² > 0,9, p-value très basse. On pourrait croire à un modèle béton. Il n'en est rien. Le poids d'une personne ne dépend pas de sa consommation de cacahuètes durant la semaine précédent la mesure.

On le verrait autant avec valreg()

library(KefiR)

valreg(reg)

5. Indépendance des résidus (test de Durbin-Watson)

Et si il y avait des catégories cachées ? Qui nous donnent l'illusion d'une relation statistique ?

Partons d'un exemple (idiot)

Je fais un modèle qui met en relation la consommation de burger et la longueur des poils et ô stupeur, je trouve une corrélation très significative : plus on mange de burgers, plus les poils sont longs !

En réalité, j'ai oublié qu'il y a deux catégories, les filles et les garçons qui n'ont ni la même longueur de poils ni la même consommation de burger (on va l'admettre pour l'exemple) : cela créer une illustion, mais l'analyse des corrélations burger-poils pour chaque sexe confirmera l'absence de lien.

sexe <- rep(c("F","M"),each=40)

burger <- rep(NA,length(sexe))

burger[sexe=="F"] <- round(abs(rnorm(40,5,3)),0)

burger[sexe=="M"] <- round(abs(rnorm(40,10,3)),0)

poils <- rep(NA,length(sexe))

poils [sexe=="F"] <- round(abs(rnorm(40,1,0.5)),1)

poils [sexe=="M"] <- round(abs(rnorm(40,2,0.8)),1)

plot(poils, burger, col=c("blue","red")[as.factor(sexe)],pch=16,cex=2)

abline(lm(burger~poils),col="black",lty=3,lwd=2)

abline(lm(burger[sexe=="M"]~poils[sexe=="M"]),col="red",lty=3,lwd=2)

abline(lm(burger[sexe=="F"]~poils[sexe=="F"]),col="blue",lty=3,lwd=2)

cor.test(poils, burger)$p.value

Un test de Durbin-Watson en triant les valeurs par sexe montrera un bon entre les résidus du des filles et les résidus des filles : on peut donc prédire la taille des résidus en fonction du sexe : le modèle doit être invalidé.

On arrivera aussi à la même démarche avec valreg() de {KefiR} :

library(KefiR)

valreg(lm(burger~poils))