*!version 1.2.1 : 2-3-95  Joseph Hilbe -  revised: 6-7-96, 7-1-96
* cloglog regression, binomial/bernoulli : OIM-based
 
program define cloglog
  version 4.0
  local options "Level(integer $S_level) EForm"
  if "`1'"!="" & substr("`1'",1,1)!="," {
  local varlist "req ex"
  #delimit ;
  local options "`options' LTolera(real -1) ITerate(integer 0) Offset(string)
     noCONStant NOLog Group RESidual Expec(real 2) DISP SCale(string)";
  #delimit cr
  local weight "fweight aweight"
  local in "opt"
  local if "opt"
  parse "`*'"
  parse "`varlist'",parse(" ")
  tempvar typ
  qui gen `typ'=.
  local vtype : type `typ'
  drop `typ'
  set type double
  qui  {
    tempvar w wt mu eta m ym touse u ll dev oldev n z ofs cdisp chi g1 g11
    tempvar v v1 wo ll0
  local y "`1'"
  mac shift
  if "`group'" != "" {
    local m = "`1'"
    mac shift
  }
  else { local m = 1 }
  global S_E_cmd
  count
  local nobs = _result(1)
  if `ltolera'<0  { local ltolera 1e-6 }
  if `iterate'<=0 { local iterate 50 }
  if `level'<10 | `level'>99  { local level 95 }
  if "`constan'"  !=""  { local cons "nocons" }
  if "`eform'"    !=""  { local eform "eform(EF)" }
  if "`offset'"   !=""  { gen `ofs' = `offset' }
     else { gen `ofs' = 0 }
  gen byte `touse'=1 `in' `if'
  replace `touse'=0 if `y'==. | `touse'==.
  if "`exp'" !="" {
     gen float `wt' `exp' if `touse'
     replace `wt'=. if `wt'<=0
     replace `touse'=0 if `wt'==.
     if ("`weight'"=="aweight")  {
        sum `wt'
        replace `wt' = `wt'/_result(3)
     }
  }
  else gen float `wt' = `touse' if `touse'
  qui reg `y' `*' if `touse'
  mac list
  local nobs=_result(1)
 
  replace `ofs'=`ofs' if `touse'
  gen `ym'=`y'/`m' if `touse'
  gen `w' = 1 if `touse'
  gen `ll' = 0 if `touse'
  gen `ll0'=0 if `touse'
 
* BINOMIAL INITIALIZATION
  local i = 1
  local ddev = 1
  gen `z'    = 0 if `touse'
  gen `dev'  = 0 if `touse'
  gen `oldev'= 1 if `touse'
  gen `u'    = 0 if `touse'
  gen `g1'   = 0 if `touse'
  gen `g11'  = 0 if `touse'
  gen `v'    = 0 if `touse'
  gen `v1'   = 0 if `touse'
  gen `wo'   = 0 if `touse'
 
* MODULE TO CALCULATE LL0
  gen `mu' =`m'* (`y'+0.5)/(`m'+1) if `touse'
  gen `eta' = ln(-ln(1-`mu'/`m')) if `touse'
  while ((abs(`ddev')>`ltolera' & `i'<`iterate') | (`i'==0))  {
      if `i' <= `expec'  {
      replace `u' = -(`y'-`mu')/((`m'-`mu')*ln((`m'-`mu')/`m'))
      replace `w'= ((`m'-`mu')*(ln((`m'-`mu')/`m')))^2 /(`mu'*(1-`mu'/`m'))
      replace `z' = `eta' + `u' - `ofs'
      regress `z' [iw=`w'*`wt'], mse1 dep(`y') `cons' dof(100000)
    }
    else {
      replace `u' = -(`y'-`mu')/((`m'-`mu')*ln((`m'-`mu')/`m'))
      replace `w'= ((`m'-`mu')*(ln((`m'-`mu')/`m')))^2 /(`mu'*(1-`mu'/`m'))
      replace `g1' = -1/((`m'-`mu')*ln((`m'-`mu')/`m'))
      replace `g11' = (-1-ln((`m'-`mu')/`m'))/((`mu'-`m')^2* /*
         */ (ln((`m'-`mu')/`m'))^2)
      replace `v' = `mu'*((`m'-`mu')/`m')
      replace `v1' = (`m'-2*`mu')/`m'
      replace `wo' = `w'+(`y'-`mu')*(`v'*`g11'+`v1'*`g1')/( (`v'^2 * `g1'^3))
      replace `z' = `eta' + ((1/`wo')*`w'*`u') -`ofs'
      regress `z' [iw=`wo'*`wt'], mse1 dep(`y') `cons' dof(100000)
    }
    cap drop `eta'
    predict `eta'
    replace `eta' = `eta' + `ofs'
    replace `mu' = `m'*(1-exp(-exp(`eta')))  /* inverse link */
    replace `oldev' = `dev'
    #delimit ;
    replace `dev' = `m'*ln((`m'-`y')/(`m'-`mu')) if (`y'==0);
    replace `dev' = `y'*ln(`y'/`mu') if (`y'==`m');
    replace `dev' = `y'*ln(`y'/`mu')+(`m'-`y')*ln((`m'-`y')/(`m'-`mu')) if
                    (`y'!=0 & `y'!=`m');
    #delimit cr
    replace `dev' = sum(`dev'*`wt')
    replace `dev' = 2 * `dev'[_N]
    replace `ll0'  = `y'*ln(`mu'/`m')+(`m'-`y')*ln(1-(`mu'/`m'))
    replace `ll0'  = sum(`ll0'*`wt')
    replace `ll0'  = `ll0'[_N]
    local ddev = `dev' - `oldev'
    local i = `i' + 1
}
 
* MAIN MODULE
  replace `mu' =`m'* (`y'+0.5)/(`m'+1) if `touse'
  replace `eta' = ln(-ln(1-`mu'/`m')) if `touse'
  local i=1
  local ddev=1
  while ((abs(`ddev')>`ltolera' & `i'<`iterate') | (`i'==0))  {
      if `i' <= `expec'  {
      replace `u' = -(`y'-`mu')/((`m'-`mu')*ln((`m'-`mu')/`m'))
      replace `w'= ((`m'-`mu')*(ln((`m'-`mu')/`m')))^2 /(`mu'*(1-`mu'/`m'))
      replace `z' = `eta' + `u' - `ofs'
      regress `z' `*' [iw=`w'*`wt'], mse1 dep(`y') `cons' dof(100000)
    }
    else {
      replace `u' = -(`y'-`mu')/((`m'-`mu')*ln((`m'-`mu')/`m'))
      replace `w'= ((`m'-`mu')*(ln((`m'-`mu')/`m')))^2 /(`mu'*(1-`mu'/`m'))
      replace `g1' = -1/((`m'-`mu')*ln((`m'-`mu')/`m'))
      replace `g11' = (-1-ln((`m'-`mu')/`m'))/((`mu'-`m')^2* /*
         */ (ln((`m'-`mu')/`m'))^2)
      replace `v' = `mu'*((`m'-`mu')/`m')
      replace `v1' = (`m'-2*`mu')/`m'
      replace `wo' = `w'+(`y'-`mu')*(`v'*`g11'+`v1'*`g1')/( (`v'^2 * `g1'^3))
      replace `z' = `eta' + ((1/`wo')*`w'*`u') -`ofs'
      regress `z' `*' [iw=`wo'*`wt'], mse1 dep(`y') `cons' dof(100000)
    }
    cap drop `eta'
    predict `eta'
    replace `eta' = `eta' + `ofs'
    replace `mu' = `m'*(1-exp(-exp(`eta')))  /* inverse link */
    replace `oldev' = `dev'
    #delimit ;
    replace `dev' = `m'*ln((`m'-`y')/(`m'-`mu')) if (`y'==0);
    replace `dev' = `y'*ln(`y'/`mu') if (`y'==`m');
    replace `dev' = `y'*ln(`y'/`mu')+(`m'-`y')*ln((`m'-`y')/(`m'-`mu')) if
                    (`y'!=0 & `y'!=`m');
    #delimit cr
    replace `dev' = sum(`dev'*`wt')
    replace `dev' = 2 * `dev'[_N]
    replace `ll'  = `y'*ln(`mu'/`m')+(`m'-`y')*ln(1-(`mu'/`m'))
    replace `ll'  = sum(`ll'*`wt')
    replace `ll'  = `ll'[_N]
    if "`nolog'"==""  {
       noi di in gr "Iteration  " `i'-1  in gr ": Log Likelihood = " /*
         */ in ye %12.6f `ll'
    }
    local ddev = `dev' - `oldev'
    local i = `i' + 1
}
local pdev = `dev'
gen `chi' = (`y'-`mu')^2 / (`mu'*(1-`mu'/`m')) if `touse'
replace `chi' = sum(`chi'*`wt')
replace `chi' = `chi'[_N]
 
