スコア関数ベクトル = $\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 公開