Kengo Nagashima



総当たり法によるモデル選択

回帰モデル等の解析におけるモデル選択を考えます。 例えば、説明変数が $x_1$ と $x_2$ だった場合、$y = x_1$, $y = x_2$, $y = x_1 + x_2$ の 3 通りのモデルのうち、どのモデルが一番良いモデルなのかを検討する状況を想定します。

上の例のように 2 変数程度であれば、コピー & ペーストして変数を修正したものを作れば充分です。 しかし、$p$ 変数になった場合、$2^p - 1$ 通りの組み合わせを生成する必要があり、5 変数もあれば コピー & ペーストは効率的ではないでしょう。

$p$ と $2^p - 1$ の関係

$p$$2^p - 1$
531
663
7127
8255
9511
101023
$p$$2^p - 1$
1532767
201048575
2533554431
301073741823
401099511627775
501125899906842620

%allvarcomb() マクロ

今回は $p$ 変数の全ての組み合わせを列挙する %allvarcomb() マクロを作成しました。
引数は、var: 変数リスト (スペース区切り)、out: 出力データセット名、info: 組み合わせ数と変数の数を記録するデータセット名、です。

%macro bincombinations(p, out);
data &out.;
  array v[*] v1-v%eval(2**&p.);
  retain f 0;
  do j = 1 to &p.;
    pn = 2**&p./2**j;
    do i = 1 to 2**&p.;
      v[i] = f;
      if mod(i, pn) = 0 then f = mod(f+1, 2);
    end;
    output;
  end;
  keep v1-v%eval(2**&p.);
proc transpose data = &out.
  out = &out.(drop=_NAME_) prefix = v;
run;
%mend bincombinations;

%macro allvarcomb(var, out, info);
data &out.;
  string = "&var.";
  do until(position <= 0);
    count + 1;
    call scan(string, count, position, length, ' ');
    word = substrn(string, position, length);
    if word ^= "" then output;
  end;
  keep word;
data _null_; set &out.;
  call symput("vnum", compress(_n_));
proc transpose data = &out.
  out = &out.(drop=_NAME_) prefix = vn;
  var word;
%bincombinations(&vnum., vs);
data &out.; set vs;
  if _n_ = 1 then set &out.;
data &out.; set &out.;
  array avn[*] vn1-vn&vnum.;
  array av[*] v1-v&vnum.;
  do i = 1 to &vnum.;
    if av[i] = 0 then avn[i] = "";
  end;
  if _n_ ^= 1 then output;
  call symput("cnum", compress(_n_)-1);
  keep vn1-vn&vnum.;
data &info.;
  cnum = &cnum.; vnum = &vnum.;
proc datasets lib = work; delete vs; run; quit;
%mend allvarcomb;

モデル選択基準を計算する部分のマクロです。
真ん中の部分を解析目的に応じて変更して利用します。
以下の例では、GENMOD Procedure で正規線型モデルを当てはめて AIC を求めています。

%macro modelselect(y, x, data, out, class=);
proc datasets lib = work; delete &out.;
%if       "&class." = "1" %then %let class = class &x.;
%else %if "&class." = ""  %then %let class = ;
%else                           %let class = class &class.;
%allvarcomb(&x., allvarcomb, info);
data _null_; set info;
  call symput("cnum", cnum);
  call symput("vnum", vnum);
run;
ods listing close;
ods results off;
%do i = 1 %to &cnum.;
data _null_;
  set allvarcomb(firstobs=&i. obs=&i.);
  %do j = 1 %to &vnum.;
    call symput("vn&j.", vn&j.);
  %end;
proc datasets lib = work; delete temp;
/************************************
modify for any purpose
************************************/
proc genmod data = &data.;
&class.;
model &y. =
  %do j = 1 %to &vnum.;
    &&vn&j. 
  %end; /
  link = id dist = normal;
  ods output ModelFit = temp;
run;
data temp; set temp(obs=7 firstobs=7); keep value;
/************************************/
data &out.;
  set %if "&i." ^= "1" %then &out.; temp;
%end;
ods results on;
ods listing;
data &out.; merge &out. allvarcomb;
proc datasets lib = work; delete allvarcomb info temp;
run; quit;
%mend modelselect;

$7$ 変数のテストデータを生成し、$AIC$ によるモデル選択を行います。
$7$ 変数なので、$2^7 - 1 = 127$ 通りのモデルを全て当てはめて、$AIC$ を求めます。

/* テストデータの生成 */
data testdata;
  do var = 0, 1; do uvar = 0, 1; do xbar = 0, 1;
  do avar = 0, 1; do x1 = 0, 1; do x2 = 0, 1; do xe = 0, 1;
    y = rand("Normal", 0*var + 1*uvar + 5*xbar + 5*avar + 1*x1 + 0*x2 + 1*xe, 1);
    output;
  end; end; end; end; end; end; end;
run;

filename fdummy dummy;
proc printto log = fdummy;
run;
%modelselect(y, var uvar xbar avar x1 x2 xe, testdata, aic, class = 1);
proc printto log = log;
run;

proc rank data = aic out = aic;
var value; ranks aicrank;
proc sort; by aicrank;
proc print;
  where aicrank <= 10;
run;
   Value    vn1    vn2     vn3     vn4     vn5    vn6    vn7    aicrank

341.2267           uvar    xbar    avar    x1            xe         1
343.0875           uvar    xbar    avar    x1     x2     xe         2
343.2253    var    uvar    xbar    avar    x1            xe         3
345.0861    var    uvar    xbar    avar    x1     x2     xe         4
374.1257           uvar    xbar    avar    x1                       5
376.0197           uvar    xbar    avar    x1     x2                6
376.1247    var    uvar    xbar    avar    x1                       7
378.0187    var    uvar    xbar    avar    x1     x2                8
380.3804           uvar    xbar    avar                  xe         9
382.2795           uvar    xbar    avar           x2     xe        10

モデル y = uvar + xbar + avar + x1 + xe が選択されました。
これはテストデータで生成したモデルと一致しているので、うまく選択できているようです。

また、%allvarcomb() でどのようなデータが生成されているのかを見たい場合は、以下を実行します。

%allvarcomb(var uvar xbar avar x1 x2 xe, allvarcomb, info, class = 1);
proc print data = allvarcomb; run;

ダウンロード

modelselect.sas

履歴