Rによる計量分析:データ解析と可視化 - 回帰分析の理解

55 

Loading....

Loading....

Loading....

Loading....

Loading....

全文

(1)

R

による計量分析:データ解析と可視化

回帰分析の理解 伊藤 岳 富山大学 経済学部 2017 年度後期 Email: gito@eco.u-toyama.ac.jp January 15, 2018 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 1 / 46

(2)

Agenda

1 線形回帰分析 2 回帰分析の解釈 3 回帰係数の検定 4 回帰分析の前提 5 内生変数とバイアス

(3)
(4)

確認:従属変数と独立変数

従属変数 (dependent variable) y: 説明される変数 (結果) ▶ 応答変数,被説明変数,結果変数 ▶ 独立変数 (independent variable) x: 説明する変数 (原因) ▶ 説明変数,予測変数 ▶ コントロール (調整/統制) 変数 ▶ 注目する独立変数以外で,結果に影響を与える変数 ▶ 統計学的/数学的には単にx ▶ 従属変数と独立変数は,分析者が (なんらかの理論的根拠に基づき) 想定・ 仮定する ▶ 「賃金(従属変数y)を教育年数(独立変数x)に回帰する」 ▶ 回帰分析のみによって,因果関係を示すことができる訳ではない

(5)

回帰式と誤差項

yi= xiβ + ϵi (1)

= β1xi1+· · · + βkxik+ ϵi, for i = 1, . . . , n,

回帰係数 βk: 独立変数と従属変数の関係を示す

▶ 第一種過誤・第二種過誤のβと混同しないように

(予測) 誤差項 (error term): ϵi= yi− ˆyi

▶ 従属変数の変化は,独立変数の変化によって完全に説明できる訳ではない

▶ 観察不可能な (モデルに取り込まれていない)「偶然」のような変数の影響

も受ける

(6)

回帰分析の記法

線形回帰モデルの数式表現

▶ 線形回帰モデルは,次のように表現できる

yi= xiβ + ϵi (2)

= β1xi1+· · · + βkxik+ ϵi, for i = 1, . . . , n,

▶ また,次のようにも表現できる

yi∼ N (xiβ, σ2), for i = 1, . . . , n. (3)

▶ あるいは,

y∼ N (xβ, σ2I), (4)

(7)

回帰分析の記法

x y o ˆ y = β0+ β1x 1 2 β0 β1 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 6 / 46

(8)

回帰分析の記法

x y o ˆ y = β0+ β1x yi ˆ yi ϵi

(9)

推定量と推定値

▶ 推定量 (estimator): 推定の結果を,標本{(yi, xi) : i = 1, . . . , n} の関数とし

て表現したもの

▶ OLS推定量βˆ1=

n

i=1(xi− ¯xn)(yi− ¯yn) n

i=1(xi− ¯xn)

2 (次スライド)

▶ 推定値 (estimate): (推定量に) 実際に標本を代入して計算した数値のこと

(10)

最小二乗法

▶ どのようにして (上の図の) 直線を引けば (回帰直線を推定すれば) よい?

▶ 最小二乗法 (ordinary least squares, OLS): 予測誤差 (残差 residual) の平方和 (二乗の和) を最小化することで,事で,観測値と予測値のズレを最小化す る方法

予測誤差 e は観測値 (点) y と予測値 (直線) ˆy のズレなので,そのズレが最

も小さくなる回帰直線を見つける

▶ 平均二乗誤差 (mean squared error, MSE):

MSE = 1 n(e 2 1, . . . , e 2 n) = 1 n ni=1 (yi− ˆyi)2 ▶ 以下は単回帰を想定した説明 (重回帰でも発想は同じ)

(11)

最小二乗法

▶ 回帰式yi= β0+ β1xi+ ϵiを考える.予測誤差eiについて, 1 ne 2 i = 1 n(yi− ˆyi) 2 = 1 n(yi− ˆβ0− ˆβ1xi) 2 (5) 1 n ni=1 e2i = 1 n ni=1 (yi− ˆβ0− ˆβ1xi)2 (6) ▶ これを最小化する回帰係数の値を,最小二乗推定量(OLS推定量)と呼ぶ ▶ MSEを最小化するには,回帰係数で偏微分して「= 0」と置いた連立方程式を解け ばよい(最適化の1階条件) ∂MSE ∂ ˆβ0 =2 n ni=1 (yi− ˆβ0− ˆβ1xi) = 0 (7) ∂MSE ∂ ˆβ1 =2 n ni=1 xi(yi− ˆβ0− ˆβ1xi) = 0 (8) ▶ この連立方程式を解けば,OLS推定量が得られる ˆ β0= ¯yn− ˆβ1x¯n (9) ˆ β1= ∑n

