留一交叉驗(yàn)證及SAS代碼

在數(shù)據(jù)量很少，用什么模型？我們總結(jié)過當(dāng)數(shù)據(jù)量很少時如何選擇模型和方法，以使得數(shù)據(jù)能夠最大限度的得到利用。

其中有一個方法就是做交叉驗(yàn)證。

我有備選的模型G(x1, x2, x3)， G(x1, x5, x6), F(x1, x2, x3)，想知道哪一個預(yù)測的效果好。不能做樣本內(nèi)預(yù)測（就是用樣本訓(xùn)練出模型，再用同樣的樣本代到模型中看準(zhǔn)確度），樣本量太少，再分成訓(xùn)練集和測試集就更少的可憐了，怎么辦？

K折交叉驗(yàn)證可以充分利用少樣本的信息。

K折交叉驗(yàn)證是將樣本分成K個子樣本集，拿出其中的K-1個子樣本集來訓(xùn)練模型，用剩下的1個子樣本集來對模型進(jìn)行驗(yàn)證；再拿出K-1個訓(xùn)練模型，留下另外1個（與上一步的不同）子樣本集進(jìn)行驗(yàn)證......，如此交叉驗(yàn)證K次，每個子樣本集驗(yàn)證1次，平均K次的結(jié)果作為一個模型的預(yù)測效果。

而本文想說的留一交叉驗(yàn)證（Leave-one-out cross validation, LOOCV）就是這種方法的極端情況：

假設(shè)只有10個樣本（真的很小啊），每次拿出其中9個來訓(xùn)練模型，用剩下一個進(jìn)行測試，得到一個測試結(jié)果（真實(shí)值與預(yù)測值的差異）；再拿出另外9個進(jìn)行訓(xùn)練，留下另外一個進(jìn)行測試......如此驗(yàn)證10次（每個樣本都能輪到一次驗(yàn)證樣本），將10次的預(yù)測效果平均，就可以評價(jià)這個模型的好壞。

留一交叉驗(yàn)證就是留下1個單樣本，將其他所有樣本拿來做訓(xùn)練。可以充分利用小樣本的信息。

下面分享一下數(shù)說君留一交叉驗(yàn)證的SAS代碼，樣本量假設(shè)為30：

*樣本量30;

%let K=30;

 

*為數(shù)據(jù)增加一個變量：index，標(biāo)識出觀測值的ID（從1到30）;

data sample;

set sample;

index = _n_;

run;

 

  *用全30個樣本建模看一下;

  proc reg data=sample;

  model y= x1 x2 x3;

  run;

 

data sample_all;

set sample;

selected = .;

replicate = .;

run;

 

*每次模型將一個樣本留作測試，其他用來訓(xùn)練樣本，重復(fù)30次，那我們就建立30個數(shù)據(jù)集，并將這30個數(shù)據(jù)集合在一起;

%macro generateData;

 

%do i = 1%to &K;

 

*每次選擇一個觀測值，其selected=0，意為測試樣本，其他29個均為1，意為訓(xùn)練樣本。

data temp;

set sample;

if index = &i thenselected = 0;

else selected = 1;

replicate =&i;

run;

 

data sample_all;

set sample_all temp;

run;

 

%end;

 

data sampleOut;

set sample_all;

where selected ^= .;

run;

 

%mend;

 

*運(yùn)行宏;

%generateData;

 

 

*slelected=0的樣本意為一個數(shù)據(jù)集中的測試樣本，我們看一下是否每個觀測值都輪到一次測試;

  proc print data=sampleOut;

  where Selected=0;

  var Selected id;

  run;

 

  data sampleOut;

  set sampleOut;

   if selected then new_y=y;

  run;

 

*計(jì)算selected=0的樣本、也就是測試樣本的預(yù)測值;

  proc reg data=sampleOut;

    model new_y=x1 x2 x3;

   by replicate;

   outputout=out1(where=(new_y=.)) predicted=y_hat;

  run;

 

  data out2;

  set out1;

   d=y-y_hat;

   absd=abs(d);

  run;

 

*畫出預(yù)測值與真實(shí)值的散點(diǎn)圖;

proc gplot data=out2;

plot y*y_hat;

run;

 

  proc summary data=out2;

  var d absd;

  outputout=out3std(d)=rmse mean(absd)=mae sum(d)=sumd;

  run;

 

 *計(jì)算y與y_hat的相關(guān)系數(shù)，以及y=y_hat的R方（這個常被用于評價(jià)模型的擬合好壞）;

 proc corr data=out2 pearson out=corr(where=(_TYPE_='CORR'));

  var y ;

  with y_hat;

 run;

 

 data corr;

  set corr;

  Rsqrd=y**2;

 run;