local pred = _result(3)
local df = `nobs'-`pred'-("`cons'"=="")
if ("`weight'" == "fweight")  {
   tempvar fwt
   egen `fwt' = sum(`wt')
   local fobs = `fwt'
   local df = `fobs' - `pred'-1
   if "`constan'"!=""  {
      local df = `fobs - `pred'
   }
 }
local df = `nobs'-`pred'-("`cons'"=="")
gen `cdisp' = `chi'/`df' if `touse'
local pdisp = `cdisp'
local pr2 = 1-(`ll'[_N]/`ll0'[_N])
local ch2 = -2*(`ll0'[_N]-`ll'[_N])
local ll = `ll'[_N]
local ll0 = `ll0'[_N]
local cdisp = `cdisp'
}
 
qui regress,  noheader `eform' level(`level')
tempname b V
mat `b' = get(_b)
mat `V' = get(VCE)
mat post `b' `V', depname(`y') obs(`nobs')
global S_E_depv "`y'"
global S_E_ll `ll'
global S_E_ll0 `ll0'
global S_E_pr2 `pr2'
global S_E_chi2 `ch2'
global S_E_nobs `nobs'
global S_E_pred `pred'
global S_E_cdis `cdisp'
global S_E_ddis `dev'/`df'
global S_E_cmd "cloglog"
 
}
 
