第10章 回帰分析の基礎
渡邊 翔吾
2019/2/28
浅野・矢内(2018)に基づき、「第10章 回帰分析の基礎」をまとめていきます。
コマンドやデータは、オーム社のHPからDLしたものを元にしています。
浅野・矢内『Rによる計量政治学』(2018年,オーム社)
第10章〜第11章 回帰分析
データの読み込みから相関分析まで
一応、一連の前処理
- パッケージの読み込み
library(tidyverse) - データの読み込み
HR <- readr::read_rds("data/hr96-17.Rds") #データの読み込み- データの整形
HR2012 <- dplyr::filter(HR, year == 2012) #2012年時点の選挙に限定して、データを整形
HR2 <- dplyr::filter(HR2012, party %in% c("LDP", "DPJ", "CGP", "SDP", "JCP")) #主要政党に限定してデータを整形。今後は、 HR2 を使用する
今回の目的
変数関係の関係性を視覚化するために、plotを練習してきた。
次のステップとして、複数の変数の関係性を示すことにしよう。
相関関係
変数間の関係性を表す統計量として、相関関係が第一に挙げられる
相関(correlation) / 相関関係(correlation coefficient)：
二変量関係を要約する際に最もよく使われる統計量で、二変数がどの程度関係しているかを測る尺度である。
cor()関数
Rでは、cor()コマンドを使って相関関係を計算できる。
- サンプル①の作成
a <- c(1,2,3)
b <- c(3,6,9)- corの実施
cor(a,b) #当たり前だが完全に正の相関を示す## [1] 1cor(b,a) #相関なので、順番が逆でも結果は変わらない## [1] 1- サンプル②の作成・実施
##サンプル②
c <- c(9,6,3)
cor(a,c) #完全に負の相関を示す## [1] -1d <- c(3.4, 1.9, 7.8) #適当
cor(a,d) #やや強い正の相関を示すようだ## [1] 0.7174337cor()の練習
ではここで、衆院選データを使って、任意の変数を分析してみよう。
ピアソンの積率相関係数を求める
「use」欠損値サンプルの除去を指定（リストワイズ）
「method」でP(ピアソン)を指定
with(HR2, cor(voteshare, expm,
use = "complete.obs", method = "p")) ## [1] 0.6649999スピアマンの順位相関係数を求める
「use」欠損値サンプルの除去を指定（リストワイズ）
「method」でS(スピアマン)を指定
with(HR2, cor(HR2$voteshare, HR2$expm,
use = "complete.obs", method = "s"))## [1] 0.7473294どうやら、ピアソンでもスピアマンでも正の相関を示しているようであることがわかる。
二変数間の関係性を、-1 ~ +1 の値で読み解くことができる。
しかしながら、単に相関関係を見るだけでは不十分である。
疑似相関の可能性
相関関係を分析している際には、常に擬似相関の可能性があることを意識するべきである。
疑似相関とは、見かけ上の相関のことである。
数値上、相関があることを示すも、理論上は異なっていたりすることである。
擬似相関の様々な例を見て、確認してみよう。
最小二乗法による回帰分析
最小二乗法による回帰分析にはlm(linear model)コマンドを使う
lmコマンドのシンタックスは以下の通り
結果を格納するobject <- lm( DV ~ IV, data = データフレーム)
lmコマンドの実践
まずは、得票率と選挙費用(exp / expm)を使って回帰分析してみる
fit_1 <- lm(voteshare ~ expm, data = HR2) #得票率 ⇐ 選挙費用(100万円)
#「fit_1」 という名前で、lm()関数の結果を格納。（結果は表示されないものの）fit_1に、lm結果が格納された。
名前には、自分で覚えやすいものを使うこと。
また、スクリプトには、#で分析内容を簡単にコメントしておくと良い。
回帰分析の結果
次に、中身を確認していく。
- 推定結果全体を確認
summary(fit_1) #推定結果のすべてを表示##
## Call:
## lm(formula = voteshare ~ expm, data = HR2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -66.289 -7.623 -3.628 7.031 66.442
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.4538 0.7874 12.01 <2e-16 ***
## expm 2.7359 0.1041 26.29 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.81 on 872 degrees of freedom
## (10 observations deleted due to missingness)
## Multiple R-squared: 0.4422, Adjusted R-squared: 0.4416
## F-statistic: 691.4 on 1 and 872 DF, p-value: < 2.2e-16- 係数のみの確認
coef(fit_1) #係数のみを表示。見易い。## (Intercept) expm
## 9.453805 2.735949- 決定係数のみを確認
summary(fit_1)$r.squared #R-squaredのみを表示## [1] 0.4422249- 自由度修正済決定係数のみを確認
summary(fit_1)$adj.r.squared #自由度修正済決定係数のみを表示## [1] 0.4415852以上から、「得票率= 9.4 + 2.7*選挙費用(100万円) 」の推定式が得られる。
R-Squared は 44%
F検定のP値は、2.2e-16 と極小
選挙費用変数の推定値に関しては、SEが0.1、t値が26、p値が2e-16 であることから、1%水準で有意であると判断できる。
重回帰分析の実践
次に、IVを複数投入する重回帰分析を実践してみる。 重回帰分析のシンタックスは、lmの単回帰分析に、** + IV2, + IV3, +・・・**と付け足すだけ。
- 気になる変数をとりあえず投入してみよう。
fit_2 <- lm(voteshare ~ expm + age, data = HR2) #得票率 ⇐ 選挙費用(100万円) , 年齢
fit_3 <- lm(voteshare ~ expm + age + previous, data = HR2) #得票率 ⇐ 選挙費用(100万円) 、 年齢、 当選回数- summaryで結果を確認
summary(fit_2)##
## Call:
## lm(formula = voteshare ~ expm + age, data = HR2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -66.936 -7.839 -3.502 7.391 66.622
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.54064 2.29832 2.411 0.0161 *
## expm 2.73250 0.10394 26.290 <2e-16 ***
## age 0.07629 0.04211 1.812 0.0703 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.79 on 871 degrees of freedom
## (10 observations deleted due to missingness)
## Multiple R-squared: 0.4443, Adjusted R-squared: 0.443
## F-statistic: 348.2 on 2 and 871 DF, p-value: < 2.2e-16結果から、expmは1%水準で正に有意だが、ageは有意ではないことがわかる。
細かな検討はここでは割愛
summary(fit_3)##
## Call:
## lm(formula = voteshare ~ expm + age + previous, data = HR2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.504 -6.329 -1.786 5.693 49.004
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.18332 2.09693 10.102 < 2e-16 ***
## expm 1.67405 0.10340 16.190 < 2e-16 ***
## age -0.23709 0.03899 -6.081 1.78e-09 ***
## previous 3.72999 0.19531 19.098 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.58 on 870 degrees of freedom
## (10 observations deleted due to missingness)
## Multiple R-squared: 0.6085, Adjusted R-squared: 0.6071
## F-statistic: 450.7 on 3 and 870 DF, p-value: < 2.2e-16結果から、expm、age、previousそれぞれが1%水準で有意であることがわかる。
先程、有意でなかったageが有意になっていることに着目。
重回帰分析では、複数の変数を同時に投入するが、それぞれ単独で回帰分析ををしている訳ではない
他の変数が一定である場合にその変数がどれほどの影響を与えるかを示している。
細かな検討はここでは割愛
しかし、何とも出力結果が分かりづらい。。。
推定結果の綺麗な出力
ここで、閑話休題
はっきりいって、結果出力が最も重要。
論文に載せる表はなんとしても綺麗にしたい。
にもかかわらず、浅野・矢内(2018)はそれに対応していない。
そこで調べてみた結果、綺麗な推定表を出力するパッケージを幾つか発見したので、ぜひ使ってみてほしい。
出力用のパッケージ
ともかく、パッケージをインストールしてみよう。
# install.packages("texreg") #パッケージのインストール
library(texreg)では、インストールしたterregコマンドを使っていく。
せっかくなので、任意でいくつかの回帰分析を行っておこう。
5つほどOLSの結果を格納しておこう。
fit_4 <- lm(voteshare ~ expm + age + previous + nocand, data = HR2)
fit_5 <- lm(voteshare ~ exppv + age + previous, data = HR2)screenregによるConsole画面への出力
まずは、シンプルに、しかし綺麗に出力してみよう。
screenreg(list(fit_1, fit_2, fit_3, fit_4, fit_5)) #最もシンプルに出力##
## =======================================================================
## Model 1 Model 2 Model 3 Model 4 Model 5
## -----------------------------------------------------------------------
## (Intercept) 9.45 *** 5.54 * 21.18 *** 36.83 *** 21.97 ***
## (0.79) (2.30) (2.10) (2.62) (2.01)
## expm 2.74 *** 2.73 *** 1.67 *** 1.66 ***
## (0.10) (0.10) (0.10) (0.10)
## age 0.08 -0.24 *** -0.25 *** -0.25 ***
## (0.04) (0.04) (0.04) (0.04)
## previous 3.73 *** 3.69 *** 3.75 ***
## (0.20) (0.19) (0.19)
## nocand -3.40 ***
## (0.37)
## exppv 0.54 ***
## (0.03)
## -----------------------------------------------------------------------
## R^2 0.44 0.44 0.61 0.64 0.63
## Adj. R^2 0.44 0.44 0.61 0.64 0.63
## Num. obs. 874 874 874 874 874
## RMSE 13.81 13.79 11.58 11.06 11.32
## =======================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05何回も出力結果を表示させるなら、オブジェクトで一連のfitを格納しておくと便利
out <- list(fit_1, fit_2, fit_3, fit_4, fit_5)
screenreg(out) #出力は割愛するテキストデータへの出力
次に、論文にそのまま掲載できるように、テキスト形式で出力してみよう。
ここで、シンタックスとオプションの細かな説明は省くが、とりあえずは大丈夫。
とにかく綺麗に出してみよう。
- tex形式
texreg(list(fit_1, fit_2, fit_3, fit_4, fit_5),file="texreg.tex", caption="texreg", digits=3, booktabs=T, dcolumn=T, center=T, use.packages=F, caption.above=T, custom.model.names=c("Model1", "Model2","Model3","Model4","Model5"))- html形式
htmlreg(list(fit_1, fit_2, fit_3, fit_4, fit_5),file="htmlreg.html", caption="htmlreg", digits=3, booktabs=T, dcolumn=T, center=T, use.packages=F, caption.above=T, custom.model.names=c("Model1", "Model2", "Model3","Model4","Model5"))それぞれ、出力結果が、Filesに追加されているはずだ。
各々、確認してみよう。
ここでのhtml出力結果は以下のようになっている。
Rからだと、即座にhtml形式で出力できるので便利。
論文に掲載可能な綺麗な推定結果を出力することができた
しかし、texregも限界がある。
reg(regression)とあるように、あくまでlm()による回帰分析にしか使うことができないようだ。
他の分析手法を使った際や、記述統計について出力するには他のパッケージが必要。 インストールし、使い込んで慣れていこう。
#install.packages("psych", "xtable", "Hmisc", "stargazer")
#library(psych)
#library(xtable)
#library(Hmisc)
#library(stargazer)
##今後、これらパッケージについても慣れていこうとりあえずはここまで。