/* DATE: 6/20/91   */
/*               LOGISTIC REGRESSION EXTENSIONS 
   extlogit.ado
   Wald statistics & partial correlations; model GOF statistics;
   Creation of diagnostic variables including Pearson residual,
         deviance, hat diagonal, modified Cook, deltaChi, deltaDev

   Joseph Hilbe, STB 
*/
cap program drop extlogit
program define extlogit
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 "-> extlogit %_varlist %_if %_in" _n
   }
   else  {
      mac def _options "*"
      parse "%_*"
   }
   parse "%_varlist",parse(" ")
   quietly logit %_*  %_if %_in
     cap drop _num
     cap drop _depvar
     cap drop logpred
     cap drop logindex
     cap drop sepred
     cap drop resid
     cap drop stresid
     cap drop _good
     cap drop llo
     cap drop _one
     cap drop hat
     cap drop dev
     cap drop cmod
     cap drop deltax
     cap drop deltad
           
   quietly {
      gen _depvar=%_1 %_if %_in
      gen _one=sum(cond(_depvar,1,0))
        replace _one=_one[_N]
          gen _num=[_n-1]+1 
      gen llo=(_one*(log(_one)))+((_num[_N]-_one)*(log(_num[_N]-_one)))-(_num[_N]*(log(_num[_N])))
      predict logpred
      predict logindex, index
      predict sepred,stdp
      gen resid=_depvar-logpred
      gen _good=sum((resid^2)/(logpred*(1-logpred)))
        replace _good=_good[_N]
      gen stresid=((resid)/(sqrt(logpred*(1-logpred))))
      gen hat=(sepred^2)*(1-logpred)*logpred
        sort logpred
        by logpred: replace hat=sum(hat)
        by logpred: replace hat=hat[_N]
      gen dev=sqrt(-2*log(logpred)) if _depvar==1
        replace dev=-sqrt(-2*log(1-logpred)) if _depvar==0
      gen cmod   = (stresid^2*hat)/(1-hat)^2
      gen deltad = dev^2 + ((stresid^2*hat)/(1-hat))
      gen deltax = ((stresid^2*hat)/(1-hat))/hat
        drop _depvar
        drop _num
     }
   cap drop _cof
   cap drop _secof
   cap drop _wald
   cap drop _parcor
   
   quietly gen _cof = .
   quietly gen _secof = .
   quietly gen _wald = .
   quietly gen _parcor = .
   mac def _i= 0
   mac shift

   di in gr "                                           "
   di in gr "                                           "
   di in wh "MODEL GOODNESS OF FIT STATISTICS"
   di in gr "** ChiSq>.05 fails to reject hypothesis that model fits"
   di in gr "                                           "
   di in gr "H-L Goodness of Fit (Mod) => " %9.4f in ye _good
   di in gr "   ChiSq sign. (df:N-k)   => " %9.4f in ye chiprob((_result(1)-_result(3)-1),_good)
   di in gr "-2 LL(Model)              => " %9.4f in ye 2*abs(_result(2))
   di in gr "   ChiSq sign. (df:N-k)   => " %9.4f in ye chiprob((_result(1)-_result(3)-1),2*abs(_result(2)))
   di in gr "-2 LL(Intercept)          => " %9.4f in ye 2*abs(llo)
   di in gr "   ChiSq sign. (df:N-k)   => " %9.4f in ye chiprob((_result(1)-_result(3)-1),2*abs(llo))
   drop _good
   drop _one
   di in gr "                                           "
   di in gr "                                           "
   di in wh "WALD STATISTICS & PARTIAL CORRELATIONS"
   di in gr _n _col(1) "No  Var         Wald   Prob(Chi)    Partial Corr"
   di in gr _dup(48) "="
   while "%_1"!=""  {
      mac def _i = %_i + 1
      quietly {
         replace _cof = _b[%_1] in %_i
         cap test %_1
         replace _secof = abs(_b[%_1])/sqrt(_result(6)) in %_i
         replace _wald  = (_cof/_secof)^2 in %_i
         replace _parcor = sqrt((_wald-2)/(2*abs(llo)))
         replace _parcor = 0 if _wald<2  
         di in gr %_i %9.4f in gr _col(5) "%_1" %9.4f in ye _col(13) _wald[%_i] %9.4f _col(22) chiprob(1,_wald[%_i]) %9.4f _col(40) _parcor[%_i] %9.4f
         mac shift
        }
}
   drop _cof
   drop _wald
   drop _secof
   drop _parcor
   drop llo
   di in gr "                                           "
   di in gr "                                           "
   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 "    logpred  =  "in ye "Probability of success (1)"
   di in gr "    resid    =  "in ye "Residual"
   di in gr "    stresid  =  "in ye "Standardized residual (Pearson)"
   di in gr "    hat      =  "in ye "Hat matrix diagonal"
   di in gr "    dev      =  "in ye "Deviance"
   di in gr "    cmod     =  "in ye "Influence on est parameter values"
   di in gr "    deltad   =  "in ye "Change in Deviance"
   di in gr "    deltax   =  "in ye "Change in Pearson chi-square"
   
end