スコア信頼区間の求め方を復習。
積二項分布モデルの尤度から、リスク比 $\phi=p_1/p_2$ のスコア信頼区間を求める。
スコア関数と Fisher 情報行列から $\mathbf{S}^t(\phi)\mathbf{I}^{-1}\mathbf{S}(\phi)=z^2$ をみたす $\phi_{L}$ と $\phi_{U}$ を求める、nuisance parameter には $\phi$ を既知として求めた最尤推定量をぶち込む。
今回は最適な数値計算法などは深く考えずに、二分法 (bisection method) で解く。
本当のところを言うと、意外にも IML の非線形最適化は最大化と最小化しか指定できないようなので、手っ取り早く二分法を使ったんですけどね。
$\mathbf{S}^t(\phi) \mathbf{I}^{-1} \mathbf{S}(\phi)-z^2=0$ を満たす解を探索します。
proc iml;
/* search range & step & criterion */
MIN = 0.01; MAX = 10; STEP = 0.01; EPS = 1E-8;
/* bisection method */
start bisection(a, b) global(EPS);
if obj_f(a) = 0 then return(a);
if obj_f(b) = 0 then return(b);
do while(1);
c = a + (b - a) / 2;
if b - a < EPS then go to break;
if -EPS <= obj_f(c) & obj_f(c) <= EPS then go to break;
if obj_f(a) * obj_f(c) < 0 then b = c;
else if obj_f(b) * obj_f(c) < 0 then a = c;
end;
break:
ret = c || obj_f(c);
return(ret);
finish;
start bisection_search(Ans) global(MIN, MAX, STEP);
AnsN = 0;
do a = MIN to MAX by STEP;
b = a + STEP;
if obj_f(a) * obj_f(b) < 0 then do;
if AnsN = 0 then Ans = bisection(a, b);
else Ans = Ans // bisection(a, b);
AnsN = AnsN + 1;
end;
end;
finish;
/* objective function */
start obj_f(phy);
n = 20; m = 20;
x = 7; y = 13;
z = Quantile('Normal', 0.975);
p2 =
(m + x + phy * (n + y) -
sqrt( (m + x + phy * (n + y)) ** 2 -
4 * phy * (n + m) * (x + y))) /
(2 * phy * (n + m));
p1 = phy * p2;
Tphy =
(x - n * p1) ** 2 / (1 - p1) ** 2 *
( (1 - p1) / (n * p1) + (1 - p2) / (m * p2) )
- z ** 2;
return(Tphy);
finish;
Ans = J(1, 2, 0); names = {Estimate Error};
call bisection_search(Ans);
print Ans [colname=names];
quit;
\[
\begin{array}{ll}
L(\phi, p_2 \mid x, y)&=
{m \choose x}p_1^x(1-p_1)^{m-x}{n \choose y}p_2^y(1-p_2)^{n-y}\\&=
{m \choose x}(\phi p_2)^x(1-\phi p_2)^{m-x}{n \choose y}p_2^y(1-p_2)^{n-y},
\end{array}
\] \[
\log L(\phi, p_2 \mid x, y)=
x\log(\phi p_2)+(m-x)\log(1-\phi p_2)+y\log(p_2)+(n-y)\log(1-p_2)+c,
\] \[
\frac{\partial \log L}{\partial \phi}=
\frac{x}{\phi}-\frac{(m-x)p_2}{1-\phi p_2}=
\frac{x-m \phi p_2}{\phi(1-\phi p_2)}
\] \[
\frac{\partial^2 \log L}{\partial \phi^2}=-\frac{x}{\phi^2}-\frac{(m-x)p_2^2}{(1-\phi p_2)^2}
\] \[
\mathrm{E}\left(-\frac{\partial^2 \log L}{\partial \phi^2}\right)
=\frac{\mathrm{E}(x)}{\phi^2}+\frac{(m-\mathrm{E}(x))p_2^2}{(1-\phi p_2)^2}=
\frac{m p_2}{\phi}+\frac{m p_2^2}{1-\phi p_2}=\frac{m p_2}{\phi(1-\phi p_2)}
\]
$\phi$ を与えたときの最尤推定値
\[
\left.\frac{\partial \log L}{\partial p_2}\right|_{p_2=\hat{p}_2}=
\frac{x}{\hat{p}_2}-\frac{(m-x)\phi}{1-\phi \hat{p}_2}+\frac{y}{\hat{p}_2}-\frac{n-y}{1-\hat{p}_2}=0
\] \[
\hat{p}_2=
\frac{n+x+(m+y)\phi-\sqrt{\{n+x+(m+y)\phi\}^2-4(x+y)(m+n)\phi}}{2(m+n)\phi}
\] $\hat{p}_2$ の分子のルートにつく符号が $+$ となる解は、$\hat{p}_2$ が $1$ を超える場合がある (これは証明可能だが省略)。
2008/11/22に書いたものを手直し?したものです
