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) |
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") |
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") |
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) |
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) |
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") |
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 |
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") |
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 |
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 |
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"])
} |
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 |
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") |
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.
