Ｒで重回帰分析
佐藤慧
2019年10月14日
なんだかんだいってまだやったことがなかったので、諸々参考にしつつ、標準データセットでやってみました。
今回使用したのは、swissデータセット。1888年のスイスのフランス語圏の47自治体の社会経済的変数（出生率など六つ）が収められています。詳しくはhelp(swiss)で。それにしても『自殺論』チックですね～。
記述統計
とりあえず記述統計を確認しましょう。
summary(swiss)## Fertility Agriculture Examination Education
## Min. :35.00 Min. : 1.20 Min. : 3.00 Min. : 1.00
## 1st Qu.:64.70 1st Qu.:35.90 1st Qu.:12.00 1st Qu.: 6.00
## Median :70.40 Median :54.10 Median :16.00 Median : 8.00
## Mean :70.14 Mean :50.66 Mean :16.49 Mean :10.98
## 3rd Qu.:78.45 3rd Qu.:67.65 3rd Qu.:22.00 3rd Qu.:12.00
## Max. :92.50 Max. :89.70 Max. :37.00 Max. :53.00
## Catholic Infant.Mortality
## Min. : 2.150 Min. :10.80
## 1st Qu.: 5.195 1st Qu.:18.15
## Median : 15.140 Median :20.00
## Mean : 41.144 Mean :19.94
## 3rd Qu.: 93.125 3rd Qu.:21.70
## Max. :100.000 Max. :26.60まあ各変数そこそこいい感じでばらついてるんじゃないでしょうか。
ヒストグラムとかも確認したいですが、今回調べて知ったのが、psych::pairs.panels関数。ヒストグラム、散布図行列、相関行列をすべてひとつの図の中に収めてくれます。
psych::pairs.panels(swiss)
わかりやすい！
回帰分析
とりあえず今回は慣れるのが目的なので仮説も何もなしに適当にやります。
従属変数はExamination（軍の試験で最高レベルで合格したその地域の出身者の数）がよさげです。
独立変数としてはAgriculture（男性労働者のうち農業従事者の割合）とEducation（中等教育以上を受けた人たちの教育年数の合計）にしてみましょう。たぶん、農業人口の多いところでは高い教育は必要がなく、それゆえ試験の値も低い、といった関係があるでしょう。
そこでモデル1として農業人口だけ投入、次にモデル2で教育程度を加え、標準偏回帰係数の値を比べたいと思います。
分析の前に
しかしこのままだと、予想される係数の正負が、農業人口はマイナスですが教育程度はプラスになってしまいます。これだと結果が直感的にわかりにくいので、農業人口を反転（最小値が最大値になるような感じで）させた新しい変数を作って分析に使いたいと思います。
数学だめなので少し考えたのですがわかりました。最小値に、最大値から当該の値を引いた数を足せばいいんですね。例としては以下のような感じです。
## ag ag2
## 1 1 9
## 2 1 9
## 3 7 3
## 4 8 2
## 5 9 1たとえば\[ag2[4] = min(ag2)+(max(ag)-ag[4]) = min(ag)+max(ag)-ag[4]\]と表せます。
この塩梅でAgricultureを反転させた変数は以下のように作りました。
for (i in 1:nrow(swiss)) {
swiss$Agriculture.reverse[i] <- min(swiss$Agriculture) + max(swiss$Agriculture) - swiss$Agriculture[i]
}summary()で比較すると、ちゃんと反転できてるみたいです。
summary(swiss)## Fertility Agriculture Examination Education
## Min. :35.00 Min. : 1.20 Min. : 3.00 Min. : 1.00
## 1st Qu.:64.70 1st Qu.:35.90 1st Qu.:12.00 1st Qu.: 6.00
## Median :70.40 Median :54.10 Median :16.00 Median : 8.00
## Mean :70.14 Mean :50.66 Mean :16.49 Mean :10.98
## 3rd Qu.:78.45 3rd Qu.:67.65 3rd Qu.:22.00 3rd Qu.:12.00
## Max. :92.50 Max. :89.70 Max. :37.00 Max. :53.00
## Catholic Infant.Mortality Agriculture.reverse
## Min. : 2.150 Min. :10.80 Min. : 1.20
## 1st Qu.: 5.195 1st Qu.:18.15 1st Qu.:23.25
## Median : 15.140 Median :20.00 Median :36.80
## Mean : 41.144 Mean :19.94 Mean :40.24
## 3rd Qu.: 93.125 3rd Qu.:21.70 3rd Qu.:55.00
## Max. :100.000 Max. :26.60 Max. :89.70実際の分析
では実際に分析します。lm()関数を使います。
モデル1
lm1 <- lm(data = swiss, Examination ~ Agriculture.reverse)
summary(lm1)##
## Call:
## lm(formula = Examination ~ Agriculture.reverse, data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.1324 -3.0931 -0.0264 4.2817 10.5378
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.78477 1.75476 3.866 0.000352 ***
## Agriculture.reverse 0.24117 0.03807 6.334 9.95e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.865 on 45 degrees of freedom
## Multiple R-squared: 0.4713, Adjusted R-squared: 0.4596
## F-statistic: 40.12 on 1 and 45 DF, p-value: 9.952e-08係数は0.1%水準で有意です。農業人口が1%上がると軍の試験がよくできる人が0.24人増えるらしい。決定係数も0.45くらいはありますね。
しかしこの効果は、教育程度に媒介されているはずです。果たして農業人口の効果は残るのか。
モデル2
lm2 <- lm(data = swiss, Examination ~ Agriculture.reverse + Education)
summary(lm2)##
## Call:
## lm(formula = Examination ~ Agriculture.reverse + Education, data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.8751 -2.9823 -0.3439 4.4592 10.3453
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.75446 1.57203 4.297 9.43e-05 ***
## Agriculture.reverse 0.14258 0.04437 3.214 0.00245 **
## Education 0.36410 0.10479 3.474 0.00116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.254 on 44 degrees of freedom
## Multiple R-squared: 0.5852, Adjusted R-squared: 0.5663
## F-statistic: 31.03 on 2 and 44 DF, p-value: 3.922e-09結果は1%水準ではありますが有意です。予想通り農業人口の標準偏回帰係数は下がっており、教育程度のほうが大きな規定力を持っています。また決定係数は0.56へ上がっています。
ただし農業人口の効果も残っています。教育程度が同程度の地域でも農業人口が少ないほうが軍の試験の成績がいい人が多いようです。
この場合の農業人口は、工業化の度合いに基づく高学歴化の背景というよりは、都鄙度を表すものと解釈するのが適切ではないでしょうか。少しでも都会に近いほうが、たとえば試験勉強用の特別な教育を受ける機会も多い、というような関係が考えられます。
結果まとめ
memisc::mtable()関数で簡単にモデルを並列させた結果の表が作れるようです。ミーミスクか。
library(memisc)
mtable(lm1, lm2)##
## Calls:
## lm1: lm(formula = Examination ~ Agriculture.reverse, data = swiss)
## lm2: lm(formula = Examination ~ Agriculture.reverse + Education, data = swiss)
##
## =============================================
## lm1 lm2
## ---------------------------------------------
## (Intercept) 6.785*** 6.754***
## (1.755) (1.572)
## Agriculture.reverse 0.241*** 0.143**
## (0.038) (0.044)
## Education 0.364**
## (0.105)
## ---------------------------------------------
## R-squared 0.471 0.585
## N 47 47
## =============================================
## Significance: *** = p < 0.001;
## ** = p < 0.01; * = p < 0.05まあこんなところでしょうか。
今回の復習
- 記述統計量の把握に便利な
psych::pairs.panels()関数 - 回帰分析は
lm(data = data, data$dependent ~ data$independent)で - 回帰分析の結果を表にするのに便利な
memisc::mtable()関数