『Stataによる計量政治学』（浅野正彦, 矢内勇生)：Rによる分析
Yuki Yanai
September 19, 2014 (rev. June 6, 2015)
第11章 回帰分析による統計的推定
準備として、作業ディレクトリを指定する。 作業ディレクトリの指定にはsetwd()関数を利用する。 ディレクトリへのパス（引用符の中身）は、各自で設定する。
setwd("~/projects/Stata-R/")この章で使うパッケージをあらかじめ読み込んでおく。
library("foreign")
library("dplyr")##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("ggplot2")
## Windows ユーザは次の行をコメントアウト
theme_set(theme_gray(base_size=12, base_family="HiraKakuProN-W3"))11.1 単回帰による統計的推定
11.1.1 単回帰モデル
まず、気温とビールの出荷量のデータセット beer2010.dta を読み込む。 データは、作業ディレクトリ内のdataというフォルダ（ディレクトリ）に保存されているとする。 読み込んだデータセットはBeerと名付ける。
Beer <- read.dta("data/beer2010.dta")10章と同じように、lm() を使って回帰分析を行おう。
fit.beer <- lm(beer ~ temp, data = Beer)
summary(fit.beer)##
## Call:
## lm(formula = beer ~ temp, data = Beer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.860 -20.242 0.882 11.037 76.939
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.546 21.735 4.028 0.00241
## temp 2.104 1.165 1.805 0.10117
##
## Residual standard error: 31.98 on 10 degrees of freedom
## Multiple R-squared: 0.2458, Adjusted R-squared: 0.1704
## F-statistic: 3.259 on 1 and 10 DF, p-value: 0.1012これで、表11.1と同じ結果が得られた。表の読み方については教科書 (p.200) を参照してほしい。 Stataの結果と比較すると、統計量の並び方が異なることと、Stataでは_consと表示されていた切片が Intercept と表示されていることがわかる。
11.1.2 信頼区間と仮説検定
回帰係数の信頼区間
Rで回帰係数の信頼区間を求めるには、lm() の結果に対して、confint() 関数を使う。
confint(fit.beer, "temp", level = 0.95)## 2.5 % 97.5 %
## temp -0.4928032 4.70094これで、表11.2 (p.202) と同じ95パーセント信頼区間が得られた。 confint() の中の level の値を変えることによって、様々なパーセンテージの信頼区間を求めることができる。
この信頼区間を図示してみよう。
scat.beer.1 <- ggplot(Beer, aes(temp, beer)) +
geom_point(size = 4) +
labs(x = "気温（℃）", y = "ビールの出荷量 (1000kl)",
title = "ビール出荷量と平均気温の散布図に信頼区間を上書きする")
scat.beer.1 + geom_smooth(method = "lm")
このように、図11.2 (p.203) とほぼ同じ図ができる。
おまけ：シミュレーションによって信頼区間を求める
信頼区間は、シミュレーションによって求めることもできる。 ここでは、arm パッケージのsim() 関数を使ってシミュレーションを行ってみよう。
## 必要なら arm パッケージをインストールする
# install.packages('arm')
library("arm")
fit.beer.sim <- sim(fit.beer, n.sims = 10000) ## 10,000回のシミュレーションこれで、fit.beer.sim にシミュレーションによって得られた1万組の回帰係数が保存された。 ここから回帰係数を取り出すには、coef(fit.beer.sim) とする。 coef(fit.beer.sim)は10000行（シミュレーションの回数分）\(\times\) 2列（係数の数：ここでは切片と気温temp）の行列である。 したがって、気温（第2列）の回帰係数だけを取り出すには、coef(fit.beer.sim)[, 2] とする。
シミュレーションの結果を上で行った回帰分析の結果と比べてみよう。
mean(coef(fit.beer.sim)[, 2]) ## sim() で得られた回帰係数の平均値## [1] 2.101922この値は、fit.beerの回帰係数2.1040682 に十分近い値であり、シミュレーションは問題なく行われたようである（シミュレーションの値は毎回異なる。毎回同じ値を得るには、set.seed(数値) で乱数の種を固定する）。
このシミュレーション結果から、気温の回帰係数の95パーセント信頼区間を求めよう。 10000個ある回帰係数のうち、中心にある95パーセント分の範囲を求めれば、それが私たちが知りたい95パーセント信頼区間となる。 95パーセント信頼区間の下限と上限はそれぞれ最小値から2.5%, 97.5% の位置にあるの、quantile() 関数を使って次のように入力する。
quantile(coef(fit.beer.sim)[,2], c(.025, .975))## 2.5% 97.5%
## -0.5275604 4.7046837これで、シミュレーションによる95％信頼区間が得られた。 上で計算したものとは少し違う値（になるはず）であるが、まったく違うというわけでもない。
最後に、このシミュレーション結果を図示する。
plot(NULL, ## 図の枠を用意する
xlim = c(min(Beer$temp), max(Beer$temp)),
ylim = c(min(Beer$beer), max(Beer$beer)),
xlab = "気温", ylab = "ビールの出荷量 (1000kl)",
main = "ビール出荷量と気温の関係：シミュレーション")
## for ループでシミュレートした回帰直線を図示する
## 描画に時間がかかるときは、10000を小さい数に変えてもよい
for (i in 1:10000) {
curve(coef(fit.beer.sim)[i, 1] + coef(fit.beer.sim)[i ,2] * x,
add = TRUE, col = "gray")
}
## 回帰直線を赤で描く
curve(coef(fit.beer)[1] + coef(fit.beer)[2]*x, add = TRUE,
col = "red", lwd = 2)
## 観測値のプロットを上書きする
points(beer ~ temp, data = Beer, pch = 16)
この図では、赤い直線が回帰直線、その周りにある灰色の部分（実際には、灰色の直線の集合）が推定の不確実性を表現している。 したがって、上で描いた信頼区間付きの図と同じような役割を果たすことができる（ここではシミュレーションで得られた回帰係数をすべて使って図を描いており、灰色部分は95%信頼区間より広い）。
仮説検定
仮説検定で利用する値は既に得られている。もう一度表示すると、
summary(fit.beer)##
## Call:
## lm(formula = beer ~ temp, data = Beer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.860 -20.242 0.882 11.037 76.939
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.546 21.735 4.028 0.00241 **
## temp 2.104 1.165 1.805 0.10117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 31.98 on 10 degrees of freedom
## Multiple R-squared: 0.2458, Adjusted R-squared: 0.1704
## F-statistic: 3.259 on 1 and 10 DF, p-value: 0.1012となっている。Stata のときと同じように、\(t\) 値と\(p\)値と利用して仮説検定の行えばよい（表11.3 とpp. 203-206の説明を参照）。 Rの初期設定では、一般的に使われる有意水準で帰無仮説を棄却できる \(p\) 値の右側に記号（* など）が付いている。この例では、切片の回帰係数の\(p\) 値の横に記号が付いており、1%の有意水準で「切片は0である」という帰無仮説を棄却できることが分かる（この例では、切片に関する仮説を棄却する意味はないが）。 この星が鬱陶しいときは、options() で記号が出ないようにする。
options(show.signif.stars = FALSE)
summary(fit.beer)##
## Call:
## lm(formula = beer ~ temp, data = Beer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.860 -20.242 0.882 11.037 76.939
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 87.546 21.735 4.028 0.00241
## temp 2.104 1.165 1.805 0.10117
##
## Residual standard error: 31.98 on 10 degrees of freedom
## Multiple R-squared: 0.2458, Adjusted R-squared: 0.1704
## F-statistic: 3.259 on 1 and 10 DF, p-value: 0.1012記号が表示されると、記号が付いているかどうかだけに気をとられ、結果を注意深く読まなくなる（ような気がする）ので、非表示にするのがお勧めである。毎回オプションを設定するのが面倒な場合は、.Rprofile に書き加えておけばよい（.Rprofile については ココ や ココ などを参照）。
11.2 重回帰分析による統計的推定
11.2.1 重回帰モデル
分析用のデータ (hr96-09.dta) を読み込み、内容を確認する （第5章のRによる分析で作った hr96-09labeled.dat を使ってもよい）。
HR <- read.dta("data/hr96-09.dta")
## hr96-09labeled.dat を使う場合は、
# load("data/hr96-09labeled.dat")
## データが適切に読み込まれたどうか確認したいなら、
# head(HR)今回の分析で利用する2009年の自民党議員のデータのみ抜き出し、HR09LDPとして保存する。
HR09LDP <- filter(HR, year == 2009, party == "LDP")重回帰分析を行い、結果と信頼区間を表示しよう。
fit.ldp.2009 <- lm(voteshare ~ I(exp / 1000000) + previous, data = HR09LDP)
summary(fit.ldp.2009)##
## Call:
## lm(formula = voteshare ~ I(exp/1e+06) + previous, data = HR09LDP)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.378 -4.671 -0.242 4.246 30.139
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.3824 1.3094 26.258 < 2e-16
## I(exp/1e+06) 0.1064 0.1084 0.982 0.327
## previous 1.3748 0.1595 8.619 4.62e-16
##
## Residual standard error: 8.078 on 287 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.229, Adjusted R-squared: 0.2236
## F-statistic: 42.63 on 2 and 287 DF, p-value: < 2.2e-16confint(fit.ldp.2009)## 2.5 % 97.5 %
## (Intercept) 31.8052069 36.9596769
## I(exp/1e+06) -0.1069589 0.3198126
## previous 1.0608636 1.6887899これで、表11.4 (p.207) と同じ結果が得られた。 表の読み方、仮説検定の方法については、教科書の「11.2.2 信頼区間と仮説検定」(pp.208-212) を参照。