第10章　回帰分析の基礎

準備として、作業ディレクトリを指定する。 作業ディレクトリの指定には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, union
library("ggplot2")
## Windows ユーザは次の行をコメントアウト
theme_set(theme_gray(base_size=12, base_family="HiraKakuProN-W3"))


10.1 線形回帰：散布図への直線の当てはめ

まず、気温とビールの出荷量のデータセット（beer2010.dta）を読み込む。 データは、作業ディレクトリ内のdataというフォルダ（ディレクトリ）に保存されているとする。 Stata形式 (.dta) で保存されたデータを読むために、foreign パッケージを使う。 読み込んだデータセットはBeerと名付ける。

Beer <- read.dta("data/beer2010.dta")

データが適切に読み込めたかどうか確認する。

Beer    ## データをすべて表示する
##    month temp    beer
## 1      1  7.0  59.415
## 2      2  6.5  81.349
## 3      3  9.1 109.654
## 4      4 12.4 123.673
## 5      5 19.0 109.576
## 6      6 23.6 156.601
## 7      7 28.0 160.500
## 8      8 29.6 148.631
## 9      9 25.1 115.968
## 10    10 18.9 105.968
## 11    11 13.5 120.190
## 12    12  9.9 185.316

month, temp, beer という3つの変数が含まれていることがわかる。 このデータセットは観測数が12しかないのでこの方法でかまわないが、観測数が大きいデータセットの場合、この方法ではデータ全体が表示されてしまうので、あまり嬉しくない。 データ全体を表示させたくないときは、

names(Beer)    ## データセットに含まれる変数名を表示する
## [1] "month" "temp"  "beer"
head(Beer)     ## データセットの一部（最初の数行）を表示する
##   month temp    beer
## 1     1  7.0  59.415
## 2     2  6.5  81.349
## 3     3  9.1 109.654
## 4     4 12.4 123.673
## 5     5 19.0 109.576
## 6     6 23.6 156.601
dim(Beer)      ## データセットに含まれる観測数と変数の数を確認する (dimはdimension)
## [1] 12  3

などが利用できる。


続いて、散布図を描く。 最も簡単な方法は、plot() 関数を使う方法である。

plot(beer ~ temp, data = Beer,
     ylab = "ビール出荷量（1000kl）", xlab = "気温（℃）", 
     main="月別のビール出荷量と平均気温の散布図")

これで十分目的は果たしているが、より見やすい図を作りたいときはggplot2 を利用する。

scat.beer <- ggplot(Beer, aes(temp, beer)) +
    geom_point(size = 4) +
    labs(x = "気温（℃）", y = "ビールの出荷量 (1000kl)",
         title = "月別のビール出荷量と平均気温の散布図")
scat.beer


次に、cor() 関数を使って相関係数を確かめる。

with(Beer, cor(beer, temp))

または

cor(Beer$beer, Beer$temp)

とする。その結果、beer とtemp の相関係数は0.4957856であることがわかる。


2変数間の直線的関係を散布図上に示すには、次のようにする。

scat.beer + geom_smooth(method = "lm", se = FALSE)


10.2 最小二乗法

R で線形回帰分析を行うには、lm() を使う。 使い方は、lm(応答変数 ~ 説明変数1 + 説明変数2 + ・・・, data = データの名前)である（応答変数と説明変数の間はハイフンではなくチルダ）。 ビールの出荷量 (beer) を応答変数、気温 (temp) を説明変数とする単回帰を model1 と呼ぶことにすると、

model.1 <- lm(beer ~ temp, data = Beer)

というコマンドで、分析結果がmodel1に保存される。


分析結果の表示は、次のように行う。

summary(model.1)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 87.5464 21.7352 4.03 0.0024
temp 2.1041 1.1655 1.81 0.1012


結果の読み方はStataとほとんど同じである（p.184の表10.2を参照）。 この結果から、切片の推定値は87.5463973、傾きの推定値は2.1040682であることがわかる。


10.3 単回帰と重回帰

10.3.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)
##   year    ku kun party              name age     status nocand   wl rank
## 1 1996 aichi   1   NFP KAWAMURA, TAKASHI  47  incumbent      7  win    1
## 2 1996 aichi   1   LDP     IMAEDA, NORIO  72       moto      7 lose    2
## 3 1996 aichi   1   DPJ     SATO, TAISUKE  53  incumbent      7 lose    3
## 4 1996 aichi   1   JCP   IWANAKA, MIHOKO  43 challenger      7 lose    4
## 5 1996 aichi   1 bunka       ITO, MASAKO  51 challenger      7 lose    5
## 6 1996 aichi   1    NP  YAMADA, HIROSHIB  51 challenger      7 lose    6
##   previous  vote voteshare eligible turnout     exp
## 1        2 66876      40.0   346774   49.22 9828097
## 2        3 42969      25.7   346774   49.22 9311555
## 3        2 33503      20.1   346774   49.22 9231284
## 4        0 22209      13.3   346774   49.22 2177203
## 5        0   616       0.4   346774   49.22      NA
## 6        0   566       0.3   346774   49.22      NA

今回の分析で利用する2009年の自民党議員のデータのみ抜き出し、HR09LDPとして保存する。

HR09LDP <- filter(HR, year == 2009, party == "LDP")


データの準備ができたので、重回帰分析を行う。 ここでは、応答変数を得票率 (voteshare)、説明変数を選挙費用 (exp) と当選回数 (previous) とする重回帰分析を行う（model2と名付ける）。

model.2 <- lm(voteshare ~ exp + previous, data = HR09LDP)

分析結果は次の通りである。

summary(model.2)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.3824 1.3094 26.26 0.0000
exp 0.0000 0.0000 0.98 0.3271
previous 1.3748 0.1595 8.62 0.0000

結果の読み方は、Stataの場合と同様である（p.187, 表10.3 を参照）。


選挙費用 (exp) の推定結果がほぼ0になっているが、これは選挙費用の測定単位（1円）が小さすぎるためである。 そこで、選挙費用の単位を100万円に変えて再分析してみる（model3と名付ける）。

model.3 <- lm(voteshare ~ I(exp / 1000000) + previous, data = HR09LDP)

Stataで同様の分析を行う場合、expを1,000,000で割ったexpmという変数を用意する必要があったが、RではI()関数を使ってlm() の中で直接計算することができる。 分析用のフォーミュラ (formula) の中で通常の演算を行いたいときはこの I() 関数を利用する（これを使わないと、通常の演算と認識してくれない。例えば、model3 のlm()の中には+があるが、これは説明変数を並べるための記号であり、加算とは認識されていない）。

分析結果は次の通りである。

summary(model.3)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.3824 1.3094 26.26 0.0000
I(exp/1e+06) 0.1064 0.1084 0.98 0.3271
previous 1.3748 0.1595 8.62 0.0000

この結果の読み方は、表10.4（p.188）の読み方と同様である。


10.3.2 単回帰と重回帰の違い

Rで残差を求めるには、lm() の返り値を residuals() の引数として与えればよい。 たとえば、上で分析したmodel3 の残差を求めるには

res3 <- residuals(model.3)

とすればよい。

『Stataによる計量政治学』のページに戻る