様々な変数をとれるということは様々な定義域をとるということ

　重回帰では様々な変数を解析に組み込めることを紹介したが、それは同時に様々な定義域の変数が解析に含まれうることを意味する。説明変数のオーダーが違えども、問題なく推定は行われるが、その時算出された係数の解釈は注意が必要である。例えば、植物の高さcmに影響を与える要因として、与えた水量、土壌窒素含量、日照量を説明変数に組み込むことを考える。与えた水量はL（リットル）、土壌窒素含量はmg、日射量はJ/m^2で測ったとしよう。すると、まあ、実験設計によるが、水量は0 ～ 数十L、窒素含量は0 ～ 数千mg、日射量は0 ～ 数百万J/m^2くらいのオーダーになるだろう。もし、水量の係数が0.5、日射量の係数が0.01と推定されたら、このとき、水量のほうが影響が強い、と判断できるだろうか。実はそんな簡単には判断できない。係数が、各説明変数が1増加した時に増加する被説明変数の量であるので、水量は1L増加すれば0.5 cm増加し、日射量は1J/m^2増加すれば0.01cm増加することを意味する。見方を変えれば、水量は数十L分の1で0.5程度の影響しかなく、日射量は数百万J/m^2の1で0.01も影響がある、と考えることができる。このように、係数をそのまま解釈するだけでは、その要因の影響力を測ることはできない。

標準化（スケーリング）で変数間の影響量を比較する

　そこで、1単位の影響量を比較できるように、変数を標準化（スケーリング）することがある。方法は様々だが、その中でも最も一般的なものは変数を平均0、標準偏差1に統一する方法である。これは変数から変数の平均を引いたのち、もとの変数の標準偏差で割れば簡単に計算できる。R上では関数一つですぐに計算できる。

では、下記で具体的なデータをもとに解析してみよう。

スケーリングの具体例

　下記のようなデータを解析することを考える。Rではscale関数を用いることで、上記で紹介したスケーリングを実行してくれる。データフレームd内には、元の変数とスケーリングした変数、そしてx1とx2を掛け算し、スケーリングした変数を格納している。これは交互作用項を模したものだが、これの必要性は後述する。

------------------------------------------------------

library(plotn)

x1 <- c(29, 12, 22, 20, 14, 24, 17, 19, 24, 14, 

        19, 14, 23, 24, 21, 23, 13, 11, 27, 25, 

        24, 15, 11, 11, 13, 29, 18, 29, 17, 27, 

        25, 13, 24, 19, 23, 16, 24, 22, 24, 29, 

        28, 23, 14, 13, 22, 22, 12, 29, 19, 28, 

        18, 27, 14, 11, 28, 29, 26, 20, 24, 21)

x2 <- c(1923, 2914, 1260, 1654, 2478, 

        1279, 2764, 1605, 2461, 2094, 

        1865, 1094, 2597, 2474, 1032, 

        1105, 1130, 2963, 2090, 2656, 

        2304, 1638, 1974, 1407, 2190, 

        1885, 2479, 2521, 1622, 2318, 

        1259, 1893, 2795, 1773, 2498, 

        2013, 1536, 1109, 1195, 2869, 

        2288, 2407, 2754, 1509, 2080, 

        1239, 1309, 1086, 2344, 2518, 

        2397, 1561, 1812, 1405, 2684, 

        2639, 1227, 1951, 1583, 2649)

y <- c(494, 419, 357, 374, 399, 376, 448, 363, 496, 348, 

       375, 257, 492, 509, 315, 346, 257, 409, 490, 521, 

       485, 323, 311, 261, 340, 486, 436, 549, 331, 499, 

       380, 324, 535, 371, 486, 363, 402, 337, 363, 580, 

       514, 482, 425, 288, 430, 344, 255, 406, 427, 539,

       426, 430, 329, 252, 547, 556, 391, 395, 400, 484)

d <- data.frame(x1 = x1, x2 = x2, y = y, sx1 = scale(x1), sx2 = scale(x2), sx12 = scale(x1*x2), sy = scale(y))

head(d)

##   x1   x2   y        sx1         sx2       sx12         sy

## 1 29 1923 494  1.4535437 -0.07978025  0.8349930  0.9879591

