Tag Archives: NA

Traitement des valeurs manquantes : remplacer les NA par une constante ?

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.