『統計学実践ワークブック』第19章の勉強メモです。
トービットモデル
分析の対象となる被説明変数のデータ範囲が制限されている場合、制限従属変数（LDV）と呼ぶ。これらは3つの場合に分けられる。
- 切断：被説明変数および説明変数がともにかけているデータ
- 打ち切り：被説明変数はかけているが、説明変数はかけていないデータ
- セレクションバイアス：データを集める段階で、特定のサンプルについてデータを完全に集められていないがために欠損が生じるデータ
このようなデータに対して分析を行った際に、データ範囲外の推定結果が出たり、誤差項が正規分布でないため偏りのある推定結果が出てしまう。そうした分布のゆがみを補正して推定するのがトービットモデルとなる。
トービットモデルには大きく２種類あり、切断データ・打ち切りデータに対してはタイプIトービットモデルを、セレクションバイアスに対してはタイプIIトービットモデルを適用する。
タイプIトービットモデル
切断回帰モデル
$y_i$を実際の観測値、$y_i^*$を潜在変数とする。潜在変数$y^*_i$、説明変数$x_i$に打ち切りが伴う。
\begin{array}{}
y_i = y_i^* & (y_i^* > 0)
\end{array}
切断正規分布の密度関数で尤度関数を構成し、最尤法によってパラメータを求める。
打ち切り回帰モデル
$y_i$を実際の観測値、$y_i^*$を潜在変数とする。潜在変数$y^*_i$に打ち切りが伴うが、説明変数$x_i$は常に観測可能である。
観測値$y_i$と潜在変数$y_i^*$の関係は以下の3つのケースに分けられる。ただし、$L,U$は定数。
-
左打ち切りが伴う場合
y_i= \begin{cases} y_i^* & (L < y_i^*) \\ L & (y_i^* \leq L) \end{cases} -
右打ち切りが伴う場合
y_i = \begin{cases} y_i^* & (y_i^* < U) \\ U & (U \leq y_i^*) \end{cases} -
両側打ち切りが伴う場合
y_i= \begin{cases} L & (y_i^* \leq L) \\ y_i^* & (L < y_i^* < U) \\ U & (U \leq y_i^*) \end{cases}
$y$の密度と、$y$が$0$になる確率の積を尤度関数とし、最尤法によってパラメータを求める。
尤度関数
打ち切り回帰モデルを考える。ただし、$\beta$はパラメータからなるベクトルで、$\epsilon_i$は互いに独立に正規分布$N(0,\sigma^2)$に従うとする。
y_i^* = x_i^\top \beta + \epsilon_i, i = 1,\cdots,n
左打ち切りの場合の尤度関数について考える。
次のように定義関数を定義する。
I(y_i^*)=
\begin{cases}
0 & (y_i \leq L) \\
1 & (L < y_i)
\end{cases}
$\psi(・)$、$\Phi(・)$をそれぞれ標準正規分布の確率密度関数、累積分布関数とする。
$L<y_i$における$y_i$の確率密度関数は、切断された部分の確率を除いて基準化するので、
f(y_i|y_i>L)=\frac{1}{\sigma}\psi(\frac{y_i-x_i^\top \beta}{\sigma})/[1-\Phi(\frac{L-x_i^\top \beta}{\sigma})]
となる。よって$L<y_i$となる確率$f(L<y_i)$は
P(L<y_i)=f(L<y_i)f(y_i|y_i>L)=\frac{1}{\sigma}\psi(\frac{y_i-x_i^\top \beta}{\sigma})
となる。
打ち切られる確率$P(y_i= L)$は
P(y_i= L)=P(y_i^*\leq L)=\Phi(\frac{L-x_i^\top\beta}{\sigma})
となる。
$\beta,\sigma$に関する尤度関数は次のようになる。
\begin{multline}
\begin{split}
L(\beta,\sigma) &= \prod^n_{i=1}P(y_i)\\
&= \prod^n_{i=1}P(L<y_i)^{I(y_i)}P(y_i=L)^{1-I(y_i)}\\
&= \prod^n_{i=1}\left[\frac{1}{\sigma}\psi(\frac{y_i-x_i^\top \beta}{\sigma})\right]^{I(y_i)}\left[\Phi(\frac{L-x_i^\top\beta}{\sigma})\right]^{1-I(y_i)}\\
&= \prod_{i:y_i^*>L}\frac{1}{\sigma}\psi(\frac{y_i-x_i^\top \beta}{\sigma})\prod_{i:y_i^*\leq L}\Phi(\frac{L-x_i^\top\beta}{\sigma})
\end{split}
\end{multline}
タイプIIトービットモデル
$y_i$を実際の観測値とし、属性が2つある場合、属性１の観測値を$y_{1i}$、属性２の観測値を$y_{2i}$とする。説明変数$x_i$に打ち切りが伴い、属性１の説明変数を$x_{1i}$、属性２の説明変数を$x_{2i}$。
\begin{multline}
\begin{split}
y_{1i} &= x_{1i}^\top\beta_1+\epsilon_1 \\
y_{2i} &= x_{2i}^\top\beta_2+\epsilon_2 \\
y_i &= T_iy_{1i}+(1-T_i)y_{2i} \\
T_i &= I(Z_i\gamma+\epsilon_{0i})
\end{split}
\end{multline}
推定方法は、Heckmanの２段推定（Heckit）と呼ばれる手順を用いることが多い。
（第1段階）全サンプルを使って労働供給に関するProbit分析を行いγの推定値を求め、γ
の推定値から逆ミルズ比を求める。
（第2段階）労働供給を行った（z＊＞0）サンプルを使って、賃金（y）を説明変数（x）と逆
ミルズ比にOLSで回帰する。
生存時間解析
ある時点から興味のあるイベントが観察されるまでの時間$T(\geq 0)$を生存時間という。
生存時間解析の主要な目的は生存関数の推定である。
S(t)=P(T > t)=\int^\infty_t f(x|\theta)dx
ただし、$f(x|\theta)$は$T$の確率密度関数で、$\theta$は未知のパラメータである。生存関数$S(t)$はイベントが発生するまでの時間が$t$を超える確率を表している。
生存時間解析のモデル構築は生存関数モデルの構築であるが、ハザード関数からアプローチするのが一般的である。
ハザード関数の例
比例ハザードモデル
予後因子などの複数の共変量$x$を考慮する必要があるが、一般化線形モデルにおいては、連結関数$g(・)$を用いて、$E(T|x)=g(x^\top\beta)$の仮定を置き、反復重み付き最小二乗法をもちいて$\beta$の推定を行う。生存時間解析において、$E(T|x)$の代わりにハザード関数を用いて、
h(t;x)=\exp(\alpha + x^\top\beta)=h_0e^{x^\top\beta},(h_0=e^\alpha >0)
を仮定する。このモデルにおけるハザードは時間に依存せず、生存時間$T$の分布は指数分布に限られる。確率密度関数$f(t)=\lambda\exp(-\lambda t)$をもつ指数分布のハザードは$\lambda$なので、$\lambda=h_0e^{x^\top\beta}$として最尤法で$\alpha,\beta$で推定できる。
指数分布を仮定した回帰モデルにおける$h_0$を$h_0(t)$で置き換えて得られるのがCox比例ハザードモデルである。
h(t;x) = h_0(t)e^{x^\top\beta},(h_0(t)>0)
$h_0(t)=h(t;x=0)$なので、$h_0(t)$は基準ハザードと呼ばれる。$h_0(t)$は$t$の関数であるが、その形を全く指定しないため、比例ハザードモデルはセミパラメトリックモデルの一種である。生存時間$T$の分布がどんな分布であっても、比例ハザードモデルに基づく解析はロバストな結果を保証する。
比例ハザードモデルの下で
\frac{h(t;x)}{h(t;x^*)}=\exp[(x-x^*)^\top\beta]
となるので、ハザードの比は時間に依存しないことがわかる。この性質を比例ハザード性という。
比例ハザードモデルを適用するときに、比例ハザード性のチェックが必要である。比例ハザードモデルの下で
\log(-\log{S(t;x)})=x^\top\beta+\log{H_0(t)}
と導ける。ただし、$H_0(t)=\int^t_0 h_0(u)du$は基準累積ハザードである。これを利用して比例ハザード性を検証できる。
カプラン・マイヤー推定量
パラメトリックモデルの適用が難しい場合、生存関数の推定量を次のように構成できる。大きさ$n$の無作為標本$t_1,\cdots,t_n$に対して、打ち切りがない場合を考える。$S(t)=1-F(t)$により、経験分布$F_n(t)$を用いて
\hat{S(t)}=1-F(t)=\frac{1}{n}\sum^n_{i=1}I(t_i>t)
で推定できる。ここで、$I(・)$は定義関数。
推定量$\hat{S(t)}$はカプラン・マイヤー推定量の特殊な場合である。
時刻$t_j$におけるカプラン・マイヤー推定量は、直前の時刻$t_{j-1}$における推定量と$t_j$まで生きていた条件のもとで$t_j$を乗り越える条件付確率の推定量との積で表すことができる。すなわち、$\hat{S(t_j)}=\hat{t_{j-1}}\times\hat{P}[T>t_j|T\geq t_j]$。
これより、カプラン・マイヤー推定量は
\hat{S(t_j)}=\prod^j_{i=1}\hat{P}[T>t_i|T\geq t_i]
とできる。
ニューラルネットワークモデル
説明変数ベクトル$x=(x_1,\cdots,x_p)^\top$を入力、被説明変数ベクトル$y=(y_1,\cdots,y_p)^\top$を出力とする。入力層と出力層のみからなるニューラルネットワークモデルでは$y$の第$i$要素$y_i$が、$x$の要素の線形結合$\omega_i^\top x$を活性化関数と呼ばれる関数$f$で変換して得られるものとする。
y_i = f(\omega_i^\top x)
活性化関数としてはロジスティック関数$f(x)=\frac{1}{1+e^{-x}}$やReLU関数$f(x)=\max\{0,x\}$が用いられる。シグモイド関数とも呼ばれる。$\omega_i$を重みベクトルとよぶ。また、$\omega_i^\top x$に定数項$\omega_{i0}$を付け加えて$\omega_{i0}+\omega_i^\top x$とすることもあり、$\omega_{i0}$をバイアス項と呼ぶ。以上の関係をまとめて$y=f_\theta(x)$と表す。ここで、$\theta$は重みベクトルやバイアス項からなるパラメータベクトルである。
$f_\theta$の形の関数を何回か合成し、合成の途中の段階を中間層、隠れ層とよぶ。特に中間層の数の多いニューラルネットワークモデルを深層ニューラルネットワークという。
例題
問１９．１
[1]
\begin{multline}
\begin{split}
L(\beta,\sigma) &= \prod_{i;y_i>L}\frac{1}{\sigma}\psi(\frac{y_i-x_i^\top\beta}{\sigma})\prod_{i;y_i\leq L}\Phi(\frac{L-x_i^\top\beta}{\sigma}) \\
&= \prod_{i;y_i>0}\frac{1}{\sigma}\psi(\frac{y_i-\beta_0-\beta_1x_{i1}-\beta_2x_{i2}}{\sigma})\prod_{i;y_i\leq 0}\Phi(\frac{-\beta_0-\beta_1x_{i1}-\beta_2x_{i2}}{\sigma})
\end{split}
\end{multline}
[2]
パラメータ数を$k$とすると
AIC = -2\times{最大対数尤度}+2k
となる。AICはモデルと真の分布の乖離を計るもので、小さい方が良い。
したがって、「日照時間＋平均気温」のモデルが最適といえる。
問１９．２
生存関数$S(t)$は
S(t)=P(T>t)=\int^\infty_t\lambda e^{-\lambda t}dt=e^{-\lambda t}
ハザード関数は
h(t)=(-\log{S(t)})'=-\frac{S(t)'}{S(t)}=\frac{\lambda e^{-\lambda t}}{e^{-\lambda t}}=\lambda
問１９．３
$h(t;x)=(-\log{S(t;x)})'$より、$S(t;x)=\exp(-\int^t_0 h(u;x)du)$となる。
よって、
\begin{multline}
\begin{split}
S(t;x) &= \exp(-\int^t_0 h_0(u)\exp(x^\top\beta)du)\\
&= \exp(-\exp(x^\top\beta)\int^t_0h_0(u)du) \\
&= \exp(-\exp(x^\top\beta)H_0(t))
\end{split}
\end{multline}
両辺に$\log(-\log{・})$をとると
\log(-\log{S(t;x)}) = x^\top\beta + \log{H_0(t)}
問１９．４
[1]
比例ハザード性$\log{(-\log{S(t;x)})}=x^\top \beta+\log{H_0(t)}$において、$H_0(t)$は対照群となる。すなわち、$H_0(t)=\int^t_0(-\log{S(u;x=0)})'du=-\log{S(t;0)}$となる。
治療群$x=1$について、$\log{(-\log{S(t;x=1)})}=\beta+\log{(-\log{S(t;0)})}$となる。
これより、治療群と対照群は平行になることが期待され、グラフは大方平行に見える。よって、比例ハザードモデルを適用するのは妥当である。
[2]
[1]より、距離は$\beta$にあたり、
\begin{multline}
\begin{split}
\beta &= \log{(-\log{S(t;x=1)})}-\log{(-\log{S(t;x=0)})}=\log{\frac{\log{S(t;x=1)})}{\log{S(t;x=0)})}}\\
&= \log{\frac{h(t;x=1)}{h(t;x=0)}}
\end{split}
\end{multline}
となる。よって、距離は治療効果の大きさを表す