General Statistics Examples

Example 8.8: Logistic and Probit Regression for Binary Response Models

A binary response Y is fit to a linear model according to

where is some smooth probability distribution function. The normal and logistic distribution functions are supported. The method is maximum likelihood via iteratively reweighted least squares (described by Charnes, Frome, and Yu 1976; Jenrich and Moore 1975; and Nelder and Wedderburn 1972). The row scaling is done by the derivative of the distribution (density). The weighting is done by , where has the counts or other weights. The following program calculates logistic and probit regression for binary response models.

  
    /* routine for estimating binary response models             */ 
    /* y is the binary response, x are regressors,               */ 
    /* wgt are count weights,                                    */ 
    /* model is choice of logit probit,                          */ 
    /* parm has the names of the parameters                      */ 
  
    proc iml ; 
  
    start binest; 
       b=repeat(0,ncol(x),1); 
       oldb=b+1;                              /* starting values */ 
       do iter=1 to 20 while(max(abs(b-oldb))>1e-8); 
          oldb=b; 
          z=x*b; 
          run f; 
          loglik=sum(((y=1)#log(p) + (y=0)#log(1-p))#wgt); 
          btransp=b`; 
          print iter loglik btransp; 
          w=wgt/(p#(1-p)); 
          xx=f#x; 
          xpxi=inv(xx`*(w#xx)); 
          b=b + xpxi*(xx`*(w#(y-p))); 
       end; 
       p0=sum((y=1)#wgt)/sum(wgt);           /* average response */ 
       loglik0=sum(((y=1)#log(p0) + (y=0)#log(1-p0))#wgt); 
       chisq=(2#(loglik-loglik0)); 
       df=ncol(x)-1; 
       prob=1-probchi(chisq,df); 
       print , 
          'Likelihood Ratio with Intercept-only Model' chisq df prob,; 
       stderr=sqrt(vecdiag(xpxi)); 
       tratio=b/stderr; 
       print parm b stderr tratio,,; 
    finish;
 

  
       /*---routine to yield distribution function and density---*/ 
    start f; 
       if model='LOGIT' then 
       do; 
          p=1/(1+exp(-z)); 
          f=p#p#exp(-z); 
       end; 
       if model='PROBIT' then 
       do; 
          p=probnorm(z); 
          f=exp(-z#z/2)/sqrt(8*atan(1)); 
       end; 
    finish; 
  
         /* Ingot data from COX (1970, pp. 67-68)*/ 
    data={ 7 1.0 0 10, 14 1.0 0 31, 27 1.0 1 56, 51 1.0 3 13, 
           7 1.7 0 17, 14 1.7 0 43, 27 1.7 4 44, 51 1.7 0 1, 
           7 2.2 0 7, 14 2.2 2 33, 27 2.2 0 21, 51 2.2 0 1, 
           7 2.8 0 12, 14 2.8 0 31, 27 2.8 1 22, 
           7 4.0 0 9, 14 4.0 0 19, 27 4.0 1 16, 51 4.0 0 1}; 
    nready=data[,3]; 
    ntotal=data[,4]; 
    n=nrow(data); 
    x=repeat(1,n,1)||(data[,{1 2}]);    /* intercept, heat, soak */ 
    x=x//x;                                        /* regressors */ 
    y=repeat(1,n,1)//repeat(0,n,1);           /* binary response */ 
    wgt=nready//(ntotal-nready);                  /* row weights */ 
    parm={intercept, heat, soak};         /* names of regressors */ 
  
    model={logit}; 
    run binest;                               /* run logit model */ 
  
    model={probit}; 
    run binest;                              /* run probit model */
 

The results are shown in Output 8.8.1.

Output 8.8.1: Logistic and Probit Regression: Results

ITER LOGLIK 1 -268.248 BTRANSP 0 0 0 ITER LOGLIK 2 -76.29481 BTRANSP -2.159406 0.0138784 0.0037327 ITER LOGLIK 3 -53.38033 BTRANSP -3.53344 0.0363154 0.0119734 ITER LOGLIK 4 -48.34609 BTRANSP -4.748899 0.0640013 0.0299201 ITER LOGLIK 5 -47.69191 BTRANSP -5.413817 0.0790272 0.04982 ITER LOGLIK 6 -47.67283 BTRANSP -5.553931 0.0819276 0.0564395 ITER LOGLIK 7 -47.67281 BTRANSP -5.55916 0.0820307 0.0567708 ITER LOGLIK 8 -47.67281 BTRANSP -5.559166 0.0820308 0.0567713 CHISQ DF PROB 11.64282 2 0.0029634 PARM B STDERR TRATIO INTERCEPT -5.559166 1.1196947 -4.964895 HEAT 0.0820308 0.0237345 3.4561866 SOAK 0.0567713 0.3312131 0.1714042

Output 8.8.1: (continued)

ITER	LOGLIK
1	-268.248

BTRANSP
0	0	0

ITER	LOGLIK
2	-71.71043

BTRANSP
-1.353207	0.008697	0.0023391

ITER	LOGLIK
3	-51.64122

BTRANSP
-2.053504	0.0202739	0.0073888

ITER	LOGLIK
4	-47.88947

BTRANSP
-2.581302	0.032626	0.018503

ITER	LOGLIK
5	-47.48924

BTRANSP
-2.838938	0.0387625	0.0309099

ITER	LOGLIK
6	-47.47997

BTRANSP
-2.890129	0.0398894	0.0356507

ITER	LOGLIK
7	-47.47995

BTRANSP
-2.89327	0.0399529	0.0362166

ITER	LOGLIK
8	-47.47995

BTRANSP
-2.893408	0.0399553	0.0362518

ITER	LOGLIK
9	-47.47995

BTRANSP
-2.893415	0.0399554	0.0362537

ITER	LOGLIK
10	-47.47995

BTRANSP
-2.893415	0.0399555	0.0362538

ITER	LOGLIK
11	-47.47995

BTRANSP
-2.893415	0.0399555	0.0362538

CHISQ	DF	PROB
12.028543	2	0.0024436

PARM	B	STDERR	TRATIO
INTERCEPT	-2.893415	0.5006009	-5.779884
HEAT	0.0399555	0.0118466	3.3727357
SOAK	0.0362538	0.1467431	0.2470561