プロビット回帰

ある種の生物 m 個体を生まれてから n 年間追跡調査したとき、各年齢において死亡率がどのように変化したかをモデルする例を示す。

i 歳の個体数を mi とし、そのうち i + 1 歳になれなかった個体数を死亡数として yi とする。i 歳のときの死亡率を pi としたとき、実際の死亡数が観測される確率分布は以下のように計算できる。

\[ P(Y = y_{i}) = \begin{pmatrix}m_{i} \\ y_{i} \end{pmatrix}p_{i}^{y_{i}}(1-p_{i})^{m_{i}-y_{i}} \]

そこで、全年齢 i = 1, 2, ..., n について考えた時、その同時確率は以下のように書ける。

\[ P(Y_{1}=y_{i}, \cdots, Y_{n}=y_{n}) = \prod_{i=1}^{n}\begin{pmatrix}m_{i} \\ y_{i} \end{pmatrix}p_{i}^{y_{i}}(1-p_{i})^{m_{i}-y_{i}} \]

次に、各年齢における死亡数を死亡率 πi としてとらえた時、E(Yi)/ mi = πi により、モデルは以下のように構築できる。

\[ E[\mathbf{Y}] = \mathbf{\pi}\] \[g(\mathbf{\pi}) = \mathbf{X}\mathbf{\beta} \]

2 つ目の式について、抗生物質の濃度が変化することによって 1 よりも大きくなると考えられる。そこで、リンク関数 g の逆関数は、その右辺を 0 から 1 の範囲に収める関数でなければならない。そのような関数にはロジット関数やプロビット関数がある。ロジット関数を利用すると以下のように書ける。

\[ g(\mathbf{\pi}) = \Phi^{-1}(\mathbf{\pi}) = \mathbf{X}\mathbf{\beta} = \begin{pmatrix}1 & x\end{pmatrix}\begin{pmatrix}\beta_{1} \\ \beta_{2}\end{pmatrix} \]

Φ-1 は標準正規分布の累積分布関数の逆関数を表す。

ここまでで、誤差構造は二項分布、リンク関数はロジット関数と利用することを確認した。

x <- 1:14
m <- c(1000, 1000, 999, 997, 996, 993, 899, 895, 885, 841, 739, 343, 65, 2)
y <- c(   0,    1,   2,   1,   3,   4,   4,  10,  44, 102, 396, 278, 63, 2)

これらのデータを解析しやすいように整理して、xx を年齢、yy をその年齢に対してその個体が生きていれば 0 を、死亡した場合 1 とする。

xx <- 0
yy <- 0
for (i in 1:length(y)) {
  xx <- c(xx, rep(i, length = m[i]))
  yy <- c(yy, rep(0, length = m[i] - y[i]), rep(1, length = y[i]))
}
xx <- xx[-1]
yy <- yy[-1]

glm 関数を利用してモデルを作成する。

probit.model <- glm(yy ~ xx, family = binomial("probit"))

モデルのパラメーターを書くにする。

summary(probit.model)
## Call:
## glm(formula = yy ~ xx, family = binomial("probit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1201  -0.1867  -0.0060   0.0000   6.5829  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -7.49460    0.21714  -34.52   <2e-16 ***
## xx           0.67265    0.02101   32.02   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6217.6  on 10653  degrees of freedom
## Residual deviance: 2883.1  on 10652  degrees of freedom
## AIC: 2887.1
## 
## Number of Fisher Scoring iterations: 8

Coefficients 項を見ると β1 = -7.49460、β2 = 0.67265 にそれぞれ対応していることが確認できる。

モデルの式から、π と年齢の間は以下の関係が成り立つ。

\[ \pi = \Phi (\beta_{1} + \beta_{2}x) \]

この式を利用すれば実際のデータの散布図上にモデル化された曲線を描くことができる。

plot(x, y / m, xlim = c(0, 14), ylim = c(0, 1),
     xlab = "age", ylab = "mortality")
par(new = T)
curve(pnorm(-7.49460 + 0.67265 * x),
      xlim = c(0, 14), ylim = c(0, 1), xlab = "", ylab = "")