Example 9.8 Logistic and Probit Regression for Binary Response Models

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

$\displaystyle  \Pr (Y=1)  $
$\displaystyle  =  $
$\displaystyle  F(\mb {X} \beta )  $
$\displaystyle \Pr (Y=0)  $
$\displaystyle  =  $
$\displaystyle  1 - F(\mb {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 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, 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 9.8.1.

Output 9.8.1: Logistic and Probit Regression: Results

iter loglik btransp  
1 -268.248 0 0 0

iter loglik btransp  
2 -76.29481 -2.159406 0.0138784 0.0037327

iter loglik btransp  
3 -53.38033 -3.53344 0.0363154 0.0119734

iter loglik btransp  
4 -48.34609 -4.748899 0.0640013 0.0299201

iter loglik btransp  
5 -47.69191 -5.413817 0.0790272 0.04982

iter loglik btransp  
6 -47.67283 -5.553931 0.0819276 0.0564395

iter loglik btransp  
7 -47.67281 -5.55916 0.0820307 0.0567708

iter loglik btransp  
8 -47.67281 -5.559166 0.0820308 0.0567713

  chisq df prob
Likelihood Ratio, Intercept-only Model 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 btransp  
1 -268.248 0 0 0

iter loglik btransp  
2 -71.71043 -1.353207 0.008697 0.0023391

iter loglik btransp  
3 -51.64122 -2.053504 0.0202739 0.0073888

iter loglik btransp  
4 -47.88947 -2.581302 0.032626 0.018503

iter loglik btransp  
5 -47.48924 -2.838938 0.0387625 0.0309099

iter loglik btransp  
6 -47.47997 -2.890129 0.0398894 0.0356507

iter loglik btransp  
7 -47.47995 -2.89327 0.0399529 0.0362166

iter loglik btransp  
8 -47.47995 -2.893408 0.0399553 0.0362518

iter loglik btransp  
9 -47.47995 -2.893415 0.0399554 0.0362537

iter loglik btransp  
10 -47.47995 -2.893415 0.0399555 0.0362538

iter loglik btransp  
11 -47.47995 -2.893415 0.0399555 0.0362538

  chisq df prob
Likelihood Ratio, Intercept-only Model 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