else  {
     if "$S_E_cmd"!="cloglog" {
          error 301
     }
     parse "`*'"
}
if `level'<10 | `level'>99 {
    local level 95
}
if "`eform'" !=""  { local eform "eform(EF)" }
 
noi di in gr _n "Cloglog Estimates"  _col(56) "Number of obs = " /*
   */ in ye %7.0g $S_E_nobs
noi di in gr _col(56) "Chi2(" in ye $S_E_pred in gr ") " _col(69) " = " /*
   */ in ye %7.2f $S_E_chi2
noi di in gr _col(56) "Prob>Chi2     = " /*
   */ in ye %7.4f chiprob($S_E_pred,$S_E_chi2)
noi di in gr "Log Likelihood = " %9.3g /*
   */ in ye $S_E_ll in gr _col(56) "Pseudo R2     = " in ye %7.4f $S_E_pr2 _n
 
mat mlout, level(`level') `eform'
if "`disp'"!="" {
    noi di in gr "chi2 dispersion     = " in ye %9.4f $S_E_cdis
    noi di in gr "deviance dispersion = " in ye %9.4f $S_E_ddis
}
 
if"`residual'"!="" {
  set type float
  qui {
  tempvar yhat stdp Hat Presid Dresid Aresid Spear Sdev pp
  gen `Dresid'=1
  gen `yhat' = `mu'
  gen `pp' = `yhat'/`m'
  predict `stdp',stdp
  gen `Hat'=(`stdp'^2 * `yhat'*(1-`pp'))
  gen `Presid' = (`y'-`yhat')/sqrt(`yhat'*(1-`pp'))
  replace `Dresid' = `m'*log((`m'-`y')/(`m'-`yhat')) if (`y'==0)
  replace `Dresid' = `y'*log(`y'/`yhat') if (`y'==`m')
  #delimit ;
  replace `Dresid' = `y'*log(`y'/`yhat') + (`m'-`y') *
          log((`m'-`y')/(`m'-`yhat')) if (`y'~=0 & `y'~=`m');
  #delimit cr
  replace `Dresid' = sqrt(2*`Dresid')
  replace `Dresid' = -`Dresid' if (`y'-`yhat')<0
  gen `Aresid' = (2.0533902*ibeta(2/3,2/3,(`y'/`m')))- /*
     */ (2.0533902*ibeta(2/3,2/3,(`pp')))
  replace `Aresid' = `Aresid'/(((`pp'*(1-`pp'))^(1/6))*sqrt((1-`Hat')/`m'))
  gen `Spear' = `Presid'/sqrt(1-`Hat')
  gen `Sdev'  = `Dresid'/sqrt(1-`Hat')
  gen _mu = `mu'
  gen _lp = `eta'
  gen _Pear = `Presid'
  gen _Dev  = `Dresid'
  gen _Hat  = `Hat'
  gen _Anscom = `Aresid'
  gen _StandP = `Spear'
  gen _StandD = `Sdev'
  noi di in gr _n "Variables created:  "
  noi di in ye " _mu (fit),                  _lp (linear predictor)"
  noi di in ye " _Pear (Pearson resid),      _Dev (deviance resid)"
  noi di in ye " _Anscom (Anscombe resid),   _Hat (hat matrix diag.)"
  noi di in ye " _StandP (standard Pearson), _StandD (standard deviance)"
  set type double
  }
}
set type `vtype'
end