allvarcomb

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

回帰モデル等の解析におけるモデル選択を考えます。 例えば、説明変数が $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$ の関係

%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

履歴

• 2014/09/06 マクロの改良 (結果変数 y を指定できる様に変更; 連続型共変量への対応; Resultsのリストを表示しないように変更)
• 結果変数 y を指定できる様に変更
• 離散型、連続型共変量が混ざっている場合に対応 (引数 class を追加; class = 1: 全部離散型, class = 変数名リスト: 指定変数だけ離散型, class = (ブランク): 全部連続型)
• Resultsのリストを表示しないように変更
• 2010/06/25 公開