*! version 1.0.0 : 9-3-93  J. Hilbe
* GLM: Poisson based censored exponential & Weibull regressions 

program define glmweib
  version 3.1
  local varlist "req ex"
  local options "`options' LTolera(real -1) ITerate(integer 0) Cen(string) Time Model(string)"
  local weight "fweight aweight"
  parse "`*'"
  parse "`varlist'",parse(" ")

  qui  {

tempvar mu eta z dev oldev cens logt alpha a1 devw oldevw ty z0 at beta t_adj
  
* EXTRACT TIME VARIABLE; LOG(T)
  local t "`1'"                /* survival time variable */
  mac shift
  gen `logt' = log(`t')        /* log time */
  gen `ty' = `t'               /* duplicate time variable */
  
* DEFINE STATUS OF CENSOR VARIABLE
  if ("`cen'"!="") { gen byte `cens'=`cen' }   /* censor variable */
     else { gen byte `cens' = 1 }
  count if `cens'==1
  local ncens=_result(1)          /* number of non-censored cases */
  count  
  local nobs=_result(1)

  if `ltolera'<0  { local ltolera 0.0001  }
  if `iterate'<=0 { local iterate 50 }
  
*INITIALIZATIONS
  local i  = 1
  gen `z'  = 0
  gen `dev' = 0
  gen `oldev' = 1
  gen `alpha' = 1
  gen `beta' = 1
  local ddev = 1
  local j = 1
  gen `oldevw' = 1  
  local ddevw = 1
  gen `devw' = 0
  gen `a1'=0
  gen `at'=0
  
  gen `z0' = sum(`cens'*`logt')
  replace `z0'=`z0'[_N]
  replace `z0'=2*`z0'
  
  if ("`model'"~="e") {    
       noi di in gr _col(11) "Deviance" _col(26) "Shape"
  }
  
* MAIN LOOP
  while ((abs(`ddevw') > `ltolera' & `j'<`iterate') | (`j'==0))  { 
     cap drop `mu'
     cap drop `eta'
     egen `mu' = mean(`cens')           /* Poisson initial values */
     gen `eta' = ln(`mu')               /* Poisson link */
     local i = 1                        /* Reset */
     local ddev = 1
     replace `dev' = 0
     
* POISSON LOOP
     while ((abs(`ddev') > `ltolera' & `i'<`iterate') | (`i'==0))  {
         replace `z' = (`eta'+(`cens'-`mu')/`mu') - (`logt'*`alpha') 
         regress `z' `*' [iw=`mu'], mse1  dep(`t')
         cap drop `eta'
         predict `eta' 
         replace `eta' = `eta' + (`logt'*`alpha')
         replace `mu' = exp(`eta')
         replace `oldev' = `dev'
         replace `dev' = (`cens'*log(`cens'/`mu') - (`cens'-`mu'))
         replace `dev' = sum(`dev')
         replace `dev' = 2*`dev'[_N]
         local ddev = `dev' - `oldev'
         if "`model'"=="e" {
            noi di in gr "Iter `i' :  " in ye %9.4f `dev'   
         }   
         local i = `i'+1
     }
         
* CALCULATION OF WEIBULL DEVIANCE
         replace `devw' = `cens'*log(`alpha'*`mu') - `mu'
         replace `devw' = sum(`devw')
         replace `devw' = `z0'- 2*`devw'[_N]
         local ddevw = `devw' - `oldevw'
         if "`model'"=="e" { noi di in gr _n "Generalized Deviance" }
         noi di in gr "Iter `j' : " in ye %9.4f `devw' " : "  %9.4f `alpha'
         if `j'==1 { local edev = `devw'}
         local j = `j'+1
         replace `oldevw' = `devw'

* RECALCULATION OF SHAPE PARAMETER
         replace `at' = `alpha'
         replace `a1' = sum((`mu'-`cens')*`logt'/`ncens')
         replace `a1' = `a1'[_N]                    
         replace `a1' = 1/`a1'
         replace `alpha'=(`alpha'+`a1')/2
         if "`model'"=="e" { local ddevw = 0.000001 }   /* expon: loop out */
 }

* LOG EXPECTED TIME FORM OPTION   
  if ("`time'"!="" & "`model'"=="e") {
  if (`cens'==1)  { replace `z' = -(`eta'+(`cens'-`mu')/`mu')+(`logt'*`at') }  
      else { replace `z' = 1-(`logt'*`at') }
    reg `z' `*' [iw=`mu'], mse1 dep(`t') 
  }
  else if ("`time'"!="" & "`model'"!="e") {
  if (`cens'==1)  { replace `z' =(-(`eta'+(`cens'-`mu')/`mu')+ /*
    */ (`logt'*`at'))/`at' }  
      else { replace `z' = 1-`eta'/log(`at') }
    reg `z' `*' [iw=`mu'*`at'^2], mse1 dep(`t') 
  }
local sescale=(1/`alpha')/sqrt(`nobs'+10)

}
* DISPLAY
if "`model'"~="e" { 
    di in gr _col(46)  "Observations         = " in ye %9.0f `nobs'
di in gr "Shape      = " in ye %9.4f `alpha'    in gr _col(46)  "Censored observ.     = " in ye %9.0f  `nobs'-`ncens' 
di in gr "Scale      = " in ye %9.4f 1/`alpha'  in gr _col(46)  "Exponential Deviance = " in ye %9.4f `edev'
di in gr "S.E. Scale = " in ye %9.4f `sescale'  in gr _col(46)  "Scaled Deviance      = " in ye %9.4f `devw'
    local devdiff = `edev'- `devw'
    di in gr _col(46)  "Prob > Chi2          = " in ye %9.4f chiprob(1,`devdiff')
}

if "`model'"=="e" { di in gr _n "Exponential regression" }
  else { di in gr  "Weibull regression" }
noi regress, noheader 
end