Dobson本 chap. 4

4.1 はじめに

• 本章では、最尤法を用いたGLMのパラメータの点推定と区間推定の方法を説明する。
  
• 多くの場合、数値的な解法が必要となるが、それらの方法は反復的で、ニュートン・ラプソンアルゴリズムに基づくことが多い。

4.2 例；圧力釜の故障時間

• 以下のデータは、圧力レベル70%の圧力釜の故障までの時間である

atsuryoku <- c(1051, 1337, 1389, 1921, 1942, 2322, 3629, 4006, 4012, 4063, 4921, 5445, 5620, 5817,
               5905, 5956, 6068, 6121, 6473, 7501, 7886, 8108, 8546, 8666, 8831, 9106, 9711, 9806,
               10205, 10396, 10861, 11026, 11214, 11362, 11604, 11608, 11745, 11762, 11895, 12044,
               13520, 13670, 14110, 14496, 15395, 16179, 17092, 17568, 17568)

• このデータの分布は以下の通りである。

hist(atsuryoku)

---

対数尤度関数とスコア関数の導出

• 故障までの時間（生存時間）の分布として、よく利用される分の１つにワイブル分布があり、確率密度関数は以下の通りである。 \[f(y; \lambda, \theta) = \frac{\lambda y^{\lambda-1}}{\theta^{\lambda}}exp[-(\frac{y}{\theta})^\lambda] \tag{4.2.1}\] ここで\(y>0\)は故障までの時間、\(\lambda\)は分布の形状を決めるパラメータ、\(\theta\)は尺度パラメータである。
• 式(4.1.1)を変形すると、指数型分布族であることがよくわかる。 \[f(y; \theta) = exp[log\lambda+(\lambda-1)logy-\lambda log\theta-(y/\theta)^{\lambda}] \tag{4.2.2}\]
• \(Y_1, ..., Y_N\)は\(N=49\)のデータであり、\(Y_i\)は独立な確率変数で、全て同じパラメータのワイブル分布に従うとすると、 \[f(y_1, ..., y_N; \theta, \lambda) = \prod^N_{i=1}\frac{\lambda y_i^{\lambda-1}}{\theta^{\lambda}}exp\biggl[-\biggl(\frac{y_i}{\theta}\biggr)^{\lambda}\biggr] \tag{4.2.3}\]
• 式(4.1.3)の対数尤度関数は、 \[l = logf(\theta; y_1,...,y_N, \lambda) = \sum^N_{i=1}\Biggl[\biggl[(\lambda-1)logy_i+log\lambda-\lambda log\theta \biggr]-\biggl(\frac{y_i}{\theta}\biggr)^2\Biggr] \tag{4.2.4}\]
• 式(4.1.4)を\(\theta\)について微分すると、以下のスコア関数が得られる。 \[\frac{dl}{d\theta} = U = \sum^N_{i=1}\biggl[\frac{-\lambda}{\theta}+\frac{\lambda y_i^{\lambda}}{\theta^{\lambda+1}}\biggr] \tag{4.2.5}\]
• 最尤推定値\(\hat{\theta}\)は、式(4.1.5)のスコア関数\(U(\theta) = 0\)の解となる。
  \(\lambda\)を既知の定数とすれば、\(\hat{\theta}\)の明示的な解は簡単に見つかる。

---

ニュートン・ラプソン近似によるパラメータ推定

• ここでは、ニュートン・ラプソン近似を使って数値解を求める。

• 関数\(t\)がx軸と交わる\(t(x) = 0\)となる\(x\)を見つけたい。
  値\(x^{(m-1)}\)での\(t\)の傾きは、距離\(x^{(m)}-x^{(m-1)}\)が小さければ、 \[\biggl[\frac{dt}{dx}\biggr]_{x = x^{(m-1)}} = t'(x^{(m-1)}) = \frac{t(x^{(m)})-t(x^{(m-1)})}{x^{(m)}-x^{(m-1)}}\] もし、\(x^{(m)}\)が\(t(x)=0\)を満たすとすれば、 \[x^{(m)} = x^{(m-1)}-\frac{t(x^{(m-1)})}{t'(x^{(m-1)]})} \tag{4.2.6}\] であり、式(4.2.6)はニュートン・ラプソンの公式である。

• 初期推定値\(x^{(1)}\)から始めて、式(4.2.6)を用いて、逐次堅持して、収束するまで反復する。
  
