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

(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推定量_βˆ₁₌

∑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 n ∑ i=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 n ∑ i=1 e2i = 1 n n ∑ i=1 (yi− ˆβ0− ˆβ1xi)2 (6) ▶ これを最小化する回帰係数の値を，最小二乗推定量(OLS推定量)と呼ぶ ▶ MSEを最小化するには，回帰係数で偏微分して「= 0」と置いた連立方程式を解けばよい(最適化の1階条件) ∂MSE ∂ ˆβ0 =−2 n n ∑ i=1 (yi− ˆβ0− ˆβ1xi) = 0 (7) ∂MSE ∂ ˆβ1 =−2 n n ∑ i=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 推定量は ˆβ = (x⊤x)−1x⊤y と書ける ▶ 重回帰分析の場合はこうした行列・ベクトルを用いた表記が一般的 ▶ 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) のうち，何 %が予測値の変動 (回帰変動∑n_i=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) あるいは，残差変動∑n_i=1(yi− ˆyi)2が全変動に対してどれだけ小さいか ▶ R2= 1 なら回帰分析によって従属変数を完全に説明でき (残差変動が 0)， R2= 0 なら全く説明できない (残差変動が 1) ▶ ただし，R2_{は独立変数を追加すれば，それが何であれ (無意味なものでも)} 必ず増加する ▶ 独立変数の数 k の違いを調整するには，自由度調整済み決定係数 R2 aを用いる (副読本・田中本第 6 章) R2_a= 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 eﬀect)という ▶ 回帰係数と同じに思えるが，より高度な一般化線形モデル(本講義では扱わない)の理解や，後述の交互作用項を含むモデルの理解には限界効果が非常に重要

(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 eﬀect)

▶ 理論的議論：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 eﬀect” が生じ

るため，さらに大変(「他の全ての独立変数」の値に依存する; 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 架空データを発生させる

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

error_sigma)

14

15 ## OLS estimate

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)

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

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

1 >cor(simulated_x, simulated_error)

2 [1] 0.6467905

(46)

R

コード：内生性バイアス

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

回帰分析の結果を出力する．回帰係数を正しく推定できている？ 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

15

16 ## OLS estimate

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 Eﬀects 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 Eﬀects 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 Eﬀects: A Product Term Is Essential When Using Logistic Regression to Test for Interaction. Political Science Research and Methods forth.