i=1(xi− ¯xn)(yi− ¯yn) n

i=1(xi− ¯xn)

2 (10)

(12)

最小二乗法

行列を使えば,OLS 推定量は ˆβ = (xx)−1xy と書ける ▶ 重回帰分析の場合はこうした行列・ベクトルを用いた表記が一般的 ▶ OLS 推定量と相関係数 ▶ βˆ1= ∑n

i=1(xi− ¯xn)(yi− ¯yn) n i=1(xi− ¯xn) 2 ▶ rxy= ∑n i=1(xi− xn)(yi− yn) √∑n i=1(xi− xn)2 √∑n i=1(yi− yn)2 ▶ 違いは分母:独立変数の分散 (回帰係数) か両者の標準偏差 (相関係数) か ▶ 「回帰分析のみで実証できること」は因果関係ではない(とは限らない) ▶ 回帰係数にせよ相関係数にせよ,2 つの変数の「変化 (variation) の関係」を 見るもの ▶ 「変化の関係」をみるのだから,両者に十分な変化がなければならない! ▶ 副読本森田本の例:市区町村の失業者数が,1つだけ500でそれ以外は1,000 の場合?

(13)
(14)

回帰直線の「当てはまり」具合

決定係数 R2:観察された従属変数の全変動 (∑n i=1(yi− y)2) のうち,何 %が予測値の変動 (回帰変動∑ni=1( ˆyi− y)2) で説明できるかを示す値 R2= ∑n i=1( ˆyi− y)2 ∑n i=1(yi− y)2 = 1 ∑n i=1(yi− ˆyi)2 ∑n i=1(yi− y)2 (11) あるいは,残差変動∑ni=1(yi− ˆyi)2が全変動に対してどれだけ小さいか ▶ R2= 1 なら回帰分析によって従属変数を完全に説明でき (残差変動が 0), R2= 0 なら全く説明できない (残差変動が 1) ▶ ただし,R2は独立変数を追加すれば,それが何であれ (無意味なものでも) 必ず増加する ▶ 独立変数の数 k の違いを調整するには,自由度調整済み決定係数 R2 aを用い る (副読本・田中本第 6 章) R2a= 1 ∑n

i=1(yi− ˆyi)

2/(n− 1 − k)

∑n

i=1(yi− y)2/(n− 1)

(15)

単回帰分析と重回帰分析

▶ 単回帰分析 ▶ 従属変数の原因を単一の独立変数に求める回帰分析 例 y = βˆ 0+ β1x1+ ϵ ▶ 重回帰分析 ▶ 従属変数の原因を複数の独立変数に求める回帰分析 例 y = βˆ 0+ β1x1+ β2x2+ β3x3+ ϵ ▶ 重回帰分析は,複数の単回帰分析を行なっている訳ではない y′= β0+ β1x1+ ϵ y′′= β0+ β2x2+ ϵ ではない ▶ 重回帰分析は,他の独立変数を一定に保ったときに,ある独立変数が従属変 数に与える影響を求めている ▶ ceteris paribus(他の条件が一定とする) ▶ 「他の独立変数のうち他の独立変数とは関係のない部分が,従属変数のうち他 の独立変数とは関係のない部分に与える影響」(浅野・矢内本, p.192) 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 13 / 46

(16)

回帰係数と限界効果

線形回帰における (偏) 回帰係数の解釈 ˆ y = β0+ β1x1+ β2x2+ β3x3+ ϵ (13) という回帰式を考える ▶ この回帰式において,独立変数 xkの偏回帰係数 βkは,他の独立変数の値 を一定に保ったとき,xkが 1 (単位) 増加したときの従属変数 (の予測値) ˆy の増加量 ▶ 重回帰分析における「偏回帰係数」の「偏」は偏微分の「偏 (partial)」と同 じ意味 ▶ つまり,「偏」は「他の (独立) 変数の値を一定に保ったとき」の意味 ▶ 「xkが 1 (単位) 増加したときの従属変数 (の予測値) ˆy の増加量 (影響の大 きさ)」のことを,限界効果 (marginal effect)という ▶ 回帰係数と同じに思えるが,より高度な一般化線形モデル(本講義では扱わな い)の理解や,後述の交互作用項を含むモデルの理解には限界効果が非常に重要

