SAS—計(jì)算K-S值及畫圖

近來，時(shí)于夜半下雨，也常在夢(mèng)里被雨擾醒。究其原因，也是因?yàn)榻鼇砉ぷ鳠o趣，本身也只喜歡寫寫SAS或Python代碼，做模型，可近來卻連數(shù)據(jù)權(quán)限都沒有，萬灰俱滅。雨聲輕輕，卻也容易愁悶得睡不著。想著要去外包，卻因自己的猶豫不決也沒有去成。
好了，不說廢話了。昨天有人在微信上問KS的計(jì)算方式。今天介紹一下KS值吧。
先看一段程序：
data logistic;

   input accident age vision drive;

datalines;

1 17 1 1

1 44 0 0

1 48 1 0

1 55 0 0

1 75 1 1

0 35 0 1

0 42 1 1

0 57 0 0

0 28 0 1

0 20 0 1

0 38 1 0

0 45 0 1

0 47 1 1

0 52 0 0

0 55 0 1

1 68 1 0

1 18 1 0

1 68 0 0

1 48 1 1

1 17 0 0

1 70 1 1

1 72 1 0

1 35 0 1

1 19 1 0

1 62 1 0

0 39 1 1

0 40 1 1

0 55 0 0

0 68 0 1

0 25 1 0

0 17 0 0

0 45 0 1

0 44 0 1

0 67 0 0

0 55 0 1

1 61 1 0

1 19 1 0

1 69 0 0

1 23 1 1

1 19 0 0

1 72 1 1

1 74 1 0

1 31 0 1

1 16 1 0

1 61 1 0

;

proc logistic data=logistic  outest=model;

model accident(event='1')=age vision drive/selection=stepwise

sle=0.05 sls=0.05 outroc=roc;

output out= pred   p=p1 ;

run;

proc npar1way data=pred  noprint;

    class accident;

    var p1;

    output out=ks;

run;
以上是一段常規(guī)的計(jì)算KS值。
首先要了解什么是KS值，K-S值(kolmogorov-smirnov curve)將總體進(jìn)行n等分組并按照違約概率降序排列，計(jì)算每一等份中違約與正常百分比的累積分布，繪制出兩者之間的差值就是K-S曲線。其中K-S曲線中的最大值即為K-S統(tǒng)計(jì)量，取值范圍在0到1之間。
以下是計(jì)算K-S值及畫圖的宏程序：
/**data：邏輯回歸后包含驗(yàn)證結(jié)果的數(shù)據(jù)集，var：違約概率變量，status：分類變量，data1：切分后的變量，Mks：最大ks值,M：分組組數(shù)**/

%macro KS(data, var, status, data1, Mks,M);

proc sort data=&data;

by &var;

run;

proc sql noprint;

 select sum(&status) into:P from &data;

 select count(*) into :Ntot from &data;

 quit;

 %let N=%eval(&Ntot-&P);

data &data1;

set &data nobs=NN;

by &var;

retain tile 1  totP  0 totN 0;

Tile_size=ceil(NN/&M);

if &status=1 then totP=totP+&status;

else totN=totN+1;

Pper=totP/&P;

Nper=totN/&N;

if _N_ = Tile*Tile_Size then

  do;

  output;

   if Tile <&M then  

       do;

         Tile=Tile+1;

         SumResp=0;

       end;

  end;    

keep Tile Pper Nper;

run;

data temp;

     Tile=0;

     Pper=0;

     NPer=0;

run;

Data &data1;

  set temp &data1;

run;

data &data1;

    set &data1;

    Tile=Tile/&M;

    label Pper='Percent of Positives';

    label NPer ='Percent of Negatives';

    label Tile ='Percent of population';

    KS=NPer-PPer;

run;

proc sql noprint;

 select max(KS) into :&Mks from &data1;

run; quit;

proc datasets library=work nodetails nolist;

 delete temp ;

run;

quit;

%mend;

%macro PlotKS(data1);

 symbol1 value=dot color=red   interpol=join  height=1;

 legend1 position=top;

 symbol2 value=dot color=blue  interpol=join  height=1;

 symbol3 value=dot color=green interpol=join  height=1;

proc gplot data=&data1;

  plot( NPer PPer KS)*Tile / overlay legend=legend1;

 run;

quit;

    goptions reset=all;

%mend;