要約
次のようなサンプル・セレクションがある場合で
ヘックマンの2段階推定はサンプル・セレクション・バイアスを除去できました。
ヘキット分析
（1）ヘキット・モデル
ヘキット・モデルはサンプル・セレクション・バイアスを考慮するモデルです。例えば、母集団では
$$Y=\beta_0+\beta_1X+u$$
でも、潜在変数S*
$$S^*= \delta_0+\delta_1X +\delta_2 Z+v$$
が0以上でないと観測されないモデルです。
ただし
$$E(u|X)=E(v|X)=0$$
$$v 〜 N(0,1)$$
$$E(u|v)=\gamma v$$
が仮定されています。
（2）ヘックマンの2段階推定
ヘキット分析では、ヘックマンの2段階推定か最尤法を用います。この記事では、ヘックマンの2段階推定を解説します。ヘックマンの2段階推定は、回帰分析を2回する方法です。
1回目はセレクション式
$$S^*= \delta_0+\delta_1X +\delta_2 Z+v$$
をプロビット分析し、δ0、δ1、δ2を推定します。推定値には「^」をつけましょう。次に、
$$\widehat{c}= \widehat{\delta_0}+\widehat{\delta_1}X+\widehat{\delta_2} Z$$
として、それぞれの標本について
$$逆ミルズ比＝E(v|v ≧\widehat{c})$$
を推定します。
最後に、構造式に逆ミルズ比の推定値を入れて
$$Y_{観測有}=\beta_0+\beta_1X+\gamma \widehat{逆ミルズ比}+\epsilon$$
について最小二乗法を行います。これがヘックマンの2段階推定です。
最小二乗法とヘックマンの2段階推定
ヘキット分析では、どのようなサンプル・セレクション・バイアスに対応できるのかを、シミュレーションで確かめてみます。
（1）母集団と標本
母集団では、構造式
$$Y=1.5+0.3X+$$
でも、潜在変数S*
$$S^*= -0.1-0.2X +0.1 Z+v$$
が0以上でないと観測されないとします（セレクション式）。そして、
$$E(u_i|v_i)=-5 v_i$$
とします。つまり、逆ミルズ比の係数は-5です。
他にもいくつかの項目を設定して、出力すると以下の結果になりました。左が母集団に従う分布で、右がYが観測されなかったサンプルを除いた標本の分布です。
（2）最小二乗推定
右図の標本を最小二乗法推定した結果が、下図です。
目的変数Xの母回帰係数が「0.3」であるのに、最小二乗推定では「-0.21」となりました。プラスがマイナスとされたのですから、深刻なサンプル・セレクション・バイアスが発生しています。しかも、1%有意でありタチが悪いです。
（3）ヘックマンの2段階推定
次にヘックマンの2段階推定によって、推定します。その結果が下です。
目的変数Xの母回帰係数が「0.3」である中、ヘックマンの2段階推定では「0.18」となりました。10%有意にすぎず、誤差は存在しますが、最小二乗法よりは改善されています。
サンプル・セレクション・バイアスの影響
（1）分布を調べる意義
さきほどの結果のみでヘックマンの2段階推定が優れていると結論することはできません。偶然、母回帰係数に近かっただけかもしれないからです。
そこで、同じモデルで乱数発生させ推定する試行を、1000回繰り返します。それぞれで得た推定値の分布を比較すれば、平均的にどちらの推定法がよいのかがわかります。
（2）推定量の分布
1000回推定した分布が下図です。
上図をみると、最小二乗法ではサンプル・セレクション・バイアスを受けており、ヘックマンの2段階推定ではサンプル・セレクション・バイアスを除去できていることがわかります。
ヘキット分析のRコード
#パッケージのインストールと呼び出し
install.packages("sampleSelection") #PCにつき1回やればよい
library(sampleSelection)
#ヘックマンの2段階推定
lm0 <- heckit(s ~ x + z, #セレクション式
y ~ x, #構造式
method = "2step", #2段階推定
data=?) #データ「？」を使う
summary(lm00)
#最尤法
lm00 <- heckit(s ~ x + z, #セレクション式
y ~ x, #構造式
method = "mls", #最尤法
data=?) #データ「？」を使う
summary(lm00)
#（おまけ）最小二乗法（←ヘキット分析ではない）
lm000 <-lm(YS ~ x,data=？)
summary(lm000)シミュレーションの方法
（1）セッティング
シミュレーションでは次のモデルを考えていました。
- n＝1500
- 乱数発生
$$説明変数：X 〜 N(-0.5,5)$$
$$Zの誤差項：e 〜 N(0,2)$$
$$構造式の誤差項の誤差項：ee 〜 N(0,0.5)$$
$$セレクション式の誤差項：uu 〜 N(0,1)$$
- セレクション式の変数作成
$$除外変数Z：z=-0.05-0.1X+e$$
$$プロビットの潜在変数S^*：ss=-0.1 – 0.2x + 0.1z+ uu$$
$$観測有無：ss≧0ならs=1、ss＜0ならs=0$$
- 構造式の変数作成
$$構造式の誤差項：u=-5uu+ee$$
$$目的変数：Y=1.5+0.3X+u$$
- サンプル・セレクション
$$観測された目的変数（※観測なしは0）：$$
$$ys = s (1.5 + 0.3x + u)$$
$$観測された目的変数（※観測なしはNA）：$$
$$ys ≠0ならYS =ys、 ys ＝0ならYS =（NA）$$
（2）シミュレーションのRコード
#サンプル数
sample_number <- 1500
#乱数発生
x <- rnorm(n=sample_number, mean =-0.5, sd=5)
e <- rnorm(n=sample_number, mean =0, sd=2)
ee <- rnorm(n=sample_number, mean =0, sd=0.5)
uu <- rnorm(n=sample_number, mean =0, sd=1)
#セレクション式
z <- -0.05 - 0.1*x + e
ss <- -0.1 - 0.2*x + 0.1 *z+ uu
s <- ifelse(ss>=0,1,0)
#構造式
u <- -5*uu + ee
y <- 1.5 + 0.3*x + u
#サンプル・セレクション
ys <- s*(1.5 + 0.3*x + u)
YS <- ifelse(ys==0,NA,ys)（3）散布図の作成
#母集団
plot(x,y)
#標本
plot(x,YS)（4）1000回の試行
1000回の試行は、処理に数分かかります。気長に待ちましょう。
#サンプル数sample_number
sample_number <- 1500
#試行回数trial_number
trial_number <- 1000
#データなし（NA）がtrial_number個並んだ箱をつくり、beta_1と名付けます
beta_1_heck_2step <- rep(NA, trial_number)
beta_1_heck_mls <- rep(NA, trial_number)
beta_1_ols <- rep(NA, trial_number)
#1回目からtrial_number回目まで繰り返し処理します
for(i in 1:trial_number){
#サンプル数
sample_number <- 1500
#乱数発生
x <- rnorm(n=sample_number, mean =-0.5, sd=5)
e <- rnorm(n=sample_number, mean =0, sd=2)
ee <- rnorm(n=sample_number, mean =0, sd=0.5)
uu <- rnorm(n=sample_number, mean =0, sd=1)
#セレクション式
z <- -0.05 - 0.1*x + e
ss <- -0.1 - 0.2*x + 0.1 *z+ uu
s <- ifelse(ss>=0,1,0)
#構造式
u <- -5*uu + ee
y <- 1.5 + 0.3*x + u
#サンプル・セレクション
ys <- s*(1.5 + 0.3*x + u)
YS <- ifelse(ys==0,NA,ys)
#ヘックマンの2段階推定
lm0 <- heckit(s ~ x + z,
YS ~ x,
method = "2step")
beta_1_heck_2step[i] <- lm0$coef[5] #coefがポイント
#最尤法
lm00 <- heckit(s ~ x + z,
YS ~ x,
method = "mls")
beta_1_heck_mls[i] <- lm00$estimate[5] #estimateがポイント
#最小二乗推定量
lm000 <-lm(YS ~ x)
beta_1_ols[i] <- lm000$coef[2]
}
#最小二乗推定量の分布
hist(beta_1_ols,xlim=c(-0.5,0.8),breaks=30)
#ヘックマンの2段階推定の分布
hist(beta_1_heck_2step,xlim=c(-0.5,0.8),breaks=30)
#最尤推定量の分布
hist(beta_1_heck_mls,xlim=c(-0.5,0.8),breaks=30)
カテゴリー