NLMIXED Procedureを使った最尤推定 (一般化線型モデルの場合)

NLMIXED Procedureで一般化線型モデルの最尤推定を行ってみよう。
今回はロジスティック回帰分析。
説明変数xのオッズ比が大体2になるようにダミーデータを発生させ、NLMIXED Procedureでロジスティック回帰分析を行うプログラムを以下に示した。

data dummy;
  OR = 2; px1 = 0.4;
  px0 = px1 / (OR * (1 - px1) + px1);
  do i = 0 to 300;
    y = 1;
    x = rand('bernoulli', px1);
    output;
  end;
  do i = 0 to 300;
    y = 0;
    x = rand('bernoulli', px0);
    output;
  end;
  drop OR px1 px0;
proc nlmixed data = dummy;
  parm alpha = 1 beta1 = 1.5;
  logit = alpha + beta1 * x; /* 線型モデルの部分 */
  p = 1 / (1 + exp(-logit)); /* リンク関数 */
  likelihood = comb(1, y) * p ** y * (1 - p) ** (1 - y);
  logl = log(likelihood);
  model i ~ general(logl);
  estimate 'exp(beta1)' exp(beta1);
run;

ダミーデータの発生は一手間かかる。
オッズ比が2になるように発生させるには上のように、オッズ比を固定して解いてあげればよい。
あとはロジットが線形になるように書いて、尤度はベルヌーイ分布を書く。
estimate statementで、オッズ比に戻してやると2に近い推定結果が得られる。
これに変量効果をぶち込んだり、少し複雑なことをやろうとするときの出発点としていい例かなと。