$K$層の積二項分布モデルの無条件の尤度関数から導かれるスコア方程式
$k$番目の層のCase群の曝露ありの人数が二項分布$X_{11k} \sim Bin(n_{1+k}, p_{1k})$に従い、$k$番目の層のControl群の曝露ありの人数が二項分布$X_{01k} \sim Bin(n_{0+k}, p_{0k})$に従い、$k$番目の層のオッズ比を$\psi_k=\frac{p_{1k}/(1-p_{1k})}{p_{0k}/(1-p_{0k})}$で表したもとで、共通オッズ比 ($\psi_1=\ldots=\psi_K=\psi$) の仮定をおいたときの積二項分布モデルの対数尤度関数は、$\mathbf{p}_{0}=(p_{01},\ldots,p_{0K})^T$とおくと、
\[
\begin{aligned}
l(\psi, \mathbf{p}_{0})=
\sum_{k=1}^K &
\log \binom{n_{1+k}}{x_{11k}}+
x_{11k} \log \psi p_{0k} -
x_{11k} \log \{1-p_{0k}(1-\psi)\} +\\&
(n_{1+k}-x_{11k}) \log (1-p_{0k})-
(n_{1+k}-x_{11k}) \log \{1-p_{0k}(1-\psi)\}+\\&
\log \binom{n_{0+k}}{x_{01k}}+
x_{01k}\log p_{0k}+
(n_{0+k}-x_{01k})\log (1-p_{0k})
\end{aligned}
\] となります。
対数尤度関数を$\psi$で偏微分すると、
\[
\frac{\partial l(\psi, \mathbf{p}_{0})}{\partial \psi}=
\sum_{k=1}^K
\frac{x_{11k}}{\psi}-
n_{1+k}\frac{p_{0k}}{1-p_{0k}(1-\psi)}
\] が得られ、これを$=0$と置いて解くことで
\[
\sum_{k=1}^K x_{11k}=
\sum_{k=1}^K n_{1+k}\frac{\hat{\psi} \hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}
\] という関係式が得られます。
また、定義から$\frac{\hat{\psi} \hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}=\hat{p}_{1k}$なので、$\sum_{k=1}^K x_{11k}=\sum_{k=1}^K n_{1+k} \hat{p}_{1k}=\sum_{k=1}^K E^A[X_{11k}; \psi]$と書くこともできます。
次に、対数尤度関数を$p_{0k}$で偏微分すると、
\[
\begin{aligned}
\frac{\partial l(\psi, \mathbf{p}_{0})}{\partial p_{0k}}&=
\frac{x_{11k}}{p_{0k}}+
\frac{x_{11k}(1-\psi)}{1-p_{0k}(1-\psi)}-
\frac{n_{1+k}-x_{11k}}{1-p_{0k}}+
\frac{(n_{1+k}-x_{11k})(1-\psi)}{1-p_{0k}(1-\psi)}+
\frac{x_{01k}}{p_{0k}}-
\frac{n_{0+k}-x_{01k}}{1-p_{0k}}\\&=
\frac{x_{11k}+x_{01k}}{p_{0k}}-
\frac{n_{1+k}+n_{0+k}-x_{11k}-x_{01k}}{1-p_{0k}}+
\frac{n_{1+k}(1-\psi)}{1-p_{0k}(1-\psi)}\\&=
\frac{1}{1-p_{0k}}\left(
\frac{x_{11k}+x_{01k}}{p_{0k}}-
\frac{n_{1+k}\psi}{1-p_{0k}(1-\psi)}-
n_{0+k}
\right)
\end{aligned}
\] が得られます。
1:1マッチングの場合の$\hat{p}_{0k}$
1:1マッチングの場合、$n_{1+k}=n_{0+k}=1$となり、$x_{11k}+x_{01k}=0, 1, 2$の3パターンだけになります。
下で使うため、3パターンの場合について対数尤度関数を$p_{0k}$で偏微分して$=0$とおいて解き、それぞれの場合の$\hat{p}_{0k}$を導出しておきます。
$x_{11k}+x_{01k}=0$のとき
\[
\frac{0}{\hat{p}_{0k}}-
\frac{\hat{\psi}}{1-\hat{p}_{0k}(1-\hat{\psi})}-
1=0
\iff
\hat{p}_{0k}=-\frac{\hat{\psi}\hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}=-\hat{p}_{1k}
\] $\hat{p}_{0k}=-\hat{p}_{1k}$が成り立つのは$\hat{p}_{0k}=0$の時しかないため、$\hat{p}_{0k}=0$が得られます。
$x_{11k}+x_{01k}=1$のとき
\[
\frac{1}{\hat{p}_{0k}}-
\frac{\hat{\psi}}{1-\hat{p}_{0k}(1-\hat{\psi})}-
1=0
\iff
\hat{p}_{0k}=\frac{1}{1 \pm \sqrt{\hat{\psi}}}
\] $\hat{p}_{0k}=\frac{1}{1-\sqrt{\hat{\psi}}}$の方は$\hat{p}_{1k}$が負になるので、$\hat{p}_{0k}=\frac{1}{1+\sqrt{\hat{\psi}}}$が得られます。
$x_{11k}+x_{01k}=2$のとき
\[
\frac{2}{\hat{p}_{0k}}-
\frac{\hat{\psi}}{1-\hat{p}_{0k}(1-\hat{\psi})}-
1=0
\iff
\hat{p}_{0k}=1,\frac{2}{1-\hat{\psi}}
\] $\hat{p}_{0k}=\frac{2}{1-\hat{\psi}}$の方は$\hat{p}_{0k}>1$か$\hat{p}_{0k}<0$になってしまうので、$\hat{p}_{0k}=1$が得られます。
1:1マッチングの場合の$\hat{\psi}$
上で求めた対数尤度関数を$\psi$で偏微分して$=0$とおいて解いた式について、1:1マッチングの場合を考えることで$\hat{\psi}$を得ることができます。
全体で$K$個の層がありますが、
- $x_{11k}=1, x_{01k}=1$となっている層の数を$n_{11}$
- $x_{11k}=1, x_{01k}=0$となっている層の数を$n_{10}$
- $x_{11k}=0, x_{01k}=1$となっている層の数を$n_{01}$
- $x_{11k}=0, x_{01k}=0$となっている層の数を$n_{00}$
と書くことにします ($K=n_{11}+n_{10}+n_{01}+n_{00}$です)。
$\sum_{k=1}^K x_{11k}=\sum_{k=1}^K n_{1+k}\frac{\hat{\psi} \hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}$に対して、今度は4パターンについて、上で導出した$\hat{p}_{0k}$も代入してそれぞれ当てはめると、
\[
\sum_{k=1}^K x_{11k}=\sum_{k=1}^K n_{1+k}\frac{\hat{\psi} \hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}
\iff
n_{11}+n_{10}+0+0=
n_{11}+n_{10}\frac{\sqrt{\hat{\psi}}}{1+\sqrt{\hat{\psi}}}+n_{01}\frac{\sqrt{\hat{\psi}}}{1+\sqrt{\hat{\psi}}}+0
\] が得られます。
これを解くと、
\[
\hat{\psi}=\left(\frac{n_{10}}{n_{01}}\right)^2
\] が得られ、1:1マッチングの時の無条件の最尤推定量が条件付き最尤推定量$\hat{\psi}_C=\frac{n_{10}}{n_{01}}$を二乗したものになっていることがわかります。
参考文献
- Breslow N. Odds ratio estimators when the data are sparse. Biometrika 1981; 68(1): 73–84.
ロジスティック回帰版
ironman さんが解いています。
1:1 Matched Case ControlでUnconditionalな分析をした場合のMLEのバイアスの導出について気になる方はぜひご覧ください。 #iron勉強メモ
— ironman (@hommedefer3) 2022年6月27日
1) @triadsou さんに教えていただいた方法https://t.co/TFS9fQl6DX
2) 1)を参考にしてロジスティック回帰としてやってみた方法https://t.co/T3TLPfF3Cs https://t.co/304NTDdGGJ
