Kengo Nagashima



NLMIXED Procedure を使ってスコア関数ベクトル = $\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 を使った方が一番良いでしょう) はそれを利用した方が良いでしょう。

履歴