ゴール
- Modified poisson regressionについての理解を深めること
やること
- GLMとその解の漸近性質についての理解を深める
- 分散のsandwich estimatorについての理解を深める
- 数式を整理した後、Rでの実装を行う（→長くなりすぎたので別の記事に回す）
参考文献
- Agresti本：Agresti, A. (2002) Categorical Data Analysis. 2nd Edition, John Wiley & Sons, Inc., New York, 320-332. http://dx.doi.org/10.1002/0471249688
- Dobson本：Dobson, A.J., & Barnett, A.G. (2018). An Introduction to Generalized Linear Models (4th ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781315182780
- Geyer. (2013) 5601 Notes: The Sandwich Estimator. https://www.stat.umn.edu/geyer/5601/notes/sand.pdf
- Freedman. On The So-Called “Huber Sandwich Estimator” and “Robust Standard Errors”. https://www.stat.berkeley.edu/~census/mlesan.pdf
- Zou G. A modified poisson regression approach to prospective studies with binary data. Am J Epidemiol. 2004 Apr 1;159(7):702-6. doi: 10.1093/aje/kwh090. PMID: 15033648.
- Fuyama K, Hagiwara Y, Matsuyama Y. A simulation study of regression approaches for estimating risk ratios in the presence of multiple confounders. Emerg Themes Epidemiol. 2021 Dec 11;18(1):18. doi: 10.1186/s12982-021-00107-2. PMID: 34895270; PMCID: PMC8665581.
Exponential family
下式の様に表現できる確率分布をexponential familyという。
p(y; \theta, \phi) = \exp \left(\frac{y\theta - b(\theta)}{a(\phi)} + c(y, \phi)\right)
\thetaが回帰係数に関係するところなので対数尤度を想定してそこだけを取ると、下記となる。
l = \frac{y\theta - b(\theta)}{a(\phi)}
ポアソン分布で考える
例えばポアソン分布を考える場合、 a(\phi) = 1, b(\theta) = e^\theta, c(y, \phi) = -\log(y!)とすると、
p(y; \theta) = \exp \left(y\theta - e^\theta - \log(y!)\right)
ここで、\theta = \log(\mu) (\muはポアソン分布に従う変数Yの平均値)とすると、
p(y; \theta, \phi) = \exp \left(y\log(\mu) - \mu - \log(y!)\right)
= \frac{\mu^y e^{-\mu}}{y!}
\theta = \log(\mu)はポアソン回帰の一般的なリンク関数(\eta = \log(\mu) = X\beta)である。このように\theta = \etaとなるようなg(\cdot)をcanonicalリンク関数という。二項分布の場合b(\theta) = log[1 + exp(\theta)]で後述の通り\mu = b'(\theta)なので\mu = \frac{exp(\theta)}{1+exp(\theta)}となる。よって、\eta = logit(\mu) = \thetaとなるため、ロジット関数は二項分布においてcanonicalである。
注1：exponential familyの記法にはいくつか種類がある。また、canonicalリンク関数の定義についても表現が異なる場合がある。今回はAgrestiを参考にした。
注2：a(\cdot)は正直よくわかっていないが、\phiがdispersionパラメータと呼ばれるものである。a(\cdot)はポアソン回帰の場合は1、二項分布の場合1/n、正規分布の場合は\phiそのもので分散に相当する。
重要な関係式の整理
微分と積分の交換が可能な元で、ある個体iの対数尤度の1階微分の期待値は、E[\frac{\partial l_i}{\partial \theta_i}] = 0より
E[y_i] = \frac{\partial b(\theta_i)}{\partial \theta_i} = \mu_i
である (\muはYの平均値パラメータ)。
ここで\muにiの添字があるのは個人的にはまだ釈然とはしておらずE[Y_i] = \muではという気もするが、\mu_iと書かれていることが多い。あるモデルを考えたあとならば、\muが個体iの共変量の値によって変わるので、E[\frac{\partial l_i}{\partial \theta_i}|\boldsymbol{X_i}] = 0を考えて、E[Y_i|\boldsymbol{X_i}]= \mu_iとするのがしっくり来る。
また、微分と積分の交換が可能な元で、対数尤度の2階微分の期待値は1階微分と次のような関係にある：
E[\frac{\partial^2 l_i}{\partial \theta_i^2}] = E[ -(\frac{\partial l_i}{\partial \theta_i})^2]
尤度関数の1階微分の期待値が0であることから：
Var[\frac{\partial l_i}{\partial \theta_i}] = E[(\frac{\partial l_i}{\partial \theta_i})^2]
となる。
この両辺にそれぞれ、- \frac{1}{a(\phi)} \cdot \frac{\partial^2 b(\theta_i)}{\partial \theta_i^2}と(\frac{y_i - \frac{\partial b(\theta)}{\partial \theta_i}}{a(\phi)})^2を代入すると：
Var[y_i] = a(\phi) \cdot \frac{\partial^2 b(\theta_i)}{\partial \theta_i^2}
を得る。
また、これらの証明においては、尤度で仮定している確率分布と積分する確率分布が同じであることが前提になっている。すなわち分布の仮定が正しいという条件の下での性質ということになると思う。
最後に、ここで対数尤度の和ではなくある個体iの対数尤度(データの確率分布で積分するので確率変数Yがある)の関係を見ていることを強調しておく。これが後々の大数の法則につながる。
推定方程式
推定方程式とは
最尤推定を行う際には対数尤度の1階微分が0になるような推定値を探す。この対数尤度の1階微分をchain ruleで分けて考えると、後々都合がいい。なお、\beta_jはj番目の回帰係数として、\eta_i = X_i^T \betaである。個体iの対数尤度の1階微分は次のように書ける。
\frac{\partial l_i}{\partial \beta_j} = \frac{\partial l_i}{\partial \theta_i} \frac{\partial \theta_i}{\partial u_i} \frac{\partial u_i}{\partial \eta_i} \frac{\partial \eta_i}{\partial \beta_j}
この式に、expotential familyの式を当てはめると：
\frac{\partial l_i}{\partial \theta_i} = \frac{y_i - \frac{\partial b(\theta_i)}{\partial \theta_i}}{a(\phi)}
\frac{\partial \theta_i}{\partial u_i} = \frac{1}{\frac{\partial u_i}{\partial \theta_i}} = \frac{1}{\frac{\partial^2 b(\theta_i)}{\partial \theta_i^2}}
\frac{\partial \eta_i}{\partial \beta_j} = x_{ij}
個人的に混乱したポイントは、個体iの共変量X_iを考えるときはこれが縦ベクトルになっているので転置を取ってX_i^T \betaにするが、Xを行列全体で考えるときはX \betaで良い?ことである。また、3つ目の式については\beta_jでの微分なのでx_iではなくx_{ij}になっている点に注意。
数式がギチギチになってきたので\frac{\partial b(\theta_i)}{\partial \theta_i} = b'(\theta_i)、\frac{\partial^2 b(\theta_i)}{\partial \theta_i^2} = b''(\theta_i)として
\frac{\partial l_i}{\partial \beta_j} = \frac{y_i - b'(\theta_i)}{a(\phi) b''(\theta_i)} \frac{\partial \mu_i}{\partial \eta_i} x_{ij}
と書ける。
ここで、b''(\theta_i)はvariance functionV(\mu_i)と呼ばれ(\theta_iが\mu_iの関数なので\mu_iでこのvariance functionを表現できる)、a(\phi) \cdot V(\mu_i) = Var(y_i) となる。また、b'(\theta_i) = \mu_iである。よって、
\frac{\partial l_i}{\partial \beta_j} = \frac{y_i - \mu_i}{Var(y_i)} \frac{\partial \mu_i}{\partial \eta_i} x_{ij} := U_{ij}
と書ける。
対数尤度の1階微分をスコア関数ともいうのでU_{ij}とい別の記号を導入した(Dobson本に倣った)。このU_{ij}=0をサンプルサイズnで足し合わせ：
U_j = \sum_{i=1}^n U_{ij} = 0
K個のパラメータ\beta_jがあるとき：
U = \begin{bmatrix} U_1 \\ U_2 \\ \vdots \\ U_K \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix}
となる (以降は0のベクトルであっても単に0と書く)。
このU_{ij}=0をestimating equationといい、U_{j}=0をestimating equationsという(単変量の問題であってもestimating equationsと複数形という説を見かけた)。スコア関数=0が推定方程式となっているというふうに考えれば良いだろうか。また、このような推定方程式の解としてパラメータを推定することをM推定ともいう。
推定方程式の解の漸近正規性
推定方程式の解\widetilde{\beta}において、スコア関数を真のパラメータ\beta_0の周りでの一次テーラー展開で近似することを考える。なお、\frac{\partial l }{\partial \beta}はベクトル\betaでの微分であるのでK次元の縦ベクトルであり、\frac{\partial^2 l }{ \partial \beta \partial \beta'} はK×Kの行列であることに注意する。
\begin{align*}
U(\widetilde{\beta}) &= \sum_{i=1}^n \frac{\partial l_i }{ \partial \beta}(\widetilde{\beta})\\
&= \frac{\partial l }{ \partial \beta}(\beta_0) +
\frac{\partial^2 l }{ \partial \beta \partial \beta'} (\beta_0) \cdot (\widetilde{\beta} - \beta_0) \equiv 0
\end{align*}
これを整理すると、
\widetilde{\beta} - \beta_0 = \lbrack -\frac{\partial^2 l }{ \partial \beta \partial \beta'} (\beta_0) \rbrack^{-1} \lbrack \frac{\partial l }{ \partial \beta}(\beta_0) \rbrack
\implies \sqrt{n}(\widetilde{\beta} - \beta_0) = \lbrack -
\frac{\partial^2 l }{ \partial \beta \partial \beta'} (\beta_0) \rbrack^{-1} \lbrack \sqrt{n}\frac{\partial l }{\partial \beta}(\beta_0) \rbrack
となる。
2つ目のベクトル\sqrt{n}\frac{\partial l }{\partial \beta}(\beta_0)は、中心極限定理によって、平均が0のベクトルで分散行列がnI_n(\beta_0)の多変量正規分布に分布収束する。
平均0については、前述した対数尤度の期待値の性質からE[\frac{\partial l_i }{\partial \beta}] = 0であるため、大数の法則によって2つ目のベクトルが0のベクトルに確率収束することから言える(l=\sum_{i=1}^{n} l_iである)。
スコア関数の分散をフィッシャー情報量といい、サンプルサイズnの時I_n(\beta_0)と書く。スコア関数の期待値が0であることから、スコア関数の分散はスコア関数の直積I_n(\beta_0) = E[(\frac{\partial l }{\partial \beta})(\frac{\partial l }{\partial \beta})^T]で表せる。
上の2つ目のベクトルでは\sqrt{n}がついているので分散にnが掛けられている：
1つ目の逆行列の中身-\frac{\partial^2 l }{ \partial \beta \partial \beta'} (\beta_0)については、上述した単変量での関係を多変量に拡張して
Var[\frac{\partial l_i}{\partial \beta}] = E \lbrack (\frac{\partial l_i }{\partial \beta})(\frac{\partial l_i }{\partial \beta})^T \rbrack = - E \lbrack
\frac{\partial^2 l_i}{\partial \beta \beta'} \rbrack
と考えることができれば、大数の法則でnI_n(\beta_0)に確率収束する：
\begin{align*}
-\frac{\partial^2 l }{ \partial \beta \partial \beta'} (\beta_0) &= -n \cdot \frac{1}{n} \frac{\partial^2 l }{ \partial \beta \partial \beta'} (\beta_0) \\
&\rightarrow -n \cdot E \lbrack \frac{\partial^2 l_i }{ \partial \beta \partial \beta'} (\beta_0) \rbrack \\
&= n \cdot E \lbrack (\frac{\partial l_i }{\partial \beta})(\frac{\partial l_i }{\partial \beta})^T \rbrack \\
&= n \cdot Var[\frac{\partial l_i}{\partial \beta}] \\
&= n \cdot I_1(\beta_0) \\
&= nI_n(\beta_0)
\end{align*}
なお、上記ではデータがIIDと仮定してI_1がn個あるという計算になっている。
あとはスラツキーの定理によって、\sqrt{n}(\widetilde{\beta} - \beta_0)が平均0、分散がnI_n^{-1}の多変量正規分布に分布収束するので、
\widetilde{\beta} \rightarrow N(\beta_0, I_n(\beta_0)^{-1})
が推定方程式の解の漸近分布となる。ここで、I_n^{-1} I_n I_n^{-1} = I_n^{-1}となっていることに注意する。
フィッシャー情報量
さて、フィッシャー情報量をスコア関数の直積で表した場合、
\begin{align*}
I_n &= E \lbrack U U^T \rbrack \\
&= E \lbrack \lbrace \sum_{i=1}^n \frac{\partial l_i}{\partial \beta} \rbrace \lbrace \sum_{i=1}^n \frac{\partial l_i}{\partial \beta} \rbrace ^T \rbrack \\
&= E \lbrack \lbrace \sum_{i=1}^n \begin{bmatrix} \frac{y_i - \mu_i}{\text{Var}(y_i)} \frac{\partial \mu_i}{\partial \eta_i} x_{i1}
\\
\frac{y_i - \mu_i}{\text{Var}(y_i)} \frac{\partial \mu_i}{\partial \eta_i} x_{i2}
\\
\vdots
\\
\frac{y_i - \mu_i}{\text{Var}(y_i)} \frac{\partial \mu_i}{\partial \eta_i} x_{ik}
\end{bmatrix} \rbrace \lbrace \sum_{i=1}^n \begin{bmatrix} \frac{y_i - \mu_i}{\text{Var}(y_i)} \frac{\partial \mu_i}{\partial \eta_i} x_{i1}
\\
\frac{y_i - \mu_i}{\text{Var}(y_i)} \frac{\partial \mu_i}{\partial \eta_i} x_{i2}
\\
\vdots
\\
\frac{y_i - \mu_i}{\text{Var}(y_i)} \frac{\partial \mu_i}{\partial \eta_i} x_{ik}
\end{bmatrix} \rbrace ^T \rbrack \\
&= E
\lbrack \sum_{i=1}^n
\begin{bmatrix}
\left(\frac{y_i - \mu_i}{\text{Var}(y_i)}\right)^2 \left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2 x_{i1}^2 & \cdots & \left(\frac{y_i - \mu_i}{\text{Var}(y_i)}\right)^2 \left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2 x_{i1} x_{ik} \\
\vdots & \ddots & \vdots \\
\left(\frac{y_i - \mu_i}{\text{Var}(y_i)}\right)^2 \left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2 x_{ik} x_{i1} & \cdots & \left(\frac{y_i - \mu_i}{\text{Var}(y_i)}\right)^2 \left(\frac{\partial \mu_i}{\partial \eta_i}\right)^2 x_{ik}^2
\end{bmatrix}
\rbrack
\end{align*}
となる。
(j,k)成分を見ていくと次のようになる。Xは所与で定数でYについての期待値を取っていることを利用すると：
\begin{align*}
E \lbrack U_{j}U_{k} \rbrack &= E \lbrack \frac{\partial l_i}{\partial \beta_j} \frac{\partial l_i}{\partial \beta_k} \rbrack \\
&= E \lbrack \frac{(y_i - \mu_i)^2}{Var(y_i)^2} (\frac{\partial \mu_i}{\partial \eta_i})^2 x_{ij} x_{ik}\rbrack \\
&= \frac{E \lbrack(y_i - \mu_i)^2 \rbrack}{Var(y_i)^2} (\frac{\partial \mu_i}{\partial \eta_i})^2 x_{ij} x_{ik} \\
&= \frac{1}{Var(y_i)} (\frac{\partial \mu_i}{\partial \eta_i})^2 x_{ij} x_{ik}
\end{align*}
と簡単になる。ただし、ここでは\mu_iのモデル化が正しくできていない場合、E \lbrack(y_i - \mu_i)^2 \rbrack \not = Var(y_i)となるはずなので注意。
ここでw_i = \frac{1}{Var(y_i)} (\frac{\partial \mu_i}{\partial \eta_i})^2と置くと、二次形式のシグマ表現になっていることがわかるので行列表現にするとフィッシャー情報量は：
I_n = E \lbrack X^T W X \rbrack
と書ける。
Canonicalリンク関数を用いた場合の重みの簡略化表記
ところで、Canonicalリンク関数を使うとw_i = \frac{1}{Var(y_i)} (\frac{\partial \mu_i}{\partial \eta_i})^2をより簡単な形で表すことができる。
まず、Canonicalリンクを採用することで\theta_i = \eta_iとなることとb'(\theta_i) = \mu_iであることから
\begin{align*}
a(\phi)V(\mu_i) &= a(\phi)b''(\theta) \\
&= a(\phi)\frac{\partial \mu_i}{\partial \theta_i} \\
&= a(\phi)\frac{\partial \mu_i}{\partial \eta_i}
\end{align*}
なので、
\begin{align*}
w_i &= \frac{1}{Var(y_i)} (\frac{\partial \mu_i}{\partial \eta_i})^2 \\
&= \frac{1}{a(\phi)} \frac{\partial \mu_i}{\partial \eta_i}
\end{align*}
となる。さらに、Canonicalリンク関数を用いる場合、逆関数の微分によって
\begin{align*}
\frac{\partial \mu_i}{\partial \eta_i} &= \frac{\partial g^{-1}(\eta_i)}{\partial \eta_i} \\
&= \frac{1}{\frac{\partial g(\mu_i)} {\partial \mu_i}}
\end{align*}
となる。この量は、
\begin{align*}
\frac{1}{\frac{\partial g(\mu_i)}{\partial \mu_i}} =
\begin{cases}
\mu_i = V(\mu_i) & \text{if } \text{ ポアソン分布-logリンク関数}\\
\mu_i(1-\mu_i) = V(\mu_i) & \text{if } \text{ 二項分布分布-logitリンク関数 }
\end{cases}
\end{align*}
となる。これらをまとめると、
\begin{align*}
w_i &= \frac{1}{Var(y_i)} (\frac{\partial \mu_i}{\partial \eta_i})^2 \\
&= \frac{1}{a(\phi)} \frac{\partial \mu_i}{\partial \eta_i} \\
&= \frac{1}{a(\phi)} V(\mu_i)
\end{align*}
となる。また、同じ関係をスコア関数に対して代入すると、
\begin{align*}
\frac{\partial l_i}{\partial \beta_j}
&= \frac{(y_i - \mu_i) x_{ij} }{a(\phi) V(\mu_i)} \frac{\partial \mu_i}{\partial \eta_i} \\
&= \frac{(y_i - \mu_i) x_{ij} }{a(\phi) V(\mu_i)} V(\mu_i) \\
&= \frac{(y_i - \mu_i) x_{ij} }{a(\phi)}
\end{align*}
なので、
\frac{\partial^2 l_i}{\partial \beta_j \beta_k} = -\frac{x_{ij}}{a(\phi)} \frac{\partial \mu_i}{\partial \beta_k}
となる。これはy_iに依存しない量になっていることがポイントである。
Dispersion parameter
これまでw_i = \frac{1}{Var(y_i)} (\frac{\partial \mu_i}{\partial \eta_i})^2と置いていた重みをw_i = \frac{1}{V(u_i)} (\frac{\partial \mu_i}{\partial \eta_i})^2としてa(\phi)を外した形を考えると、フィッシャー情報量は\frac{1}{a(\phi)} E[X^TWX]と括り出すことができる。そして、少しsloppyな感じでa(\phi)を\phiに置き換えた上で\betaの分散共分散を考えると\phi E[(X^TWX)^{-1}]を得る。
Dispersion parameterをデータから推定したい場合には、平均が正しくモデル化されていれば、\frac{\sum_{i=1}^n (y_i - \mu_i)^2}{V(\mu_i)}が自由度(n-K)のカイ二乗分布に従うので（Kはパラメータの数）、IIDを仮定してMethod of Momentsアプローチで、
\begin{align*}
\frac{1}{\phi} E[\frac{\sum_{i=1}^n (y_i - \mu_i)^2}{V(\mu_i)}] = n-k \\
\implies \widehat{\phi} = \frac{\sum_{i=1}^n \lbrack (y_i - \mu_i)^2/V(\mu_i) \rbrack}{n-k}
\end{align*}
という推定量を得る。この値が何かしらのExpotential familyの分布を仮定した時点で何かしらの値に固定されてしまうため、データのばらつきをうまく表現できていないモデルとなってしまうことがある。モデルの想定より大きい場合をoverdispersion小さい場合をunderdispersionといい、ポアソンの場合はそれぞれ\phi>1と\phi<1に相当する。
Quasi-likelihood estimation
GLMにおいては、Var(y_i) = a(\phi)V(\mu_i)のように分散も平均に依存する関係(mean-variane relationship)であるため、\mu_iのモデリングが重要になってくる。実際、\phiは\betaに依存しないので、\betaの推定には関係がない。さらに、その平均のモデリングさえうまくいっていれば、\phiを推定することも可能であることを既に見た。そこで、y_iの分布自体を仮定するのではなく、分散とvariance functionの関係Var(y_i) = \phi V(\mu_i)だけを仮定する方法をQuasi-likelihood estimationという。通常のGLMの枠組みでは、y_iの分布を仮定するとV(\mu_i)と\phiが決まってしまうため、そこを自由にするわけである。
実際には、V(\mu_i)としては想定したy_iの分布から決まるvariance functionを用いつつ、\phiをデータから推定するという風に使われる。ポアソン分布を仮定しつつもover-/under-dispersionを考慮したいならば、1)ポアソン回帰を実施→2)\phiをデータから推定→3)共分散行列に\hat{\phi}を掛けることで、quasi-poisson regression modelを実行したことになる。
ロバストな共分散行列
Sandwich estimator
ヘシアン行列と対数尤度関数の直積の期待値が一致しない、すなわち、E \lbrack (\frac{\partial l_i }{\partial \beta})(\frac{\partial l_i }{\partial \beta})^T \rbrack \not = - E \lbrack \frac{\partial^2 l_i}{\partial \beta \beta'} \rbrackの場合には、ヘシアン行列の期待値がスコア関数の分散と一致しない。尤度関数の確率分布と実際のデータの生成されている確率分布が異なるときは、この関係が成り立たない場面の一つと思われる。
I_n^{-1} I_n I_n^{-1}の左と右に出現していたフィッシャー情報量の逆行列は、確率収束先であるヘシアンの期待値-E\lbrack \frac{\partial^2}{\partial \beta \partial \beta'} \rbrackをスコア関数の積の期待値に近似していたので、別の量に変えるべくA^{-1}と書くことにしてヘシアンのまま用いることにする：
A^{-1} = - E\lbrack \frac{\partial^2 l }{ \partial \beta \partial \beta'} (\beta_0) \rbrack ^{-1}
次に真ん中のI_nは、スコア関数の直積をE \lbrack(y_i - \mu_i)^2 \rbrack = Var(y_i)によって簡略化していたので、直接スコア関数の直積を用いてBと書くことにする：
B = E \lbrack \lbrace \sum_{i=1}^n \frac{y_i - \mu_i}{\text{Var}(y_i)} \frac{\partial \mu_i}{\partial \eta_i} X_{i} \rbrace \lbrace \frac{y_i - \mu_i}{\text{Var}(y_i)} \frac{\partial \mu_i}{\partial \eta_i} X_{i} \rbrace^T \rbrack
こうしてできたA^{-1} B A^{-1}がHuber sandwich estimatorというロバスト分散推定量であり、この分散行列の分散の平方根をとったものが、Huber-White standard errorsと呼ばれる推定量である。
Observed vs Expected フィッシャー情報量(FIM)
上述した通り、ヘシアンをベースにした量をフィッシャー情報量(FIM)として使うことが多いが、それでも期待値計算には積分が必要であるため、より簡易的に計算できる代替量が欲しい。そこでよく使われるのがObserved FIM：{J_n(\beta_0)} = - \sum_{i=1}^n \frac{\partial l_i}{\partial \beta \partial \beta'}(\beta_0)である。ただし、\beta_0は分からないので代わりに最尤推定量(MLE)である\widetilde{\beta}を代入したJ_{n}(\widetilde{\beta})で手元のデータから推定することになる。
もう一つ、Expected FIMというものがある。これは-\sum_{i=1}^n E \lbrack \frac{\partial l_i}{\partial \beta \partial \beta'}(\beta_0) \rbrackと定義されるものである。こちらも真のパラメータの値は分からないのでMLEの推定値を代入する。大学院のノートを見返してみたところ、どうやら仮説検定などにおいて、仮定されている確率分布から理論上のFIMを負のヘシアン行列として導出し、そこに帰無仮説のパラメータ値を代入することでExpected FIMの値として検定統計量に用いるらしい（スコア統計量; スコア検定）。
また重要なポイントとして、どちらのFIMもCanonicalリンク関数を使ってMLEの推定値を代入する場合一致する。なぜなら、前述したようにCanonicalリンク関数の場合、ヘシアンがy_iに依存しないためy_iの期待値を取る操作をしても変化がないからである。
推定アルゴリズム
Iteratively Weighted Least Square (IWLS)
上記推定方程式を解けばパラメータを推定することができるが、GLMが当てはめる確率分布は正規分布を除いて分散が平均に依存するために、一度に推定方程式を解くことができない。Var(y_i) = a(\phi) V(\mu_i)だからである（正規分布の場合、V(\mu_i)=1なのでこの問題が生じない）。よって、推定を繰り返していく最適化が必要になる。
スコア関数の一次テーラー展開を考えて、(m)回目の推定値から(m+1)回目の推定値の関係を近似すると
\begin{align*}
U(\beta^{(m+1)}) &= \sum_{i=1}^n \frac{\partial l_i}{\partial \beta}(\beta^{(m+1)})\\
&= \frac{\partial l}{\partial \beta}(\beta^{(m)}) + \frac{\partial^2 l }{ \partial \beta \partial \beta'} (\beta^{(m)}) \cdot (\beta^{(m+1)} - \beta^{(m)}) \equiv 0
\end{align*}
これを整理して、
\begin{align*}
\beta^{(m+1)} = \beta^{(m)} - \lbrack \frac{\partial^2 l }{ \partial \beta \partial \beta'} (\beta^{(m)}) \rbrack^{-1}\frac{\partial l}{\partial \beta}(\beta^{(m)})
\end{align*}
前述したようにヘシアンをスコアの直積の形で表した上で重み行列Wを(m)回目の推定値\beta^{(m)}を用いることとし、Yと\muをそれぞれy_iと\mu_iの縦ベクトルとすると、
\begin{align*}
\lbrack X^TW^{(m)}X \rbrack \beta^{(m+1)} &= \lbrack X^TW^{(m)}X \rbrack \beta^{(m)} - X^T W^{(m)}(Y-\mu) \frac{\partial \eta}{\partial \mu} \\
&= X^TW^{(m)} \lbrace X\beta^{(m)}-(Y-\mu)\frac{\partial \eta}{\partial \mu} \rbrace
\end{align*}
となる。Z=X\beta^{(m)}-(Y-\mu)\frac{\partial \eta}{\partial \mu}をworking response variableなどともいう。この最終的なX^TW^{(m)}X \beta = X^TW^{(m)}ZはZをアウトカムとしたWeighted Least Squareの形をしているので、いわゆる線形回帰の繰り返しと同等である。
Newton Raphson
IWLSの章はヘシアンをスコア関数の直積で近似した方法であったが、ヘシアンを用いて更新していく方法も考えられる。このとき、\beta^{(m+1)}の更新式の中にあるヘシアンに対してObserved FIMを使う方法がNewton Raphsonである。
Fisher Scoring
一方で、ヘシアンに対してExpected FIMを使う方法がFisher scoringである。ただし、上述したようにCanonicalリンク関数を用いる場合、Newton RaphsonとFisher scoringは一致する。
Modified poisson regression model
さて、ようやくここからModified poisson regression modelについてである。Zou (2004) によると、ポアソン回帰を当てはめた上で分散をSandwich estimatorで推定するというものなので、パラメータの分散についてのみ考えれば良さそうである。
Fuyama et al. (2021)によると、ポアソン回帰モデルはlog-binomialモデルのスコア関数がポアソン回帰とequivalentであることから、log-binomialモデルのリスク比を推定していることにもなるらしい。
まず尤度関数を考える。ポアソン分布の確率分布から愚直に考えていっても良いが、せっかくここまで整理してきたものがあるのでその結果を使うと、ポアソン回帰モデルのスコア関数は次のようになる：
\begin{align*}
U = \sum_{i=1}^n \frac{\partial l_i}{\partial \beta} &= \sum_{i=1}^n \frac{(y_i - \mu_i)}{a(\phi)} X_i \\
&= \lbrack y_i - exp(X_i^T \beta) \rbrack X_i \equiv 0
\end{align*}
参考までにlog-binomialモデルのスコア関数についても考える。ポアソン回帰モデルとリンク関数が同じなので、b''(\theta_i) = \mu_i (1-\mu_i)であることだけに注意していくと：
\begin{align*}
U = \sum_{i=1}^n \frac{\partial l_i}{\partial \beta} &= \sum_{i=1}^n \frac{(y_i - \mu_i)}{\frac{1}{n} \mu_i (1-\mu_i)} \mu_i X_i \\
&\approx \sum_{i=1}^n \lbrack \frac{y_i - exp(X_i^T \beta)}{1-exp(X_i^T \beta)} \rbrack X_i \equiv 0
\end{align*}
となる。なので\mu_i=exp(X_i^T \beta)が小さいようなレアイベントの場合には分母が1に近づき、ポアソン回帰モデルのスコア関数と近しくなるということだろう。
あとは、Sandwich estimatorに必要な部品であるAとBについて考えれば良い。まず、ヘシアンについては、その(j,k)成分を整理した結果を利用することで：
\begin{align*}
\frac{\partial^2 l_i}{\partial \beta_j \beta_k}
&= -\frac{x_{ij}}{a(\phi)} \frac{\partial \mu_i}{\partial \beta_k} \\
&= -\frac{x_{ij}}{a(\phi)} \frac{\partial \mu_i}{\partial \eta_i} \frac{\partial \eta_k}{\partial \beta_k} \\
&= -x_{ij} V(\mu_i)x_{ik} \\
&= -x_{ij} exp(X_i^T \beta) x_{ik}
\end{align*}
となるので、sandwich estimatorのBreadに相当するAの推定値は、最尤推定量\widetilde{\beta}を代入したObserved FIMを用いて、
\begin{align*}
\hat{A} = J_n(\widetilde{\beta}) &= \sum_{i=1}^n \frac{\partial^2 l_i}{\partial \beta \partial \beta'}(\widetilde{\beta}) \\
&= -\sum_{i=1}^n X_i \exp(X_i^T \widetilde{\beta}) X_i^T \\
&= -X^T W X \quad \text{where W is diagnoal and } w_{ii} = \exp(X_i^T \widetilde{\beta})
\end{align*}
となる。途中、iの要素で見ているのでn全体で見ているときと転置が逆になっている（X_iが縦ベクトルなので、X^Tの1列目に相当）。
同様にしてBの推定値は、
\begin{align*}
\hat{B} &= U(\widetilde{\beta}) U(\widetilde{\beta})^T \\
&= X^TGX \text{ where G is diagnoal and } g_{ii} = \lbrack y_i - exp(X_i^T \widetilde{\beta}) \rbrack^2
\end{align*}
となるので、modified poisson regression modelのsandwich estimatorは
\begin{align*}
Var(\widetilde{\beta})_{sandwich} &= \lbrack -X^T W X \rbrack^{-1} \lbrack X^TGX \rbrack \lbrack -X^T W X \rbrack^{-1} \\
&= \lbrack X^T W X \rbrack^{-1} \lbrack X^TGX \rbrack \lbrack X^T W X \rbrack^{-1}
\end{align*}
となる。
感想など
- 結局GLMをめちゃくちゃ復習することになってGWが終わった・・・。
- 推定方程式はGEEにも繋がるので、次はGEEをちゃんと勉強したい。
- Modified poisson regression modelの発想ってquasi-poisson modelと似ているので、その比較をちょっとやってみたいかもと思えてきた。
- ベクトルの微分が苦手すぎる & 行列表記に直すのが苦手すぎる。
- Pearson residualsが自由度(n-K)のカイ二乗分布に従う件、実はちゃんと証明できないので誰か助けて・・・。
- 次にRでのスクラッチ実装の記事を書く予定だけれど、Nが小さい場合うまくsandwich::sandwich()と一致しないレベルの実装である。より正確な実装ができるようになるにはまだ時間がかかりそう・・・。
Discussion