統計解析R演習2
演習の概要
- 重回帰分析
注意事項
- 慶應義塾大学SFCで開講している「統計解析」の授業履修者向けの演習用ページです。
- 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.6.0を使っています。
- Rの更新などにより、Rコードなどを予告無しに変更する場合があります。
データ分析の準備
データのダウンロード
- 演習で用いるデータはここからダウンロードしてください。
作業用フォルダの指定
- ダウンロードしたzipファイルを解凍し、作成されたフォルダを作業用フォルダとして指定してください。
setwd()関数を使うなどして、作業用フォルダを指定できます。例えば、HDの直下にDSBdataフォルダを作成した場合は以下のようになります。
setwd("~/DSBdata")データの読み込み
- 『家計調査』2017年9月、第２−６表「年間収入階級別１世帯当たり１か月間の収入と支出」から作成した所得階級と支出のデータ
kakei201709_l18.csvを用います。 - ここでは、
read.csv()関数を使って、CSVファイルを読み込みます。
kakei <- read.csv("kakei201709_l18.csv", header=TRUE, sep=",")- 重回帰分析で用いる３つの変数
INC1,CONS1,WORKの散布図をプロットします。
plot(kakei[,4:6])
回帰分析
重回帰分析
lm()関数を使って、線形回帰モデル（単回帰モデル）を最小二乗推定します。- まず、
INC1のみを説明変数とする単回帰モデルを推定します。
kakei.lm1 <- lm(CONS1~INC1, data=kakei)
summary(kakei.lm1)
Call:
lm(formula = CONS1 ~ INC1, data = kakei)
Residuals:
Min 1Q Median 3Q Max
-3.5039 -1.0488 0.3011 1.2763 2.3413
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.263033 0.768529 18.56 3.02e-12 ***
INC1 0.216059 0.009384 23.02 1.08e-13 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.774 on 16 degrees of freedom
Multiple R-squared: 0.9707, Adjusted R-squared: 0.9689
F-statistic: 530.1 on 1 and 16 DF, p-value: 1.08e-13- 次に、
INC1とWORKを説明変数とする重回帰モデルを推定します。
kakei.lm2 <- lm(CONS1~INC1+WORK, data=kakei)
summary(kakei.lm2)
Call:
lm(formula = CONS1 ~ INC1 + WORK, data = kakei)
Residuals:
Min 1Q Median 3Q Max
-2.8671 -0.9996 -0.1032 1.3300 1.9492
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.03768 1.11306 10.815 1.77e-08 ***
INC1 0.18704 0.01419 13.181 1.19e-09 ***
WORK 2.98212 1.19404 2.498 0.0246 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.539 on 15 degrees of freedom
Multiple R-squared: 0.9793, Adjusted R-squared: 0.9765
F-statistic: 354.9 on 2 and 15 DF, p-value: 2.337e-13- 残差プロット、QQプロット、S-Lプロット、梃子比とクックの距離（R-Lプロット）を以下のようにして図示します。
plot(kakei.lm2)
- 比較のために、
WORKのみを説明変数とする単回帰モデルを推定します。
kakei.lm3 <- lm(CONS1~WORK, data=kakei)
summary(kakei.lm3)
Call:
lm(formula = CONS1 ~ WORK, data = kakei)
Residuals:
Min 1Q Median 3Q Max
-6.5383 -3.6378 -0.5404 2.5427 14.3834
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.656 3.556 1.871 0.0797 .
WORK 15.870 2.354 6.742 4.73e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.287 on 16 degrees of freedom
Multiple R-squared: 0.7396, Adjusted R-squared: 0.7234
F-statistic: 45.45 on 1 and 16 DF, p-value: 4.733e-06標準化したデータを用いた重回帰分析
scale()関数を用いると、データを標準化できます。
z.kakei <- as.data.frame(scale(kakei))
kakei.lm4 <- lm(z.kakei[,5]~z.kakei[,4]+z.kakei[,6])
summary(kakei.lm4)
Call:
lm(formula = z.kakei[, 5] ~ z.kakei[, 4] + z.kakei[, 6])
Residuals:
Min 1Q Median 3Q Max
-0.28523 -0.09944 -0.01026 0.13231 0.19391
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.102e-16 3.610e-02 0.000 1.0000
z.kakei[, 4] 8.529e-01 6.471e-02 13.181 1.19e-09 ***
z.kakei[, 6] 1.616e-01 6.471e-02 2.498 0.0246 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1531 on 15 degrees of freedom
Multiple R-squared: 0.9793, Adjusted R-squared: 0.9765
F-statistic: 354.9 on 2 and 15 DF, p-value: 2.337e-13AICによるモデル間の比較
AIC()関数を用いて、モデル推定結果を比較できます。
AIC(kakei.lm1, kakei.lm2, kakei.lm3)BICによるモデル間の比較
BIC()関数を用いて、モデル推定結果を比較できます。
BIC(kakei.lm1, kakei.lm2, kakei.lm3)多重共線性
carパッケージのvif()関数を用いて、多重共線性をチェックできます。install.packages("car")でパッケージをインストールしてください。
library(car)
car::vif(kakei.lm2) INC1 WORK
3.035161 3.035161 - VIFが10以下のとき多重共線性がないと判断されることから、この２変数の間に多重共線性があるとはいえない。
以上で今回の演習は終了です。
LS0tCnRpdGxlOiAi57Wx6KiI6Kej5p6QUua8lOe/kjIiCiNhdXRob3I6ICJUb21veXVraSBGdXJ1dGFuaSIKI2RhdGU6ICIyMDE5LzgvOSIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKICBodG1sX2RvY3VtZW50OiBkZWZhdWx0Ci0tLQoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkKa25pdHI6Om9wdHNfY2h1bmskc2V0KGNhY2hlID0gVFJVRSkKa25pdHI6Om9wdHNfY2h1bmskc2V0KGVycm9yID0gVFJVRSkKYGBgCgojIyDmvJTnv5Ljga7mpoLopoEKLSDph43lm57luLDliIbmnpAKCiMjIyDms6jmhI/kuovpoIUKLSDmhbbmh4nnvqnlob7lpKflraZTRkPjgafplovorJvjgZfjgabjgYTjgovjgIzntbHoqIjop6PmnpDjgI3jga7mjojmpa3lsaXkv67ogIXlkJHjgZHjga7mvJTnv5LnlKjjg5rjg7zjgrjjgafjgZnjgIIKLSDlv4XjgZrjgZfjgoLlhajjgabjga7jg5Djg7zjgrjjg6fjg7Pjga5S44KET1Pjgafli5XkvZznorroqo3jgpLooYzjgaPjgabjgYTjgb7jgZvjgpPjgILjgZPjga7mvJTnv5LnlKjjg5rjg7zjgrjjgpLkvZzmiJDjgZfjgabjgYTjgovmrrXpmo7jgafjga/jgIFSMy42LjDjgpLkvb/jgaPjgabjgYTjgb7jgZnjgIIKLSBS44Gu5pu05paw44Gq44Gp44Gr44KI44KK44CBUuOCs+ODvOODieOBquOBqeOCkuS6iOWRiueEoeOBl+OBq+WkieabtOOBmeOCi+WgtOWQiOOBjOOBguOCiuOBvuOBmeOAgiAKCiMjIOODh+ODvOOCv+WIhuaekOOBrua6luWCmQojIyMg44OH44O844K/44Gu44OA44Km44Oz44Ot44O844OJCi0g5ryU57+S44Gn55So44GE44KL44OH44O844K/44GvW+OBk+OBk10oaHR0cDovL3dlYi5zZmMua2Vpby5hYy5qcC9+bWF1bnovRFNCMTkvRFNCZGF0YS56aXAp44GL44KJ44OA44Km44Oz44Ot44O844OJ44GX44Gm44GP44Gg44GV44GE44CCCgojIyMg5L2c5qWt55So44OV44Kp44Or44OA44Gu5oyH5a6aCi0g44OA44Km44Oz44Ot44O844OJ44GX44Gfemlw44OV44Kh44Kk44Or44KS6Kej5YeN44GX44CB5L2c5oiQ44GV44KM44Gf44OV44Kp44Or44OA44KS5L2c5qWt55So44OV44Kp44Or44OA44Go44GX44Gm5oyH5a6a44GX44Gm44GP44Gg44GV44GE44CCYHNldHdkKClg6Zai5pWw44KS5L2/44GG44Gq44Gp44GX44Gm44CB5L2c5qWt55So44OV44Kp44Or44OA44KS5oyH5a6a44Gn44GN44G+44GZ44CC5L6L44GI44Gw44CBSETjga7nm7TkuIvjgatgRFNCZGF0YWDjg5Xjgqnjg6vjg4DjgpLkvZzmiJDjgZfjgZ/loLTlkIjjga/ku6XkuIvjga7jgojjgYbjgavjgarjgorjgb7jgZnjgIIKCmBgYHtyIHNldGVudjEsIGVjaG89VFJVRSwgZXZhbD1GQUxTRX0Kc2V0d2QoIn4vRFNCZGF0YSIpCmBgYAoKYGBge3Igc2V0ZW52MiwgZWNobz1UUlVFLCBldmFsPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQpzZXR3ZCgifi9Hb29nbGUgRHJpdmUgRmlsZSBTdHJlYW0v44Oe44Kk44OI44KZ44Op44Kk44OV44KZL0xlY3R1cmUgJiBFeGNlY2lzZS8yMDE5XzIwNF/np4vjg7vntbHoqIjop6PmnpAv5o6I5qWt6LOH5paZUm1kXzIwMTkwODA5L0RTQmRhdGEiKQpgYGAKCiMjIyDjg4fjg7zjgr/jga7oqq3jgb/ovrzjgb8KLSDjgI7lrrboqIjoqr/mn7vjgI8yMDE35bm0OeaciOOAgeesrO+8kuKIku+8luihqOOAjOW5tOmWk+WPjuWFpemajue0muWIpe+8keS4luW4r+W9k+OBn+OCiu+8keOBi+aciOmWk+OBruWPjuWFpeOBqOaUr+WHuuOAjeOBi+OCieS9nOaIkOOBl+OBn+aJgOW+l+majue0muOBqOaUr+WHuuOBruODh+ODvOOCv2BrYWtlaTIwMTcwOV9sMTguY3N2YOOCkueUqOOBhOOBvuOBmeOAggotIOOBk+OBk+OBp+OBr+OAgWByZWFkLmNzdigpYOmWouaVsOOCkuS9v+OBo+OBpuOAgUNTVuODleOCoeOCpOODq+OCkuiqreOBv+i+vOOBv+OBvuOBmeOAggoKYGBge3IgcmVhZGRhdGExLCBlY2hvPVRSVUUsIGV2YWw9RkFMU0UsIG1lc3NhZ2U9VFJVRSwgd2FybmluZz1UUlVFfQprYWtlaSA8LSByZWFkLmNzdigia2FrZWkyMDE3MDlfbDE4LmNzdiIsIGhlYWRlcj1UUlVFLCBzZXA9IiwiKQpgYGAKCmBgYHtyIHJlYWRkYXRhMiwgZWNobz1GQUxTRSwgZXZhbD1UUlVFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQprYWtlaSA8LSByZWFkLmNzdigifi9Hb29nbGUgRHJpdmUgRmlsZSBTdHJlYW0v44Oe44Kk44OI44KZ44Op44Kk44OV44KZL0xlY3R1cmUgJiBFeGNlY2lzZS8yMDE5XzIwNF/np4vjg7vntbHoqIjop6PmnpAv5o6I5qWt6LOH5paZUm1kXzIwMTkwODA5L0RTQmRhdGEva2FrZWkyMDE3MDlfbDE4LmNzdiIsIGhlYWRlcj1UUlVFLCBzZXA9IiwiKQpgYGAKCi0g6YeN5Zue5biw5YiG5p6Q44Gn55So44GE44KL77yT44Gk44Gu5aSJ5pWwYElOQzFgLCBgQ09OUzFgLCBgV09SS2Djga7mlaPluIPlm7PjgpLjg5fjg63jg4Pjg4jjgZfjgb7jgZnjgIIKYGBge3IgcGxvdDIsIGZpZy5oZWlnaHQ9MS41LCBmaWcud2lkdGg9MS41fQpwbG90KGtha2VpWyw0OjZdKQpgYGAKCiMjIOWbnuW4sOWIhuaekAojIyMg6YeN5Zue5biw5YiG5p6QCi0gYGxtKClg6Zai5pWw44KS5L2/44Gj44Gm44CB57ea5b2i5Zue5biw44Oi44OH44Or77yI5Y2Y5Zue5biw44Oi44OH44Or77yJ44KS5pyA5bCP5LqM5LmX5o6o5a6a44GX44G+44GZ44CCCi0g44G+44Ga44CBYElOQzFg44Gu44G/44KS6Kqs5piO5aSJ5pWw44Go44GZ44KL5Y2Y5Zue5biw44Oi44OH44Or44KS5o6o5a6a44GX44G+44GZ44CCCmBgYHtyIGxtMSwgZWNobz1UUlVFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPVRSVUV9Cmtha2VpLmxtMSA8LSBsbShDT05TMX5JTkMxLCBkYXRhPWtha2VpKQpzdW1tYXJ5KGtha2VpLmxtMSkKYGBgCgotIOasoeOBq+OAgWBJTkMxYOOBqGBXT1JLYOOCkuiqrOaYjuWkieaVsOOBqOOBmeOCi+mHjeWbnuW4sOODouODh+ODq+OCkuaOqOWumuOBl+OBvuOBmeOAggpgYGB7ciBsbTIsIGVjaG89VFJVRSwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1UUlVFfQprYWtlaS5sbTIgPC0gbG0oQ09OUzF+SU5DMStXT1JLLCBkYXRhPWtha2VpKQpzdW1tYXJ5KGtha2VpLmxtMikKYGBgCgotIOaui+W3ruODl+ODreODg+ODiOOAgVFR44OX44Ot44OD44OI44CBUy1M44OX44Ot44OD44OI44CB5qKD5a2Q5q+U44Go44Kv44OD44Kv44Gu6Led6Zui77yIUi1M44OX44Ot44OD44OI77yJ44KS5Lul5LiL44Gu44KI44GG44Gr44GX44Gm5Zuz56S644GX44G+44GZ44CCCmBgYHtyIHJlc2lkcGxvdCwgZmlnLmhlaWdodD0xLjUsIGZpZy53aWR0aD0xLjV9CnBsb3Qoa2FrZWkubG0yKQpgYGAKCi0g5q+U6LyD44Gu44Gf44KB44Gr44CBYFdPUktg44Gu44G/44KS6Kqs5piO5aSJ5pWw44Go44GZ44KL5Y2Y5Zue5biw44Oi44OH44Or44KS5o6o5a6a44GX44G+44GZ44CCCmBgYHtyIGxtMywgZWNobz1UUlVFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPVRSVUV9Cmtha2VpLmxtMyA8LSBsbShDT05TMX5XT1JLLCBkYXRhPWtha2VpKQpzdW1tYXJ5KGtha2VpLmxtMykKYGBgCgojIyMg5qiZ5rqW5YyW44GX44Gf44OH44O844K/44KS55So44GE44Gf6YeN5Zue5biw5YiG5p6QCi0gYHNjYWxlKClg6Zai5pWw44KS55So44GE44KL44Go44CB44OH44O844K/44KS5qiZ5rqW5YyW44Gn44GN44G+44GZ44CCCmBgYHtyIHNjYWxlbG0sIGVjaG89VFJVRX0Kei5rYWtlaSA8LSBhcy5kYXRhLmZyYW1lKHNjYWxlKGtha2VpKSkKa2FrZWkubG00IDwtIGxtKHoua2FrZWlbLDVdfnoua2FrZWlbLDRdK3oua2FrZWlbLDZdKQpzdW1tYXJ5KGtha2VpLmxtNCkKYGBgCgojIyMgQUlD44Gr44KI44KL44Oi44OH44Or6ZaT44Gu5q+U6LyDCi0gYEFJQygpYOmWouaVsOOCkueUqOOBhOOBpuOAgeODouODh+ODq+aOqOWumue1kOaenOOCkuavlOi8g+OBp+OBjeOBvuOBmeOAggpgYGB7ciBBSUMsIGVjaG89VFJVRX0KQUlDKGtha2VpLmxtMSwga2FrZWkubG0yLCBrYWtlaS5sbTMpCmBgYAoKIyMjIEJJQ+OBq+OCiOOCi+ODouODh+ODq+mWk+OBruavlOi8gwotIGBCSUMoKWDplqLmlbDjgpLnlKjjgYTjgabjgIHjg6Ljg4fjg6vmjqjlrprntZDmnpzjgpLmr5TovIPjgafjgY3jgb7jgZnjgIIKYGBge3IgQklDLCBlY2hvPVRSVUV9CkJJQyhrYWtlaS5sbTEsIGtha2VpLmxtMiwga2FrZWkubG0zKQpgYGAKCiMjIyDlpJrph43lhbHnt5rmgKcKLSBgY2FyYOODkeODg+OCseODvOOCuOOBrmB2aWYoKWDplqLmlbDjgpLnlKjjgYTjgabjgIHlpJrph43lhbHnt5rmgKfjgpLjg4Hjgqfjg4Pjgq/jgafjgY3jgb7jgZnjgIJgaW5zdGFsbC5wYWNrYWdlcygiY2FyIilg44Gn44OR44OD44Kx44O844K444KS44Kk44Oz44K544OI44O844Or44GX44Gm44GP44Gg44GV44GE44CCCmBgYHtyIFZJRiwgZWNobz1UUlVFLCB3YXJuaW5nPVRSVUV9CmxpYnJhcnkoY2FyKQpjYXI6OnZpZihrYWtlaS5sbTIpCmBgYAotIFZJRuOBjDEw5Lul5LiL44Gu44Go44GN5aSa6YeN5YWx57ea5oCn44GM44Gq44GE44Go5Yik5pat44GV44KM44KL44GT44Go44GL44KJ44CB44GT44Gu77yS5aSJ5pWw44Gu6ZaT44Gr5aSa6YeN5YWx57ea5oCn44GM44GC44KL44Go44Gv44GE44GI44Gq44GE44CCCgrku6XkuIrjgafku4rlm57jga7mvJTnv5Ljga/ntYLkuobjgafjgZnjgIIKCgo=