SAS Logistic回归：一个完整例子

 SAS Logistic回归:代码及输出报告详解

这篇将作为五一后一个讲稿的阅读材料之一，先整出来就搁这。如果没有耐心读下去，你可以立即转到以下的参考资料，该篇所有的知识都来自它们：

Cody, R.F. and Smith, J.K. Applied Statistics and the SAS Programming
Language,4th
ed..NJ: 
Prentice-Hall,1997.这书已经出第五版
了，北大图书馆只有这第四版。非常容易上手的一本书，前半部分用input和datalines让读者专心做统计，后半部分从导入导出数据开始阐述SAS的通用编程语言。这本书用的是SAS8.这里我们只关注它第九章Multiple-Regression
Analysis的最后Logistic Regression部分。我这篇的例子即来于此，有简化；

SAS
OnlineDoc V8
, 或者SAS OnlineDoc V9
,
是要花功夫熟悉它们的结构了。以前我四处下载了数G的电子书，现在才发觉还是它们好使。体例上V8和V9一样，你找到SAS/STAT-->SAS/STAT
User's Guide
-->The LOGISTIC Procedure,
就可以跟着学习了，文字都非常简明。

Logistic回归处理因变量是分类型变量如“0、1”的情形。一下就假设你至少对它模模糊糊有些印象，比如说我们用p表示正例（如输出变量为“1”）的概率，那么p/(1-p)就被称作odds
ratio，对p做logit变换记做logit(p),它等于log(p/(1-p),我们回归方程的形式就如logit(p)=log(p/(1-p)=a+bx,你可以把它理解成向量形式。

假设我们有一个数据，45个观测值，四个变量，包括：

age（年龄，数值型）；

vision（视力状况，分类型，1表示好，0表示有问题）；

drive（驾车教育，分类型，1表示参加过驾车教育，0表示没有） 和

一个分类型输出变量accident（去年是否出过事故，1表示出过事故，0表示没有）。我们的目的就是要考察前三个变量与发生事故的关系。

为了保持程序的可读性，以下我就直接在代码中一行一行敲入数据，在实际工作中，超过20个观测值的数据就要写程序导入了。另外，之间的语句都是注释，SAS运行时会把它们忽略掉：

 

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
descending

;

     model
accident=age vision
drive; run
;

 

proc
logistic
data=logistic
descending;

   model accident=age vision
drive  /

               
selection
=forward

;run;

 

 

 

 

 

 

 

 

 

 

----------------------------------------------------------------------------------------------------------------------------

用上面这个过程步乙替代过程步甲，再运行一遍。这两个过程步的输出结果大同小异，只是过程步乙多了个forward选项。以下用的是过程步乙的输出结果，其中黑体字是输出结果本身，我做的注释语句以红笔描出,包括数字编号。 

                 (1
)
 
The SAS System 10:47 Tuesday, May 4, 2007 1

                            
The LOGISTIC Procedure

                            
  Model Information

Data
Set                            
WORK.LOGISTIC

Response
Variable                         
accident

Number of Response
Levels           
2

Number of
Observations                
45

Model                                             
binary logit

Optimization
Technique                 Fisher's
scoring

 

                            
Response
Profile
                                     

Ordered 
Value    
accident          Total
Frequency

  
1              
1                   
25

      
2              
0                   
20

Probability modeled is accident=1.

(1)
给出了本模型的基本信息，意思大多自明。需要注意的是Response
Profile
中，accident=1排在首位。前面我们说过，SAS的Logistic回归方程log(odds)默认的形式是处理那个变量值比较小的，加上descending选项后，accident=1就排在首位了

(2)
  Forward
Selection Procedure

Step 0. Intercept entered:

                         
Model Convergence Status

Convergence criterion (GCONV=1E-8) satisfied.

                         
Residual Chi-Square Test

Chi-Square              DF               
Pr > ChiSq

 
10.7057                  
3                   
 0.0134

 

Step 1
. Effect vision
entered:

                             
Model Convergence Status

Convergence criterion (GCONV=1E-8)
satisfied. 

                              
Model Fit
Statistics                                               

Criterion     Intercept Only     Intercept and
Covariates
AIC
                 63.827                   59.244

SC
                   65.633                  
62.857

-2 Log
L
           61.827                 
55.244

                      
Testing Global Null Hypothesis: BETA=0

Test                      
Chi-Square   
DF    Pr
> ChiSq
Likelihood
Ratio
     
6.5830        
1        
0.0103

Score
                      
6.4209        
1        
0.0113

Wald                       
6.0756         1        
0.0137

                      
Residual Chi-Square Test

Chi-Square           
DF     
Pr > ChiSq

     
4.9818             
2        
0.0828

Step 2
. Effect drive
entered:

                    
Model Convergence Status

                  Convergence
criterion (GCONV=1E-8)
satisfied.                      

                        
Model Fit Statistics

Criterion      Intercept Only   
 Intercept and Covariates

AIC                  63.827                    
56.287

SC
                  
65.633                     61.707

-2 Log
L
          
61.827                    
50.287

 

                
Testing Global Null Hypothesis: BETA=0

Test                       
Chi-Square     DF  
Pr > ChiSq
Likelihood
Ratio
        
11.5391        2     
0.0031

Score
                         
10.5976       
2     
0.0050

Wald                           
8.5949        
2     
0.0136

                
Residual Chi-Square Test

Chi-Square     DF   Pr
> ChiSq

    
0.1293        1    
0.7191

NOTE: No (additional) effects met the 0.05 significance level
for entry into the model.

(2)
给出了自变量进入模型的次序。先是截距项Step
0了，不管它。Step 1
vision第一个进入模型，附带了很多评估它对因变量预测能力的指标。-2 Log
L
 和Score
 用来检测自变量是否显著。-2 Log L
中的L就是Likelihood Ratio，它的p值是0.0103
，Score
  的p值是0.0113
，都小于0.05，故vision是一个很显著的解释变量。还有，AIC
(Akaike
Information Criterion)和SC
(Schwarz
Criterion)两个信息量标准用来比较不同的模型，它们数值越小，模型变现就越好，这个接下来我们看step2 drive
变量进入模型后的情况。
我们可以看到模型的表现变好了，因为这是AIC
和SC
的值变小了，-2
Log L
 和Score
对应的
p值也更小。

                (3)
  Summary of Forward Selection

Step    Effect
Entered       
DF    
 Number In
    
Score
Chi-Square       
Pr > ChiSq

1           
vision                   
1          
1                     
6.4209                       
0.0113

2           
drive                     1           
2                    
4.8680                         0.0274

 (3)
总结了我们模型使用的前向选择方法，包括自变量进入模型的次序，以及每个自变量的卡方值和p值。

                  
 (4)
Analysis of
Maximum Likelihood Estimates

Parameter   DF  
Estimate   
Standard Error  Wald
Chi-Square   Pr >
ChiSq

Intercept    
1     
0.1110
              
0.5457          
0.0414                
0.8389

vision          
1     
1.7137
              
0.7049          
5.9113                
0.0150

drive            1    
-1.5000
              
0.7037           4.5440                
0.0330

(4)
给出了模型参数的估计，据此可
以写出改回归方程的形式是log(odds of having an
accident)=log(p/(1-p))=0.1110+1.7137*vision-1.5000*drive。我们知道，odds=p/(1-
p),有p=odds/(1+odds)。假设有个哥们，视力没问题但没有受过驾车教育（vision=0,drive=0），代入方程，有
log(odds)=0.1110，再odds=exp(0.110)=1.1174
，p=1.1174/2.1174=0.5277，即我们说这人发生事故的概率为0.5277；又另一个，视力有问题同样没受过驾车教育（vision=1,drive=0），同样的步骤得
log(odds)=1.8249，odds=exp(1.8249)=6.2022
，p=6.2022/7.2022=0.8612，即这人发生事故的概率为0.8612（视力多重要）。

                 (5)
Odds Ratio Estimates

Effect  
    Point
Estimate
             
95% Wald Confidence Limits
vision
          
 5.550
                               1.394          
22.093

drive             0.223                              
0.056           
0.886

 (5)
是对比率Odds
Ratio的估计。这里对vision的odds ratio的点估是 5.550
 ，这个比率由odds(having
an accident|vision=1,drive=0)比odds(having an
accident|vision=0,drive=0)得来，根据(4)
的计算，就是6.2022
/
1.1174=5.550

.还有，对vision来说，95%
的置信区间不包括1，说明vision是一个非常显著的解释变量（比率的置信区间不包括1，就跟p值小于0.05一样是一个规则）。

                
(6)
Association of Predicted
Probabilities and Observed Responses
Percent
Concordant
            
67.2              
Somers' D 0.532
Percent
Discordant
             
14.0               Gamma
0.655

Percent
Tied                       
18.8               
Tau-a 0.269

Pairs                                    
500                  c
0.766

(6)
这个东西就有些复杂,大概说些预测概率与观测道德因变量间的关联性，我们看到一致性比率Percent
Concordant 为67.2 %，不一致性比率Percent
Discordant 为14.0%，说明预测值与观测值在现有水平上有较强的关联性,回归模型有很强的预测能力。

最后注意到，上面我们用过程步乙得出的输出结果没有age这个自变量，用过程步甲得出的输出有，这不是因为年龄在这个预测模型中不重要，而是上面以数值型面目出现的年龄在是否出现事故的两组人中分布不均匀，为了解决这个问题，把年龄分组就是，不多说了。