（この記事のもともとの場所はこちらのウェブサイトになります）

従属変数が，すべてのケースについては観察されない場合がある．これがランダム（Missing Completely At Random, MCAR）に生じていればよいが，たいていはそうではなく，もし系統的な偏りがある場合には，何らかの補正が必要になってくる．その際に用いられるのが，サンプルセレクションモデルである．

サンプルセレクションモデルで頻出する概念が，逆ミルズ比（inverse mills ratio）である．この逆ミルズ比というのが，なかなか直感的に分かりにくい値で，何度も誤解してしまっていたので，なるべく分かりやすく考えるため，逆ミルズ比を図示しながら考えてみる．

サンプルセレクションモデルは，従属変数が観察されるかどうかを考慮したモデルを作るという点で，打ち切り・切断データの分析手法であるトービットモデルと極めて似た特徴を持っている．なので，打ち切り・切断データの分析手法であるType-1 Tobit modelを事例に逆ミルズ比について考えてみよう．

Type-1 Tobit model

以下ではType-1 Tobit modelについて述べる．詳細については Maddala(1983), Breen(1996), Long(1997)などを参照されたい．

\(Y_i\)の潜在変数\(Y_i^*\)を考える．ここで回帰式\(Y_i^* = X_i\beta + u_i, \ u_i\sim N(0,\sigma^2)\)はOLSの仮定を満たすものとする．\(Y_i\)，閾値\(c\)を超える場合には観察でき，\(c\)を下回る場合にはすべてcとなっている（センサーされている）ものとする．このとき\(Y_i\)は以下のように表せる．

