読者です 読者をやめる 読者になる 読者になる

スコア信頼区間

スコア信頼区間の求め方を復習。

積二項分布モデルの尤度から、リスク比 $\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に書いたものを手直し?したものです