Un rapide billet pour répondre à une question posée à la fin du cours de ce matin (ST5100), par Jean-Pierre Liégeois, jeune lecteur du Var (pour préserver un peu d’anonymat)
Dans mon stage, quand on avait des valeurs manquantes, on me disait de remplacer par -1, puis de rajouter une indicatrice comme quoi la variable vaut -1. Ça permet de ne supprimer ni variable, ni observations. On peut faire ça ?
Si je formalise un peu, on va simuler ici des données, disons x_1 et x_2, on génère ensuite des données suivant un modèle, de la forme y=\beta_0+\beta_1x_1+\beta_2x_2+\varepsilon. Une proportion \alpha de x_1 sera transformée en NA. Ce que suggérais Jean-Pierre, c’est de remplacer les valeurs manquantes par -1, puis d’ajuster un modèle y=\beta_0+\beta_1x_1+\beta_{-1}\mathbf{1}(x_1=-1)+\beta_2x_2+\varepsilon. Côté code, c’est assez simple. Par défaut, la stratégie de R est de supprimer les valeurs manquantes. Si 50% des données de x_1 sont manquantes, la moitié des lignes sont supprimées
n=1000 x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.05 indice=sample(1:n,size=round(n*alpha)) base=data.frame(y=y,x1=x1) base$x1[indice]=NA reg=lm(y~x1+x2,data=base) |
Au lieu de générer un unique échantillon, on va en simuler 10,000, et regarder la distribution de \widehat{\beta}_1,
m=10000 B=rep(NA,m) for(s in 1:m){ x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.5 indice=sample(1:n,size=round(n*alpha)) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=NA reg=lm(y~x1+x2,data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white",xlab="missing values = 50%") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") |
Bien sur, avec un taux de valeurs manquantes plus faibles – disons \alpha=5\% – on perd moins d’observations, et donc l’estimateur a une variance plus faible.
Tentons maintenant la stratégie consistant à remplacer les valeurs manquantes par des valeurs numériques fixes, et de rajouter une indicatrice,
B=rep(NA,m) for(s in 1:m){ x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.5 indice=sample(1:n,size=round(n*alpha)) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=-1 reg=lm(y~x1+x2+I(x1==(-1)),data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") |
Ce qui ne change pas grand chose, on en conviendra… y compris si le taux de valeurs manquantes passe à 5%,
On peut se demander ce qui se passe si le shift n’est plus de 1 mais de 10 (a priori, c’est arbitraire, la variable x_1 pouvant être plus ou moins dispersée… -1 pour une variable entre 0 et 1, ou entre 0 et 1000, ça n’est pas tout à fait pareil). Mais non, par exemple avec toujours 5% de valeurs manquantes, on a
Si on regarde notre échantillon, en particulier le nuage de points (x_1,y), on observe
ici, les valeurs manquantes sont choisies au hasard, de manière totalement indépendante,
x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.3333333 indice=sample(1:n,size=round(n*alpha)) clr=rep("black",n) clr[indice]="red" plot(x1,y,col=clr) |
(ici avec 1/3 de valeurs manquantes, en rouge). Mais on pourrait supposer que les valeurs manquantes sont les plus grandes valeurs de x_1, par exemple,
x1=runif(n) x2=runif(n) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.3333333 indice=sample(1:n,size=round(n*alpha),prob = x1^3) clr=rep("black",n) clr[indice]="red" plot(x1,y,col=clr) |
On peut se demander ce que ça donnerait sur l’estimateur \widehat{\beta}_1
Ça ne change pas grand chose, mais on a plus de variance, si on regarde bien. Dernier essai : que se passe-t-il si les variables x_1 et x_2 sont maintenant corrélées,
B=rep(NA,m) library(mnormt) r=.8 S = matrix(c(1,r,r,1),2,2) for(s in 1:m){ x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=-1 reg=lm(y~x1+x2+I(x1==(-1)),data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") |
Cette fois, on a un estimateur biaisé (de l’ordre de 10% sur cet exemple numérique). Manifestement, cette technique n’est pas très concluante…
Je pourrais ajouter que cette méthode ne revient pas à la première, même si la distribution des estimateurs est proche
set.seed(1) x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=-1 reg1=lm(y~x1+x2+I(x1==(-1)),data=base) coefficients(reg1) (Intercept) x1 x2 I(x1 == (-1))TRUE 1.0988005 1.7454385 -0.5149477 3.1000668 base$x1[indice]=NA reg2=lm(y~x1+x2,data=base) coefficients(reg2) (Intercept) x1 x2 1.1123953 1.8612882 -0.6548206 |
Comme je le disais (lors de la discussion qui a suivi le cours) une méthode plus prometteuse est l’imputation. L’idée est de prédire une valeur pour les x_1 manquants. On pourrait être tenté de prendre \overline{x}_1, la moyenne des x_1 observés. Mais on sait que les valeurs manquantes sont justement pour les grandes valeurs de x_1, ici, donc on doit pouvoir faire mieux ! On sait aussi que x_1 et x_2 sont corrélés ici. Corrélés positivement, en plus. Autrement dit, si x_2 est grand, on sait que le x_1 (non observé) devait être grand. Le plus simple est de faire un modèle linéaire, x_1=\alpha_0+\alpha_2x_2+\eta_i, calibré sur les valeurs non-manquantes. Puis on utilise \widehat{x}_1=\widehat{\alpha}_0+\widehat{\alpha}_2x_2 pour les valeurs manquantes. C’est simpliste, mais pourquoi pas ? On estime alors le modèle sur cette nouvelle base.
for(s in 1:m){ x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=NA reg0=lm(x1~x2,data=base[-indice,]) base$x1[indice]=predict(reg0,newdata=base[indice,]) reg=lm(y~x1+x2,data=base) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") |
sur l’exemple numérique, on obtient
base$x1[indice]=NA reg0=lm(x1~x2,data=base[-indice,]) base$x1[indice]=predict(reg0,newdata=base[indice,]) reg3=lm(y~x1+x2,data=base) coefficients(reg3) (Intercept) x1 x2 1.1593298 1.8612882 -0.6320339 |
Cette méthode a au moins pu corriger du biais…
Après, si on regarde attentivement, on a exactement la même valeur qu’avec la première méthode qui consiste à supprimer les lignes avec des valeurs manquantes !
summary(reg3) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.15933 0.06649 17.435 < 2e-16 *** x1 1.86129 0.21967 8.473 < 2e-16 *** x2 -0.63203 0.20148 -3.137 0.00176 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.051 on 997 degrees of freedom Multiple R-squared: 0.1094, Adjusted R-squared: 0.1076 F-statistic: 61.23 on 2 and 997 DF, p-value: < 2.2e-16 summary(reg2) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.11240 0.06878 16.173 < 2e-16 *** x1 1.86129 0.21666 8.591 < 2e-16 *** x2 -0.65482 0.20820 -3.145 0.00172 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.037 on 797 degrees of freedom (200 observations deleted due to missingness) Multiple R-squared: 0.1223, Adjusted R-squared: 0.12 F-statistic: 55.5 on 2 and 797 DF, p-value: < 2.2e-16 |
Au lieu de faire une régression linéaire, on peut utiliser une autre méthode d’imputation. Par exemple prendre la moyenne des k valeurs de x_1 (observées) pour les k individus ayant des x_2 les plus proches du x_2 de l’individu ayant x_1 manquant :
kpp=function(i,basena,k=5){ x2=basena$x2[i] sb=basena[!is.na(basena$x1),] idx=rank(abs(sb$x2-x2)) mean(sb[which(idx<=k),"x1"]) } |
Sur notre base simulée on obtient
base$x1[indice]=NA base0=base for(j in indice) base0$x1[j]=kpp(j,base0,k=5) reg4=lm(y~x1+x2,data=base) coefficients(reg4) (Intercept) x1 x2 1.197944 1.804220 -0.806766 |
Si on regarde ce que ça donne sur nos 10,000 simulations, on a (c’est un peu long, car j’ai codé ça très rapidement, et pas du tout de manière optimale)
for(s in 1:m){ x=rmnorm(n,varcov = S) x1=pnorm(x[,1]) x2=pnorm(x[,2]) e=rnorm(n,.2) y=1+2*x1-x2+e alpha=.2 indice=sample(1:n,size=round(n*alpha),prob = x1^3) base=data.frame(y=y,x1=x1,x2=x2) base$x1[indice]=NA base0=base for(j in indice) base0$x1[j]=kpp(j,base0,k=5) reg=lm(y~x1+x2,data=base0) B[s]=coefficients(reg)[2] } hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white") lines(density(B),lwd=2,col="blue") abline(v=2,lty=2,col="red") |
Le biais semble ici plus faible que sans imputation… autrement dit, les méthodes d’imputation me semblent plus robustes que cette stratégie visant à remplacer des NA par une valeur arbitraire, et de rajouter une indicatrice dans la régression.