(17)

回帰係数と限界効果

線形回帰における (偏) 回帰係数の解釈 ▶ 解析では,統計的有意性 (statistical significance)と同程度以上に,実質的有 意性 (substantial significance; 効果の大小)が重要 ▶ 実質的有意性を解釈する上では,限界効果が鍵になる ▶ 一般化線形モデル(identity link以外)では線形回帰のような直感的な解釈が困 難なため,限界効果の計算・作図がより重要になる 限界効果と (偏) 回帰係数 ▶ 「偏」の字から直感的にわかる通り,xkの限界効果は (13) 式を xkで偏微 分すれば求められる.(13) 式を x1偏微分すれば, ∂ ˆy ∂xk = βk (14) ▶ 具体的な解釈は,xkと y の「単位」(“1” が何を意味するか) に依存する 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 15 / 46

(18)

回帰係数と限界効果

限界効果,実質的有意性の重要さ ▶ よく聞く表現:「この変数はトウケイテキニユーイだった.大事だ」 ▶ 気をつけるべき点:統計的に有意だからといって,実質的に有意とは限らな い (「統計的に有意でも,大して影響のない」変数もある) 「他の独立変数が一定のとき」の意味をよく考える ▶ 常に「一定に保つ」ことができる訳ではない ▶ 例:二乗項や交互作用項/交差項(interaction term)が含まれる回帰式 ▶ 交互作用(interaction):ある独立変数の効果が,他の独立変数の値によって異 なること ▶ いずれの場合でも,ある独立変数 xk「だけ」を動かして,その限界効果を 提示することは不可能 (掛け算してるから!) ▶ 観察データ (observational data) では,往々にして独立変数間に相関がある

(19)

回帰分析の解釈:ダミー変数

▶ ダミー変数 (dummy variable):特定の場合に 1 をとり,他の場合には 0 をと る二値の変数 例 大学院卒 (大学院卒なら 1,その他なら 0) ▶ 0 の値をとるグループを,基準グループ (reference group) と呼ぶ x y o ˆ y = β0+ β1x ˆ y = β0+ β1x + β2 β2 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 17 / 46

(20)

回帰分析の解釈:変数変換

▶ 収益率や弾力性等を計算したい場合,変数を対数変換する ▶ 対数変換した場合,以下のように (近似的に) 解釈できる 例 コブ・ダグラス型生産関数,貿易の重力モデル (対数線形モデル) 従属変数 独立変数 回帰係数 β の解釈 y x x が 1 単位増加したとき,y が β 単位増 加する ln y x x が 1 単位増加したとき,y が 100× β% 増加する (半弾力性 semi-elasticity) β = ∆y/y ∆x y ln x x が 1% 増加したとき,y が β/100 単位 増加する ln y ln x x が 1% 増加したとき,y が β% 増加す る (弾力性elasticity) β = ∆y/y ∆x/x

(21)

補足:対数関数の性質

ln(1) = 0 (15) ln(xr) = r ln(x) (16) ln(x1x2) = ln(x1) + ln(x2) (17) ln ( x1 x2 ) = ln(x1)− ln(x2) (18) d ln(x) dx = 1 x (19) 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 19 / 46

(22)

限界効果の計算:交互作用項がある場合

交互作用項 (interaction term) とは

x1× x2のような,「かけ算」で表現される (独立) 変数

x1× x2は “product term,” x1と x2は “constitutive term(s)” (main effect)

▶ 理論的議論:Berry et al. (2010, 2012); Brambor et al. (2006); Clark et al. (2006); Kam & Franzese (2007); Rainey (2015)

交互作用項の解釈と限界効果 ▶ よくある間違い (1):x1 (e.g., 投薬の有無) や x2(e.g., 性別) の交互作用項 x1× x2を回帰式に投入したから,x1と x2をバラバラに回帰式に投入する 必要はない! ▶ 誤:y = β0+ β3x1x2+ ϵ ▶ 正:y = β0+ β1x1+ β2x2+ β3x1x2+ ϵよくある間違い (2):交互作用項の回帰係数 β3が有意だった! だから,x1 や x2単独ではなく,その組み合わせが重要! (もうちょっと作業が必要)