## 2 12 2914 419 -1.4881519  1.62782730 -0.3315744  0.1184469

## 3 22 1260 357  0.2422573 -1.22220589 -0.7380978 -0.6003499

## 4 20 1654 374 -0.1038245 -0.54329834 -0.4374679 -0.4032604

## 5 14 2478 399 -1.1420700  0.87654890 -0.3470546 -0.1134230

## 6 24 1279 376  0.5883391 -1.18946669 -0.5711809 -0.3800735

plotn(d[,1:3]) #図1の描画

plotn(d[,4:7]) #図2の描画

------------------------------------------------------

確認すると、x1（スケーリングしたsx1）とx2（sx2）はともにyと正の相関があることがわかる。偏差の分布は均等なので正規分布を仮定した重回帰を実行すればよいだろう。まず、x1とx2を変数とした解析を行う。交互作用はいったん考慮しないことにする。

------------------------------------------------------

fit1 <- lm(y ~ x1 + x2, data = d)

fit2 <- summary(fit1)

fit2

## 

## Call:

## lm(formula = y ~ x1 + x2, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -14.1616  -3.6122   0.1291   2.9072  15.4082 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept)  3.85506    3.24947   1.186     0.24    

## x1          10.09241    0.11749  85.903   <2e-16 ***

## x2           0.10005    0.00117  85.518   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 5.192 on 57 degrees of freedom

## Multiple R-squared:  0.9965, Adjusted R-squared:  0.9964 

## F-statistic:  8114 on 2 and 57 DF,  p-value: < 2.2e-16

------------------------------------------------------

解析の結果、x1の係数の推定値は10程度、x2の推定値は0.1程度となった。ともに有意であるが、さて先の議論でこの係数をそのまま解釈して、x1のほうが影響が強いと判断できないことを述べた。x1は10～30の値、x2は1000～3000の値をとる。図1を見ると横軸の長さがそろっているため、x1が10～30を動く時と、x2が1000～3000を動く時のyの影響の程度は同程度のように見える。まさにこの図の感覚がスケーリングである。実際、今度はスケーリングした変数を説明変数として解析してみよう。

------------------------------------------------------

fit3 <- lm(sy ~ sx1 + sx2, data = d)

fit4 <- summary(fit3)

fit4

## 

## Call:

## lm(formula = sy ~ sx1 + sx2, data = d)

## 

## Residuals:

##       Min        1Q    Median        3Q       Max 

## -0.164183 -0.041878  0.001496  0.033705  0.178635 

## 

## Coefficients:

##               Estimate Std. Error t value Pr(>|t|)    

## (Intercept) -2.152e-16  7.770e-03    0.00        1    

## sx1          6.762e-01  7.871e-03   85.90   <2e-16 ***

## sx2          6.731e-01  7.871e-03   85.52   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 0.06019 on 57 degrees of freedom

## Multiple R-squared:  0.9965, Adjusted R-squared:  0.9964 

## F-statistic:  8114 on 2 and 57 DF,  p-value: < 2.2e-16

------------------------------------------------------

すると今度は、sx1とsx2の係数は0.67程度と推定され、効果の大きさが同程度であることがわかる。このようにスケーリングをすれば、定義域の異なる変数間でも効果の大きさを比較することができる。係数の推定値はスケーリングによって変わったが、t値はスケーリング前と変わらず、それゆえP値も変わらない。つまり、（後述の状況を除けば）x1やx2が有意かどうかは、x1、x2およびyのスケーリングによって影響を受けない。なお、切片だけは説明変数のスケーリングによって有意かどうか変わるので注意である。今回はさらにyもスケーリングしたので、syは平均がぴったり0になり、有意になることはない。つまり、切片に関してはスケーリング後に解釈する意味はないと言えるだろう。

　では、今度は交互作用を入れた時にどうなるかを確認してみよう。交互作用を説明変数に入れたときは注意が必要になる。まずはスケーリング前で解析する。

------------------------------------------------------

fit5 <- lm(y ~ x1 * x2, data = d)

fit6 <- summary(fit5)

fit6

## 

## Call:

## lm(formula = y ~ x1 * x2, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -14.0597  -3.7751   0.0214   3.0516  15.3303 

