拓端tecdat|R語言缺失值的處理：線性回歸模型插補

原文鏈接： http://tecdat.cn/?p=14528

 

 

在當我們缺少值時，系統會告訴我用-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")

 

當然，丟失值的比率較低-丟失的觀測值較少，因此估計量的方差較小。

 

現在讓我們嘗試以下策略：用固定的數值替換缺失的值，並添加一個指標，

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")

 

 

不會有太大變化，遺漏值的比率下降到5％，

 

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

 

如果我們查看樣本，尤其是未定義的點，則會觀察到

 

缺失值是完全獨立地隨機選擇的，

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

 

 

（此處缺失值的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)

 

有人可能想知道，估計量會給出什么？

 

它變化不大，但是如果仔細觀察，我們會有更多差異。如果未定義變量會發生什么，

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.  
     

 

這次，我們有一個有偏差的估計量。

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")

 

在數字示例中，我們得到

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")

 

這里的偏差似乎比沒有插補時要弱一些，換句話說，在我看來，插補方法似乎比旨在用任意值替換NA並在回歸中添加指標的策略更強大。

 

參考文獻

1.用SPSS估計HLM層次線性模型模型

2.R語言線性判別分析（LDA），二次判別分析（QDA）和正則判別分析（RDA）

3.基於R語言的lmer混合線性回歸模型

4.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析

5.在r語言中使用GAM（廣義相加模型）進行電力負荷時間序列分析

6.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM

7.R語言中的嶺回歸、套索回歸、主成分回歸：線性模型選擇和正則化

8.R語言用線性回歸模型預測空氣質量臭氧數據

9.R語言分層線性模型案例