回帰モデルと内生性とR
この記事は
R言語 Advent Calendar12日目の記事です
書いていること
- 回帰分析の話
- なぜなら私が文字に起こす上で一番楽しいのが回帰分析だから。
- その中でも「内生性」に関する話をする
- Rで内生性回避のための操作変数法を実装する
- 二段階最小二乗法
- 一般化モーメント法
- [Advanced] 内部操作変数法
実装
回帰分析の設定
回帰モデルと内生性
末石(2015)に従って、教育年数と賃金の関係に関する回帰モデルを考えます。
を賃金、
を教育年数、
を個人の能力とした場合に、以下のようなモデル(賃金方程式)を設定したとしましょう。
$$
\begin{equation}
\log{W_i} = \beta_0 + \beta_1S_i + \beta_2 A_i+ u_i
\end{equation}
$$
ただし、各観測で各変数は独立同分布に従い、
は誤差項、
誤差項に関する条件期待値は0、すなわち = 0]を仮定します。
経済学や社会学の領域では教育年数と賃金の関係については数多くの先行研究が蓄積されていますし、
直感的にも「教育年数が長いほど、賃金が高い」という関係性の仮説は想像できると思います。
同様に、個人の能力と教育年数も何らかの関係性を持つことが考えられます。
教育年数が長いほど、個人の能力が高い、という仮説は想定可能だからです。
その関係性をこの回帰モデルで考察する場合、の値が具体的にどうであるかが関心の対象となります。
通常、このモデルは最小二乗法(Ordinary Least Square; OLS)で推定され、OLS推定量は不偏性、一致性、漸近正規性の3つの性質を持ちます *1。
ここで、この個人の能力は観測できていないと仮定し、誤差項と個人の能力の影響が区別できない状態を考えます。
すなわち誤差項をと置いた以下のモデルを考えます。
$$
\begin{equation}
\log{W_i} = \beta_0 + \beta_1S_i + e_i
\end{equation}
$$
上記のモデルは、
- 誤差項として観測できていない「個人の能力」が含まれる
- 個人の能力は教育年数と相関関係があると考えられる
ので、教育年数と誤差項との間に相関が生じている構造になります。
結果として、このモデルが本来満たすべき=0]を満たしません。
この状態を内生性(Endogeneity)と呼び、内生性を持つモデルではOLSによる推定は一致性を持ちません
具体的には、の推定値
は
に確率収束することが示せます。
操作変数法
パラメータ推定値が不偏性を持たない場合、同じ設定でデータを取得して追試を行った場合に全く異なった値が推定される可能性を否定できず、一致性を持たない場合はデータの量を増やすと全く異なる値に推定されてしまう可能性を否定できないので、特に学術的には避けたい状況です。
こうした状態を回避する方法として操作変数法という手法が知られています。
そこで、回帰モデルに対して、 以下の条件を満たすような変数
を考えます*2。
つまり「内生変数と相関して、誤差項
と相関しないような変数」をうまく探してモデルに反映できれば良いということになります。
実装の1パターンとして2段階最小二乗法があります。これはざっくりいうと以下のような手順を実現することになります。
- 内生変数を操作変数で回帰(
)
- 1で得た推定値
を説明変数にして
に回帰
これによって得られる推定量は2段階最小二乗推定量と呼ばれ、操作変数の条件を満たしていれば、 この推定量は一致性を持つことが知られています。
Rによる確認
データ
Wooldredgeのサイトからmorzデータを取得します。
上記の設定に近いデータをえることができる点でデータを採用します。wageには欠損値がありますが、今回は欠損値を削除して用います。
今回は操作変数としてfatheduc(父親の教育年数)、motheduc(母親の教育年数)を扱います。いわゆる教育の再生産。
# use libraries library(tidyverse) library(foreign) library(AER) library(stargazer) library(sem) # 二段階最小二乗法がある library(momentfit) # 一般化モーメント法 library(REndo) # Internal Instrumental Variable Method set.seed(20231212) mroz_rawdata <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/mroz.dta") mroz_usedata <- mroz_rawdata %>% dplyr::mutate(log_wage = log(wage)) %>% dplyr::select(wage, log_wage, educ, exper, fatheduc, motheduc) %>% dplyr::filter(!is.na(wage))
通常の重回帰
# 普通の重回帰 lm_model <- lm( log(wage) ~ educ + fatheduc + motheduc, data = mroz_usedata)
興味のある変数はeducで、今回は0.12くらい、ということが確認できます。
Call:
lm(formula = log(wage) ~ educ + fatheduc + motheduc, data = mroz_usedata)
Residuals:
Min 1Q Median 3Q Max
-2.98160 -0.31189 0.06732 0.40075 2.16340
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.134268 0.186513 -0.720 0.472
educ 0.123928 0.016131 7.683 1.09e-13 ***
fatheduc -0.008041 0.011589 -0.694 0.488
motheduc -0.018082 0.012177 -1.485 0.138
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.6779 on 424 degrees of freedom
Multiple R-squared: 0.1275, Adjusted R-squared: 0.1213
F-statistic: 20.66 on 3 and 424 DF, p-value: 1.646e-12
二段階最小二乗法
Rで二段階最小二乗法を行うにはsem::tsls関数もありますが、
結局内生変数を操作変数で回帰して、その予測値を用いる点で変わりはないので、
今回はlmでゴリ推します。
tls_first <- lm(educ ~ exper + fatheduc + motheduc, data = mroz_usedata) mroz_usedata$educ_fitted <- fitted(tls_first) tls_second <- lm(log(wage) ~ educ_fitted + exper, data = mroz_usedata)
結果、educ_fittedの効果を見るのですが0.050で、普通の回帰分析と比較すると値自体は半分程度になっていますね。
# 1段階目の推定結果
Call:
lm(formula = educ ~ fatheduc + motheduc, data = mroz_usedata)
Residuals:
Min 1Q Median 3Q Max
-7.8319 -0.9413 0.1086 1.0835 6.4865
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.48014 0.32111 29.523 < 2e-16 ***
fatheduc 0.18810 0.03363 5.593 4.01e-08 ***
motheduc 0.15637 0.03582 4.365 1.60e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.039 on 425 degrees of freedom
Multiple R-squared: 0.2081, Adjusted R-squared: 0.2043
F-statistic: 55.83 on 2 and 425 DF, p-value: < 2.2e-16
# 2段階目の推定結果
Call:
lm(formula = log_wage ~ educ_fitted, data = mroz_usedata)
Residuals:
Min 1Q Median 3Q Max
-3.2293 -0.3801 0.0531 0.4138 2.0675
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.55102 0.42580 1.294 0.196
educ_fitted 0.05049 0.03352 1.506 0.133
Residual standard error: 0.7221 on 426 degrees of freedom
Multiple R-squared: 0.005297, Adjusted R-squared: 0.002962
F-statistic: 2.268 on 1 and 426 DF, p-value: 0.1328
一般化モーメント法
上記では触れていませんが一般化モーメント法も同様に操作変数法の推定量として用いられることが多いです。
その理由は二段階最小二乗推定量よりも効率的に推定ができるという点にあります。
詳細は末石(2015)を参照してください。
gmm_model <- momentfit::momentModel( log(wage) ~ educ, educ ~ fatheduc + motheduc, data = mroz_usedata ) gmm_fit <- momentfit::gmmFit(gmm_model) gmm_fit
educの係数を見ると、2段階最小二乗法で推定した結果と一致します。
Model based on moment conditions ********************************* Moment type: linear Covariance matrix: iid Number of regressors: 2 Number of moment conditions: 3 Number of Endogenous Variables: 1 Sample size: 428 Estimation: Two-Stage Least Squares coefficients: (Intercept) educ 0.55102048 0.05049048
内部操作変数法
これは最近Gui et al.(2023)などが整理・Rでの実装(REndoパッケージ)を提案している手法で、
「外部での操作変数の選定や評価は難しく、弱操作変数の問題もあるため、外部変数に頼らずに内生性を補正できないか？」という試みです。
この補正の方法はいくつか提案があるようですが、今回は内生変数が一つという設定なので、 Ebbs et al(2005)で提案されている潜在操作変数法(Latent Instrumental Variables; LIV)を使ってみます。
# latent instrumental variable IIV_model <- REndo::latentIV( log(wage) ~ educ, optimx.args = list(itnmax = 50000), data = mroz_usedata, )
結果は以下のようになります。今回興味のある効果はeducですが、上記の結果ともまた違うっぽい。
Call:
REndo::latentIV(formula = log(wage) ~ educ, data = mroz_usedata,
optimx.args = list(itnmax = 50000))
Coefficients:
Estimate Std. Error z-score Pr(>|z|)
(Intercept) 0.30935 0.34496 0.897 0.37034
educ 0.07100 0.02707 2.623 0.00903 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Further parameters estimated during model fitting:
pi1 pi2 theta5 theta6 theta7 theta8
11.7334 15.9946 0.7767 0.4621 0.1285 1.9577
(see help file for details)
Initial parameter values:
(Intercept)=-0.1852 educ=0.1086 pi1=12.6589 pi2=14.9443
theta5=0.5 theta6=1 theta7=0.5 theta8=1
The value of the log-likelihood function: 1379.178
AIC: -2742.356 , BIC: -2709.883
通常の回帰分析の結果よりは、TSLSやGMMでの推定結果に寄っている結果ですが、それでもちょっと違う値ですね。
終わりに
内生性とMMM
MMMに関する内生性への言及は安永ら(2022)のレビュー論文に以下のように記載されています。
内生性が生じる例として，独立変数がモデル内の別の変数から影響を受けるときが挙げられる． 例えば，(1) 式の独立変数である xi,t をテレビ CM の出稿量，従属変数である yi を売上とするモデルがあるとしよう．このとき，過去の売上をもとに テレビ CM の出稿量を決定するデータ生成プロセスが存在すると内生性が生じる．
これの回避のためには多変量時系列回帰モデルの提案をしているものの、その実装上にも課題が残るというふうに述べています。
さいごに
操作変数法も難しいし、内生性も難しい。
あと付け焼き刃で色々試そうとするとボロが出るので良くない！！
でも記事執筆間に合ったので良し！Enjoy！！
参考文献
Kei Sakamoto, Instrumental Variables Method (操作変数法) with R
Deepblue社テックブログ
Ebbes, P., Wedel, M., Böckenholt, U. et al. Solving and Testing for Regressor-Error (in)Dependence When no Instrumental Variables are Available: With New Evidence for the Effect of Education on Income. Quant Market Econ 3, 365–392 (2005). https://doi.org/10.1007/s11129-005-1177-6
Gui, R., Meierer, M., Schilter, P., & Algesheimer, R. (2023). REndo: Internal Instrumental Variables to Address Endogeneity. Journal of Statistical Software, 107(3), 1–43. https://doi.org/10.18637/jss.v107.i03
末石直也, 2015, 『計量経済学 ミクロデータ分析へのいざない』, 日本評論社
TJO,MMM (Media/Marketing Mix Modeling)を回すなら、まずGeorge E. P. Boxの格言を思い出そう
西山慶彦, 新谷元嗣, 川口大司, 奥井亮, 2019, 『計量経済学』, 有斐閣
