『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』
第10章 重回帰分析（Ⅱ）
SUGINO Isamu
2019-11-15 13:09:57, Build with R3.5.2
0 全体の章構成
1 重回帰モデルの複雑化
1-1 二乗項(squared term)
\[y=\alpha_{0}+\alpha_{1}(x_{1}-\gamma)^2=(\alpha_{0}+\alpha_{1}\gamma^2)-2\alpha_{1}\gamma\;x_{1}+\alpha_{1}x_{1}^2\] \((\alpha_{0}+\alpha_{1}\gamma^2)=\beta_{0}\)，\(-2\alpha_{1}\gamma=\beta_{1}\)，\(x_{1}^2=x_{2}\)，\(\alpha_{1}=\beta_{2}\)と書き換えれば， \[y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}\] と云う普通の重回帰の式になる（但し独立変数間の相関は極めて強い）。
統計分析ソフトはこの後者の様な重回帰モデルとして係数を計算するだけである。
計算結果から放物線の頂点を逆算すると， \[(\gamma,\alpha_{0})=(-\frac{\beta_{1}}{2\alpha_{1}},\beta_{0}-\alpha_{1}\gamma^2)=(-\frac{\beta_{1}}{2\beta_{2}},\beta_{0}-\alpha_{1}\frac{\beta_{1}^2}{4\alpha_{1}^2})=(-\frac{\beta_{1}}{2\beta_{2}},\beta_{0}-\frac{\beta_{1}^2}{4\beta_{2}})\] 本文の例を，サポートウェブの表記や流儀に合わせて実行してみよう。
data01 <- read.csv("practice.csv")
sum100101 <- summary(lm(income ~ age + edu + I(edu^2), data=data01))
sum100101
Call:
lm(formula = income ~ age + edu + I(edu^2), data = data01)
Residuals:
Min 1Q Median 3Q Max
-546.11 -227.55 -59.18 171.87 1958.94
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1466.440 631.237 2.323 0.020773 *
age 5.053 2.243 2.252 0.024941 *
edu -258.732 92.445 -2.799 0.005428 **
I(edu^2) 11.411 3.392 3.365 0.000856 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 320.4 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.1173, Adjusted R-squared: 0.1094
F-statistic: 14.8 on 3 and 334 DF, p-value: 4.537e-09教育年数eduの1次の項も二次の項も1%水準で有意なので，教育年数は個人年収に対して非線形の（曲線的な）関係を有していると考えられる。その下に凸な放物線の頂点のx座標を求めると，
coef100101 <- sum100101$coefficients[1:4, 1]
coef100101(Intercept) age edu I(edu^2)
1466.440211 5.052575 -258.732377 11.411289 # 放物線の頂点のx座標を求める
(-1)*coef100101[3]/(2*coef100101[4]) edu
11.33668 .edu <- c(9:18)
.y <- coef100101[1] + coef100101[2]*mean(data01$age, na.rm=T) + coef100101[3]*.edu + coef100101[4]*(.edu^2)
RorB <- c("red", "grey", "red", "red", "grey", "red", "grey", "red", "grey", "red")
plot(.edu, .y, xlim=c(range(data01$edu, na.rm=T)), pch=15, col=RorB,
bty="n", xlab="教育年数", ylab="収入", ylim=c(200, 800),
main="教育年数と年収の予測値の関係（平均年齢で計算）")
単純に直線的に上昇しないのが不思議に思われるだろうが，データの中に専修学校高等課程（教育年数11年として計算）の回答者がわずかに含まれ，この人たちの平均年収が中卒や高卒以上の人達と比べてかなり低いためにこうした結果になっている。ケース数や教育年数カテゴリ別平均年収を確認してみよう。
table(data01$edu)
9 11 12 14 16 18
15 3 110 123 113 12 (.observed <- by(data01$income, data01$edu, mean, na.rm=T))data01$edu: 9
[1] 273.3333
--------------------------------------------------------
data01$edu: 11
[1] 208.3333
--------------------------------------------------------
data01$edu: 12
[1] 287.3786
--------------------------------------------------------
data01$edu: 14
[1] 236.4865
--------------------------------------------------------
data01$edu: 16
[1] 537.5
--------------------------------------------------------
data01$edu: 18
[1] 530因みに，それぞれの教育年数グループで平均年齢が異なるのでグループごとの平均年齢を用いて予測値を求めると，317万円，208万円，244万円，313万円，476万円，690万円となる。
# 教育年数のそれぞれのカテゴリごとの平均年齢を用いて年収の予測値を計算する
v_edu <- c(9, 11 ,12, 14, 16, 18)
m_age <- by(data01$age, data01$edu, mean, na.rm=T)
.predict <- coef100101[1] + coef100101[2]*m_age + coef100101[3]*v_edu + coef100101[4]*v_edu^2
.predictdata01$edu: 9
[1] 317.4867
--------------------------------------------------------
data01$edu: 11
[1] 208.3056
--------------------------------------------------------
data01$edu: 12
[1] 243.5885
--------------------------------------------------------
data01$edu: 14
[1] 312.9305
--------------------------------------------------------
data01$edu: 16
[1] 475.557
--------------------------------------------------------
data01$edu: 18
[1] 690.092この予測値ととグループごとの実測値平均を重ねると以下の通りである。
plot(v_edu, .predict, xlim=c(range(data01$edu, na.rm=T)), pch=15, col="#ff000080",
bty="n", xlab="教育年数", ylab="収入", ylim=c(200, 800),
main=paste("年収の予測値（赤，群別平均年齢で計算）と実測値（青）\n（赤の曲線は全体平均年齢で作図）"))
par(new=T)
plot(.x0 <- c(90:180)/10, coef100101[1] + coef100101[2]*mean(data01$age, na.rm=T) + coef100101[3]*.x0 + coef100101[4]*.x0^2,
xlim=c(range(data01$edu, na.rm=T)), type="l", col="#ff000060",
bty="n", xlab="", ylab="", ylim=c(200, 800), axes=F, main="")
par(new=T)
plot(v_edu, .observed, xlim=c(range(data01$edu, na.rm=T)), pch=16, col="#0000ff80",
bty="n", xlab="", ylab="", ylim=c(200, 800), axes=F, main="")
1-2 交互作用項(interaction term)
\[\begin{alignat*}{3} y_{i} &=\beta_{0}+\beta_{1}\color{red}{d_{i}}&&+\beta_{2}x_{2i}+\beta_{3}\color{red}{d_{i}}x_{2i}&&+\beta_{4}x_{3i}+\epsilon_{i},\;\;\;\;(\color{red}{d_{i}}=0\;or\;1)\\ \color{red}{d_{i}}=0\Rightarrow y_{i} &=\beta_{0}&&+\beta_{2}x_{2i}&&+\beta_{4}x_{3i}+\epsilon_{i}\\ \color{red}{d_{i}}=1\Rightarrow y_{i} &=\beta_{0}+\beta_{1}&&+\beta_{2}x_{2i}+\beta_{3}x_{2i}&&+\beta_{4}x_{3i}+\epsilon_{i}\\ &=\beta_{0}+\beta_{1}&&+(\beta_{2}+\beta_{3})x_{2i}&&+\beta_{4}x_{3i}+\epsilon_{i}\\ \end{alignat*}\]
本文では，男性ダミーの主効果を投入して個人年収を説明する重回帰分析で，教育年数は二乗項のみを入れている。その理由も含めて以下に結果を示す。
教育年数の一次の項と二次の項を同時投入するといずれも有意にならず，分散説明率は.3467（自由度調整済みで.3389）である。教育年数の一次の項のみだと.3465(.3406)，二次の項のみだと.3467(.3408)。よって，僅かではあるが二乗項のみの投入が最善である。
data01$male <- c(1, 0)[data01$sex]
summary(lm(income ~ male + age + edu + I(edu^2), data=data01))
Call:
lm(formula = income ~ male + age + edu + I(edu^2), data = data01)
Residuals:
Min 1Q Median 3Q Max
-718.02 -169.48 -46.59 112.37 1787.38
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -421.531 571.202 -0.738 0.46105
male 349.208 32.293 10.814 < 2e-16 ***
age 5.393 1.933 2.790 0.00557 **
edu 11.373 83.474 0.136 0.89171
I(edu^2) 1.134 3.073 0.369 0.71246
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 276.1 on 333 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3467, Adjusted R-squared: 0.3389
F-statistic: 44.19 on 4 and 333 DF, p-value: < 2.2e-16summary(lm(income ~ male + age + edu, data=data01))
Call:
lm(formula = income ~ male + age + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -624.486 153.254 -4.075 5.76e-05 ***
male 352.892 30.670 11.506 < 2e-16 ***
age 5.366 1.929 2.782 0.00571 **
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16summary(reg100102a <- lm(income ~ male + age + I(edu^2), data=data01))
Call:
lm(formula = income ~ male + age + I(edu^2), data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.95 -169.43 -47.52 108.84 1787.45
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -345.3163 115.2701 -2.996 0.00294 **
male 347.8911 30.7684 11.307 < 2e-16 ***
age 5.3980 1.9297 2.797 0.00545 **
I(edu^2) 1.5503 0.2847 5.446 1e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3467, Adjusted R-squared: 0.3408
F-statistic: 59.08 on 3 and 334 DF, p-value: < 2.2e-16この結果から幾つかのパタンについて予測値を計算してみよう。
coef100102a <- summary(reg100102a)$coefficients; coef100102a Estimate Std. Error t value Pr(>|t|)
(Intercept) -345.316252 115.2700728 -2.995715 2.942732e-03
male 347.891096 30.7683839 11.306772 2.572902e-25
age 5.398030 1.9297491 2.797271 5.452397e-03
I(edu^2) 1.550349 0.2846746 5.446038 1.000328e-07sum(coef100102a[,1] * c(1, 1, 45, 12^2)) # 男性，45歳，高校卒[1] 468.7364sum(coef100102a[,1] * c(1, 0, 40, 16^2)) # 女性，40歳，四大卒[1] 267.4943次に，性別と年齢の交互作用を投入した重回帰分析の結果である。
summary(reg100102b <- lm(income ~ male * age + I(edu^2), data=data01))
Call:
lm(formula = income ~ male * age + I(edu^2), data = data01)
Residuals:
Min 1Q Median 3Q Max
-781.69 -164.71 -50.10 98.12 1728.13
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -207.4699 133.6383 -1.552 0.121
male -2.4228 176.7552 -0.014 0.989
age 2.1679 2.5034 0.866 0.387
I(edu^2) 1.6095 0.2849 5.649 3.46e-08 ***
male:age 7.6551 3.8040 2.012 0.045 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.4 on 333 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3545, Adjusted R-squared: 0.3468
F-statistic: 45.73 on 4 and 333 DF, p-value: < 2.2e-16# male * age は male + age + male:age でも同じこれについて上と同じパタンで予測値を求めてみる。ただし，定数項，男性ダミー主効果，年齢主効果は10%水準でも有意ではない（＝0でないとは言えない）ので，計算には含めない。
coef100102b <- summary(reg100102b)$coefficients; coef100102b[,1](Intercept) male age I(edu^2) male:age
-207.469909 -2.422840 2.167872 1.609500 7.655062 sum(coef100102b[4:5,1] * c(12^2, 1*45)) # 男性，45歳，高校卒[1] 576.2457sum(coef100102b[4:5,1] * c(16^2, 0*40)) # 女性，40歳，四大卒[1] 412.0319全く有意にならなかった年齢主効果，男性ダミー主効果を除去したモデルも確認してみよう。
summary(lm(income ~ male + male : age + I(edu^2), data=data01))
Call:
lm(formula = income ~ male + male:age + I(edu^2), data = data01)
Residuals:
Min 1Q Median 3Q Max
-779.65 -166.99 -50.42 92.16 1730.12
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -102.7950 56.9769 -1.804 0.072108 .
male -99.5486 136.5629 -0.729 0.466538
I(edu^2) 1.5848 0.2834 5.593 4.65e-08 ***
male:age 9.7673 2.9180 3.347 0.000909 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.3 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3531, Adjusted R-squared: 0.3473
F-statistic: 60.77 on 3 and 334 DF, p-value: < 2.2e-16summary(lm(income ~ male : age + I(edu^2), data=data01))
Call:
lm(formula = income ~ male:age + I(edu^2), data = data01)
Residuals:
Min 1Q Median 3Q Max
-751.51 -177.80 -50.67 94.53 1756.19
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -95.3141 56.0058 -1.702 0.0897 .
I(edu^2) 1.5346 0.2747 5.587 4.78e-08 ***
male:age 7.6944 0.6536 11.772 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.1 on 335 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3521, Adjusted R-squared: 0.3482
F-statistic: 91.01 on 2 and 335 DF, p-value: < 2.2e-16summary(reg100102c <- lm(income ~ age + male : age + I(edu^2), data=data01))
Call:
lm(formula = income ~ age + male:age + I(edu^2), data = data01)
Residuals:
Min 1Q Median 3Q Max
-781.30 -164.85 -50.10 98.15 1728.49
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -208.4125 114.4184 -1.821 0.0694 .
age 2.1896 1.9319 1.133 0.2579
I(edu^2) 1.6090 0.2823 5.700 2.64e-08 ***
age:male 7.6037 0.6582 11.552 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3545, Adjusted R-squared: 0.3487
F-statistic: 61.15 on 3 and 334 DF, p-value: < 2.2e-16自由度調整済み分散説明率(Adjusted R-squared)に着目すると，本文でも採用している最後のモデルが僅かに良い値である。しかし本文脚注で述べた様に，モデル選択において用いられる基準であるAIC（赤池情報量規準）やBIC（ベイズ情報量規準）を計算すると，年齢と男性ダミーの双方の主効果を除去したモデルが最も良好となる（結果は省略）。
世帯年収と年齢，性別，教育年数の分析
本文の，世帯収入を説明する重回帰分析を行う。教育年数の二乗項を入れないのは，予備的に分析した結果，二乗項を投入するよりも一次の項を投入した方が僅かに良好だった為である。
summary(reg100102d <- lm(fincome ~ age + male * edu, data=data01))
Call:
lm(formula = fincome ~ age + male * edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-913.25 -281.18 -57.49 204.88 1920.85
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1272.063 336.570 -3.779 0.000191 ***
age 9.573 3.312 2.891 0.004137 **
male 790.292 375.726 2.103 0.036297 *
edu 114.381 21.157 5.406 1.35e-07 ***
male:edu -61.381 26.874 -2.284 0.023096 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 438.5 on 289 degrees of freedom
(90 observations deleted due to missingness)
Multiple R-squared: 0.1314, Adjusted R-squared: 0.1194
F-statistic: 10.93 on 4 and 289 DF, p-value: 2.893e-08全ての係数が有意であるが，本文でも書いた通り，ダミー変数，それも交互作用を含む場合には，ダミー変数の値によって予測式を書き分けた方がきちんと理解出来る。
summary(reg100102d <- lm(fincome ~ age + male * edu, data=data01))
Call:
lm(formula = fincome ~ age + male * edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-913.25 -281.18 -57.49 204.88 1920.85
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1272.063 336.570 -3.779 0.000191 ***
age 9.573 3.312 2.891 0.004137 **
male 790.292 375.726 2.103 0.036297 *
edu 114.381 21.157 5.406 1.35e-07 ***
male:edu -61.381 26.874 -2.284 0.023096 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 438.5 on 289 degrees of freedom
(90 observations deleted due to missingness)
Multiple R-squared: 0.1314, Adjusted R-squared: 0.1194
F-statistic: 10.93 on 4 and 289 DF, p-value: 2.893e-08coef100102d <- summary(reg100102d)$coefficients
round(coef100102d[, c(1, 4)], 3) Estimate Pr(>|t|)
(Intercept) -1272.063 0.000
age 9.573 0.004
male 790.292 0.036
edu 114.381 0.000
male:edu -61.381 0.023# 幾つかのパタンについて予測値を計算する
sum(coef100102d[, 1] * c( 1, 45, 1, 12, 1*12))[1] 584.9999sum(coef100102d[, 1] * c( 1, 45, 1, 16, 1*16))[1] 796.9982sum(coef100102d[, 1] * c( 1, 40, 0, 16, 0*16))[1] 940.9438sum(coef100102d[, 1] * c( 1, 40, 0, 12, 0*12))[1] 483.4202変数のセンタリング
ここでは，結果を読み取り易くする為に，0を取り得ない値を，「平均からの偏差」に変換して分析してみよう。
現在の結果では定数項も大幅にマイナスとなっており，意味のある値を表現していない。
また，交互作用項の解釈が平易になる様に，男性ダミーではなく女性ダミーとして投入しよう。
\[\begin{alignat*}{3} &\color{red}{no-centering}&&&&\\ y_{i} &=\beta_{0}+\beta_{1}age_{i}&&+\beta_{2}female_{i}+\beta_{3}edu_{i}&&+\beta_{4}female_{i}\times edu_{i}&&+\epsilon_{i}\\ &\color{green}{centering}&&&&\\ y_{i} &=\beta_{0}'+\beta_{1}(age_{i}-\overline{age})&&+\beta_{2}female_{i}+\beta_{3}(edu_{i}-\overline{edu})&&+\beta_{4}female_{i}\times (edu_{i}-\overline{edu})&&+\epsilon_{i}\\ \end{alignat*}\]
data01$d_age <- data01$age - mean(data01$age, na.rm=T)
data01$d_edu <- data01$edu - mean(data01$edu, na.rm=T)
data01$female <- c(0, 1)[data01$sex]
mean(data01$age, na.rm=T); mean(data01$edu, na.rm=T)[1] 45.80208[1] 13.92021summary(reg100102e <- lm(fincome ~ d_age + female * d_edu, data=data01))
Call:
lm(formula = fincome ~ d_age + female * d_edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-913.25 -281.18 -57.49 204.88 1920.85
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 694.449 39.523 17.571 < 2e-16 ***
d_age 9.573 3.312 2.891 0.00414 **
female 64.149 52.156 1.230 0.21971
d_edu 53.000 16.767 3.161 0.00174 **
female:d_edu 61.381 26.874 2.284 0.02310 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 438.5 on 289 degrees of freedom
(90 observations deleted due to missingness)
Multiple R-squared: 0.1314, Adjusted R-squared: 0.1194
F-statistic: 10.93 on 4 and 289 DF, p-value: 2.893e-08coef100102e <- summary(reg100102e)$coefficients
round(coef100102e[, c(1, 4)], 3) Estimate Pr(>|t|)
(Intercept) 694.449 0.000
d_age 9.573 0.004
female 64.149 0.220
d_edu 53.000 0.002
female:d_edu 61.381 0.023confint(reg100102e) 2.5 % 97.5 %
(Intercept) 616.659842 772.23723
d_age 3.054708 16.09089
female -38.503536 166.80250
d_edu 19.998359 86.00080
female:d_edu 8.487227 114.27542定数項の値は，年齢，教育年数共に平均値に等しい男性の世帯年収平均を意味する。性別の主効果は消えてしまった。
本人の教育年数によって世帯年収平均は上昇するが，その上昇幅は女性の方が2倍余りもある（主効果と交互作用効果の比較）。
1-3 多重共線性(multi-collinearity)
\[y=\beta_{0}+\beta_{1}x_{1}+...+\beta_{j}x_{j}+...+\beta_{p}x_{p}+\epsilon\] この時， \[x_{j}=\alpha_{0}+\alpha_{1}x_{1}+...+\alpha_{j-1}x_{j-1}+\alpha_{j+1}x_{j+1}+...+\alpha_{p}x_{p}+\delta\] の分散説明率を\(R_{x_{j}}^2\)とする。この時，トレランスとVIF(Variance Inflation Factor)は以下の様に定義される。 \[\begin{align*} tolerance_{x_{j}}&=1-R_{x_{j}}^2\\ VIF_{x_{j}}&=\frac{1}{tolerance_{x_{j}}}=\frac{1}{1-R_{x_{j}}^2}\\ \end{align*}\]
トレランスは小さい程，VIFは大きい程問題とされる。目安は明確に定められないが，経験的には，\(VIF>2\)となったら少し気に留める方が良いかも知れない。
第9章で説明しているが，（標本）偏回帰係数\(\hat{\beta}_{x_{j}}=b_{j}\)の標準誤差の推定値は，下記の式で表される。 \[\hat{\sigma}_{b_{j}}=\frac{\hat{\sigma}_{\epsilon}}{\sqrt{SS_{x_j}}\sqrt{1-R_{x_{j}}^2}}=\frac{\sqrt{MS_{residual}}}{\sqrt{SS_{x_j}}\sqrt{1-R_{x_{j}}^2}}=\frac{\sqrt{MS_{residual}}}{\sqrt{SS_{x_j}}}\sqrt{VIF_{x_{j}}}\] \(H_{0}:\;\beta_{x_{j}}=0\)の検定統計量tは， \[t=\frac{b_{j}}{\hat{\sigma}_{b_{j}}}\sim\;t_{(n-p-1)}\] であるが，分母の標準誤差推定値が大きければ絶対値は小さくなりがちであり，偏回帰係数の推定値自体が大きくても統計的に有意にならない事が有り得る。
「点推定値は大きいのに，統計的に有意にならない！」と云うのがこの問題なのである。言いかえれば，点推定値が大きくても，信頼区間が非常に広くなって0を挟んでしまうのである。
本文と同様の重回帰分析を行って，toleranceやvifを計算してみよう。
library(car)Warning: package 'car' was built under R version 3.5.3Loading required package: carDatasummary(reg100103a <- lm(income ~ age + male + edu, data=data01))
Call:
lm(formula = income ~ age + male + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -624.486 153.254 -4.075 5.76e-05 ***
age 5.366 1.929 2.782 0.00571 **
male 352.892 30.670 11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16vif(reg100103a) age male edu
1.051380 1.012674 1.060861 1/vif(reg100103a) # これがトレランスの値 age male edu
0.9511306 0.9874847 0.9426309 summary(reg100103b <- lm(income ~ age + edu + I(edu^2), data=data01))
Call:
lm(formula = income ~ age + edu + I(edu^2), data = data01)
Residuals:
Min 1Q Median 3Q Max
-546.11 -227.55 -59.18 171.87 1958.94
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1466.440 631.237 2.323 0.020773 *
age 5.053 2.243 2.252 0.024941 *
edu -258.732 92.445 -2.799 0.005428 **
I(edu^2) 11.411 3.392 3.365 0.000856 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 320.4 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.1173, Adjusted R-squared: 0.1094
F-statistic: 14.8 on 3 and 334 DF, p-value: 4.537e-09summary(reg100103c<- lm(income ~ age * male + edu, data=data01))
Call:
lm(formula = income ~ age * male + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-781.90 -164.61 -50.88 98.10 1727.93
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -497.0488 164.9765 -3.013 0.00279 **
age 2.1105 2.5026 0.843 0.39965
male -0.1679 176.7227 -0.001 0.99924
edu 43.6989 7.7425 5.644 3.56e-08 ***
age:male 7.7188 3.8055 2.028 0.04332 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.4 on 333 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3544, Adjusted R-squared: 0.3467
F-statistic: 45.71 on 4 and 333 DF, p-value: < 2.2e-16vif(reg100103b) age edu I(edu^2)
1.052632 112.211490 112.378552 vif(reg100103c) age male edu age:male
1.786086 33.935555 1.072961 34.059939 モデルreg100103bについては，標準誤差が著しく大きくなってt検定が有意にならないという問題は生じていない。モデルreg100103cでは特にmaleの主効果で標準誤差が著しく大きいので，これを独立変数から除外しても良いかも知れない。
summary(reg100103c <- lm(income ~ age * male + edu, data=data01))
Call:
lm(formula = income ~ age * male + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-781.90 -164.61 -50.88 98.10 1727.93
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -497.0488 164.9765 -3.013 0.00279 **
age 2.1105 2.5026 0.843 0.39965
male -0.1679 176.7227 -0.001 0.99924
edu 43.6989 7.7425 5.644 3.56e-08 ***
age:male 7.7188 3.8055 2.028 0.04332 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.4 on 333 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3544, Adjusted R-squared: 0.3467
F-statistic: 45.71 on 4 and 333 DF, p-value: < 2.2e-16summary(reg100103d <- lm(income ~ age + edu + age:male , data=data01))
Call:
lm(formula = income ~ age + edu + age:male, data = data01)
Residuals:
Min 1Q Median 3Q Max
-781.87 -164.62 -50.88 98.11 1727.96
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -497.1084 152.3763 -3.262 0.00122 **
age 2.1120 1.9291 1.095 0.27437
edu 43.6980 7.6735 5.695 2.71e-08 ***
age:male 7.7153 0.6564 11.754 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3544, Adjusted R-squared: 0.3486
F-statistic: 61.13 on 3 and 334 DF, p-value: < 2.2e-16vif(reg100103c); vif(reg100103d) age male edu age:male
1.786086 33.935555 1.072961 34.059939 age edu age:male
1.064434 1.057094 1.016386 confint(reg100103c); confint(reg100103d) 2.5 % 97.5 %
(Intercept) -821.5763840 -172.521287
age -2.8123529 7.033366
male -347.8015019 347.465758
edu 28.4685541 58.929205
age:male 0.2328914 15.204782 2.5 % 97.5 %
(Intercept) -796.846584 -197.37017
age -1.682615 5.90665
edu 28.603501 58.79247
age:male 6.424062 9.00649AIC(reg100103c, reg100103d)maleの主効果を除外したモデルreg100103dを計算すると，先のモデルreg100103cよりも，定数項を含むそれぞれの要因の標準誤差が小さくなって信頼区間が狭くなり，自由度調整済み分散説明率は上昇し，赤池情報量基準AICは低下するというように，すべての指標においてreg100103dのほうが優れている。
定義に従ったVIFの計算
本文で触れた，carパッケイジを使わずにvifを計算する方法を紹介しておく。本項冒頭のモデルreg100103aについて自力でVIFを求めるには下記の様にする。ただしこの方法は，特に独立変数の数が多い場合には手間がかかる。
Full <- complete.cases(data01$income, data01$age, data01$male, data01$edu)
y <- data01$income[Full]
x1 <- data01$age[Full]
x2 <- data01$male[Full]
x3 <- data01$edu[Full]
summary(lm(x1 ~ x2 + x3))
Call:
lm(formula = x1 ~ x2 + x3)
Residuals:
Min 1Q Median 3Q Max
-15.5571 -6.2986 0.3625 6.4429 15.3625
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 57.9054 2.9719 19.484 < 2e-16 ***
x2 -0.4701 0.8683 -0.541 0.589
x3 -0.8624 0.2139 -4.031 6.88e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.809 on 335 degrees of freedom
Multiple R-squared: 0.04887, Adjusted R-squared: 0.04319
F-statistic: 8.606 on 2 and 335 DF, p-value: 0.0002266summary(lm(x2 ~ x1 + x3))
Call:
lm(formula = x2 ~ x1 + x3)
Residuals:
Min 1Q Median 3Q Max
-0.5429 -0.4170 -0.3568 0.5512 0.7339
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.151117 0.272882 0.554 0.5801
x1 -0.001860 0.003435 -0.541 0.5886
x3 0.024969 0.013711 1.821 0.0695 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4912 on 335 degrees of freedom
Multiple R-squared: 0.01252, Adjusted R-squared: 0.00662
F-statistic: 2.123 on 2 and 335 DF, p-value: 0.1213summary(lm(x3 ~ x1 + x2))
Call:
lm(formula = x3 ~ x1 + x2)
Residuals:
Min 1Q Median 3Q Max
-5.5845 -1.4946 0.1396 1.6300 4.3274
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.12295 0.62928 25.621 < 2e-16 ***
x1 -0.05364 0.01331 -4.031 6.88e-05 ***
x2 0.39261 0.21559 1.821 0.0695 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.948 on 335 degrees of freedom
Multiple R-squared: 0.05737, Adjusted R-squared: 0.05174
F-statistic: 10.19 on 2 and 335 DF, p-value: 5.038e-05# vif
vif(reg100103a) # carパッケイジのvif( )関数を使用した結果 age male edu
1.051380 1.012674 1.060861 1/(1 - summary(lm(x1 ~ x2 + x3))$r.squared) # age[1] 1.051381/(1 - summary(lm(x2 ~ x1 + x3))$r.squared) # male[1] 1.0126741/(1 - summary(lm(x3 ~ x1 + x2))$r.squared) # edu[1] 1.060861# tolerance
1/vif(reg100103a) # carパッケイジのvif( )関数を使用した結果 age male edu
0.9511306 0.9874847 0.9426309 1 - summary(lm(x1 ~ x2 + x3))$r.squared # age[1] 0.95113061 - summary(lm(x2 ~ x1 + x3))$r.squared # male[1] 0.98748471 - summary(lm(x3 ~ x1 + x2))$r.squared # edu[1] 0.94263091-4 変数のコントロールとは
本文で例示している予測値を求める。既に何度か紹介した方法である。
summary(reg100103a)
Call:
lm(formula = income ~ age + male + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -624.486 153.254 -4.075 5.76e-05 ***
age 5.366 1.929 2.782 0.00571 **
male 352.892 30.670 11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16coef100103a <- summary(reg100103a)$coefficients; coef100103a Estimate Std. Error t value Pr(>|t|)
(Intercept) -624.486186 153.253683 -4.074853 5.757668e-05
age 5.366085 1.928999 2.781798 5.712928e-03
male 352.892055 30.670080 11.506069 4.949757e-26
edu 42.031195 7.734500 5.434248 1.063014e-07sum(coef100103a[, 1]*c(1, 40, 0, 14))[1] 178.5939sum(coef100103a[, 1]*c(1, 40, 0, 16))[1] 262.6563sum(coef100103a[, 1]*c(1, 40, 1, 16))[1] 615.5484各カテゴリの実際の標本平均は次の様にして求める事が出来る。
# cases1，2，3は，各ケースが指定の条件に該当するか否かをTRUE(1)かFALSE(0)で表現するヴェクトル
cases1 <- data01$sex==2 & data01$age >=35 & data01$age <=44 & data01$edu==14
cases2 <- data01$sex==2 & data01$age >=35 & data01$age <=44 & data01$edu==16
cases3 <- data01$sex==1 & data01$age >=35 & data01$age <=44 & data01$edu==16
summary(data01$income[cases1]) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 25.0 125.0 149.2 200.0 800.0 2 summary(data01$income[cases2]) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 25.0 200.0 432.7 600.0 1750.0 3 summary(data01$income[cases3]) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 500.0 600.0 627.6 862.5 1125.0 5 発展1 一般線型モデル(GLM, LM)
かつては一般線型モデル\((General\;Linear\; Model)\)の事を，GLMと呼ぶ事が一般的であったが，現在では一般化線型モデル\((General\color{red}{ized}\;Linear\;Model)\)の事をGLMと略す事が多く，一般線型モデルの方は単にLMと呼ぶ事が多くなって来ている。
性別変数を3通りの方法で重回帰分析に投入したのが以下である。それぞれの結果をよく見比べて，何処が違うかを確認した上で，1番目と3番目についてはそれぞれ予測式を具体的に考えて，本質的に同じである事を確認しよう。
summary(lm(income ~ age + sex + edu, data=data01))
Call:
lm(formula = income ~ age + sex + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 81.298 163.342 0.498 0.61901
age 5.366 1.929 2.782 0.00571 **
sex -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16summary(lm(income ~ age + factor(sex) + edu, data=data01))
Call:
lm(formula = income ~ age + factor(sex) + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -271.594 155.380 -1.748 0.08139 .
age 5.366 1.929 2.782 0.00571 **
factor(sex)2 -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16data01$female <- c(0, 1)[data01$sex]
summary(lm(income ~ age + female + edu, data=data01))
Call:
lm(formula = income ~ age + female + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -271.594 155.380 -1.748 0.08139 .
age 5.366 1.929 2.782 0.00571 **
female -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16本文で言及している残りのモデルについても結果を掲示しておくので読み解いてみよう。（表示結果を簡略化する為にsummary( )ではなくlm( )の結果だけを示す。）
lm(income ~ age * sex + edu, data=data01)
Call:
lm(formula = income ~ age * sex + edu, data = data01)
Coefficients:
(Intercept) age sex edu age:sex
-497.3846 17.5482 0.1679 43.6989 -7.7188 lm(income ~ age * factor(sex) + edu, data=data01)
Call:
lm(formula = income ~ age * factor(sex) + edu, data = data01)
Coefficients:
(Intercept) age factor(sex)2 edu
-497.2167 9.8293 0.1679 43.6989
age:factor(sex)2
-7.7188 lm(income ~ age + edu + age:sex , data=data01)
Call:
lm(formula = income ~ age + edu + age:sex, data = data01)
Coefficients:
(Intercept) age edu age:sex
-497.108 17.543 43.698 -7.715 lm(income ~ age + edu + age:factor(sex) , data=data01)
Call:
lm(formula = income ~ age + edu + age:factor(sex), data = data01)
Coefficients:
(Intercept) age edu age:factor(sex)2
-497.108 9.827 43.698 -7.715 発展2 ‘良い’モデルの検討
発展2-1 変数選択
まずは本文と同様の変数減少法の結果を紹介した後，それ以外の変数選択法についても簡単に紹介する。
出力そのままではstepごとの区切りが分かりにくいかも知れないので，各ステップをピンクの枠線で囲んである。が何もしない場合であり，どれかの独立変数を除去したらAICが低下するかどうかを確かめながら，それ以上除去する独立変数がなくなった時点で終了する。最後に，採用されたモデルの偏回帰係数が出力される。
step(lm(fincome ~ age * edu * factor(sex), data=data01))Start: AIC=3585.14
fincome ~ age * edu * factor(sex)
Df Sum of Sq RSS AIC
- age:edu:factor(sex) 1 147674 55184080 3583.9
<none> 55036407 3585.1
Step: AIC=3583.93
fincome ~ age + edu + factor(sex) + age:edu + age:factor(sex) +
edu:factor(sex)
Df Sum of Sq RSS AIC
- age:factor(sex) 1 23618 55207698 3582.1
<none> 55184080 3583.9
- age:edu 1 393611 55577691 3584.0
- edu:factor(sex) 1 1086688 56270768 3587.7
Step: AIC=3582.05
fincome ~ age + edu + factor(sex) + age:edu + edu:factor(sex)
Df Sum of Sq RSS AIC
- age:edu 1 371070 55578768 3582.0
<none> 55207698 3582.1
- edu:factor(sex) 1 1063716 56271414 3585.7
Step: AIC=3582.02
fincome ~ age + edu + factor(sex) + edu:factor(sex)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
- edu:factor(sex) 1 1003253 56582021 3585.3
- age 1 1606902 57185670 3588.4
Call:
lm(formula = fincome ~ age + edu + factor(sex) + edu:factor(sex),
data = data01)
Coefficients:
(Intercept) age edu factor(sex)2
-481.771 9.573 53.000 -790.292
edu:factor(sex)2
61.381 これをst1というオブジェクトに格納してsummary( )を見ると以下の通りである。また，各ステップの要点だけを表形式に表現するには，オブジェクトの中のanova情報を表示させる。オブジェクトst1に格納する際に上記と同じ変数選択の過程が出力されるが以下では省略する。
st1 <- step(lm(fincome ~ age * edu * factor(sex), data=data01))Start: AIC=3585.14
fincome ~ age * edu * factor(sex)
Df Sum of Sq RSS AIC
- age:edu:factor(sex) 1 147674 55184080 3583.9
<none> 55036407 3585.1
Step: AIC=3583.93
fincome ~ age + edu + factor(sex) + age:edu + age:factor(sex) +
edu:factor(sex)
Df Sum of Sq RSS AIC
- age:factor(sex) 1 23618 55207698 3582.1
<none> 55184080 3583.9
- age:edu 1 393611 55577691 3584.0
- edu:factor(sex) 1 1086688 56270768 3587.7
Step: AIC=3582.05
fincome ~ age + edu + factor(sex) + age:edu + edu:factor(sex)
Df Sum of Sq RSS AIC
- age:edu 1 371070 55578768 3582.0
<none> 55207698 3582.1
- edu:factor(sex) 1 1063716 56271414 3585.7
Step: AIC=3582.02
fincome ~ age + edu + factor(sex) + edu:factor(sex)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
- edu:factor(sex) 1 1003253 56582021 3585.3
- age 1 1606902 57185670 3588.4summary(st1) # 最終的に選択されたモデルの結果
Call:
lm(formula = fincome ~ age + edu + factor(sex) + edu:factor(sex),
data = data01)
Residuals:
Min 1Q Median 3Q Max
-913.25 -281.18 -57.49 204.88 1920.85
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -481.771 304.899 -1.580 0.11518
age 9.573 3.312 2.891 0.00414 **
edu 53.000 16.767 3.161 0.00174 **
factor(sex)2 -790.292 375.726 -2.103 0.03630 *
edu:factor(sex)2 61.381 26.874 2.284 0.02310 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 438.5 on 289 degrees of freedom
(90 observations deleted due to missingness)
Multiple R-squared: 0.1314, Adjusted R-squared: 0.1194
F-statistic: 10.93 on 4 and 289 DF, p-value: 2.893e-08names(st1) # オブジェクトに格納されている情報の名前一覧 [1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "na.action" "contrasts" "xlevels" "call"
[13] "terms" "model" "anova" st1$anova # 変数選択の過程で検討されたモデルの分散分析表st1$anovaの結果は，1が出発点のモデル，そこから2，3，4と，AICが低下する限り変数を減少させていっている。
AIC(st1) # ある（単純な）計算方法におけるAIC[1] 4418.357logLik(st1) # 対数尤度，dfとして最尤推定されるパラメタの個数kが示されている。独立変数の個数4個+2になっている。'log Lik.' -2203.179 (df=6)(-2)*logLik(st1)+2*attributes(logLik(st1))$df # -2LL+2kを計算'log Lik.' 4418.357 (df=6)extractAIC(st1) # 別の方法で計算されたAIC[1] 5.000 3582.021step( )関数で使用されているのはextractAIC( )の計算法と同じである事が分かる。AICはモデル同士のAICの差だけが重要で，それについてはextraxtAIC( )で計算してもAIC( )で計算しても同じであるので気にしなくて良い。
# extractAIC(st1)の結果を再現する
log(deviance(st1)/nrow(st1$model))*nrow(st1$model)+2*st1$rank[1] 3582.021変数減少法の選択の最初に与えたモデルについても，AIC( )とextractAIC( )の結果の再現を示しておく。
# step( )の最初のfullモデルについてもextractAIC( )の値を再現しよう
model2 <- lm(fincome ~ age * edu * factor(sex), data=data01)
log(deviance(model2)/nrow(model2$model))*nrow(model2$model)+2*model2$rank[1] 3585.138extractAIC(model2)[1] 8.000 3585.138AIC(model2) # 単純な方のAIC[1] 4421.474(-2)*logLik(model2)+2*attributes(logLik(model2))$df'log Lik.' 4421.474 (df=9)いろいろな変数選択法
さて，漸く，よりカスタマイズした変数選択法の紹介である。
# 変数増加法において，モデルが複雑になるにつれて有効ケースが減少してしまうので，
# 最初から完備ケース分析の準備をしておく必要がある。
full <- complete.cases(data01$fincome, data01$age, data01$edu, data01$sex)
fincome_c <- data01$fincome[full]
age_c <- data01$age[full]
edu_c <- data01$edu[full]
sex_c <- data01$sex[full]
# 完備ケース分析用の変数はデータフレイムに格納していない点に注意
# 年齢＋教育年数＋性別の主効果だけからなるモデルを出発点(object)として，
# それらの二次の交互作用，三次の交互作用まで含むモデルを最大とする(scope)場合の
# 変数増加法(forward)。つまり，交互作用のどれを入れると良いかを検討している
step(object=lm(fincome_c ~ age_c + edu_c + factor(sex_c)),
scope=list(upper=~ age_c * edu_c * factor(sex_c)),
direction="forward")Start: AIC=3585.28
fincome_c ~ age_c + edu_c + factor(sex_c)
Df Sum of Sq RSS AIC
+ edu_c:factor(sex_c) 1 1003253 55578768 3582.0
<none> 56582021 3585.3
+ age_c:edu_c 1 310607 56271414 3585.7
+ age_c:factor(sex_c) 1 16097 56565924 3587.2
Step: AIC=3582.02
fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
+ age_c:edu_c 1 371070 55207698 3582.1
+ age_c:factor(sex_c) 1 1077 55577691 3584.0
Call:
lm(formula = fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c))
Coefficients:
(Intercept) age_c edu_c
-481.771 9.573 53.000
factor(sex_c)2 edu_c:factor(sex_c)2
-790.292 61.381 # 変数増減法
step(lm(fincome_c ~ age_c * edu_c * factor(sex_c)), direction="both")Start: AIC=3585.14
fincome_c ~ age_c * edu_c * factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:edu_c:factor(sex_c) 1 147674 55184080 3583.9
<none> 55036407 3585.1
Step: AIC=3583.93
fincome_c ~ age_c + edu_c + factor(sex_c) + age_c:edu_c + age_c:factor(sex_c) +
edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:factor(sex_c) 1 23618 55207698 3582.1
<none> 55184080 3583.9
- age_c:edu_c 1 393611 55577691 3584.0
+ age_c:edu_c:factor(sex_c) 1 147674 55036407 3585.1
- edu_c:factor(sex_c) 1 1086688 56270768 3587.7
Step: AIC=3582.05
fincome_c ~ age_c + edu_c + factor(sex_c) + age_c:edu_c + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:edu_c 1 371070 55578768 3582.0
<none> 55207698 3582.1
+ age_c:factor(sex_c) 1 23618 55184080 3583.9
- edu_c:factor(sex_c) 1 1063716 56271414 3585.7
Step: AIC=3582.02
fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
+ age_c:edu_c 1 371070 55207698 3582.1
+ age_c:factor(sex_c) 1 1077 55577691 3584.0
- edu_c:factor(sex_c) 1 1003253 56582021 3585.3
- age_c 1 1606902 57185670 3588.4
Call:
lm(formula = fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c))
Coefficients:
(Intercept) age_c edu_c
-481.771 9.573 53.000
factor(sex_c)2 edu_c:factor(sex_c)2
-790.292 61.381 # 変数減少法
step(lm(fincome_c ~ age_c * edu_c * factor(sex_c)), direction="backward")Start: AIC=3585.14
fincome_c ~ age_c * edu_c * factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:edu_c:factor(sex_c) 1 147674 55184080 3583.9
<none> 55036407 3585.1
Step: AIC=3583.93
fincome_c ~ age_c + edu_c + factor(sex_c) + age_c:edu_c + age_c:factor(sex_c) +
edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:factor(sex_c) 1 23618 55207698 3582.1
<none> 55184080 3583.9
- age_c:edu_c 1 393611 55577691 3584.0
- edu_c:factor(sex_c) 1 1086688 56270768 3587.7
Step: AIC=3582.05
fincome_c ~ age_c + edu_c + factor(sex_c) + age_c:edu_c + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
- age_c:edu_c 1 371070 55578768 3582.0
<none> 55207698 3582.1
- edu_c:factor(sex_c) 1 1063716 56271414 3585.7
Step: AIC=3582.02
fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c)
Df Sum of Sq RSS AIC
<none> 55578768 3582.0
- edu_c:factor(sex_c) 1 1003253 56582021 3585.3
- age_c 1 1606902 57185670 3588.4
Call:
lm(formula = fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c:factor(sex_c))
Coefficients:
(Intercept) age_c edu_c
-481.771 9.573 53.000
factor(sex_c)2 edu_c:factor(sex_c)2
-790.292 61.381 かなり長くなってしまったので，自分でよく見比べて異同を確認して欲しい。変数増加法では，出発点のモデルを変えると，最終結果がbothやbackwardとは異なってしまう事もある。object=lm(fincome_c ~ 1)などとして確認すると良い。
step(object=lm(fincome_c ~ 1),
scope=list(upper=~ age_c * edu_c * factor(sex_c)),
direction="forward")Start: AIC=3615.43
fincome_c ~ 1
Df Sum of Sq RSS AIC
+ edu_c 1 5333194 58652205 3591.8
+ age_c 1 867851 63117549 3613.4
<none> 63985400 3615.4
+ factor(sex_c) 1 76310 63909090 3617.1
Step: AIC=3591.85
fincome_c ~ edu_c
Df Sum of Sq RSS AIC
+ age_c 1 1818871 56833335 3584.6
<none> 58652205 3591.8
+ factor(sex_c) 1 283437 58368768 3592.4
Step: AIC=3584.58
fincome_c ~ edu_c + age_c
Df Sum of Sq RSS AIC
<none> 56833335 3584.6
+ age_c:edu_c 1 377115 56456219 3584.6
+ factor(sex_c) 1 251313 56582021 3585.3
Call:
lm(formula = fincome_c ~ edu_c + age_c)
Coefficients:
(Intercept) edu_c age_c
-785.90 75.05 10.16 発展2-2 誤差減少率
前項2-2のst1と同等のモデル（完備ケース分析の処理済みのもの）について，分散説明率と誤差減少率が一致する事を（数式から自明ではあるが）一応確認しておく。
# st1と同じモデル。但し完備ケース分析の処理をしたもの。
reg100203 <- lm(fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c*factor(sex_c))
summary(reg100203) # Multiple R-squaredを確認
Call:
lm(formula = fincome_c ~ age_c + edu_c + factor(sex_c) + edu_c *
factor(sex_c))
Residuals:
Min 1Q Median 3Q Max
-913.25 -281.18 -57.49 204.88 1920.85
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -481.771 304.899 -1.580 0.11518
age_c 9.573 3.312 2.891 0.00414 **
edu_c 53.000 16.767 3.161 0.00174 **
factor(sex_c)2 -790.292 375.726 -2.103 0.03630 *
edu_c:factor(sex_c)2 61.381 26.874 2.284 0.02310 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 438.5 on 289 degrees of freedom
Multiple R-squared: 0.1314, Adjusted R-squared: 0.1194
F-statistic: 10.93 on 4 and 289 DF, p-value: 2.893e-08SST <- var(fincome_c)*(length(fincome_c)-1); SST # 従属変数の偏差平方和[1] 63985400SSr <- var(reg100203$residuals)*(length(reg100203$residuals)-1); SSr # 残差平方和[1] 55578768# 誤差減少率の定義式に代入
1-SSr/SST # 平方和から計算[1] 0.13138361-var(reg100203$residuals)/var(fincome_c) # （不偏）分散から計算[1] 0.1313836第10章の【練習問題】の解答
1) の答え
以下の通りで，vifは非常に大きいが，ageとage^2はいずれも5%有意である。
data01$male <- c(1, 0)[data01$sex]
summary(reg10a1 <- lm(edu ~ male + age + I(age^2), data=data01))
Call:
lm(formula = edu ~ male + age + I(age^2), data = data01)
Residuals:
Min 1Q Median 3Q Max
-5.6913 -1.4822 0.2537 1.7052 4.2831
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.571921 3.322319 6.794 4.34e-11 ***
male 0.318687 0.202483 1.574 0.1164
age -0.343527 0.148608 -2.312 0.0213 *
I(age^2) 0.003216 0.001628 1.975 0.0490 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.935 on 372 degrees of freedom
(8 observations deleted due to missingness)
Multiple R-squared: 0.06148, Adjusted R-squared: 0.05392
F-statistic: 8.124 on 3 and 372 DF, p-value: 2.974e-05car::vif(reg10a1) # carパッケイジから関数vif( )だけを呼び出すやり方 male age I(age^2)
1.008407 140.787202 140.688687 age^2の係数が正なので，年齢と教育年数の関係は下に凸な放物線となる。その頂点のx座標（年齢）は，以下の計算から約53.4歳となる。データに含まれるのは30歳～59歳であるが，53歳以下の対象者では，若いほど教育年数が長い傾向があるが，50代半ば以降では（何故か）年長者の方が教育年数が長いという関係になる（50代後半では高学歴層の回答率が高かった可能性も考えられる）。
coef_a1 <- summary(reg10a1)$coefficients
coef_a1 Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.571920597 3.322318546 6.794027 4.340435e-11
male 0.318686634 0.202482838 1.573895 1.163619e-01
age -0.343527205 0.148608268 -2.311629 2.134385e-02
I(age^2) 0.003215804 0.001628121 1.975163 4.898858e-02coef_a1[3, 1]/((-2)*coef_a1[4, 1])[1] 53.412332) の答え
まず問題文の前半の結果は以下の通りである。vifの結果はどれもかなり悪いが，検定は5%水準～10%水準で有意になっている。
summary(reg10a2 <- lm(edu ~ male + age + I(age^2) + male:age, data=data01))
Call:
lm(formula = edu ~ male + age + I(age^2) + male:age, data = data01)
Residuals:
Min 1Q Median 3Q Max
-5.9341 -1.5716 0.3187 1.6643 4.4892
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.079177 3.395222 6.208 1.43e-09 ***
male 2.587744 1.169656 2.212 0.0275 *
age -0.300924 0.149609 -2.011 0.0450 *
I(age^2) 0.002997 0.001626 1.844 0.0660 .
male:age -0.049555 0.025162 -1.969 0.0496 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.927 on 371 degrees of freedom
(8 observations deleted due to missingness)
Multiple R-squared: 0.0712, Adjusted R-squared: 0.06118
F-statistic: 7.11 on 4 and 371 DF, p-value: 1.594e-05car::vif(reg10a2) male age I(age^2) male:age
33.90969 143.79339 141.34738 33.95809 問題の後半である。表示の仕方に少し工夫をしたものも付加してある。自分で調べて解読すると良い。
names(reg10a2) [1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "na.action" "xlevels" "call" "terms"
[13] "model" coef10a2 <- reg10a2$coefficients # 定数項と偏回帰係数の行列を取り出す
p.age <- 50 # 予測式で求める年齢を，調査対象の30～59歳から自由に指定。
p.male <- 1 # 男性なら1，女性なら0とする
profile <- c(1, p.male, p.age, p.age^2, p.male*p.age) # 独立変数行列
p.male; p.age[1] 1[1] 50sum(coef10a2*profile)[1] 13.63609# 一目でわかる出力を工夫してみる
sprintf("%s，%d歳の教育年数の予測値%.1f年", c("女性", "男性")[p.male+1], p.age, sum(coef10a2*profile))[1] "男性，50歳の教育年数の予測値13.6年"3) の答え
まずは分析結果を確認し，予測式を男女別に数式表現する。
summary(lm(income ~ age + sex + edu, data=data01))
Call:
lm(formula = income ~ age + sex + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 81.298 163.342 0.498 0.61901
age 5.366 1.929 2.782 0.00571 **
sex -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16 予測式は，\(income = 81.3 + 5.4*age - 353*sex + 42*edu\) となる。
ここで，age年齢は調査時満年齢，eduは教育年数であるが，sexは男性=1，女性=2という値である事に注意する。男性と女性でそれぞれ予測式を書き分けると以下の様になる。
男性の予測式：\(income = 81.3 + 5.4*age - 353*1 + 42*edu = -272 + 5.4*age + 42*edu\)
女性の予測式：\(income = 81.3 + 5.4*age - 353*2 + 42*edu = -625 + 5.4*age + 42*edu\)
summary(lm(income ~ age + factor(sex) + edu, data=data01))
Call:
lm(formula = income ~ age + factor(sex) + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -271.594 155.380 -1.748 0.08139 .
age 5.366 1.929 2.782 0.00571 **
factor(sex)2 -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16factor(sex)2というのは，factor(sex)が2の値を取る場合に，という意味であり，factor(sex)が1の値（男性）の場合は関係しない。
予測式は，\(income = -272 + 5.4*age - 353*factor(sex)2 + 42*edu\) となる。
factor(sex)2の項の意味に注意して予測式を書き分ける。
男性の予測式：\(income = -272 + 5.4*age + 42*edu\)
女性の予測式：\(income = -272 + 5.4*age - 353 + 42*edu = -625 + 5.4*age + 42*edu\)
となり，上の結果と一致する。factor(sex)を投入するのは，女性ダミー変数を投入するのと全く同じ事になる。
data01$female <- c(0, 1)[data01$sex]
summary(lm(income ~ age + female + edu, data=data01))
Call:
lm(formula = income ~ age + female + edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -271.594 155.380 -1.748 0.08139 .
age 5.366 1.929 2.782 0.00571 **
female -352.892 30.670 -11.506 < 2e-16 ***
edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16定数項に意味を持たせるには
因みに，年齢や教育年数を平均からの偏差に変換すると，定数項は平均的なプロファイルの人の年収平均を意味するので分かりやすい。
data01$d_age <- data01$age - mean(data01$age, na.rm=T)
data01$d_edu <- data01$edu - mean(data01$edu, na.rm=T)
mean(data01$age, na.rm=T); mean(data01$edu, na.rm=T)[1] 45.80208[1] 13.92021summary(lm(income ~ d_age + female + d_edu, data=data01))
Call:
lm(formula = income ~ d_age + female + d_edu, data = data01)
Residuals:
Min 1Q Median 3Q Max
-717.50 -169.66 -46.16 113.65 1787.86
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 559.267 23.430 23.870 < 2e-16 ***
d_age 5.366 1.929 2.782 0.00571 **
female -352.892 30.670 -11.506 < 2e-16 ***
d_edu 42.031 7.735 5.434 1.06e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 275.7 on 334 degrees of freedom
(46 observations deleted due to missingness)
Multiple R-squared: 0.3465, Adjusted R-squared: 0.3406
F-statistic: 59.02 on 3 and 334 DF, p-value: < 2.2e-16