ブートストラップ、クロスバリデーション、ジャックナイフ法について

質問

SASにおいて、ブートストラップ法やクロスバリデーション法(ジャックナイフ法)を実行することはできますか。

回答

一つの方法として、以下のリンクに記載のサンプルマクロ%JACKおよび%BOOTが用意されています。
[Downloads]タブにマクロを定義するためのSASプログラムがあります。

Sample 24982: Jackknife and Bootstrap Analyses

http://support.sas.com/kb/24/982.html

上記のサンプルマクロでは、任意のプロシジャによって出力された推定値をもとにブートストラップ法や

ジャックナイフ法に基づいた信頼区間をいくつかの方法(例:パーセンタイル法)によって求めることができます。

本マクロでは単純無作為抽出しか実行できません。

なお、いくつかのプロシジャではクロスバリデーション法やブートストラップ法を用いています。

一例を以下に記載します。

・DISCRIMプロシジャ

CROSSVALIDATEオプションではleave-one-outクロスバリデーションを用いて誤判別率を推定しています。

・LOGISTICプロシジャ

MODELステートメントのCTABLEオプションでは、誤判別率を求める時にクロスバリデーションを

近似した方法を採用しています。

・PLSプロシジャ

CV= オプションでは、成分の数を決める際にクロスバリデーション法を実行することができます。

・MULTTESTプロシジャ

Cochran-Armitage検定やFisherの正確検定などのP値を調整するために、ブートストラップ法を用いることができます。

詳細は以下のリンクのSAS Noteをご参照ください。

Usage Note 22220: Procedures with bootstrapping, cross validation, or jackknifing capabilities

http://support.sas.com/kb/22/220.html

またSAS/STATのSURVEYSELECTプロシジャを用いることでも、ブートストラップ法を実装することができます。

以下に簡単なサンプルを記載いたします。ブートストラップ法は大きく4つのステップから構成されています。

サンプルプログラム(Compute a bootstrap confidence interval in SAS)：

data sample(keep=x);
   set Sashelp.Iris(where=(Species="Virginica") rename=(SepalWidth=x));
run;
 
/* 1. オリジナルデータの統計値を計算する: 歪度 = 0.366 */
proc means data=sample nolabels Skew;  var x;  run;

%let NumSamples = 5000;       /* ブートストラップのリサンプル数 */

/* 2. ブートストラップサンプルをたくさん生成する */
proc surveyselect data=sample NOPRINT seed=1
     out=BootSSFreq(rename=(Replicate=SampleID))
     method=urs              /*リサンプル */
     samprate=1              /* 各ブートストラップサンプルにはN個の観測値がある */
     /* 頻度を抑えるためのオプション */
     reps=&NumSamples;       /* NumSamplesブートストラップリサンプルを生成する */
run;

/* 3. 各ブートストラップサンプルの統計量を計算する */
proc means data=BootSSFreq noprint;
   by SampleID;
   freq NumberHits;
   var x;
   output out=OutStats skew=Skewness;  /* 標本分布に近似 */
run;

/* ブートストラップ分布を視覚化する */
title "Bootstrap Distribution";
%let Est = 0.366;
proc sgplot data=OutStats;
   label Skewness= ;
   histogram Skewness;
   /* オプション: 観測値で基準線と95%信頼区間を描く */
   refline &Est / axis=x lineattrs=(color=red) 
                  name="Est" legendlabel="Observed Statistic = &Est";
   refline -0.44737 0.96934  / axis=x lineattrs=(color=blue) 
                  name="CI" legendlabel="95% CI";
   keylegend "Est" "CI";
run;

proc means data=OutStats nolabels N StdDev;
   var Skewness;
run;

/* 4. 統計的推論を行うために近似した標本分布を使う */
proc univariate data=OutStats noprint;
   var Skewness;
   output out=Pctl pctlpre =CI95_
          pctlpts =2.5  97.5       /* 95%ブートストラップ信頼区間を計算する */
          pctlname=Lower Upper;
run;
 
proc print data=Pctl noobs; run;

実行結果:

分析変数 : x
歪度
0.3659491

分析変数 : Skewness
N	標準偏差
5000	0.3624289

CI95_Lower	CI95_Upper
-0.44737	0.96934