留一交叉驗證及SAS代碼

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

其中有一個方法就是做交叉驗證。

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

K折交叉驗證可以充分利用少樣本的信息。

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

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

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

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

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

*樣本量30;

%let K=30;

 

*為數(shù)據(jù)增加一個變量：index，標識出觀測值的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;

 

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

%macro generateData;

 

%do i = 1%to &K;

 

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

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;

 

*運行宏;

%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;

 

*計算selected=0的樣本、也就是測試樣本的預測值;

  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;

 

*畫出預測值與真實值的散點圖;

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;

 

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

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

  var y ;

  with y_hat;

 run;

 

 data corr;

  set corr;

  Rsqrd=y**2;

 run;