{"word_count": 960, "character_count": 9637}
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 [latex]x_1[/latex] et [latex]x_2[/latex], on génère ensuite des données suivant un modèle, de la forme [latex]y=\beta_0+\beta_1x_1+\beta_2x_2+\varepsilon[/latex]. Une proportion [latex]\alpha[/latex] de [latex]x_1[/latex] 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 [latex display="true"]y=\beta_0+\beta_1x_1+\beta_{-1}\mathbf{1}(x_1=-1)+\beta_2x_2+\varepsilon[/latex]. 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 [latex]x_1[/latex] sont manquantes, la moitié des lignes sont suppriméesn=1000x1=runif(n)x2=runif(n)e=rnorm(n,.2)y=1+2*x1-x2+ealpha=.05indice=sample(1:n,size=round(n*alpha))base=data.frame(y=y,x1=x1)base$x1[indice]=NAreg=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 [latex]\widehat{\beta}_1[/latex],m=10000B=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 [latex]\alpha=5\%[/latex] - 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 [latex]x_1[/latex]  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 aSi on regarde notre échantillon, en particulier le nuage de points  [latex](x_1,y)[/latex], on observeici, 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+ealpha=.3333333indice=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 [latex]x_1[/latex], par exemple,x1=runif(n)x2=runif(n)e=rnorm(n,.2)y=1+2*x1-x2+ealpha=.3333333indice=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 [latex]\widehat{\beta}_1[/latex]Ç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 [latex]x_1[/latex] et [latex]x_2[/latex] sont maintenant corrélées,B=rep(NA,m)library(mnormt)r=.8S = 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 procheset.seed(1)x=rmnorm(n,varcov = S)x1=pnorm(x[,1])x2=pnorm(x[,2])e=rnorm(n,.2)y=1+2*x1-x2+ealpha=.2indice=sample(1:n,size=round(n*alpha),prob = x1^3)base=data.frame(y=y,x1=x1,x2=x2)base$x1[indice]=-1reg1=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]=NAreg2=lm(y~x1+x2,data=base)coefficients(reg2)(Intercept)          x1          x2   1.1123953   1.8612882  -0.6548206Comme 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 [latex]x_1[/latex] manquants. On pourrait être tenté de prendre  [latex]\overline{x}_1[/latex], la moyenne des [latex]x_1[/latex] observés. Mais on sait que les valeurs manquantes sont justement pour les grandes valeurs de [latex]x_1[/latex], ici, donc on doit pouvoir faire mieux ! On sait aussi que [latex]x_1[/latex]  et [latex]x_2[/latex] sont corrélés ici. Corrélés positivement, en plus. Autrement dit, si [latex]x_2[/latex] est grand, on sait que le [latex]x_1[/latex] (non observé) devait être grand. Le plus simple est de faire un modèle linéaire, [latex]x_1=\alpha_0+\alpha_2x_2+\eta_i[/latex], calibré sur les valeurs non-manquantes. Puis on utilise [latex]\widehat{x}_1=\widehat{\alpha}_0+\widehat{\alpha}_2x_2[/latex] 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 obtientbase$x1[indice]=NAreg0=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.6320339Cette 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(&#62;|t|)    (Intercept)  1.15933    0.06649  17.435  &#60; 2e-16 ***x1           1.86129    0.21967   8.473  &#60; 2e-16 ***x2          -0.63203    0.20148  -3.137  0.00176 ** ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 1.051 on 997 degrees of freedomMultiple R-squared:  0.1094,Adjusted R-squared:  0.1076 F-statistic: 61.23 on 2 and 997 DF,  p-value: &#60; 2.2e-16 summary(reg2) Coefficients: Estimate Std. Error t value Pr(&#62;|t|)    (Intercept)  1.11240    0.06878  16.173  &#60; 2e-16 ***x1           1.86129    0.21666   8.591  &#60; 2e-16 ***x2          -0.65482    0.20820  -3.145  0.00172 ** ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual 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: &#60; 2.2e-16Au lieu de faire une régression linéaire, on peut utiliser une autre méthode d'imputation. Par exemple prendre la moyenne des [latex]k[/latex] valeurs de [latex]x_1[/latex] (observées) pour les [latex]k[/latex] individus ayant des [latex]x_2[/latex] les plus proches du [latex]x_2[/latex] de l'individu ayant [latex]x_1[/latex] 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&#60;=k),"x1"])}Sur notre base simulée on obtientbase$x1[indice]=NAbase0=basefor(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.806766Si 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.