Chapter Contents
Chapter Contents
Previous
Previous
Next
Next
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

\Pr(Y=1) & = & F(X \beta) \\Pr(Y=0) & = & 1 - F(X \beta)
where F 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; Jennrich 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 w/p(1 - p), where w has the counts or other weights. The following code 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 below.

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


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

Chapter Contents
Chapter Contents
Previous
Previous
Next
Next
Top
Top

Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.