ダメ出し：データに NA が含まれるときには注意

Regression Analysis on R　kosugitti

標準化偏回帰係数を求めるためにもとのデータを変数ごとに標準化したものを使って重回帰分析をするという手順を示している。

reg2 <- lm(kokugo ~ sansuu + eigo, data = sample)
summary(reg2)

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  22.0649     6.7036    3.29   0.0014 **
sansuu       -0.0789     0.0872   -0.91   0.3675    
eigo          0.8036     0.0340   23.64   <2e-16 ***

sample.z <- scale(sample[4:10])
sample.z <- data.frame(sample.z)
reg3 <- lm(kokugo ~ sansuu + eigo, data = sample.z)
summary(reg3)

しかし，出てきた結果がおかしい。

             Estimate Std. Error t value Pr(>|t|)    
 (Intercept)  0.00409    0.03882    0.11     0.92    
 sansuu      -0.03618    0.03996   -0.91     0.37    
 eigo         0.92324    0.03905   23.64   <2e-16 ***

どこがおかしいか分からない人もいるとは思うが，変数を標準化して重回帰分析をすると，(Intercept) は 0 にならなければならないのだ。

なぜこうなったかという理由は，元のデータ sample に欠損値があるためである。
第 1 には，scale 関数は，変数ごとに NA を自動的に除いて，平均値=0，標準偏差=1 になるように標準化する。

> (d <- data.frame(x=c(1, 3, 4, NA, 7), y=c(3, NA, 7, 3, 2)))
   x  y
1  1  3
2  3 NA
3  4  7
4 NA  3
5  7  2
> (d2 <- data.frame(scale(d)))
     x          y
1 -1.1 -0.3382407
2 -0.3         NA
3  0.1  1.4657098
4   NA -0.3382407
5  1.3 -0.7892283
> apply(d2, 2, mean, na.rm=TRUE)
            x             y
-6.938894e-18  0.000000e+00
> apply(d2, 2, sd, na.rm=TRUE)
x y
1 1

第 2 には，lm 関数は分析に関与する変数（独立変数と従属変数）のどれかが NA のものは分析に使われない。
従って，lm(y ~ x, data=d2) には，2 行目と 4 行目のデータは使われない。
2 行目と 4 行目のデータを除いて d2 の平均値と標準偏差を見るとどうなるか。

> apply(d2[-c(2, 4),], 2, mean, na.rm=TRUE)
        x         y
0.1000000 0.1127469
> apply(d2[-c(2, 4),], 2, sd, na.rm=TRUE)
       x        y
1.200000 1.193201

分析に使われるデータは「標準化されたものではない」。

正しくやるためにはどうするか。
まず，実際に分析に使用する変数だけを含むデータフレームを作り，na.omit なりを使って欠損値を含まないデータフレームにする。その後，scale で標準化する。そのデータフレームを lm で分析する。

> (d3 <- data.frame(scale(na.omit(d2))))
              x          y
1 -1.000000e+00 -0.3779645
3  1.156482e-17  1.1338934
5  1.000000e+00 -0.7559289
> apply(d3, 2, mean, na.rm=TRUE)
            x             y
-7.016595e-17 -1.848565e-17
> apply(d3, 2, sd, na.rm=TRUE)
x y
1 1
> summary(lm(y ~ x, data=d3))

Call:
lm(formula = y ~ x, data = d3)

Residuals:
      1       3       5
-0.5669  1.1339 -0.5669

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.205e-17  8.018e-01   0.000    1.000
x           -1.890e-01  9.820e-01  -0.192    0.879

Residual standard error: 1.389 on 1 degrees of freedom
Multiple R-squared: 0.03571,    Adjusted R-squared: -0.9286
F-statistic: 0.03704 on 1 and 1 DF,  p-value: 0.879

さらには，普通は標準化偏回帰係数はこんな風にしては求めない（2 回も重回帰分析をする必要はない）。