スコア関数ベクトル = $\mathbf{0}$ の解を求める方法

本ページでは、推定方程式の解を求めることを想定して、NLMIXED Procedure を用いてスコア関数ベクトル = $\mathbf{0}$ の解を求めたい場合、どのようなプログラムを実行すればよいのかを考えます。

スコア関数ベクトル = 0 の解

スコア関数ベクトル = $\mathbf{0}$ の解を求める計算は、非線形連立方程式を解く事に相当します。 したがって、$p$ 個の非線形連立方程式 $F_1(y)=F_2(y)=\ldots=F_p(y)=0$ を解くという事で、$\{F_1(y)\}^2+\{F_2(y)\}^2+\ldots+\{F_p(y)\}^2$ を最小化してあげればよいという事が分かります。

例えば、ロジスティック回帰分析の場合、スコア関数ベクトルは、

\[ S(\mathbf{\beta})=\frac{\partial p}{\partial \mathbf{\beta}}\left\{\frac{y-p}{p(1-p)}\right\}, \]

\[ \frac{1}{1+\exp(-\mathrm{logit})}, \]

\[ \mathrm{logit}=\mu+x_1\beta_1+x_2\beta_2+\ldots, \]

であり、

\[ -\{S(\mathbf{\beta})\}^{\mathrm{T}}\{S(\mathbf{\beta})\}, \]

を最大化すればよいことになります。

プログラム

以下のプログラムでは、data ステップでダミーデータを生成し、NLMIXED Procedure を使ったスコア関数ベクトル = $\mathbf{0}$ の解、NLMIXED Procedure を使った最尤推定、LOGISTIC Procedure を使った最尤推定を行っています。


  data x;
    call streaminit(20120427);
    n1 = 100; n2 = 100; n3 = 100;
    mu = 0;
    beta1 = 0.5;
    beta2 = 1.0;
    do i = 1 to n1;
      x1 = 0; x2 = 0;
      logit = mu + x1 * beta1 + x2 * beta2;
      p = 1 / (1 + exp(-logit));
      y = rand('Bernoulli', p);
      output;
    end;
    do i = n1 + 1 to n1 + n2;
      x1 = 1; x2 = 0;
      logit = mu + x1 * beta1 + x2 * beta2;
      p = 1 / (1 + exp(-logit));
      y = rand('Bernoulli', p);
      output;
    end;
    do i = n1 + n2 + 1 to n1 + n2 + n3;
      x1 = 0; x2 = 1;
      logit = mu + x1 * beta1 + x2 * beta2;
      p = 1 / (1 + exp(-logit));
      y = rand('Bernoulli', p);
      output;
    end;
  run;

  /* スコア関数ベクトル = 0 の解 */
  proc nlmixed data = x(keep = y x1 x2);
    parms mu = 0 beta1 = 0 beta2 = 0;
    logit = mu + x1 * beta1 + x2 * beta2;
    p = 1 / (1 + exp(-logit));
    dp = exp(-logit) / (1 + exp(-logit))**2;
    Score = -(
      ( 1 * dp * (y - p) / (p * (1 - p))) ** 2 +
      (x1 * dp * (y - p) / (p * (1 - p))) ** 2 +
      (x2 * dp * (y - p) / (p * (1 - p))) ** 2
    );
    model y ~ general(Score);
  run;

  /* 対数尤度関数を最大化 */
  proc nlmixed data = x(keep = y x1 x2);
    parms mu = 0 beta1 = 0 beta2 = 0;
    logit = mu + x1 * beta1 + x2 * beta2;
    p = 1 / (1 + exp(-logit));
    model y ~ binary(p);
  run;

  /* LOGISTIC Procedure */
  proc logistic data = x;
    class x1 x2 / param = ref ref = first;
    model y(event='1') = x1 x2;
  run;
  

実行すると以下の結果が得られました。


  /* スコア関数ベクトル = 0 の解 */
  Parameter  Estimate        SE

  mu           0.5322    0.3033
  beta1       -0.1683    0.3671
  beta2        0.2679    0.3830


  /* 対数尤度関数を最大化 */
  Parameter  Estimate        SE

  mu           0.5322    0.2071
  beta1       -0.1683    0.2902
  beta2        0.2679    0.2994


  /* LOGISTIC Procedure */
  パラメータ      推定値    標準誤差

  Intercept       0.5322      0.2071
  x1             -0.1683      0.2902
  x2              0.2679      0.2994
  

推定値自体は全て一致しました。 また、スコア関数ベクトル = $\mathbf{0}$ の解を求めた場合、当然ですが標準誤差の値をそのまま使うことができませんので、自力で計算する必要があります。 本ページでは上記の方法でも正しいパラメータ推定値を得ることができる事をしめしましたが、信頼できる方法がある場合 (本ページの例では LOGISTIC Procedure を使った方が一番良いでしょう) はそれを利用した方が良いでしょう。

履歴

  • 2012/04/27 公開