R語言缺失值的處理：線性回歸模型插補(bǔ)
原文鏈接：?http://tecdat.cn/?p=14528
?
在當(dāng)我們?nèi)鄙僦禃r，系統(tǒng)會告訴我用-1代替，然后添加一個指示符，該變量等于-1。這樣就可以不刪除變量或觀測值。
視頻
缺失值的處理：線性回歸模型插補(bǔ)
我們在這里模擬數(shù)據(jù)，然后根據(jù)模型生成數(shù)據(jù)。未定義將轉(zhuǎn)換為NA。一般建議是將缺失值替換為-1，然后擬合未定義的模型。默認(rèn)情況下，R的策略是刪除缺失值。如果未定義50％，則缺少數(shù)據(jù)，將刪除一半的行
n=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)
?
我們模擬10,000，然后看看未定義的分布，
m=10000B=rep(NA,m)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")
?
?
當(dāng)然，丟失值的比率較低-丟失的觀測值較少，因此估計(jì)量的方差較小。
?
現(xiàn)在讓我們嘗試以下策略：用固定的數(shù)值替換缺失的值，并添加一個指標(biāo)，
B=rep(NA,m)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")
?
?
不會有太大變化，遺漏值的比率下降到5％，
?
例如仍有5％的缺失值，我們有
?
如果我們查看樣本，尤其是未定義的點(diǎn)，則會觀察到
?
缺失值是完全獨(dú)立地隨機(jī)選擇的，
x1=runif(n)plot(x1,y,col=clr)
?
?
（此處缺失值的1/3為紅色）。但可以假設(shè)缺失值的最大值，例如，
x1=runif(n)clr=rep("black",n)clr[indice]="red"plot(x1,y,col=clr)
?
?
有人可能想知道，估計(jì)量會給出什么？
?
它變化不大，但是如果仔細(xì)觀察，我們會有更多差異。如果未定義變量會發(fā)生什么，
for(s in 1:m){base$x1[indice]=-1reg=lm(y~x1+x2+I(x1==(-1)),data=base)B[s]=coefficients(reg)[2]}
?
?
這次，我們有一個有偏差的估計(jì)量。
set.seed(1)indice=sample(1:n,size=round(n*alpha),prob = x1^3)base$x1[indice]=-1coefficients(reg1)(Intercept) ? ? ? ? ? ? ? ?x1 ? ? ? ? ? ? ? ?x2 I(x1 == (-1))TRUE1.0988005 ? ? ? ? 1.7454385 ? ? ? ?-0.5149477 ? ? ? ? 3.1000668base$x1[indice]=NAcoefficients(reg2)(Intercept) ? ? ? ? ?x1 ? ? ? ? ?x21.1123953 ? 1.8612882 ?-0.6548206
?
正如我所說的，一種更好的方法是推算。這個想法是為未定義的缺失預(yù)測值預(yù)測。最簡單的方法是創(chuàng)建一個線性模型，并根據(jù)非缺失值進(jìn)行校準(zhǔn)。然后在此新基礎(chǔ)上估算模型。
for(s in 1:m){base$x1[indice]=NAreg0=lm(x1~x2,data=base[-indice,])base$x1[indice]=predict(reg0,newdata=base[indice,])reg=lm(y~x1+x2,data=base)}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")
?
?
在數(shù)字示例中，我們得到
base$x1[indice]=NAcoefficients(reg3)(Intercept) ? ? ? ? ?x1 ? ? ? ? ?x21.1593298 ? 1.8612882 ?-0.6320339
?
這種方法至少能夠糾正偏差
然后，如果仔細(xì)觀察，我們獲得與第一種方法完全相同的值，該方法包括刪除缺少值的行。
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 ‘ ’ 1Residual standard error: 1.051 on 997 degrees of freedomMultiple R-squared: ?0.1094, Adjusted R-squared: ?0.1076F-statistic: 61.23 on 2 and 997 DF, ?p-value: < 2.2e-16Coefficients: 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 ‘ ’ 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.12F-statistic: ?55.5 on 2 and 797 DF, ?p-value: < 2.2e-16
?
除了進(jìn)行線性回歸外，還可以使用另一種插補(bǔ)方法。
在模擬的基礎(chǔ)上，我們獲得
for(j in indice) base0$x1[j]=kpp(j,base0,k=5)reg4=lm(y~x1+x2,data=base)coefficients(reg4)(Intercept) ? ? ? ? ?x1 ? ? ? ? ?x21.197944 ? ?1.804220 ? -0.806766
?
如果我們看一下10,000個模擬中的樣子，就會發(fā)現(xiàn)
for(s in 1:m){base0=basefor(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")
?
?
這里的偏差似乎比沒有插補(bǔ)時要弱一些，換句話說，在我看來，插補(bǔ)方法似乎比旨在用任意值替換NA并在回歸中添加指標(biāo)的策略更強(qiáng)大。
?
參考文獻(xiàn)
1.用SPSS估計(jì)HLM層次線性模型模型
2.R語言線性判別分析（LDA），二次判別分析（QDA）和正則判別分析（RDA）
3.基于R語言的lmer混合線性回歸模型
4.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
5.在r語言中使用GAM（廣義相加模型）進(jìn)行電力負(fù)荷時間序列分析
6.使用SAS，Stata，HLM，R，SPSS和Mplus的分層線性模型HLM
7.R語言中的嶺回歸、套索回歸、主成分回歸：線性模型選擇和正則化
8.R語言用線性回歸模型預(yù)測空氣質(zhì)量臭氧數(shù)據(jù)
9.R語言分層線性模型案例
