計量経済学応用
回帰分析：交差項の利用と複合仮説の検定
矢内 勇生
2019-05-288
準備
今回使うパッケージを読み込む。インストール済みでなパッケージは、install.packages() でインストールしてから使う。
Macユーザのみ次の行も実行。
交差項の利用
浅野正彦・矢内勇生 (2018)『Rによる計量政治学』第14章の内容に沿って説明する。
データを手に入れる
浅野・矢内 (2018) で用いられているデータhr-data.RdsのCSV版 をダウンロードし、プロジェクト内のdataディレクトリに保存する。
dir.create("data")
download.file(url = "https://raw.githubusercontent.com/yukiyanai/quant-methods-R/master/data/hr-data.csv",
destfile = "data/hr-data.csv")ダウンロードできたら、データを読み込む。
データの中身を確認する。
## Observations: 8,803
## Variables: 22
## $ year <dbl> 1996, 1996, 1996, 1996, 1996, 1996, 1996, 1996, 1996,…
## $ ku <chr> "aichi", "aichi", "aichi", "aichi", "aichi", "aichi",…
## $ kun <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3,…
## $ status <chr> "現職", "元職", "現職", "新人", "新人", "新人", "新人", "現職", "元職",…
## $ name <chr> "KAWAMURA, TAKASHI", "IMAEDA, NORIO", "SATO, TAISUKE"…
## $ party <chr> "NFP", "LDP", "DPJ", "JCP", "others", "kokuminto", "i…
## $ party_code <dbl> 8, 1, 3, 2, 100, 22, 99, 8, 1, 3, 2, 10, 100, 99, 22,…
## $ previous <dbl> 2, 3, 2, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 3, 1,…
## $ wl <chr> "当選", "落選", "落選", "落選", "落選", "落選", "落選", "当選", "落選",…
## $ voteshare <dbl> 40.0, 25.7, 20.1, 13.3, 0.4, 0.3, 0.2, 32.9, 26.4, 25…
## $ age <dbl> 47, 72, 53, 43, 51, 51, 45, 51, 71, 30, 31, 44, 61, 4…
## $ nocand <dbl> 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 7, 7, 7,…
## $ rank <dbl> 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3,…
## $ vote <dbl> 66876, 42969, 33503, 22209, 616, 566, 312, 56101, 449…
## $ eligible <dbl> 346774, 346774, 346774, 346774, 346774, 346774, 34677…
## $ turnout <dbl> 49.2, 49.2, 49.2, 49.2, 49.2, 49.2, 49.2, 51.8, 51.8,…
## $ exp <dbl> 9828097, 9311555, 9231284, 2177203, NA, NA, NA, 12940…
## $ expm <dbl> 9.828097, 9.311555, 9.231284, 2.177203, NA, NA, NA, 1…
## $ vs <dbl> 0.400, 0.257, 0.201, 0.133, 0.004, 0.003, 0.002, 0.32…
## $ exppv <dbl> 28.341505, 26.851941, 26.620462, 6.278449, NA, NA, NA…
## $ smd <chr> "当選", "落選", "落選", "落選", "落選", "落選", "落選", "当選", "落選",…
## $ party_jpn <chr> "新進党", "自民党", "民主党", "共産党", "その他", "国民党", "無所属", "新進党…2005年の衆院選について、voteshare, exppv, eligibleの3変数だけ残す。select() という名前の関数は、dplyrパッケージ (tidyverseで読み込んだ) と interplot パッケージの両者に含まれるので、dplyr::select() と書いて、dplyrのselect() を使う。
## 記述統計
2005年の衆院選データについて、上で選んだ3変数の基本的な統計量を確認する。
## voteshare exppv eligible
## Min. : 0.60 Min. : 0.148 Min. :214235
## 1st Qu.: 8.80 1st Qu.: 8.352 1st Qu.:297385
## Median :34.80 Median :22.837 Median :347866
## Mean :30.33 Mean :24.627 Mean :344654
## 3rd Qu.:46.60 3rd Qu.:35.269 3rd Qu.:397210
## Max. :73.60 Max. :89.332 Max. :465181
## NA's :4exppv に欠測値が4つあることがわかる。
3変数の標準偏差を求める。
## voteshare exppv eligible
## 19.23023 17.90740 63898.23016得票率 voteshare と 一人当たり選挙費用 exppv の関係を、散布図で確認する。
plt_vs_ex <- ggplot(HR05, aes(x = exppv, y = voteshare)) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "有権者一人当たり選挙費用", y = "得票率")
plot(plt_vs_ex)## Warning: Removed 4 rows containing non-finite values (stat_smooth).## Warning: Removed 4 rows containing missing values (geom_point).
有権者数 elgible と exppv の関係を、散布図を使って確認する。
plt_vs_el <- ggplot(HR05, aes(x = eligible, y = voteshare)) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "有権者数", y = "得票率")
plot(plt_vs_el)
交差項を使った回帰分析
有権者一人あたり選挙費用 exppv が得票率 voteshare に与える影響は、有権者数 eligible によって変わるだろうか。これを確かめるために、exppv と eligible の交差項 (interaction term) を回帰分析に投入し、以下のモデルを推定する。
\[ 得票率_i = \beta_0 + \beta_1 選挙費用_i + \beta_2 有権者数_i + \beta_3 選挙費用_i \times 有権者数_i + 誤差 \]
このモデルを`lm()で推定するためには、次のコードを使う。
つまり、変数1 * 変数2 という書き方をすると、重回帰分析に 変数1、変数2、変数1と変数2の交差項という3つの項が加えられる。
結果を確認してみよう。
##
## Call:
## lm(formula = voteshare ~ exppv * eligible, data = HR05)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.290 -8.158 -3.295 8.472 44.886
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.638e+00 3.793e+00 2.541 0.0112
## exppv 1.917e-02 1.140e-01 0.168 0.8665
## eligible -1.483e-06 1.073e-05 -0.138 0.8901
## exppv:eligible 2.549e-06 3.497e-07 7.289 6.44e-13
##
## Residual standard error: 12.58 on 981 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.572, Adjusted R-squared: 0.5707
## F-statistic: 437 on 3 and 981 DF, p-value: < 2.2e-16この結果は、どのように解釈できるだろうか。例えば、exppv の係数の推定値 0.0192 はどういう意味を持っているだおるか。通常の回帰分析と同様に考えると、他の変数の値を一定に保ったとき、 exppvが1単位増えると、得票率が0.0192 ポイント上がるという意味である。
しかし、他の変数を一定に保つことは可能だろうか？このモデルでは、特殊な状況を除き、それは不可能である。exppv の値を変えれば、exppv \(\times\) eligibleの値も変わってしまうからである。、他の変数を一定に保つことができるのは、eligble = 0 のときだけある。しかし、有権者数がゼロの選挙区は存在しないので、0.0192という数字をそのまま解釈することはできない。
この例からわかる通り、交差項を含む回帰分析の結果を、表から読み取ることは一般的に非常に難しい。これは、表を提示するだけでは、分析結果の報告として不十分だということを意味する。
先ほどのモデルでは説明変数をそのまま使っていた。次に、説明変数を中心化して、同様のモデルを推定してみよう。
HR05 <- HR05 %>%
mutate(exppv_c = exppv - mean(exppv, na.rm = TRUE),
eligible_c = eligible - mean(eligible))
model_2 <- lm(voteshare ~ exppv_c * eligible_c, data = HR05)
summary(model_2)##
## Call:
## lm(formula = voteshare ~ exppv_c * eligible_c, data = HR05)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.290 -8.158 -3.295 8.472 44.886
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.123e+01 4.187e-01 74.595 < 2e-16
## exppv_c 8.977e-01 2.514e-02 35.702 < 2e-16
## eligible_c 6.129e-05 6.576e-06 9.320 < 2e-16
## exppv_c:eligible_c 2.549e-06 3.497e-07 7.289 6.44e-13
##
## Residual standard error: 12.58 on 981 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.572, Adjusted R-squared: 0.5707
## F-statistic: 437 on 3 and 981 DF, p-value: < 2.2e-16この結果は、どのように解釈できるだろうか。例えば、exppv の係数の推定値 0.8977 はどういう意味を持っているだおるか。通常の回帰分析と同様に考えると、他の変数の値を一定に保ったとき、 exppv_cが1単位増えると、得票率が0.8977 ポイント上がるという意味である。
では、他の変数を一定に保つことは可能だろうか？このモデルでも、特殊な状況を除き、それは不可能である。exppv_c の値を変えれば、exppv_c \(\times\) eligible_c の値も変わってしまうからである。、他の変数を一定に保つことができるのは、eligble_c = 0 のときだけある。eligible_c は有権者数 eligible を中心化した変数なので、eligble_c = 0というのは、「有権者数が平均値のとき」という意味である。したがって、「有権者数が平均値のとき」という特殊な場合に限り、0.8977 という数字の意味を解釈することができる。
このように、説明変数を中心化することにより、少しだけ結果が解釈しやすくなった。しかし、有権者数が平気値のとき以外については、解釈が不可能である。よって、もう少し工夫が必要である。そのような工夫の一つが、図を使った方法である。
交差項の影響を可視化する
散布図に複数の回帰直線を描く
交差項の効果を可視化する方法を色々考えられるが、最も単純な（簡単とは限らない）方法は、調整変数である「有権者数」の値をいくつかの代表的な値に設定し、それらの値ごとに異なる回帰直線（特定の有権数の下で、得票率を有権者一人当たり選挙費用に回帰した場合の回帰直線）を求め、散布図上に回帰直線を上書きするという方法である。実際にやってみよう。
まず、有権者数の平均値と標準偏差を求める。
## [1] 344654.3## [1] 63898.23有権者の数が「平均 \(\pm\) 標準偏差」の場合について、回帰直線を求めることにする。 上で推定した model_1 の結果から、 \[ \widehat{得票率} = 9.6382872 + 0.0191673 選挙費用 + -1.4832785\times 10^{-6} 有権者数 + 2.5489967\times 10^{-6} 選挙費用 \times 有権者数 \] である。したがって、有権差数が平均 \(-\) 標準偏差の場合の回帰直線の切片は、
## (Intercept)
## 9.221848であり、傾きは、
## exppv
## 0.7348136同様に、有権差数が平均 \(+\) 標準偏差の場合の回帰直線の切片は、
## (Intercept)
## 9.03229であり、傾きは、
## exppv
## 1.060566これらを図示してみよう。
Mac用のコード
## Macの場合
plt_int <- ggplot(HR05, aes(x = exppv, y = voteshare)) +
geom_point(pch = 16) +
geom_abline(intercept = intercept1, slope = slope1,
linetype = "dashed") +
geom_abline(intercept = intercept2, slope = slope2) +
ylim(0, 100) +
labs(x = "選挙費用（有権者一人当たり：円）", y = "得票率 (%)") +
geom_text(label = "得票率 = 9.22 + 0.72・選挙費用\n(有権者数：平均 - 標準偏差)",
x = 65, y = 8, family = "HiraginoSans-W3") +
geom_text(label = "得票率 = 9.03 + 1.06・選挙費用\n(有権者数：平均 + 標準偏差)",
x = 40, y = 90, family = "HiraginoSans-W3")Windows用のコード
## Windowsの場合
plt_int <- ggplot(HR05, aes(x = exppv, y = voteshare)) +
geom_point(pch = 16) +
geom_abline(intercept = intercept1, slope = slope1,
linetype = "dashed") +
geom_abline(intercept = intercept2, slope = slope2) +
ylim(0, 100) +
labs(x = "選挙費用（有権者一人当たり：円）", y = "得票率 (%)") +
geom_text(label = "得票率 = 9.22 + 0.72・選挙費用\n(有権者数：平均 - 標準偏差)",
x = 65, y = 8) +
geom_text(label = "得票率 = 9.03 + 1.06・選挙費用\n(有権者数：平均 + 標準偏差)",
x = 40, y = 90)## Warning: Removed 4 rows containing missing values (geom_point).
限界効果を可視化する
次に、有権者数の変化によって、「選挙費用が得票率に与える影響」がどのように変化するかを可視化ししよう。調整変数（この例の有権者数）が特定の値のとき、説明変数1単位の変化が結果変数に与える影響のことを限界効果 (marginal effect) と呼ぶ。（交差項がない重回帰分析の場合、係数の推定値そのものが限界効果である。） 交差項がある場合の限界効果は、interplot::interplot() で図示できる。
int_1 <- interplot(m = model_1, var1 = "exppv", var2 = "eligible") +
labs(x = "有権者数",
y = "選挙費用が得票率に与える影響")
plot(int_1)
この図から、選挙費用が得票率に与える影響（縦軸）は、有権者数とともに大きくなることがわかる。言い換えると、選挙費用の影響は、有権者が少ない場合には相対的に小さく、有権者が多い場合には相対的に大きくなることがわかる。
このように、交差項を含む回帰分析を行う場合には、限界効果を図示することで、結果が解釈しやすくなる。
複合仮説の検定
データの準備
教科書（田中 [2015]）第6章の例を使って説明する。使用するデータは、6_1_income.csv である。このデータをサポートページ からデータのzipファイル をダウンロードして展開し、6_4_minshu.csv をプロジェクト内の data フォルダに保存する。
保存できたら、データを読み込む。
## Parsed with column specification:
## cols(
## exper = col_double(),
## exper2 = col_double(),
## yeduc = col_double(),
## income = col_double(),
## lincome = col_double()
## )回帰モデルと推定
以下のミンサー方程式を推定する。 \[ \ln(賃金) = \beta_0 + \beta_1 修学年数 + \beta_2 就業可能年数 + \beta_3 (修行可能年数)^2 + U \] ただし、Uは誤差項である。
このモデルの係数を、Rで推定する。
##
## Call:
## lm(formula = lincome ~ yeduc + exper + exper2, data = myd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.0389 -0.3214 0.1681 0.5124 2.1326
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.4855019 0.1107823 22.44 <2e-16
## yeduc 0.1175467 0.0070603 16.65 <2e-16
## exper 0.1961736 0.0074935 26.18 <2e-16
## exper2 -0.0063811 0.0003162 -20.18 <2e-16
##
## Residual standard error: 0.7983 on 4295 degrees of freedom
## Multiple R-squared: 0.2066, Adjusted R-squared: 0.206
## F-statistic: 372.8 on 3 and 4295 DF, p-value: < 2.2e-16これまでは、「修学年数は賃金に影響するか」のような疑問に答えるため、次のような仮説を検定してきた。
- 帰無仮説: \(\beta_1 = 0\)
- 対立仮説: \(\beta_1 \neq 0\)
この例では、修学年数 yeduc の係数の\(p\)値が \(2\times 10^{-16}\) 未満なので、有意水準5パーセントで上の帰無仮説は棄却される。
複合仮説とF検定
これに対し、「ミンサー方程式に就業可能年数とその二乗という2つの変数を加えることに意味はあるか」という疑問をもったとしよう。この疑問に答えるためには、次のような仮説が必要になる。
- 帰無仮説: \(\beta_2 = \beta_3 = 0\)
- 対立仮説: \(\beta2\) と \(\beta_3\) のうち少なくとも一方は0でない
このように複数の母数（パラメタ）について同時に考える仮説を複合仮説 と呼ぶ。複合仮説の検定には\(F\)分布の特徴を利用した\(F\)検定が使われる。
帰無仮説が正しい場合、 \[ \ln(賃金) = \beta_0 + \beta_1 修学年数 + U \] となるはずである。こんモデル（null モデル）を推定する。
##
## Call:
## lm(formula = lincome ~ yeduc, data = myd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6801 -0.2897 0.2211 0.5496 2.2060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.309553 0.100755 42.772 <2e-16
## yeduc 0.070772 0.007204 9.824 <2e-16
##
## Residual standard error: 0.8861 on 4297 degrees of freedom
## Multiple R-squared: 0.02197, Adjusted R-squared: 0.02174
## F-statistic: 96.52 on 1 and 4297 DF, p-value: < 2.2e-16このnull モデルと先ほど推定したモデル（model_mincer; ミンサー方程式）に違いがあるかどうか検定したい。そのために\(F\)検定を行う。
検定に利用する値（検定統計量 \(F\)）は、次の式によって求められる。 \[ F = \frac{(\mathrm{SSR}_{\mathrm{null}} - \mathrm{SSR}_{\mathrm{alt}}) / q }{\mathrm{SSR}_{\mathrm{alt}} / (n - k - 1)} \] ここで、SSR (squared sum of residuals) は残差平方和のことである（本によっては、RSS [residual sum of squared] とも表記される）。\(\mathrm{SSR}_{\mathrm{null}}\)は null モデルのSSR、\(\mathrm{SSR}_{\mathrm{alt}}\)は検証対象となる重回帰モデル（この例では、ミンサー方程式）のSSRである。また、\(q\)は帰無仮説に登場する等号の数（i.e., null モデルと検証対象となるモデルの右辺の項数の差）である。この例の場合は\(q = 2\)である。\(k\)は検証対象となる重回帰モデルで推定されるパラメタの数である。ミンサー方程式では、\(\beta_0\)、\(\beta_1\)、\(\beta_2\)、\(\beta_3\)の4つのパラメタを推定するので、\(k=4\)である。 \(n\)はサンプルサイズ（観測数）であり、この例では
## [1] 4299である。 よって、\(n - k - 1 = 4299 - 4 - 1 = 4295\)である。この値は、次のように取り出すこともできる。
## [1] 4295上の式を使って、\(F\)統計量を計算してみよう。
SSR_null <- sum((fit_mincer_null$residuals)^2)
SSR_alt <- sum((fit_mincer$residuals)^2)
q <- 2
(F_mincer <- ((SSR_null - SSR_alt) / q) / (SSR_alt / fit_mincer$df.residual))## [1] 499.7542検定に使う\(F\)統計量は、499.75 である。
帰無仮説が正しいとすると、この値は、自由度 (\(q\), \(n-k-1\)) の\(F\)分布にしたがうことが知られている。この特徴を利用して、検定を実施する。\(t\)分布とは異なり、\(F\)分布は0以上の値しかとらない。よって、片側検定を実施する。自由度が (2, 4295) の\(F\)分布は、以下の確率密度をもつ。
df <- tibble(x = seq(from = 0, to = 10, length.out = 1000)) %>%
mutate(dens = df(x, df1 = 2, df2 = 4295))
fdens <- ggplot(df, aes(x = x, y = dens)) +
geom_line() +
labs(x = "", y = "確率密度", title = "F(2, 4295)")
plot(fdens)
この分布で、上で計算した検定統計量の数値以上の値を取る確率があらかじめ決めた有意水準を下回るとき、帰無仮説を棄却する。有意水準を5%にする場合、この分布で\(F\)値が499.75以上になる確率が 0.05 未満なら、帰無仮説を棄却する。
上の図から、\(F\)値が499.75 以上になる確率はほぼゼロになると思われるが、念のために計算してみよう。
## [1] 7.461626e-196この値はほぼゼロであり、0.05よりは小さい。よって、5%の有意水準で帰無仮説は棄却される。
同様に、\(F\)検定の臨界値と検定統計量を比較することによって検定を行うこともできる。自由度 (2, 4295) の\(F\)分布の95パーセンタイル地点は、
## [1] 2.997823である。上で計算した検定統計量は 499.75 であり、この臨界値よりも大きいので、5%の有意水準で帰無仮説は棄却される。
以上の結果から、重回帰分析に就業可能年数とその二乗という二つの項を加えることは統計的には意味があることがわかる。
演習
課題3で分析した教科書の練習問題7-Bと7-Cについて、
- 交差項の効果を可視化してみよう
- \(F\)検定を実行してみよう