在当我们缺少值时，系统会告诉我用-1代替，然后添加一个指示符，该变量等于-1。这样就可以不删除变量或观测值。

我们在这里模拟数据，然后根据模型生成数据。未定义将转换为NA。一般建议是将缺失值替换为-1，然后拟合未定义的模型。默认情况下，R的策略是删除缺失值。如果未定义50％，则缺少数据，将删除一半的行

  1.  n=1000
  2.  x1=runif(n)
  3.  x2=runif(n)
  4.  e=rnorm(n,.2)
  5.  y=1+2*x1-x2+e
  6.  alpha=.05
  7.  indice=sample(1:n,size=round(n*alpha))
  8.  base=data.frame(y=y,x1=x1)
  9.  base$x1[indice]=NA
  10.  reg=lm(y~x1+x2,data=base)

我们模拟10,000，然后看看未定义的分布，

  1.  m=10000
  2.  B=rep(NA,m)
  3.   
  4.   
  5.   
  6.   
  7.   
  8.   
  9.   
  10.  hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white",xlab="missing values = 50%")
  11.  lines(density(B),lwd=2,col="blue")
  12.  abline(v=2,lty=2,col="red")

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_R语言

 

当然，丢失值的比率较低-丢失的观测值较少，因此估计量的方差较小。

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_R语言_02

 

现在让我们尝试以下策略：用固定的数值替换缺失的值，并添加一个指标，

  1.  B=rep(NA,m)
  2.   
  3.   
  4.   
  5.   
  6.   
  7.   
  8.  hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
  9.  lines(density(B),lwd=2,col="blue")
  10.  abline(v=2,lty=2,col="red")

 

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_R语言_03

 

不会有太大变化，遗漏值的比率下降到5％，

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_R语言_04

 

例如仍有5％的缺失值，我们有

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_缺失值_05

 

如果我们查看样本，尤其是未定义的点，则会观察到

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_ci_06

 

缺失值是完全独立地随机选择的，

  1.  x1=runif(n)
  2.   
  3.   
  4.   
  5.   
  6.   
  7.   
  8.  plot(x1,y,col=clr)

 

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_未定义_07

 

（此处缺失值的1/3为红色）。但可以假设缺失值的最大值，例如，

  1.  x1=runif(n)
  2.   
  3.   
  4.   
  5.   
  6.   
  7.   
  8.  clr=rep("black",n)
  9.  clr[indice]="red"
  10.  plot(x1,y,col=clr)

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_ci_08

 

有人可能想知道，估计量会给出什么？

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_ci_09

 

它变化不大，但是如果仔细观察，我们会有更多差异。如果未定义变量会发生什么，

  1.   
  2.   
  3.  for(s in 1:m){
  4.   
  5.   
  6.   
  7.   
  8.   
  9.   
  10.   
  11.   
  12.   
  13.   
  14.  base$x1[indice]=-1
  15.  reg=lm(y~x1+x2+I(x1==(-1)),data=base)
  16.  B[s]=coefficients(reg)[2]
  17.  }
  18.   

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_R语言_10

 

这次，我们有一个有偏差的估计量。

  1.  set.seed(1)
  2.   
  3.   
  4.   
  5.   
  6.   
  7.   
  8.  indice=sample(1:n,size=round(n*alpha),prob = x1^3)
  9.   
  10.  base$x1[indice]=-1
  11.   
  12.   
  13.   
  14.   
  15.  coefficients(reg1)
  16.  (Intercept) x1 x2 I(x1 == (-1))TRUE
  17.  1.0988005 1.7454385 -0.5149477 3.1000668
  18.  base$x1[indice]=NA
  19.   
  20.   
  21.  coefficients(reg2)
  22.  (Intercept) x1 x2
  23.  1.1123953 1.8612882 -0.6548206

正如我所说的，一种更好的方法是推算。这个想法是为未定义的缺失预测值预测。最简单的方法是创建一个线性模型，并根据非缺失值进行校准。然后在此新基础上估算模型。

  1.  for(s in 1:m){
  2.   
  3.   
  4.   
  5.   
  6.   
  7.   
  8.  base$x1[indice]=NA
  9.  reg0=lm(x1~x2,data=base[-indice,])
  10.  base$x1[indice]=predict(reg0,newdata=base[indice,])
  11.  reg=lm(y~x1+x2,data=base)
  12.   
  13.   
  14.   
  15.   
  16.  }
  17.  hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
  18.  lines(density(B),lwd=2,col="blue")
  19.  abline(v=2,lty=2,col="red")

 

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_R语言_11

在数字示例中，我们得到

  1.  base$x1[indice]=NA
  2.   
  3.   
  4.   
  5.   
  6.   
  7.  coefficients(reg3)
  8.  (Intercept) x1 x2
  9.  1.1593298 1.8612882 -0.6320339

这种方法至少能够纠正偏差

然后，如果仔细观察，我们获得与第一种方法完全相同的值，该方法包括删除缺少值的行。

  1.   
  2.   
  3.  Coefficients:
  4.  Estimate Std. Error t value Pr(>|t|)
  5.  (Intercept) 1.15933 0.06649 17.435 < 2e-16 ***
  6.  x1 1.86129 0.21967 8.473 < 2e-16 ***
  7.  x2 -0.63203 0.20148 -3.137 0.00176 **
  8.  ---
  9.  Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  10.   
  11.  Residual standard error: 1.051 on 997 degrees of freedom
  12.  Multiple R-squared: 0.1094, Adjusted R-squared: 0.1076
  13.  F-statistic: 61.23 on 2 and 997 DF, p-value: < 2.2e-16
  14.   
  15.   
  16.   
  17.  Coefficients: Estimate Std. Error t value Pr(>|t|)
  18.  (Intercept) 1.11240 0.06878 16.173 < 2e-16 ***
  19.  x1 1.86129 0.21666 8.591 < 2e-16 ***
  20.  x2 -0.65482 0.20820 -3.145 0.00172 **
  21.  ---
  22.  Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  23.   
  24.  Residual standard error: 1.037 on 797 degrees of freedom
  25.  (200 observations deleted due to missingness)
  26.  Multiple R-squared: 0.1223, Adjusted R-squared: 0.12
  27.  F-statistic: 55.5 on 2 and 797 DF, p-value: < 2.2e-16

除了进行线性回归外，还可以使用另一种插补方法。

在模拟的基础上，我们获得

  1.   
  2.   
  3.  for(j in indice) base0$x1[j]=kpp(j,base0,k=5)
  4.  reg4=lm(y~x1+x2,data=base)
  5.  coefficients(reg4)
  6.  (Intercept) x1 x2
  7.  1.197944 1.804220 -0.806766

如果我们看一下10,000个模拟中的样子，就会发现

  1.  for(s in 1:m){
  2.   
  3.   
  4.   
  5.   
  6.   
  7.   
  8.   
  9.  base0=base
  10.  for(j in indice) base0$x1[j]=kpp(j,base0,k=5)
  11.  reg=lm(y~x1+x2,data=base0)
  12.  B[s]=coefficients(reg)[2]
  13.  }
  14.  hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
  15.  lines(density(B),lwd=2,col="blue")
  16.  abline(v=2,lty=2,col="red")

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_缺失值_12

 

这里的偏差似乎比没有插补时要弱一些，换句话说，在我看来，插补方法似乎比旨在用任意值替换NA并在回归中添加指标的策略更强大。

 

拓端tecdat|R语言编程指导缺失值的处理：线性回归模型插补_缺失值_13