(23)

限界効果の計算:交互作用項がある場合

▶ 交互作用項がある場合,product term も constitutive term も,単一の回帰

係数のみでは解釈できない ▶ たとえば,y = β0+ β1x1+ β2x2+ β3x1x2+ ϵ の回帰式を考えるx1の限界効果を計算するには回帰式を x1で偏微分すればよいので, ∂y ∂x1 = β1+ β3x2 (20) ▶ すなわち,x1の限界効果はβ1ではなく,β1+ β3x2 ▶ x2の限界効果も同様に求められる ∂y ∂x2 = β2+ β3x1 (21) ▶ x2の限界効果はβ2ではなく,β2+ β3x1 ▶ y = β0+ β1x1+ β2x21+ ϵ のような,二乗項 (“2” である必要はない.3 以上 でもよい) を含む回帰式についても同様 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 21 / 46

(24)

限界効果の計算:交互作用項がある場合

含意 (1) x1と x2の「掛け算」が入るのだから,x1を「だけ」を変化させ,「他の変 数を一定に」することはできない (2) x2= 0 (あるいは β3= 0) でない限り,x1の限界効果は (β1だけでなく) x2 の水準に依存する (2’) x2の水準を固定しなければx1の限界効果を計算できないが,x2の水準によっ てx1の限界効果が変化する ▶ 実際問題,「様々なx2の水準におけるx1の限界効果を示す図」がなければ解釈 できない! ▶ 典型的には,x2をx軸に,x2がある値のときのx1の限界効果をy軸にとっ たグラフを使う

(3) 一般化線形モデル (identity link 以外) の場合,“compression effect” が生じ

