/* DATE: 4\24\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