/* 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 */