るため,さらに大変(「他の全ての独立変数」の値に依存する; Berry et al., 2010;

(25)
(26)

回帰係数の

t

検定

tk = ˆ βk− βk SE( ˆβk) は,自由度 n− k − 1 の t 分布に従うβkは真の回帰係数,βˆkはそのOLS推定値 ▶ 中心極限定理の応用:nが十分大きければ,標準正規分布近似が使える ▶ OLS推定量の漸近正規性(副読本・鹿野本第11章) ▶ OLS推定量の分散については,副読本・鹿野本第5, 11章 ▶ 検定統計量 tkと自由度 n− k − 1 の t 分布を使った t 検定が行なえる ▶ 帰無仮説H0:βk= 0 ▶ 対立仮説H1:βk̸= 0

(27)

回帰係数の

t

検定

平均値の t 検定と同様の手順で,βkの仮説検定ができる ▶ 帰無仮説H0:βk= 0 ▶ 対立仮説H1:βk̸= 0 (両側検定) ▶ 有意水準 α = 0.05 なら,平均値の検定でも出てきた「1.96SE」が目安H0は βk= 0 なので, ▶ tk= ˆ βk− βk SE( ˆβk) = βˆk− 0 SE( ˆβk) = βˆk SE( ˆβk) ▶ |tk| > 1.96 (あるいは約 2)なら,その回帰係数の効果は統計的に有意 (statistically significant) ▶ 「 ˆβk± 1.96SE( ˆβk) (あるいは約 2SE( ˆβk)) の範囲 = 95%信頼区間に 0 が含 まれなければ」でも同様 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 24 / 46

(28)

R

での出力例:先週配布の演習資料

(29)

R

コード:線形回帰

演習資料 ▶ R での線形回帰の実行方法と,OLS 推定量の漸近正規性を理解する ▶ 「演習資料」ページの「回帰分析のシミュレーション (1):simulating_ols.r」 ▶ URL: http://cfes-project.eco.u-toyama.ac.jp/wp-content/ uploads/simulating_ols.r 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 26 / 46

(30)

R

コード:線形回帰

架空データを生成するためのパラメータを設定する 1 ## Set sample size

2 sample_size = 10^3 3

4 ## Set regression parameters 5 beta_0 = 1 6 beta_1 = 2 7 x_mu = 2 8 x_sigma = 1 9 error_mu = 0 10 error_sigma = 1 架空データを発生させる 1 set.seed(123456)

2 simulated_x =rnorm(sample_size,mean = x_mu, sd= x_sigma)

3 simulated_y =beta_0 +beta_1*simulated_x + rnorm(sample_size,mean = error_mu, sd=

error_sigma)

(31)

R

コード:線形回帰

架空データを生成するためのパラメータを設定する 1 ## Set sample size

2 sample_size = 10^3 3

4 ## Set regression parameters 5 beta_0 = 1 6 beta_1 = 2 7 x_mu = 2 8 x_sigma = 1 9 error_mu = 0 10 error_sigma = 1 架空データを発生させる

1 simulated_x =rnorm(sample_size,mean = x_mu, sd= x_sigma)

2 simulated_y =beta_0 +beta_1*simulated_x + rnorm(sample_size,mean = error_mu, sd=

error_sigma)

3 simulated_dat =data_frame(simulated_y = simulated_y, simulated_x = simulated_x)

(32)

R

コード:線形回帰

架空データの内容を確認する

1 simulated_dat

架空データについて,OLS 推定を行なう 1 ## OLS estimate

2 simulated_ols =lm(simulated_y~simulated_x, data= simulated_dat)

(33)

R

コード:線形回帰

真の回帰係数 β1= 2 に近い推定結果が得られている?

1 >summary(simulated_ols) 2

3 Call:

4 lm(formula= simulated_y~simulated_x,data = simulated_dat) 5

6 Residuals:

7 Min 1QMedian 3Q Max

8 -2.8559 -0.7363 -0.0107 0.7402 3.6052

9

10 Coefficients:

11 Estimate Std. Errortvalue Pr(>|t|)

12 (Intercept) 0.98336 0.07287 13.50 <2e-16 ***

13 simulated_x 2.00431 0.03250 61.67 <2e-16***

14

---15 Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘ ’. 0.1 ‘ ’1

16

17 Residual standard error: 1.019on 998 degrees of freedom

18 Multiple R-squared: 0.7921, Adjusted R-squared: 0.7919

19 F-statistic: 3803 on 1 and 998 DF, p-value: < 2.2e-16

(34)

R

コード:

OLS

推定量の漸近分布

以上の架空データの生成と回帰分析を 10,000 回繰り返し,β1の OLS 推定量の分

布を確認する

1 set.seed(123456)

2 n_simulationz = 10^4

3 sim_betaz = numeric(n_simulationz)

4

5 ## Generate and estimate parameters

6 for(i in 1:n_simulationz) {

7 ## Print progress

8 cat(str_c(i,"␣out␣of␣", n_simulationz,"␣done.␣\r")) 9

10 ## Sample data

11 simulated_x =rnorm(sample_size,mean = x_mu, sd= x_sigma)

12 simulated_y =beta_0 +beta_1*simulated_x + rnorm(sample_size,mean = error_mu, sd=

error_sigma)

13 simulated_dat =data_frame(simulated_y = simulated_y, simulated_x = simulated_x)

14

15 ## OLS estimate

16 simulated_ols =lm(simulated_y~simulated_x, data= simulated_dat)

17

18 ## Store result

19 sim_betaz[i] = simulated_ols$coefficients[2]

(35)

R

コード:

OLS

推定量の漸近分布

1.90 1.95 2.00 2.05 2.10 0 2 4 6 8 10 12

Sample OLS estimates

D

en

si

ty

(36)
(37)

回帰分析の根源的仮定

線形回帰モデル yi = β0+ β1xi+ ϵi, i = 1, . . . , n を考える (yiも xiも確率 変数) ▶ OLS 推定量は, ˆ β0= ¯yn− ˆβ1x¯n ˆ β1= ∑n

i=1∑n(xi− ¯xn)(yi− ¯yn) i=1(xi− ¯xn)2 = sxy sxxβˆ0と ˆβ1が β0と β1の不偏推定量・一致推定量になるための必要条件は, 仮定1 外生性:E[ϵi|xi] = 0 仮定2 独立標本:異なる観察(xi, yi)と(xj, yj)は互いに独立 ▶ 教科書・浅野・矢内本第12章,副読本・鹿野本第10章,副読本・田中本第 5–6章も参照のこと 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 33 / 46

(38)

回帰分析の根源的仮定

仮定 1:外生性 ▶ E[ϵi|xi] = 0 を満たす独立変数 xiを,外生変数 (xiは外生的) という ▶ xiが外生変数ならば, ▶ E[yi|xi] = β0+ β1xi+ E[ϵi|xi] =0 = β0+ β1xi(線形回帰モデルの整合性) ▶ E[ϵi] = 0, E[xiϵi] = 0 (xiϵiの直交) ▶ Cov(xi, ϵi) = 0 (xiϵiの無相関) が成り立つ 仮定 2:独立標本 ▶ 標本が互いに独立ならば,xi以外の独立変数は ϵiについてなんら情報をも たない

ϵi= yi− β0− β1xiとすれば,E[ϵi|xi] = E[ϵi|x1, . . . , xn]

(39)

OLS

推定量の不偏性と一致性

強い外生性 (strong exogeneity)

E[ϵi|x1, . . . , xn] = 0, i = 1, . . . , n (22) E[ϵi|xi] = 0 (仮定 1) かつ E[ϵi|xi] = E[ϵi|x1, . . . , xn] (仮定 2)

OLS 推定量の不偏性 (unbiasedness) E[ ˆβ1] = β1 (23) OLS 推定量の期待値が母回帰係数になる OLS 推定量の一致性 (consistency) plim ˆβ1= β1+ Cov(xi, ϵi) Var(xi) = β1 (24) n→ ∞ で,OLS 推定量が母回帰係数に確率収束する ▶ 大数の法則と中心極限定理を思い出す ▶ 仮定 1 より,Cov(xi, ϵi) = 0 に注意 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 35 / 46

(40)

OLS

推定量の不偏性と一致性

出所:副読本・田中本 p.80 ▶ 小標本のデータ環境では,推定量に不偏性E[ˆθ] = θ が最低限求められる ▶ 大標本の (漸近的な) データ環境では,推定量に一致性plim ˆθ = θ が最低限 求められる ▶ 大数の法則と中心極限定理は,大標本理論(漸近理論)の二大定理 ▶ 実際の分析では,一致性しか実現できないことも多い

(41)
(42)

内生性バイアス:独立変数と誤差項の相関

OLS 推定量の一致性 (再掲) plim ˆβ1= β1+ Cov(xi, ϵi) Var(xi) = β1 ▶ Cov(xi, ϵi) = 0 (xiが外生変数) のとき,OLS 推定量は母回帰係数 β1の一 致推定量になる

Cov(xi, ϵi)̸= 0 (xiが内生変数) のとき,OLS 推定量は母回帰係数 β1の一

致推定量にならない! plim ˆβ1= β1+ Cov(xi, ϵi) Var(xi) = β1+bias( ˆβ1)̸= β1 (25) ▶ 外生変数(exogenous variable):誤差項と相関しない独立変数 ▶ 内生変数(endogenous variable):誤差項と相関する独立変数

(43)

内生性バイアス:独立変数と誤差項の相関

xiが内生変数のとき,OLS 推定量は一致性をもたず,バイアス (bias)を もつ ▶ バイアス(偏り):推定量(の期待値)と真のパラメータとの乖離,ズレ ▶ 仮定が満たされていれば,標本サイズが十分大きいとき,OLS推定量を真のパ ラメータとみなしてよい(推定量が真のパラメータに確率収束する) ▶ 推定量が一致性を持たない=標本サイズを大きくしても,「真の回帰係数」に 到達できない ▶ xiが内生変数のとき, ˆβ1は bias( ˆβ1) の分だけ,真の β から乖離した値を とってしまう ▶ (25)式のbias( ˆβ1) = Cov(xi, ϵi) Var(xi) を,内生性バイアスと呼ぶ ▶ bias( ˆβ1) > 0なら,過大評価 ▶ bias( ˆβ1) < 0なら,過小評価 ▶ 上のスライドで用いた R コード (架空データの回帰分析) では,xiが外生変 数だったので一致推定できていた 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 38 / 46

(44)

R

コード:内生性バイアス

演習資料 ▶ 内生性バイアスが生じる架空データをシミュレーションする ▶ 「演習資料」ページの「内生変数バイアスのシミュレーション: simulating_bias.r」 ▶ URL: http://cfes-project.eco.u-toyama.ac.jp/wp-content/ uploads/simulating_bias.r

(45)

R

コード:内生性バイアス

架空データを生成するためのパラメータを設定する (先ほどと同じ) 1 sample_size = 10^3 2 beta_0 = 1 3 beta_1 = 2 4 x_mu = 2 5 x_sigma = 1 6 error_mu = 0 7 error_sigma = 1 x と ϵ が相関する架空データを発生させる (何が変わったか注意) 1 ## Set bias 2 bias = 0.5 3 4 ## Sample data

5 simulated_error =rnorm(sample_size, mean= error_mu,sd = error_sigma)

6 simulated_x =rnorm(sample_size,mean = x_mu, sd= x_sigma)

7 simulated_x = (1 + bias*simulated_error)*simulated_x ## x is correlated with epsilon

8 simulated_y =beta_0 +beta_1*simulated_x + simulated_error

9 simulated_dat =data_frame(simulated_y = simulated_y, simulated_x = simulated_x)

x と ϵ が (想定通り) 相関しているか確認する

1 >cor(simulated_x, simulated_error)

2 [1] 0.6467905

(46)

R

コード:内生性バイアス

架空データを回帰分析する

1 simulated_ols =lm(simulated_y~simulated_x, data= simulated_dat)

回帰分析の結果を出力する.回帰係数を正しく推定できている? 1 >summary(simulated_ols)

2

3 Call:

4 lm(formula= simulated_y~simulated_x,data = simulated_dat) 5

6 Residuals:

7 Min 1QMedian 3Q Max

8 -2.74299 -0.53724 -0.09216 0.43467 2.93483

9

10 Coefficients:

11 Estimate Std. Errortvalue Pr(>|t|)

12 (Intercept) 0.14678 0.04016 3.655 0.000271 ***

13 simulated_x 2.43096 0.01609 151.123 < 2e-16 ***

14

---15 Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘ ’. 0.1 ‘ ’1

16

17 Residual standard error: 0.7569 on 998 degrees of freedom

18 Multiple R-squared: 0.9581, Adjusted R-squared: 0.9581

(47)

R

コード:内生性バイアス

母回帰係数,推定結果,バイアスを出力する 1 >beta_1## true beta

2 [1] 2

3 > simulated_ols$coefficients[2]## estimated beta

4 simulated_x

5 2.430964

6 >cov(simulated_x, simulated_error)/var(simulated_x) ## bias

7 [1] 0.4309637

8 >beta_1 +cov(simulated_x, simulated_error)/var(simulated_x)

9 [1] 2.430964

(25) 式 plim ˆβ1= β1+

Cov(xi, ϵi)

Var(xi) = β1+bias( ˆβ1)̸= β1を思い出す

(48)

R

コード:内生性バイアス

以上のシミュレーションを,先ほどと同じ要領で 1,000 回繰り返す

1 set.seed(123456)

2 n_simulationz = 10^3

3 sim_betaz = numeric(n_simulationz)

4

5 ## Generate and estimate parameters

6 for(i in 1:n_simulationz) {

7 ## Print progress

8 cat(str_c(i,"␣out␣of␣", n_simulationz,"␣done.␣\r")) 9

10 ## Sample data

11 simulated_error =rnorm(sample_size, mean= error_mu,sd = error_sigma)

12 simulated_x = bias*simulated_error*runif(sample_size)## x is correlated with

epsilon

13 simulated_y =beta_0 +beta_1*simulated_x + simulated_error

14 simulated_dat =data_frame(simulated_y = simulated_y, simulated_x = simulated_x)

15

16 ## OLS estimate

17 simulated_ols =lm(simulated_y~simulated_x, data= simulated_dat)

18

19 ## Store result

20 sim_betaz[i] = simulated_ols$coefficients[2]

(49)

R

コード:内生性バイアス

同じ要領でヒストグラムを描くと・・・ 2.0 2.1 2.2 2.3 2.4 2.5 0 5 10 15 20 beta estimates D en si ty 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 44 / 46

(50)

なぜ

Cov(x

i

, ϵ

i

)

̸= 0

が生じるのか?

代表的な原因 ▶ 除外 (省略) 変数 (omitted variable) ▶ 前回配布の演習資料はこのケースをシミュレーション ▶ 測定誤差 (measurement error) ▶ 同時性 (simultaneity)

(51)

次回講義と課題

▶ 次回講義 ▶ 今日の続き(回帰分析とバイアス) ▶ データの自動取得(API,ウェブスクレイピング) ▶ 文献と課題 必須 星野・田中『Rによる実証分析』第4–6章(教科書)

推奨 Gelman & Hill. Data analysis. Chaps. 3–4 (教科書)

推奨 浅野・矢内『Stataによる計量政治学』第7–12章(副読本) 推奨 石田ほか『Rによるスクレイピング入門』(副読本) 推奨 森田『実証分析入門』第4–8章(副読本) 課題 (1)講義資料「R言語の基礎,オブジェクトとその要素へのアクセス」と「R によるデータの読み込みと書き出し」をRStudioで練習しておくこと.(2)シ ミュレーションを再度行なってみること ▶ 特に,(例示したコードとは異なり)負の内生変数バイアスが生じるシミュレー ションを行なうこと(配布したコードを一部変更すればできる) 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 46 / 46

(52)

再掲:統計的仮説検定への導入

大きさ n の標本について,以下の検定統計量 t 値 (t value) は,自由度 v = n− 1 の t 分布 t(v) に従う(X, s2はそれぞれ標本平均と不偏分散) t =X− µ s2/n = X− µ SE (26) ▶ 中心極限定理で出てきた z = X−µ σ2/n ∼ N (0, 1) と同じく,t 値もX− µ に着している (分母が違う:z の σ2は厳密にいうと母分散)n→ ∞ で t 分布は標準正規分布 N (0, 1) に近付くので,n が大きければ標 準正規分布で近似できる (t が z に近付く)nや正規分布近似について回りくどい言い方をしていた理由の1つ(n→ ∞) ▶ 教科書が正規分布近似を使う際,「nが十分に大きな数であるならば」と断って いることに注意(e.g.,第5章) ▶ 検定統計量によって従う分布は異なるが(t分布とは限らない),以下で説明す る検定の基本的手続きは同じ(e.g., F分布, χ2分布)

(53)

再掲:統計的仮説検定

t 分布,t 値と p 値t 分布に従う確率変数について,cL=−tn−1,0.025 (下側臨界値) から cU = tn−1,0.025 (上側臨界値) の区間 [cL, cU] に,データの約 95%が収まる (α = 0.05 を考えているので,α/2 = 0.025)大きさ n の標本から得た t 値は,自由度 n− 1 の t 分布に従うcL, cUの値は,既出のRコードを使えば求められる ▶ nが十分大きければ,N (0, 1)と区間[−1.96, 1.96]を用いてもよい ▶ t 値による検定では α を事前に定めた上で,t > cU (t < cL) かを判定する 0 (X = µ0) cU cL α 2 α 2 棄却域 [−∞, cl] 受容域 [cl, cu] 棄却域 [cu,∞] 伊藤 岳 (Toyama/NIHU) R による計量分析 January 15, 2018 2 / 4

(54)

再掲:統計的仮説検定

t 分布,t 値と p 値大きさ n の標本から得た t 値は,自由度 n− 1 の t 分布に従うp 値 (p value):帰無仮説が正しいとき,標本から得た t 値以上に (絶対値で) 分布の中心からからかけ離れた値をとる確率 (両側検定の場合) ▶ t 値が 0 から離れるほど,p 値は小さくなる:「帰無仮説が正しいとすれば, 以上な t 値が出た」(起こり得ないことが起きた) と考えるt 値による検定とは異なり,p 値による検定では「異常性」を p 値によって 確率的に示すだけで棄却はしない (「棄却の妥当性」を示す) p 2 p 2

(55)

Berry, William D; Jacqueline H R DeMeritt & Justin Esarey (2010) Testing for Interaction Effects in Binary Logit and Probit Models: Is an Interaction Term Essential? American Journal of Political Science 54(1): 248–266.

Berry, William D; Matt Golder & Daniel Milton (2012) Improving Tests of Theories Positing Interaction.

Journal of Politics 74(3): 653–671.

Brambor, Thomas; William Clark & Matt Golder (2006) Understanding interaction models: Improving empirical analyses. Political Analysis 14(1): 63–82.

Clark, William R; Michael J Gilligan & Matt Golder (2006) A simple multivariate test for asymmetric hypotheses. Political Analysis 14(3): 311–331.

Hanmer, Michael J & Kerem Ozan Kalkan (2013) Behind the Curve: Clarifying the Best Approach to Calculating Predicted Probabilities and Marginal Effects from Limited Dependent Variable Models.

American Journal of Political Science 57(1): 263–277.

Kam, Cindy D & Robert J Franzese (2007) Modeling and Interpreting Interactive Hypotheses in Regression

Analysis: A Refresher and Some Practical Advice. Ann Arbor, MI: University of Michigan Press.

Rainey, Carlisle (2015) Compression and Conditional Effects: A Product Term Is Essential When Using Logistic Regression to Test for Interaction. Political Science Research and Methods forth.

Updating...

関連した話題 :