## 

## Coefficients:

##              Estimate Std. Error t value Pr(>|t|)    

## (Intercept) 6.506e+00  8.868e+00   0.734    0.466    

## x1          9.963e+00  4.198e-01  23.731   <2e-16 ***

## x2          9.873e-02  4.272e-03  23.110   <2e-16 ***

## x1:x2       6.408e-05  1.992e-04   0.322    0.749    

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 5.233 on 56 degrees of freedom

## Multiple R-squared:  0.9965, Adjusted R-squared:  0.9963 

## F-statistic:  5325 on 3 and 56 DF,  p-value: < 2.2e-16

------------------------------------------------------

上記のように推定できた。実は交互作用を含めずにデータ生成したので、有意な効果がない状態が真のモデルである。さて、スケーリング後の変数で交互作用を考慮して解析してみよう。

------------------------------------------------------

fit7 <- lm(sy ~ sx1 * sx2, data = d)

fit8 <- summary(fit7)

fit8

## 

## Call:

## lm(formula = sy ~ sx1 * sx2, data = d)

## 

## Residuals:

##       Min        1Q    Median        3Q       Max 

## -0.163002 -0.043766  0.000248  0.035379  0.177732 

## 

## Coefficients:

##               Estimate Std. Error t value Pr(>|t|)    

## (Intercept) -0.0002319  0.0078654  -0.029    0.977    

## sx1          0.6759496  0.0079654  84.860   <2e-16 ***

## sx2          0.6731386  0.0079340  84.842   <2e-16 ***

## sx1:sx2      0.0024915  0.0077443   0.322    0.749    

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 0.06067 on 56 degrees of freedom

## Multiple R-squared:  0.9965, Adjusted R-squared:  0.9963 

## F-statistic:  5325 on 3 and 56 DF,  p-value: < 2.2e-16

------------------------------------------------------

すると、交互作用を考慮しなかったときと異なり、sx1とsx2のt値が変わっている。つまり、表示上は分からないがP値も変わっている（交互作用だけは変わらない）。このようにスケーリング時の交互作用の扱いが不適切だと、変数の効果の有意か否かが、スケーリングの前後で変わってしまう。これは、交互作用項をスケーリングしなかったことが問題である。そこで、交互作用項に相当するx1×x2をスケーリングした値をsx12としてデータフレームd内に格納している。これをもとのモデル式の代わりに説明変数に加えることにしよう。

------------------------------------------------------

fit9 <- lm(sy ~ sx1 + sx2 + sx12, data = d)

fit10 <- summary(fit9)

fit10

## 

## Call:

## lm(formula = sy ~ sx1 + sx2 + sx12, data = d)

## 

## Residuals:

##       Min        1Q    Median        3Q       Max 

## -0.163002 -0.043766  0.000248  0.035379  0.177732 

## 

## Coefficients:

##               Estimate Std. Error t value Pr(>|t|)    

## (Intercept) -2.165e-16  7.832e-03   0.000    1.000    

## sx1          6.675e-01  2.813e-02  23.731   <2e-16 ***

## sx2          6.643e-01  2.874e-02  23.110   <2e-16 ***

## sx12         1.324e-02  4.117e-02   0.322    0.749    

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 0.06067 on 56 degrees of freedom

## Multiple R-squared:  0.9965, Adjusted R-squared:  0.9963 

## F-statistic:  5325 on 3 and 56 DF,  p-value: < 2.2e-16

------------------------------------------------------

すると、上記の解析結果はfit6と同じt値のなった。このように交互作用項の考慮時は注意が必要である。上記の現象が偶然ではないことをもう一つのデータセットで確認してみよう。

------------------------------------------------------

x1 <- c(11, 16, 15, 29, 27, 28, 27, 29, 24, 27, 

        16, 29, 11, 19, 24, 25, 14, 10, 21, 10, 

        24, 27, 23, 24, 12, 21, 18, 21, 29, 13, 

        11, 13, 16, 23, 21, 29, 28, 27, 22, 24, 

        21, 25, 17, 16, 19, 14, 23, 12, 10, 14, 

        13, 10, 15, 14, 27, 22, 19, 15, 10, 28)