\[ Y_i = \left\{ \begin{array}{ll} Y_i^* = X_i \beta + u_i & \mathrm{if} \ \ Y_i^* > c \\ 0 & \mathrm{if} \ \ Y_i^* \leq c \end{array} \right. \]

たとえば，\(Y^*\sim N(100,100), \ c = 0\)のときは以下のような感じになる．\(Y^*\)が0のところで打ち切られ，それ以下の値はすべて0に貼り付いているようなかたち．

x <- seq(-200,400,length=600)
plot(x, dnorm(x,mean=100,sd=100), type="l", xlab="Y*", ylab="", yaxt="n", main="normal distribution and truncation")
abline(v = 0)
xvals <- seq(-200, 0, length=200)
dvals <- dnorm(xvals, mean=100, sd=100)
polygon(c(xvals,rev(xvals)), c(rep(0,200), rev(dvals)), col="gray")

\(X_i\)を条件づけたときの\(Y_i\)の期待値は以下のように表せる． \[ E(Y_i | X_i) = \Pr(Y_i^* > c | X_i)E(Y_i | Y_i^* > c, X_i) + [1 - \Pr(Y_i > c | X_i)]\times c \] 以降は，簡単のために\(c = 0\)とする．すると上式は以下のようになる． \[ E(Y_i | X_i) = \Pr(Y_i^* > 0 | X_i)E(Y_i | Y_i^* > 0, X_i) \] さていま，\(Y^*_i > 0\)のときの\(Y_i\)の期待値への回帰分析を行う．\(Y_i\)そのものは正規分布していないから，\(Y_i\)の値でそのまま回帰分析をした場合，係数\(\beta\)はBLUE（Best Linear Unbiased Estimator, 最良線形不偏推定量）とならない．そこで，閾値の超えやすさを考慮して推定する，というのがサンプルセレクションモデルの基本的な発想である．Type-1 Tobitの場合は以下のように書ける． \[ \begin{aligned} E(Y_i | Y_i^* > 0, X_i) &= E(Y_i | X_i\beta + u_i > 0) \\ &= E(Y_i | u_i > -X_i\beta) \\ &= X_i\beta + E(u_i | u_i > -X_i\beta) \\ &= X_i\beta + \sigma E(u_i / \sigma | u_i / \sigma > -X_i\beta / \sigma) \\ &= X_i\beta + \sigma \frac{\phi(X_i\beta / \sigma)}{\Phi(X_i\beta / \sigma)} \end{aligned} \] \(\phi(\cdot)\)は密度関数，\(\Phi(\cdot)\)は分布関数である．ここで，\(\lambda(X_i\beta / \sigma) = \phi(X_i\beta / \sigma) / \Phi(X_i\beta / \sigma)\)とおく．この\(\lambda(\cdot)\)のことを逆ミルズ比（inverse mills ratio）という．

逆ミルズ比

逆ミルズ比はどこから出てきたのだろうか．先ほどの式変形を確認しておこう．3つめの式変形について，\(E(u_i | u_i > -X_i\beta)\)とは，\(-X_i\beta\)を超えるときの\(u_i\)の期待値を意味している．4つめの式変形では，\(u_i/\sigma\)を\(u_i\)を標準偏差で割って標準化するという操作をしている． \(u_i / \sigma\)を\(-X_i\beta /\sigma\)で打ち切ったときの分布（切断正規分布）は以下のような形状になる．

#install.packages("truncnorm") #切断正規分布を作るためのパッケージ"truncnorm"をインストール．
library(truncnorm)
x <- seq(-3,3,length=1000)
plot(x, dnorm(x, mean=0, sd=1), type="l", ylim=c(0,0.5), xlab="u/σ", ylab="", main="truncated normal distribution")
lines(x, dtruncnorm(x, a=-0.8, mean=0, sd=1), type="l", col="red")
legend("topright", 
       legend=c("standard normal dist","truncated"),
       lty=c(1,1),
       col=c("black","red"),
       bty="n"
       )

標準正規分布を\(-X_i\beta /\sigma\)で打ち切ったとき，切断正規分布の条件つき期待値\(E(u_i|u_i > -X_i\beta\)は，打ち切りが起こるところ（\(-X_i\beta/\sigma\)）までの面積\(\Phi(X_i\beta / \sigma)\)で密度関数\(u_i/\sigma\)を割った値となる．（ここからのあたりがやや分かりにくく，注意すべきポイント）

この切断正規分布から，\(-X_i\beta\)を大きくするほど（\(X_i\beta\)を小さくするほど）打ち切り部分に含まれる\(\phi(u_i/\sigma)\)の面積\(\Phi(u_i/\sigma)\)は広くなっていくわけなので，打ち切りが起こりやすくなっているといえる．またその結果，条件つき期待値\(E(u_i | u_i > -X_i\beta)\)は大きくなっていくことがわかる．このときは，もともとの密度関数\(\phi(u_i/\sigma)\)からかなりずれた状態になっている．

逆に，\(-X_i\beta\)を小さくするほど（\(X_i\beta\)を大きくするほど）打ち切り部分に含まれる\(\phi(u_i/\sigma)\)の面積\(\Phi(u_i/\sigma)\)は狭くなっていくわけなので，打ち切りが起こりにくくなっているといえる．またその結果，条件つき期待値\(E(u_i | u_i > -X_i\beta)\)は小さくなっていくことがわかる．このときは，もともとの密度関数\(\phi(u_i/\sigma)\)に近い状態である．

plot(x, dnorm(x, mean=0, sd=1), type="l", ylim=c(0,0.6), xlab="u/σ", ylab="", main="relations between -Xβ and truncated normal distribution")
lines(x, dtruncnorm(x, a=-0.4, mean=0, sd=1), type="l", col="red")
lines(x, dtruncnorm(x, a=-1.5, mean=0, sd=1), type="l", lty=2, col="red")
legend("topright", 
       legend=c("standard normal dist","large -Xβ","small -Xβ"),
       lty=c(1,1,2),
       col=c("black","red","red"),
       bty="n"
       )

ここまでくれば，あとは\(X_i\beta\)の値を変化させたときに，密度関数を分布関数で除した逆ミルズ比\(\lambda(X_i\beta / \sigma)\)（すなわち，条件つき期待値\(E(u_i | u_i > -X_i\beta)\)）がどのような挙動を示すかを見ればよい．

x <- seq(-3,3,length=1000)
plot(x, dnorm(x,mean=0,sd=1), type="l", ylim=c(0,3), xlab="Xβ", ylab="", main="change of mills ratio when Xβ are changed")
lines(x, pnorm(x,mean=0,sd=1), type="l",  lty=2)
lines(x, dnorm(x,mean=0,sd=1)/pnorm(x,mean=0,sd=1), type="l", col="red")
legend("topright", 
       legend=c("pdf","cdf","invmills"),
       lty=c(1,2,1),
       col=c("black","black","red"),
       bty="n"
       )

一見してわかるとおり，\(\lambda(X_i\beta / \sigma)\)は，\(X_i\beta\)が正の方向に大きくなるほど0に近づき，負の方向に大きくなるほど0から離れていく．

先に述べたとおり，\(X_i\beta\)が大きくなるとはすなわち，打ち切りが起こりにくくなることを意味していた．ということは，逆ミルズ比\(\lambda(X_i\beta / \sigma)\)は，打ち切りの起こりやすさを表していると読むことができる．

逆ミルズ比を統制することの意味

ここまでの関係をまとめると以下のようになる．関心となる従属変数の背後に正規分布する潜在変数があり，ある閾値（今回は0とした）を超えた場合はその値が観察できるが，閾値を下回る場合はすべてその閾値でしか観察できないという状況がある．このとき，観察できる従属変数に対してそのまま回帰分析を行うと，係数にバイアスが生じる．そこで用いるのがトービットモデルであった．このモデルは，以下のように表せる． \[ \begin{aligned} E(Y_i | Y_i^* > 0, X_i) &= X_i\beta + \sigma E(u_i| u_i > -X_i\beta) \\ &= X_i\beta + \sigma E(u_i / \sigma | u_i / \sigma > -X_i\beta / \sigma) \\ &= X_i\beta + \sigma \frac{\phi(X_i\beta / \sigma)}{\Phi(X_i\beta / \sigma)} \end{aligned} \] このモデルにおいて独立変数\(X_i\)は，従属変数の値を変化させるし，従属変数が閾値を超えて観察される確率（打ち切り確率）にも影響をおよぼす変数となる．その関係は以下のようになっている．

\(X_i\beta\)が大きくなる
⇔逆ミルズ比\(\lambda(X_i\beta / \sigma)\)が0に近づく
⇔打ち切りが起こりにくくなり，その結果条件つき期待値\(E(u_i | u_i > - X_i \beta)\)が小さくなる

\(X_i\beta\)が小さくなる
⇔逆ミルズ比\(\lambda(X_i\beta / \sigma)\)が0から遠ざかる
⇔打ち切りが起こりやすくなり，その結果条件つき期待値\(E(u_i | u_i > - X_i \beta)\)が大きくなる

以上より，逆ミルズ比を統制するということはすなわち，打ち切りの起こりやすさを統制することを意味しているといえる．

Sample selection model

Type-1 Tobitでは，従属変数に影響する変数と，従属変数が観察されるかどうかに影響する変数が同一であると考えた．しかし，これらの変数は別々である可能性も当然考えられる．この場合に用いられるのが，Type-2 Tobit，あるいはHeckit（Heckman’s selection model）と呼ばれる方法である． 従属変数\(Y_i\)は，\(Z_i = 1\)のときのみ観察でき，\(Z_i = 0\)のときは観察できないとする． \[ Y_i = \left\{ \begin{array}{ll} Y_i^* = X_i \beta + u_i & \mathrm{if} \ \ Z_i = 1 \\ \cdot & \mathrm{if} \ \ Z_i = 0 \end{array} \right. \] また\(Z_i\)は以下のように表せる． \[ Z_i = \left\{ \begin{array}{ll} 1 & \mathrm{if} \ \ W_i\gamma + e_i > 0\\ 0 & \mathrm{if} \ \ W_i\gamma + e_i \leq 0 \end{array} \right. \] ただし\(u_i \sim N(0,\sigma_u^2), \ e_i \sim N(0,1), \mathrm{cov}(u_i, e_i) = \sigma_{ue}\)とする．\(\ e_i \sim N(0,1)\)としているのは，\(Z_i\)に関するモデルをプロビットモデルにしたがって推定するためで，ロジットモデルでもまったく問題ない．

このとき，\(Z = 1\)のときの\(Y_i\)の期待値への回帰分析は以下のように表せる． \[ \begin{aligned} E(Y_i | Z_i = 1, X_i) &= E(Y^*_i | W_i\gamma + e_i > 0) \\ &= E(Y^*_i | e_i > -W_i\gamma) \\ &= X_i\beta + E(u_i | e_i > -W_i\gamma) \\ &= X_i\beta + \sigma_{ue} E(e_i | e_i > -W_i\gamma) \\ &= X_i\beta + \sigma_{ue} \frac{\phi(W_i\gamma)}{\Phi(W_i\gamma)} \end{aligned} \] ここまでくれば，以降の考え方は先のType-1 Tobitと同じようになる．逆ミルズ比\(\lambda(W_i\gamma)/\Phi(W_i\gamma)\)はやはり，値が高いほど従属変数の欠損が起こりやすいということを示す．先ほどと異なるのは，従属変数が観察されるかどうかに関するモデルと，従属変数の水準に関するモデルとが別々にある，という点だけである．

参考文献

Breen, Richard. 1996. Regression Models Censored, Sample-Selected, or Truncated Data. Thousand Oaks: Sage.

Long, J.Scott. 1997. Regression Models for Categorical and Limited Dependent Variables. Thousand Oaks: Sage.

Maddala, G. S. 1983. Limited-Dependent and Qualitative Variables in Econometrics. Cambridge: Cambridge University Press.