SAS逻辑回归之二分类

数据集这里用的是australian,有14个自变量Xi,一个因变量Y，Y值只取0或1。

代码如下：

/*逻辑回归数据集australian（690个观测值，每个含14个属性，目标变量y(0、1)）*/  
    /*导入数据集australian到逻辑库work中*/  
    proc import out=aus  
        datafile="\\vmware-host\Shared Folders\桌面\SAS\\data\australian.csv"      /*文件路径*/  
        dbms=csv replace;                               /*文件类型指定*/  
        delimiter=',';  
        getnames=yes;                                   /*是否将第一列作为列名*/  
    run;  
      
    /*查看数据集*/  
    proc print data=aus;  
    run;  
      
    /****************************  使用交叉验证法选择最优模型  *****************************************/  
      
    /*利用10-折交叉验证法计算测试集上的预测准确率*/  
    %let k=10;                            /*定义宏变量-交叉验证的折数k*/  
    %let rate=%sysevalf((&k-1)/&k);       /*给出交叉验证的样本抽样比率（因为宏变量k的本质是文本，不能直接参与运算，要将其视为数字计算要用%evalf or %sysevalf）*/  
      
    /*生成交叉验证的10个样例，保存在cv中*/  
    proc surveyselect data=aus   
                  out=cv            /*生成的样例全部放在数据集cv中*/  
                  seed=158  
                  samprate=&rate    /*抽样比率设定，宏变量rate的调用要加&*/  
                  outall            /*输出全部数据*/  
                  reps=10;          /*指定样本重复的次数*/  
    run;  
      
    /*交叉验证的生成数据集中，selected列为1表示该行为训练集样本，0表示测试集样本，这里为new_y赋值，  
      若selected=1，则可获得Y的值，若为0，该行的new_y为空。接下来给出new_y为空行的预测值。*/  
    data cv;  
      set cv;  
       if selected then new_y=Y;  
      run;  
      
    /*逻辑回归主程序 - 10折交叉验证*/  
    ods output parameterestimates=paramest   /*输出交叉验证的参数估计值*/  
               association=assoc;            /*输出交叉验证的C统计量*/  
    proc logistic data=cv des;            /*des控制以Y=1来建模*/  
        /* class new_y (param=ref ref='yes');  若new_y是分类变量，则用class对其参数化处理，这里选择处理方式为ref，以“yes”作为参考水平，以便于后续odds的计算*/  
        model new_y=X1-X14 / SELECTION=STEPWISE SLE=0.1 SLS=0.1;  
        by replicate;                         /*以交叉验证的组别来分组建模*/  
        output out=out1(where=(new_y=.))    /*只给出测试集的预测结果（即new_y为空的样本）*/  
               p=y_hat;  
    run;  
    ods output close;  
      
    data out1;    
        set out1;  
        if y_hat>0.5 then pred=_LEVEL_ ;     /* PHAT为logistic方程针对每个观察体计算的属于该组别的概率，若PHAT>0.5，则属于该组别（这里level为1），否则，属于另一组别 */  
        else pred=0;                     /* 本例为二分类，概率依照level（1）计算，因此另一类为0 */  
    run;  
      
    /*汇总交叉验证的结果*/  
    /*计算预测准确率(测试集中预测准确的样本占预测总样本的概率)*/  
    data out2;  
        set out1;  
        if Y=pred then d=1;  /*d为真实值和预测值的误差,这里设无误差为1，有误差为0*/  
        else d=0;  
    run;  
      
    proc summary data=out2;  
        var d;  
        by replicate;  
        output out=out3 sum(d)=d1;   /*预测正确的个数*/  
    run;  
      
    data out3;  
        set out3;  
        acc=d1/_freq_;   /*预测准确率*/  
        keep replicate acc;  
    run;  
      
    /*结果中加入交叉验证的C统计量（度量观测值和预测值之间的一致性，越大越好）*/  
    data assoc;  
        set assoc;  
        where label2="c";  
        keep replicate cvalue2;  
    run;  
      
    /*合并交叉验证的统计结果*/  
    data cvresult;  
    merge assoc(in=ina) out3(in=inb);  
    keep replicate cvalue2 acc;  
    run;  
      
    proc print data=cvresult;  
    title'交叉验证组号、c统计量、预测准确率';  
    run;  
      
    title '交叉验证最优模型选择：组号、预测准确率';  
    ods output SQL_Results=cvparam;      /*保存最优模型结果在cvparam数据集中*/  
    proc sql ;  
        select replicate,acc from cvresult having acc=max(acc);  
    quit;  
    ods output close;  
      
      
      
    /***************** 以交叉验证的最优结果组进行建模  *************************************/  
    /*以最优组合从cv的10个样例中拿出最优样例，作为训练集和测试集*/  
    /*取出最优组号对应的selected=1的行，作为训练集train，其余的作为测试集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 - 数据集Neur - 建模方法 STEPWISE ---------------------------';  
      
    /* 逻辑回归主程序 - 通过训练集建立logistic模型*/  
    proc logistic data=train  DES                    /*根据分类值从大到小选择建模组别，此处为yes*/  
                        covout outest=Nout_step  /*输出建模参数估计值及变量间的协方差矩阵*/  
                        outmodel=model            /*输出建模结果（若想要通过已有的建模结果来预测新数据集，这里可以用inmodel实现）*/  
                        simple;                          /*输出变量的简单统计量*/   
            /* class Y (param=ref ref='yes');  若Y是分类变量，则用class对其参数化处理，这里选择处理方式为ref，以“yes”作为参考水平，以便于后续odds的计算*/  
            MODEL Y=X1-X14                             /*logistic回归模型：反应变量=自变量1 2 3...*/  
                          / SELECTION=STEPWISE           /*选择建模方式 - 逐步排除法*/    
                            SLE=0.1 SLS=0.1              /*变量在模型中的显著程度，默认为0.05*/   
                            details                      /*输出模型界定的过程，包括自变量的检定和相关系数的值*/  
                            lackfit                      /*输出HL拟合优度*/  
                            RSQ                          /*模型解释度R方*/  
                            STB                          /*输出标准化模型后的参数*/  
                            CL                           /*参数估计和置信区间*/  
                            itprint                      /*输出分析每个步骤的统计量*/  
                            corrb                        /*输出变量的相关矩阵*/  
                            covb                         /*输出变量的协方差矩阵*/  
                            ctable                       /*输出不同阈值下的二分类变量的分组情况，类似于ROC曲线上的每个点的值*/  
                            influence                    /*输出观察体中每个变量统计量，便于找出对分析结果影响力较大的观察体*/  
                            IPLOTS ;                     /*针对influence的结果画出图形，影响力过高的观察体在图形上都会显得特别突出*/  
     score data=train outroc=train_roc;            /*通过score语句得到训练集上一系列的sensitivity和specificity，画出ROC曲线*/  
     score data=test   
           out=test_pred   
           outroc=test_roc;                      /*通过score来预测测试集，结果保存在test_pred中，画出ROC曲线*/  
    OUTPUT out=train_pred                        /*保存模型预测结果在该数据集中，数据集中包含的列由以下添加的统计量给出*/  
                P=PHAT  lower=LCL upper=UCL              /*输出文件中包含每个观察体属于logistic方程预测组别的概率，用PHAT作列名，LCL和UCL为置信上下限的值*/  
                RESCHI=RESCHI  RESDEV=RESDEV             /*Pearson残差和偏差残差，找出与模型不太符合的观察体*/  
                DIFCHISQ=DIFCHISQ  DIFDEV=DIFDEV         /*检测观察体对对皮尔森卡方适合度和对偏激统计量的影响程度，越大说明与模型越不符*/  
                                                         /* 还可加入的统计量：C、CBAR、DFBETAS、H、XBETA、STDXBETA */  
               / ALPHA=0.1;                              /*界定P值的信赖度，默认为0.05，对应信赖度为95%，这里为90%*/  
    run;       
    quit;  
      
    /*   
    逻辑回归主程序 - 根据logistic模型对测试集进行预测（有需要时可使用独立的logistic过程对新数据进行预测）  
    proc logistic inmodel=model;                 
        SCORE data=test                           
              outroc=predict_roc;                
    run;        
    */  
      
    /* 训练集的预测结果中只给出了预测概率，接下来根据0.5分界将观察体归到具体的类中，加一列“pred”（预测组别）*/  
    data train_pred;    
        set train_pred;  
        if PHAT>0.5 then pred=_LEVEL_ ;     /* PHAT为logistic方程针对每个观察体计算的属于该组别的概率，若PHAT>0.5，则属于该组别（这里level为1），否则，属于另一组别 */  
        else pred=0;                      
    run;  
      
    /* 输出混淆矩阵 - 训练集*/  
    ods output CrossTabFreqs=ct_train;   /*保存混淆矩阵表（训练集）*/  
    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 '训练集上的预测准确率';  
    run;  
      
      
    /* 输出混淆矩阵及准确率等指标 - 测试集*/  
    ods output CrossTabFreqs=ct_test; /*保存混淆矩阵表（测试集）*/  
    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 '测试集上的预测准确率';  
    run;