x2 <- c(23670, 19806, 25488, 10129, 14578,

        22931, 27996, 11691, 17701, 27682, 

        14425, 20454, 11144, 26650, 21536, 

        13302, 29111, 10644, 12781, 12778, 

        22916, 21033, 21666, 17664, 23902, 

        27773, 24032, 16411, 14815, 23414, 

        17985, 23192, 28104, 14679, 18951, 

        10575, 20639, 10239, 28167, 15482,

        18474, 28102, 20795, 16267, 18332,

        20950, 26634, 13944, 18407, 22981, 

        17064, 17753, 12172, 12677, 16901, 

        23180, 16591, 14173, 22202, 13347)

y <- c(299853, 353397, 414552, 305100, 422427, 

       676352, 798798, 355171, 458245, 792914, 

       256978, 623659, 141582, 546435, 555404, 

       355238, 443693, 126937, 285541, 142941, 

       578971, 604219, 525018, 455532, 326791,

       620997, 466808, 368823, 447955, 343273,

       218216, 346470, 499396, 355470, 422177,

       319443, 602848, 293874, 658762, 399613,

       409420, 732941, 384127, 277569, 379446,

       317638, 646993, 183893, 220378, 354906, 

       240408, 203499, 196442, 190683, 479233, 

       552928, 333692, 240758, 259763, 394668)

d <- data.frame(x1 = x1, x2 = x2, y = y, sx1 = scale(x1), sx2 = scale(x2), sx12 = scale(x1*x2), sy = scale(y))

head(d)

##   x1    x2      y        sx1        sx2       sx12         sy

## 1 11 23670 299853 -1.3719333  0.8312999 -0.7216988 -0.6308081

## 2 16 19806 353397 -0.5834659  0.1256395 -0.3662197 -0.3048940

## 3 15 25488 414552 -0.7411594  1.1633110  0.0452169  0.0673471

## 4 29 10129 305100  1.4665494 -1.6416163 -0.5118362 -0.5988704

## 5 27 14578 422427  1.1511624 -0.8291206  0.1161920  0.1152810

## 6 28 22931 676352  1.3088559  0.6963405  1.6787127  1.6608834

plotn(d[,1:3]) #図3の描画

plotn(d[,4:7]) #図4の描画

------------------------------------------------------

確認すると、x1とx2はyと正の相関がありそうだ。実はこのデータは交互作用を含むように設定している。しいて言うなら、x1のデータから交互作用がありそうにも見えるが、まあ、これくらいなら交互作用無しで解析しても問題ないくらいだと思う。

------------------------------------------------------

fit1 <- lm(y ~ x1 + x2, data = d)

fit2 <- summary(fit1)

fit2

## 

## Call:

## lm(formula = y ~ x1 + x2, data = d)

## 

## Residuals:

##    Min     1Q Median     3Q    Max 

## -77161 -25191    617  20896  80856 

## 

## Coefficients:

##               Estimate Std. Error t value Pr(>|t|)    

## (Intercept) -3.641e+05  2.215e+04  -16.44   <2e-16 ***

## x1           1.774e+04  7.164e+02   24.77   <2e-16 ***

## x2           2.187e+01  8.297e-01   26.35   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 34860 on 57 degrees of freedom

## Multiple R-squared:  0.9565, Adjusted R-squared:  0.955 

## F-statistic: 626.6 on 2 and 57 DF,  p-value: < 2.2e-16

------------------------------------------------------

解析の結果、x1の係数は1.8×10^4、x2の係数は22程度とわかる。これらの効果はスケーリングするとどうなるだろうか。

------------------------------------------------------

fit3 <- lm(sy ~ sx1 + sx2, data = d)

fit4 <- summary(fit3)

fit4

## 

## Call:

## lm(formula = sy ~ sx1 + sx2, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -0.46967 -0.15333  0.00376  0.12719  0.49216 

## 

## Coefficients:

##               Estimate Std. Error t value Pr(>|t|)    

## (Intercept) -3.159e-16  2.740e-02    0.00        1    

## sx1          6.849e-01  2.765e-02   24.77   <2e-16 ***

## sx2          7.288e-01  2.765e-02   26.35   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 0.2122 on 57 degrees of freedom

