Rでトービット・モデルによる打ち切りデータの推定

金額で測った不動産の保有量を考えましょう。大半の貧乏人はゼロでありますが、一定以上の所得がある人から保有量が増加してきます。集合住宅、持ち家、そして別荘や投資物件など増加していくわけです。
こういう説明変数の値に対して、一定まではゼロ（もしくは定数）、もしくは一定以上がゼロの推定モデルを切断モデルと言います。説明変数の全体の領域に対して線形にならないので、OLSで分析するとバイアスが入ります。バイアスをコントロールするためにTyle I〜IVまでのトービット・モデル（Tobit Model）が考案されています。一番、単純なタイプIモデルで、打ち切りデータの推定をしてみましょう。

1. Tobit Model

Tobit Modelでは、ある $x_i$ における、確率 $P(y_i|y_i=0)$ で打ち切る場合と、確率 $P(y_i|y_i>0)$ で打ち切らない場合の同時確率密度を最大化するように尤度Lが定めて最尤法をかけます（詳しくは「Applied Econometrics, 2009 - Tobit モデル」などを参照）。

$L = \Pi^{n}_{i=1} \log L_i$
$y_i=0$ のとき、 $L_i = P(y_i|y_i=0)$
$y_i>1$ のとき、 $L_i = {P(y_i|y_i>0)f(y_i|y_i>0)}$

$n$ は観測数、 $f(\cdot)$ は正規分布の確率密度関数です。導出は省略しますが、 $\mu$ を誤差項として、 $P(y_i|y_i=0)=P(y_i|yi\leq0)$ と $y_i=x_i \beta + \mu$ 、 $P(y_i|y_i>0)$ が $f(y_i|y_i>0)$ の分母になる事に気付いたら、尤度関数の導出は容易です。頑張って導出すると、対数尤度は以下のようになります。

$\log L = \sum^{n}_{i=1} \log L_i$
$y_i=0$ のとき、 $\log L_i = \log (1-\Phi(\frac{x_i \beta}{\sigma}))$
$y_i>1$ のとき、 $\log L_i = \log \frac{1}{\sigma}\phi(\frac{y_i - x_i\beta}{\sigma}))$

$\sigma$ は標準偏差です。場合分けで美しくは無いですが、コンピューター的には問題ないです。良く見ると正規化された正規分布が入っているわけで、ここは正規化しなくても問題はありません。

2. optim()関数で最尤推定

対数尤度関数をRの関数で書き直して、実際に推定をしてみましょう。気象庁のページから、金沢の気温と積雪量のデータを整理しておきました。yを被説明変数、Xを説明変数の行列、パラメーター・ベクトルはp[1]が $\sigma$ 、p[2]、p[3]、p[3]、･･･をXの係数としています。これは面倒な方法です。手を抜きたい人は第3節を見てください。

# データを読み込む
df <- read.table("kanazawa_accumulation_of_snow.txt", header=TRUE, sep="\t")
# 推定式を設定
frml <- AS ~ TEMP
# 初期値を決める為にOLSで回帰
r_ols <- lm(frml, data=df)
# 初期パラメーターを設定
p <- c(sqrt(deviance(r_ols)/df.residual(r_ols)), coefficients(r_ols))
# 説明変数の行列を整理
X <- model.matrix(frml, data=df)
# 被説明変数を整理
y <- df$AS
# 最適化アルゴリズムBFGSで最尤法を行う
optim(p, function(p){
  sigma <- p[1]
  beta <- matrix(p[-1], ncol=1)
  xb <- X%*%beta
  -1*sum(
    ifelse(y<=0,
      log(1-pnorm(xb/sigma)),
      log(dnorm((y - xb)/sigma)/sigma)
    )
# -1*sum(...) の代わりに -1*sum( (y<=0)*log(1-pnorm(xb/sigma))) + (y>0)*log(dnorm((y - xb)/sigma)/sigma) ) の方がベクター演算になって速いと突っ込まれたので追記
  )
}, method="BFGS", hessian=TRUE)

最尤法については「Rと手作業で覚える最尤法」を参照してください。

3. VGAMパッケージで最尤推定

対数尤度関数を書かないといけないTobit Modelの分析が面倒だって？ ─ そうですね。optim()関数ではヘッセ行列を見ないと係数の標準偏差も分かりませんし、不便です。便利に扱うには、VGAMパッケージを使いましょう。

library(VGAM)
r_tobit <- vglm(frml, tobit(Lower=0), trace=TRUE, data=df)
summary(r_tobit)

4. Tobit vs OLS

Tobit ModelとOLSで推定したときの違いについて簡単に説明します。厳密には逆ミルズ非などの意味を調べる必要がありますが、OLSだと大量にあるゼロの値に引きずられてバイアスがかかる場合でも、Tobitでは大きくは影響されません。係数が大きく変わります。下図は2節・3節の推定結果をプロットしたものですが、バイアスが補正されています。

Tobit

	係数	標準誤差	z値
切片項	135.4718	12.84249	10.549
気温	-17.7652	2.25276	-7.886

OLS

	係数	標準誤差	z値
切片項	45.4029	4.3352	10.473
気温	-2.1489	0.2542	-8.453

係数の符号は変化がありませんが、係数の大きさや有意性が変化しているのがわかります。

5. おわりに

トービット・モデルは所得と保有資産の関係を見るときに良く使われ、OLS以外の計量手法を使わなければいけない理由を明確に図示できる分析モデルです。見た目的に説明してみようと思ったのですが、公開データで適合度が高いものを探すのに苦労しました。個票データなどで分析すると重宝するのかも知れませんが、それ以外では適応できそうなデータセットは多くは無いようです。