SAS邏輯回歸之二分類

數(shù)據(jù)集這里用的是australian,有14個(gè)自變量Xi,一個(gè)因變量Y，Y值只取0或1。

代碼如下：

/*邏輯回歸數(shù)據(jù)集australian（690個(gè)觀測(cè)值，每個(gè)含14個(gè)屬性，目標(biāo)變量y(0、1)）*/  
    /*導(dǎo)入數(shù)據(jù)集australian到邏輯庫(kù)work中*/  
    proc import out=aus  
        datafile="\\vmware-host\Shared Folders\桌面\SAS\\data\australian.csv"      /*文件路徑*/  
        dbms=csv replace;                               /*文件類型指定*/  
        delimiter=',';  
        getnames=yes;                                   /*是否將第一列作為列名*/  
    run;  
      
    /*查看數(shù)據(jù)集*/  
    proc print data=aus;  
    run;  
      
    /****************************  使用交叉驗(yàn)證法選擇最優(yōu)模型  *****************************************/  
      
    /*利用10-折交叉驗(yàn)證法計(jì)算測(cè)試集上的預(yù)測(cè)準(zhǔn)確率*/  
    %let k=10;                            /*定義宏變量-交叉驗(yàn)證的折數(shù)k*/  
    %let rate=%sysevalf((&k-1)/&k);       /*給出交叉驗(yàn)證的樣本抽樣比率（因?yàn)楹曜兞縦的本質(zhì)是文本，不能直接參與運(yùn)算，要將其視為數(shù)字計(jì)算要用%evalf or %sysevalf）*/  
      
    /*生成交叉驗(yàn)證的10個(gè)樣例，保存在cv中*/  
    proc surveyselect data=aus   
                  out=cv            /*生成的樣例全部放在數(shù)據(jù)集cv中*/  
                  seed=158  
                  samprate=&rate    /*抽樣比率設(shè)定，宏變量rate的調(diào)用要加&*/  
                  outall            /*輸出全部數(shù)據(jù)*/  
                  reps=10;          /*指定樣本重復(fù)的次數(shù)*/  
    run;  
      
    /*交叉驗(yàn)證的生成數(shù)據(jù)集中，selected列為1表示該行為訓(xùn)練集樣本，0表示測(cè)試集樣本，這里為new_y賦值，  
      若selected=1，則可獲得Y的值，若為0，該行的new_y為空。接下來給出new_y為空行的預(yù)測(cè)值。*/  
    data cv;  
      set cv;  
       if selected then new_y=Y;  
      run;  
      
    /*邏輯回歸主程序 - 10折交叉驗(yàn)證*/  
    ods output parameterestimates=paramest   /*輸出交叉驗(yàn)證的參數(shù)估計(jì)值*/  
               association=assoc;            /*輸出交叉驗(yàn)證的C統(tǒng)計(jì)量*/  
    proc logistic data=cv des;            /*des控制以Y=1來建模*/  
        /* class new_y (param=ref ref='yes');  若new_y是分類變量，則用class對(duì)其參數(shù)化處理，這里選擇處理方式為ref，以“yes”作為參考水平，以便于后續(xù)odds的計(jì)算*/  
        model new_y=X1-X14 / SELECTION=STEPWISE SLE=0.1 SLS=0.1;  
        by replicate;                         /*以交叉驗(yàn)證的組別來分組建模*/  
        output out=out1(where=(new_y=.))    /*只給出測(cè)試集的預(yù)測(cè)結(jié)果（即new_y為空的樣本）*/  
               p=y_hat;  
    run;  
    ods output close;  
      
    data out1;    
        set out1;  
        if y_hat>0.5 then pred=_LEVEL_ ;     /* PHAT為logistic方程針對(duì)每個(gè)觀察體計(jì)算的屬于該組別的概率，若PHAT>0.5，則屬于該組別（這里level為1），否則，屬于另一組別 */  
        else pred=0;                     /* 本例為二分類，概率依照level（1）計(jì)算，因此另一類為0 */  
    run;  
      
    /*匯總交叉驗(yàn)證的結(jié)果*/  
    /*計(jì)算預(yù)測(cè)準(zhǔn)確率(測(cè)試集中預(yù)測(cè)準(zhǔn)確的樣本占預(yù)測(cè)總樣本的概率)*/  
    data out2;  
        set out1;  
        if Y=pred then d=1;  /*d為真實(shí)值和預(yù)測(cè)值的誤差,這里設(shè)無誤差為1，有誤差為0*/  
        else d=0;  
    run;  
      
    proc summary data=out2;  
        var d;  
        by replicate;  
        output out=out3 sum(d)=d1;   /*預(yù)測(cè)正確的個(gè)數(shù)*/  
    run;  
      
    data out3;  
        set out3;  
        acc=d1/_freq_;   /*預(yù)測(cè)準(zhǔn)確率*/  
        keep replicate acc;  
    run;  
      
    /*結(jié)果中加入交叉驗(yàn)證的C統(tǒng)計(jì)量（度量觀測(cè)值和預(yù)測(cè)值之間的一致性，越大越好）*/  
    data assoc;  
        set assoc;  
        where label2="c";  
        keep replicate cvalue2;  
    run;  
      
    /*合并交叉驗(yàn)證的統(tǒng)計(jì)結(jié)果*/  
    data cvresult;  
    merge assoc(in=ina) out3(in=inb);  
    keep replicate cvalue2 acc;  
    run;  
      
    proc print data=cvresult;  
    title'交叉驗(yàn)證組號(hào)、c統(tǒng)計(jì)量、預(yù)測(cè)準(zhǔn)確率';  
    run;  
      
    title '交叉驗(yàn)證最優(yōu)模型選擇：組號(hào)、預(yù)測(cè)準(zhǔn)確率';  
    ods output SQL_Results=cvparam;      /*保存最優(yōu)模型結(jié)果在cvparam數(shù)據(jù)集中*/  
    proc sql ;  
        select replicate,acc from cvresult having acc=max(acc);  
    quit;  
    ods output close;  
      
      
      
    /***************** 以交叉驗(yàn)證的最優(yōu)結(jié)果組進(jìn)行建模  *************************************/  
    /*以最優(yōu)組合從cv的10個(gè)樣例中拿出最優(yōu)樣例，作為訓(xùn)練集和測(cè)試集*/  
    /*取出最優(yōu)組號(hào)對(duì)應(yīng)的selected=1的行，作為訓(xùn)練集train，其余的作為測(cè)試集test*/  
    proc sql ;  
        create table train as  
        select * from cv where replicate in (select replicate from cvparam)  
        having selected=1;  
        create table test as  
        select * from cv where replicate in (select replicate from cvparam)  
        having selected=0;  
    run;  
      
    TITLE '--------Logistic Regression - 數(shù)據(jù)集Neur - 建模方法 STEPWISE ---------------------------';  
      
    /* 邏輯回歸主程序 - 通過訓(xùn)練集建立logistic模型*/  
    proc logistic data=train  DES                    /*根據(jù)分類值從大到小選擇建模組別，此處為yes*/  
                        covout outest=Nout_step  /*輸出建模參數(shù)估計(jì)值及變量間的協(xié)方差矩陣*/  
                        outmodel=model            /*輸出建模結(jié)果（若想要通過已有的建模結(jié)果來預(yù)測(cè)新數(shù)據(jù)集，這里可以用inmodel實(shí)現(xiàn)）*/  
                        simple;                          /*輸出變量的簡(jiǎn)單統(tǒng)計(jì)量*/   
            /* class Y (param=ref ref='yes');  若Y是分類變量，則用class對(duì)其參數(shù)化處理，這里選擇處理方式為ref，以“yes”作為參考水平，以便于后續(xù)odds的計(jì)算*/  
            MODEL Y=X1-X14                             /*logistic回歸模型：反應(yīng)變量=自變量1 2 3...*/  
                          / SELECTION=STEPWISE           /*選擇建模方式 - 逐步排除法*/    
                            SLE=0.1 SLS=0.1              /*變量在模型中的顯著程度，默認(rèn)為0.05*/   
                            details                      /*輸出模型界定的過程，包括自變量的檢定和相關(guān)系數(shù)的值*/  
                            lackfit                      /*輸出HL擬合優(yōu)度*/  
                            RSQ                          /*模型解釋度R方*/  
                            STB                          /*輸出標(biāo)準(zhǔn)化模型后的參數(shù)*/  
                            CL                           /*參數(shù)估計(jì)和置信區(qū)間*/  
                            itprint                      /*輸出分析每個(gè)步驟的統(tǒng)計(jì)量*/  
                            corrb                        /*輸出變量的相關(guān)矩陣*/  
                            covb                         /*輸出變量的協(xié)方差矩陣*/  
                            ctable                       /*輸出不同閾值下的二分類變量的分組情況，類似于ROC曲線上的每個(gè)點(diǎn)的值*/  
                            influence                    /*輸出觀察體中每個(gè)變量統(tǒng)計(jì)量，便于找出對(duì)分析結(jié)果影響力較大的觀察體*/  
                            IPLOTS ;                     /*針對(duì)influence的結(jié)果畫出圖形，影響力過高的觀察體在圖形上都會(huì)顯得特別突出*/  
     score data=train outroc=train_roc;            /*通過score語句得到訓(xùn)練集上一系列的sensitivity和specificity，畫出ROC曲線*/  
     score data=test   
           out=test_pred   
           outroc=test_roc;                      /*通過score來預(yù)測(cè)測(cè)試集，結(jié)果保存在test_pred中，畫出ROC曲線*/  
    OUTPUT out=train_pred                        /*保存模型預(yù)測(cè)結(jié)果在該數(shù)據(jù)集中，數(shù)據(jù)集中包含的列由以下添加的統(tǒng)計(jì)量給出*/  
                P=PHAT  lower=LCL upper=UCL              /*輸出文件中包含每個(gè)觀察體屬于logistic方程預(yù)測(cè)組別的概率，用PHAT作列名，LCL和UCL為置信上下限的值*/  
                RESCHI=RESCHI  RESDEV=RESDEV             /*Pearson殘差和偏差殘差，找出與模型不太符合的觀察體*/  
                DIFCHISQ=DIFCHISQ  DIFDEV=DIFDEV         /*檢測(cè)觀察體對(duì)對(duì)皮爾森卡方適合度和對(duì)偏激統(tǒng)計(jì)量的影響程度，越大說明與模型越不符*/  
                                                         /* 還可加入的統(tǒng)計(jì)量：C、CBAR、DFBETAS、H、XBETA、STDXBETA */  
               / ALPHA=0.1;                              /*界定P值的信賴度，默認(rèn)為0.05，對(duì)應(yīng)信賴度為95%，這里為90%*/  
    run;       
    quit;  
      
    /*   
    邏輯回歸主程序 - 根據(jù)logistic模型對(duì)測(cè)試集進(jìn)行預(yù)測(cè)（有需要時(shí)可使用獨(dú)立的logistic過程對(duì)新數(shù)據(jù)進(jìn)行預(yù)測(cè)）  
    proc logistic inmodel=model;                 
        SCORE data=test                           
              outroc=predict_roc;                
    run;        
    */  
      
    /* 訓(xùn)練集的預(yù)測(cè)結(jié)果中只給出了預(yù)測(cè)概率，接下來根據(jù)0.5分界將觀察體歸到具體的類中，加一列“pred”（預(yù)測(cè)組別）*/  
    data train_pred;    
        set train_pred;  
        if PHAT>0.5 then pred=_LEVEL_ ;     /* PHAT為logistic方程針對(duì)每個(gè)觀察體計(jì)算的屬于該組別的概率，若PHAT>0.5，則屬于該組別（這里level為1），否則，屬于另一組別 */  
        else pred=0;                      
    run;  
      
    /* 輸出混淆矩陣 - 訓(xùn)練集*/  
    ods output CrossTabFreqs=ct_train;   /*保存混淆矩陣表（訓(xùn)練集）*/  
    ods trace on;  
    proc freq data=train_pred;  
        tables Y*pred;  
    run;  
    ods trace off;  
    ods output close;  
      
    proc sql;  
        create table acc1 as  
        select sum(percent) from ct_train where (Y=pred and Y ^=.);  
    proc print data=acc1;  
    title '訓(xùn)練集上的預(yù)測(cè)準(zhǔn)確率';  
    run;  
      
      
    /* 輸出混淆矩陣及準(zhǔn)確率等指標(biāo) - 測(cè)試集*/  
    ods output CrossTabFreqs=ct_test; /*保存混淆矩陣表（測(cè)試集）*/  
    proc freq data=test_pred;  
        tables F_Y*I_Y ;  
    run;  
    ods output close;  
      
    proc sql;  
        create table acc2 as  
        select sum(percent) from ct_test where (F_Y=I_Y and F_Y ^='');  
    proc print data=acc2;  
    title '測(cè)試集上的預(yù)測(cè)準(zhǔn)確率';  
    run;