*! version 1.0.0  7-7-92  J. Hilbe  STB-11 sg16
* GLM - binomial:bernoulli distribution; logit link
* limited options include: tolerance, if/in, offset, exposure

program define glmlogis
   version 3.0
   local varlist "req ex"
   local options "`options' LTolerance(real .0001) Offset(string) Exposur(string)"
   local in "opt"
   local if "opt"
   parse "`*'"
   parse "`varlist'", parse(" ")
   qui {
       tempvar mu eta w z dev oldev LNEX touse 
       local y "`1'"
       mac shift

*  EXPOSURE: OFFSET: IF/IN
       if "`exposur'"~=""  {
          if "`offset'"~= "" { error 198 }
            _crcunab `exposur'
            local expostr "$S_1"
            gen double `LNEX' = ln(`expostr')
            local offset `LNEX'
            local offmiss "| `offset'==."
       }
       else if "`offset'"~="" {
          _crcunab `offset'
          local offset "$S_1"
          local offmiss "| `offset'==."
          local expostr "exp($S_1)"
       }
       gen byte `touse'=1 `if' `in'       
       replace `touse'=0 if `y'==. | `touse'==. `offmiss'
       
       if "`offset'"~="" {  
           local poffset "+`offset'"
           local moffset "-`offset'"
           local offopt "offset(`offset')"
       }   

*INITIALIZATION
       sum `y' if `touse'
       local nobs = _result(1)
       gen  `mu'     = (`y' +.5)/2  if `touse'        
       gen  `eta'    = log(`mu'/(1-`mu')) if `touse'
       local i = 1
       gen `w'=0 if `touse'
       gen `z'=0  if `touse'
       gen `dev'=1  if `touse'
       gen `oldev'=1  if `touse'
       local ddev = 1  
     

* LOCAL SCORING ALGORITHM: IRLS   : simple convergence procedure
       while (abs(`ddev') > `ltolera')   {      
           replace  `w' = `mu'* (1-`mu') if `touse' 
           replace  `z' = `eta' + (`y'-`mu')/`w' `moffset' if `touse'
           regress `z' `*' [iw=`w'], mse1 dof(100000) dep(`y')
           drop `eta'
           predict `eta' if `touse'
           replace `eta' = `eta' `poffset'
           replace `mu' = 1/(1+exp(-`eta'))
           replace `oldev' = `dev' 
           replace `dev' = cond(`y',log(`mu'),log(1-`mu'))
           replace `dev' = sum(`dev') 
           replace `dev' = -2*`dev'[_N]
           local ddev = `dev' - `oldev'
           noi di in gr "Iter `i' : Dev = " in ye `dev'
           local i = `i'+1
        }
    }
local df = `nobs'-_result(3)-1

* DISPLAY
di in gr _col(57) "No obs.    = "  in ye %9.0g `nobs'
di in gr _col(57) "Deviance   = "  in ye %9.4f `dev'
di in gr _col(57) "Prob>chi2  = "  in ye %9.4f chiprob(`df',`dev')
di in gr _col(57) "Dispersion = "  in ye %9.4f 1  /* not calculated */

di in gr "GLM - Binomial:Bernoulli (logit link)" 
noi regress , noheader 

end