• 最尤推定の場合、スコア関数を使ってニュートン・ラプソンの公式を表すと、 \[\theta^{(m)} = \theta^{(m-1)}-\frac{U^{(m-1)}}{U'^{(m-1)}} \tag{4.2.7}\]
• \(\lambda = 2\)のワイブル分布の場合、 \[U = -\frac{2×N}{\theta}+\frac{2×\sum{y_i^2}}{\theta^3}\]
• \(U\)の導関数は以下の通りである。
  \[\begin{align} \frac{dU}{d\theta} &= U' = \sum_{i=1}^N \biggl[\frac{\lambda}{\theta^2}-\frac{\lambda(\lambda+1)y_i^{\lambda}}{\theta^{\lambda+2}}\biggr] \\ &= \frac{2×N}{\theta^2}-\frac{2×3×\sum{y_i^2}}{\theta^4} \tag{4.2.8} \end{align}\]
• 最尤推定では、\(E(U')\)によって\(U'\)を近似する。
  指数型分布族では、\(E(U')\)について、 \[\begin{align} E(U') &= b''(\theta) E[a(Y)] + c''(\theta) \\ &=b''(\theta)\biggl[-\frac{c'(\theta)}{b'(\theta)}\biggr]+c''(\theta) \\ &= -var(U') = -\Im \end{align}\] が成立し、情報量\(\Im\)は、 \[\begin{align} \Im &= E(-U') = E\biggl[-\sum^N_{i=1}U_i\biggr] = \sum^N_{i=1}\biggl[E(-U_i')\biggr] \\ &= \sum^N_{i=1}\biggl[\frac{b''(\theta)c'(\theta)}{b'(\theta)}-c''(\theta)\biggr] = \frac{\lambda^2N}{\theta^2} \tag{4.2.9} \end{align}\]
• 式(4.2.7)と(4.2.9)から、 \[\theta^{(m)}=\theta^{(m-1)}+\frac{U^{(m-1)}}{\Im^{(m-1)}} \tag{4.2.10}\] であり、式(4.2.10)による方法をスコア法という。

• \(\hat{\theta}\)の分散は\(\Im = E(-U')\)に反比例し、\(\hat{\theta}\)の標準誤差は、近似的に \[SE(\hat{\theta}) = \sqrt{1/\Im}\]

---

4.3 最尤推定

• GLMの性質を満たす確率変数\(Y_1, ..., Y_N\)について考える。
• \(E(Y_i) = \mu_i\)および\(g(\mu_i) = \boldsymbol{X}^T_i\boldsymbol{\beta}\)で\(\boldsymbol{\beta}\)を推定する。
  
• 各\(y_i\)に対する対数尤度関数
  \[l_i = y_ib(\theta_i)+c(\theta_i)+d(y_i) \tag{4.3.1}\]
• ここで、以下を定義しておく。 \[E(Y_i) = \mu_i = -c'(\theta_i)/b'(\theta_i) \tag{4.3.2}\] \[var(Y_i) = [b''(\theta_i)c'(\theta_i)-c''(\theta_i)b'(\theta_i)]/[b'(\theta_i)]^3 \tag{4.3.3}\] \[g(Y_i) = \boldsymbol{X}^T_i\boldsymbol{\beta} = \eta_i \tag{4.3.4}\]
• 全ての\(y_i\)に対する対数尤度関数
  \[l = \sum^N_{i = 1}l_i = \sum y_ib(\theta_i) + \sum c(\theta_i) + \sum d(y_i) \tag{4.3.5}\]
• パラメ－タ\(\beta_j\)の最尤推定量を得るために、微分の連鎖法則に基づく以下の式を利用する。 \[\frac{\delta l}{\delta \beta_j} = U_j = \sum^N_{i = 1} \biggl[\frac{\delta li}{\delta\beta_j}\biggr] = \sum^N_{i = 1}\biggl[\frac{\delta l_i}{\delta \theta_i}\frac{\delta \theta_i}{\delta \mu_i}\frac{\delta \mu_i}{\delta \beta_j}\biggr] \tag{4.3.6}\]
• 最右辺を各項に分けて考える
  \[\begin{align} \frac{\delta l_i}{\delta \theta_i} = y_ib'(\theta_i)+c'(\theta_i) &= b'(\theta_i)\Biggl(y_i + \frac{c'(\theta_i)}{b'(\theta_i)}\Biggr) \\ &= b'(\theta_i)(y_i-\mu_i) \end{align} \tag{4.3.7}\] \[\begin{align} \frac{\delta \theta_i}{\delta \mu_i} &= 1/\biggl(\frac{\delta \mu_i}{\delta \theta_i}\biggr) \\ \frac{\delta \mu_i}{\delta \theta_i} &= \frac{-c''(\theta_i)}{b'(\theta_i)} + \frac{c'(\theta_i)b''(\theta_i)}{[b'(\theta_i)]^2} \\ &= \frac{b''(\theta_i)c'(\theta_i) - c''(\theta_i)b'(\theta_i)}{[b'(\theta_i)]^2}×\frac{b'(\theta_i)}{b'(\theta_i)} \\ &= b'(\theta_i)var(Y_i) \tag{4.3.8} \end{align}\] \[\frac{\delta \mu_i}{\delta \beta_j} = \frac{\delta \mu_i}{\delta \eta_i} \frac{\delta \eta_i}{\delta \beta_j} = \frac{\delta \mu_i}{\delta \eta_i}x_{ij} \tag{4.3.9}\]
• 式(4.3.7), (4.3.8), (4.3.9)を式(4.3.6)に代入すると、 \[\begin{align} U_j &= \sum^N_{i = 1} \biggl[b'(\theta_i)(y_i-\mu_i)\frac{1}{b'(\theta_i)var(Y_i)}x_{ij}\frac{\delta \mu_i}{\delta \eta_i} \\ &= \sum^N_{i = 1}\biggl[\frac{y_i - \mu_i}{var(Y_i)} x_{ij} \frac{\delta \mu_i}{\delta \eta_i} \biggr] \tag{4.3.10} \end{align}\] となり、スコアという。
  
• 情報行列\(\boldsymbol{\Im}\)は、\(U_j, j = 1, ..., p\)の分散共分散行列に等しいが、(j, k)要素を\(\Im_{jk} = E[U_jU_k]\)
  
• \(Y_i\)が独立で\(E[(Y_i-\mu_i)(Y_l - \mu_l)] = 0, i \neq l\)が成立することを考慮すると、
  式(4.3.10)より、 \[\begin{align} \Im_{jk} &= E \Biggl[\sum^N_{i = 1} \biggl[\frac{(Y_i - \mu_i)}{var(Y_i)}x_{ij}\frac{\delta \mu_i}{\delta \eta_i}\biggr]\sum^N_{l=1}\biggl[\frac{(Y_l-\mu_l)}{var(Y_l)}x_{ik}\frac{\delta \mu_l}{\delta \eta_l}\biggr]\Biggl] \\ &= \sum^N_{i = 1}\frac{E[(Y_i-\mu_i)]x_{ij}x_{ik}}{[var(Y_i)]^2}\biggl(\frac{\delta \mu_i}{\delta \eta_i}\biggr)^2 \\ &= \sum^N_{i=1}\frac{x_{ij}x_{ik}}{var(Y_i)}\biggl( \frac{\delta \mu_i}{\delta \eta_i} \biggr)^2 \tag{4.3.11} \end{align}\]
• 単一のパラメータの場合、スコア法の推定方程式は、パラメータベクトル\(\boldsymbol{\beta}\)に対して、 \[\boldsymbol{b}^{(m)} = \boldsymbol{b}^{(m-1)}+[\boldsymbol{\Im}^{(m-1)}]^{-1}\boldsymbol{U}^{(m-1)} \tag{4.3.12}\] のように一般化できる。
  • \(\boldsymbol{b}^{(m)}\)：m回目の反復時におけるパラメータ\(\beta_1, ..., \beta_p\)の推定値ベクトル
  • \([\boldsymbol{\Im}^{(m-1)}]^{-1}\)：式(4.3.12)の要素\(\Im_{jk}\)をもつ情報行列の逆行列
  • \(\boldsymbol{U}^{(m-1)}\)：式(4.3.11)の要素をもつベクトル
• 式(4.3.12)の両辺に\(\boldsymbol{\Im}^{(m-1)}\)をかける
  \[\boldsymbol{\Im}^{(m-1)}\boldsymbol{b}^{(m)} = \boldsymbol{\Im}^{(m-1)}\boldsymbol{b}^{(m-1)}+\boldsymbol{U}^{(m-1)} \tag{4.3.13}\]
• 式(4.3.11)を行列で表現すると \[\boldsymbol{\Im} = \boldsymbol{X}^T\boldsymbol{WX} \tag{4.3.14}\]
  • \(\boldsymbol{W}\)：\(w_{ii} = \frac{1}{var(Y_i)}\biggl(\frac{\delta \mu_i}{\delta \eta_i}\biggr)^2\)を対角要素としてもつ、N × N対角行列
• 式(4.3.11)および(4.3.12)から、式(4.3.13)の右辺は、第j要素として以下をもつベクトルである。 \[\sum^P_{k=1}\sum^N_{i=1}\frac{x_{ij}x_{ik}}{var(Y_i)}\biggl(\frac{\delta \mu_i}{\delta \eta_i}\biggr)^2b_k^{(m-1)}+\sum^N_{i=1}\biggl[\frac{y_i-\mu_i}{var(Y_i)}x_{ij}\frac{\delta \mu_i}{\delta \eta_i}\biggr]\]
• 式(4.3.13)の右辺を行列で表現すると、\(\boldsymbol{X}^T\boldsymbol{Wz}\)である
  • \(\boldsymbol{z}\)：\(z_i = \sum^P_{k = 1}x_{ik}b_k^{(m-1)}+(y_i-\mu_i)\biggl(\frac{\delta \eta_i}{\delta \mu_i}\biggl)\)
• 以上のことから、反復計算式(14.3.13)は、 \[\boldsymbol{X}^T\boldsymbol{WXb}^{(m)} = \boldsymbol{X}^T\boldsymbol{Wz} \tag{4.3.15}\]
• 式(4.3.15)は\(\boldsymbol{z}\)と\(\boldsymbol{W}\)が\(\boldsymbol{b}\)に依存するため、反復的に解かないといけないことを除くと、線形モデルに対して重み付き最小二乗法を適用して得られる正規方程式と同じである。
  GLMに対する最尤推定は、反復重み付き最小二乗法の手順で得られる(Charnes et al. 1976)。

---

4.4 ポアソン反応変数に対する回帰分析の例　　

再記しておきたい式

\[w_{ii} = \frac{1}{var(Y_i)}\biggl(\frac{\delta \mu_i}{\delta \eta_i}\biggr)^2 \tag{4.4.1}\] \[z_i = \sum^P_{k=1}x_{ik}b_k^{(m-1)}+(y_i-\mu_i)\biggl(\frac{\delta \eta_i}{\delta \mu_i}\biggr) \tag{4.4.2}\] \[\boldsymbol{X}^T\boldsymbol{WXb}^{(m)} = \boldsymbol{X}^T\boldsymbol{Wz} \tag{4.4.3}\]

• 応答変数\(y_i\)がポアソン分布に従う人工データ

data <- tibble(y = c(2, 3, 6, 7, 8, 9, 10, 12, 15),
               x = c(-1, -1, 0, 0, 0, 0, 1, 1, 1))
ggplot(data, aes(x = x, y = y)) + geom_point()

• 図を見ればわかるように、xの値が大きくなるに連れて、yの平均値が大きくなり、分散も大きくなっている。
  この観察は、ポアソン分布の以下の性質を満たす。 \[E(Y_i) = var(Y_i) \tag{4.4.4}\]

• \(Y_i\)と\(x_i\)に関係を以下の直線でモデル化してみる。 \[\begin{align} E(Y_i) &= \mu_i = \beta_1 + \beta_2 x_i \\ &= \boldsymbol{x}_i^T\boldsymbol{\beta} \end{align}\]

• ここで、\(i = 1, ..., N\)で \[\boldsymbol{\beta} = \begin{bmatrix} \beta_1 \\ \beta_2 \end{bmatrix}, \;\;\; \boldsymbol{x}_i = \begin{bmatrix} 1 \\ x_i \end{bmatrix}\] であり、リンク関数は\(g(\mu_i)\)は恒等関数とする(identity)。 \[g(\mu_i) = \mu_i = \boldsymbol{x}_i^T\boldsymbol{\beta} = \eta_i\] \(\frac{\delta \mu_i}{\delta \eta_i}=1\)となり、式(4.4.1)と(4.4.2)は簡単な形になる。
  式(4.4.1)とポアソン分布の性質から、 \[w_{ii} = \frac{1}{var(Y_i)} = \frac{1}{\beta_1 + \beta_2x_i}\] となり、\(\beta\)に対する推定値\(\boldsymbol{b} = \begin{bmatrix} b_1 \\ b_2 \end{bmatrix}\)を用いると、式(4.4.2)は、 \[z_i = b_1 + b_2x_i + (y_i - b_1 - b_2x_i) = y_i\] となるので、式(4.4.3)の係数と右辺は、 \[\boldsymbol{\Im} = \boldsymbol{X}^T\boldsymbol{WX} = \begin{bmatrix} \sum^N_{i=1}\frac{1}{b_1 + b_2x_i} & \sum_{i=1}^N\frac{x_i}{b_1+b_2x_i} \\ \sum^N_{i=1}\frac{x_i}{b_1+b_2x_i} & \sum^N_{i=1}\frac{x_i^2}{b_1+b_2x_i} \end{bmatrix}\] \[\boldsymbol{X}^T\boldsymbol{Wz} = \begin{bmatrix} \sum_{i=1}^N\frac{y_i}{b_1+b_2x_i} \\ \sum_{i=1}^N \frac{x_iy_i}{b_1+b_2x_i} \end{bmatrix}\] と表され、最尤推定料は、以下の式を反復計算することによって得られる。 \[(\boldsymbol{X}^T\boldsymbol{WX})^{(m-1)}b^{(m)} = (\boldsymbol{X}^T\boldsymbol{Wz})^{(m-1)}\]

---

実データの計算

• dataのような\(N = 9\)のデータの場合、 \[\boldsymbol{y} = \boldsymbol{z} = \begin{bmatrix} 2 \\ 3 \\.\\.\\.\\15 \end{bmatrix} \;\;\; , \;\;\; \boldsymbol{X} = \begin{bmatrix} \boldsymbol{x}^t_1 \\ \boldsymbol{x}^T_2 \\ .\\.\\.\\ \boldsymbol{x}_9^T \end{bmatrix}=\begin{bmatrix} 1&-1\\1&-1\\.&.\\.&.\\.&.\\1&1\end{bmatrix}\] である。
• 初期推定値を\(b_1^{(1)}=7\)と\(b_2^{(1)} = 5\)とする。これらの値を用いて計算すると、 \[(\boldsymbol{X}^T\boldsymbol{WX})^{(1)} = \begin{bmatrix}1.821429&-0.75\\-0.75&1.25 \end{bmatrix}, \;\;\;\;\; (\boldsymbol{X}^T\boldsymbol{Wz})^{(1)} = \begin{bmatrix} 9.869048\\0.583333\end{bmatrix}\] となり、したがって \[\begin{align} \boldsymbol{b}^{(2)} &= \biggl[ (\boldsymbol{X}^T\boldsymbol{WX}^{(1)}) \biggr]^{-1}(\boldsymbol{X}^T\boldsymbol{Wz})^{(1)} \\ &= \begin{bmatrix} 0.729167&0.4375\\ 0.4375 & 1.0625\end{bmatrix}\begin{bmatrix} 9.869048\\0.583333\end{bmatrix} \\ &= \begin{bmatrix} 7.4514\\4.9375 \end{bmatrix} \end{align}\] を得る。この反復プロセスを収束するまで続ける。
• 得られた最尤推定値は、\(\hat{\beta} = 7.45163\)と\(\hat{\beta}=4.93530\)であり、これらの値における情報行列\(\boldsymbol{\Im} = \boldsymbol{X}^T\boldsymbol{WX}\)の逆行列は、 \[\boldsymbol{\Im}^{-1} = \begin{bmatrix} 0.7817&0.4166\\0.4166&1.1863\end{bmatrix}\] となる。これは\(\hat{\boldsymbol{\beta}}\)の分散共分散行列である。 \(\beta_2\)の95%CIは \[4.9353±1.96\sqrt{1.1863} = (2.80, 7.07)\]

---

Rによる行列計算

• Dobson本の人工データdataを使って、スコア法による最尤推定を手計算風にRでやってみる。

# 関数を定義しておく
XWX <- function(x, b1, b2){
  element_11 <- sum(1/(b1+(b2*x)))
  element_12 <- sum(x/(b1+(b2*x)))
  element_21 <- sum(x/(b1+(b2*x)))
  element_22 <- sum(x^2/(b1+(b2*x)))
  xwx <- matrix(c(element_11, element_12, element_21, element_22), ncol = 2, byrow = TRUE)
  return(xwx)
}

XWz <- function(x, y, b1, b2){
  upper <- sum(y/(b1+(b2*x)))
  lower <- sum((y*x)/(b1+(b2*x)))
  xwz <- matrix(c(upper, lower), ncol = 1)
  return(xwz)
}

b_new <- function(mat1, mat2){
  b_new = solve(mat1) %*% mat2
  return(b_new)
}

# b_newを収束するまで反復計算
#  初期値を与える
b1 <- 7
b2 <- 5
N <- 1000 #　適当な数字（反復計算の最大回数）

# 反復計算
for(i in 1:N){
  mat1 <- XWX(data$x, b1 = b1, b2 = b2)
  mat2 <- XWz(data$x, data$y, b1 = b1, b2 = b2)
  bnew <- b_new(mat1, mat2)
  if(abs(b1-bnew[1,]) <= 0.00001 & abs(b2-bnew[2,]) <= 0.00001) break
  b1 <- bnew[1, ]
  b2 <- bnew[2, ]
}

# 結果を見てみる
bnew

##          [,1]
## [1,] 7.451633
## [2,] 4.935300