* DATE: 10/11/91   as sqv1.3 
/*               LOGISTIC REGRESSION INFLUENCE STATISTICS 
inflogit.ado
   Creation of diagnostic variables including: 
   Pearson residual, Standardized Pearson residual,
   deviance, hat diagonal, Cook, deltaX, deltaD, and deltaB
   < Pearson X^2 and Deviance statistics also provided >
   
   Program adjusts for m-asymptotic covariate pattern distribution
   
   Joseph Hilbe, STB/CRC 
*/
cap program drop inflogit
program define inflogit
version 2.1
   mac def _varlist "req ex min(2)"
   mac def _if "opt"
   mac def _in "opt"
   if "%_*"=="" {
      parse "%LOGIODDS %LOGIFIN"
      di in gr "-> inflogit %_varlist %_if %_in" _n
   }
   else  {
      mac def _options "*"
      parse "%_*"
   }
   parse "%_varlist",parse(" ")
   qui count
   mac def _origobs=_result(1)
   capture  {
      if "%_if`%_in"~="" {
         keep %_if %_in
      }
   qui dropmiss %_varlist
   count
   mac def _obs=_result(1)
   mac def _lhs "%_1"
   mac shift
   mac def _rhs "%_*"

   quietly logit %_lhs %_rhs %_if %_in
     cap drop _depvar1
     cap drop logindex
     cap drop sepred
     cap drop presid
     cap drop hat
     cap drop dev
     cap drop cook
     cap drop deltax
     cap drop deltad
     cap drop stpresid
     cap drop _m
     cap drop _y
     cap drop pred
     cap drop mpred
     cap drop deltab
     cap drop _count
     cap drop pearson
     cap drop deviance
     cap drop _i
     
   quietly {
      gen _depvar1=%_lhs %_if %_in
      predict pred
      gen mpred=pred   /* will later take value of pred * _m  */
      gen _y=0         /* numb of positive (1) obser per c.p. */
      gen _m=0         /* numb of obser per covariate pattern */
      sort %_rhs
           by %_rhs:  replace mpred=sum(mpred)
           by %_rhs:  replace mpred=mpred[_N]
           by %_rhs:  replace _m=[_N] 
           by %_rhs:  replace _y=sum(cond(_depvar1,1,0))
           by %_rhs:  replace _y=_y[_N]
           by %_rhs:  gen _count=1 if _n==1 
           replace _count=sum(_count)
           replace _count=_count[_N]        /* numb of c.p.'s */
           
      predict logindex, index
      predict sepred,stdp
      gen presid=(_y-mpred)/sqrt((mpred*(1-pred)))
      gen hat=(sepred^2)*(1-pred)*(mpred)
      gen stpresid=presid/sqrt(1-hat)
      #delimit ;
      gen dev= sqrt((2*((_y*log(_y/mpred))+((_m-_y) * 
                  log((_m-_y)/(_m*(1-pred))))))) if _m>1;
      #delimit cr
      replace dev= sqrt(2*(_m*log(_m/(_m*(1-pred))))) if _m>1 & _y==0
      replace dev=-dev if (_y-mpred<0) & _m>1
      replace dev=sqrt(-2*log(pred)) if _depvar1==1 & _m==1
      replace dev=-sqrt(-2*log(1-pred)) if _depvar1==0 & _m==1
      
      gen cook   = (presid^2*hat)/(1-hat)
      gen deltad = dev^2 + ((presid^2*hat)/(1-hat))
      gen deltax = ((presid^2*hat)/(1-hat))/hat
      gen deltab = (presid^2*hat)/(1-hat)^2
      drop _depvar1
      }

by %_rhs: gen pearson=presid^2 if _n==1
replace pearson=sum(pearson)
replace pearson=pearson[_N]
by %_rhs: gen deviance=dev^2 if _n==1
replace deviance=sum(deviance)
replace deviance=deviance[_N]

while "%_1"~="" {    /* calculate numb of predictors */
mac def _i=%_i+1
  mac shift
}
di _n(2)
di in gr "Number of Predictors         = " in ye %9.0f %_i
di in gr "Number of Non-Missing Obs    = " in ye %9.0f %_obs
di in gr "Number of Covariate Patterns = " in ye %9.0f _count
di in gr "Pearson X^2 Statistic        = " in ye %14.4f pearson
#delimit ;
di in gr "  P>chi2(" _count-(%_i+1) ")" _col(29)" = " 
     in ye %14.4f   chiprob(_count-(%_i+1), pearson);
di in gr "Deviance                     = " in ye %14.4f deviance;
di in gr "  P>chi2(" _count-(%_i+1) ")" _col(29)" = " 
     in ye %14.4f   chiprob(_count-(%_i+1), deviance);
#delimit cr

di _n(2)
   di in wh "Additional diagnostic variables created..."
   di in gr "    logindex =  "in ye "Logit; Index value"
   di in gr "    sepred   =  "in ye "Standard error of index"
   di in gr "    pred     =  "in ye "Probability of success (1)"
   di in gr "    mpred    =  "in ye "Prob of covariate pattern success"
   di in gr "    presid   =  "in ye "Pearson Residual"
   di in gr "    stpesid  =  "in ye "Standardized Pearson Residual"
   di in gr "    hat      =  "in ye "Hat matrix diagonal"
   di in gr "    dev      =  "in ye "Deviance"
   di in gr "    cook     =  "in ye "Cook's distance"
   di in gr "    deltad   =  "in ye "Change in Deviance"
   di in gr "    deltax   =  "in ye "Change in Pearson chi-square"
   di in gr "    deltab   =  "in ye "Difference in coefficient due to"
   di in ye "                deletion of observation and others"
   di in ye "                sharing same covariate pattern"

cap drop _count   
cap drop pearson
cap drop deviance
cap drop _y
cap drop _m
cap drop _i
mac def _rc=_rc
if %_rc==0 {exit}
exit %_rc
end


program define dropmiss
version 2.1
  mac def _varlist "req ex min(2)"
  parse "%_*"
  parse "%_varlist",parse(" ")
  while "%_1"~="" {
     drop if %_1==.
     mac shift
  }
end