## Multiple R-squared:  0.9565, Adjusted R-squared:  0.955 

## F-statistic: 626.6 on 2 and 57 DF,  p-value: < 2.2e-16

------------------------------------------------------

すると、sx1の係数は0.68程度、sx2の係数は0.73程度と、なんと係数は逆転してしまった。まあ、これくらいなら同程度の強さ、ということになるだろうが、こういうこともある。t値は変わっていないので、有意かどうかには影響を与えていないことがわかるだろう

　では、今度は交互作用を入れた時にどうなるかを確認してみよう。

------------------------------------------------------

fit5 <- lm(y ~ x1 * x2, data = d)

fit6 <- summary(fit5)

fit6

## 

## Call:

## lm(formula = y ~ x1 * x2, data = d)

## 

## Residuals:

##      Min       1Q   Median       3Q      Max 

## -10235.8  -4073.6   -697.5   4135.9   9682.9 

## 

## Coefficients:

##               Estimate Std. Error t value Pr(>|t|)    

## (Intercept) -3.785e+03  8.071e+03  -0.469 0.640935    

## x1           1.875e+02  3.736e+02   0.502 0.617785    

## x2           1.760e+00  4.283e-01   4.109 0.000131 ***

## x1:x2        9.861e-01  2.008e-02  49.120  < 2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 5297 on 56 degrees of freedom

## Multiple R-squared:  0.999,  Adjusted R-squared:  0.999 

## F-statistic: 1.89e+04 on 3 and 56 DF,  p-value: < 2.2e-16

------------------------------------------------------

上記のように推定できた。交互作用を入れるとx1の効果は有意ではなくなってしまった。交互作用を正しくスケーリングせずに解析すると、

------------------------------------------------------

fit7 <- lm(sy ~ sx1 * sx2, data = d)

fit8 <- summary(fit7)

fit8

## 

## Call:

## lm(formula = sy ~ sx1 * sx2, data = d)

## 

## Residuals:

##       Min        1Q    Median        3Q       Max 

## -0.062304 -0.024796 -0.004245  0.025175  0.058938 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)    

## (Intercept) 0.008976   0.004167   2.154   0.0355 *  

## sx1         0.734957   0.004324 169.979   <2e-16 ***

## sx2         0.706160   0.004227 167.054   <2e-16 ***

## sx1:sx2     0.208431   0.004243  49.120   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 0.03224 on 56 degrees of freedom

## Multiple R-squared:  0.999,  Adjusted R-squared:  0.999 

## F-statistic: 1.89e+04 on 3 and 56 DF,  p-value: < 2.2e-16

------------------------------------------------------

このように今度はx1が有意になった。解析結果が変わってしまうことの典型例となっている。正しくスケーリングすると、

------------------------------------------------------

fit9 <- lm(sy ~ sx1 + sx2 + sx12, data = d)

fit10 <- summary(fit9)

fit10

## 

## Call:

## lm(formula = sy ~ sx1 + sx2 + sx12, data = d)

## 

## Residuals:

##       Min        1Q    Median        3Q       Max 

## -0.062304 -0.024796 -0.004245  0.025175  0.058938 

## 

## Coefficients:

##               Estimate Std. Error t value Pr(>|t|)    

## (Intercept) -1.239e-18  4.163e-03   0.000 1.000000    

## sx1          7.237e-03  1.442e-02   0.502 0.617785    

## sx2          5.866e-02  1.428e-02   4.109 0.000131 ***

## sx12         9.545e-01  1.943e-02  49.120  < 2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 0.03224 on 56 degrees of freedom

## Multiple R-squared:  0.999,  Adjusted R-squared:  0.999 

## F-statistic: 1.89e+04 on 3 and 56 DF,  p-value: < 2.2e-16

------------------------------------------------------

すると、上記の解析結果はfit6と同じt値になる。

　このようにスケーリングによって異なるオーダーの説明変数でも係数の大きさから影響の大きさを比較が可能となる。また、今回の解析では影響がないが、ベイズ推定をするときはスケーリングによって推定値が収束しやすくなることが知られている。あまりに変数間のオーダーに違いがあると収束しにくくなるらしい。覚えておいて損はないだろう。