NLMIXED Procedureを使った最尤推定

最近使わなくなってしまったNLMIXED Procedureについて少しだけ書いてみようと思う。

尤度 (対数尤度) を明示的に書けて、ヘンテコなモデル (部分的に変量効果が欠けていたり…) を想定しない場合には簡単に最尤推定して、$H_0: \theta_i = 0$ (パラメータ $\theta_i$ が $0$ かどうか) の検定もしくはパラメータの漸近信頼区間 (Wald 型) を求めることができる。
多変量正規分布の変量効果 (1サブジェクトのみ、施設と時間とかは無理) を仮定した、パラメータに対して非線形な混合効果モデルの解析ができる。
FO (テーラー展開一次近似)、Gauss 求積法 (Laplacian; Gauss求積法の積分点が1の場合)、Importance Sampling などなど、周辺尤度を求めるための数値積分、モンテカルロ積分が一応使える。
簡単な例として、2 項分布 $\mathrm{Bin}(n, p)$ に従う乱数を $100$ 個発生させて、そのデータから最尤法によりパラメータ $p$ を推定してみるプログラムを書いてみる。

/* n = 10, p = 0.2 */
data dummy;
  do i = 1 to 100;
    x = rand('binomial', 0.2, 10);
    output;
  end;
proc nlmixed data = dummy;
  parm p = 0.1; /* 初期値を適当に指定 */
  likelihood = Comb(10, x) * p ** x * (1 - p) ** (10 - x);
  logl = log(likelihood);
  model i ~ general(logl);
run;

下から 2 行目の model i の i は、データセット内の変数であれば何を指定してもいい(今のところ推定値とかは変わらない模様)。