*! version 2.0.1  8-10-93  J. Hilbe             STB-18: sg16.5
* Generalized linear models 

program define glm
   version 3.0
   local options "LEvel(integer $S_level) EForm"
   if "`1'"~="" & substr("`1'",1,1)~="," {
     local varlist "req ex"
     local options "`options' LTolerance(real -1) ITerate(integer 0) Fam(string) Link(string) Pow(real 1.0) K(real 1.0) DISP(real 1.0) SCale(string) Offset(string) EXposure(string) Group"
     local in "opt"
     local if "opt"
     local weight "fweight aweight"
   }
   parse "`*'"
   parse "`varlist'", parse(" ")

   if (("`scale'"=="") & ("`fam'"=="bin" | "`fam'"=="poi" | "`fam'"=="nb")) { local scale "n" } 
   if ("`scale'"=="") { local scale "d" }
   else { local scale = substr("`scale'",1,1) }

   if `level'<10 | `level'>99  {
     local level 95
   }
   if "`link'"=="" { local link default }
   if "`eform'"~="" {
     if ("`link'"=="l") { local eform "eform(Odds Ratio)" }
     else if ("link'"=="p" | "`link'"=="c") { local eform "eform(ExpB)" }
     else if ("`fam'"=="poi") { local eform "eform(IRR)" }
     else { local eform "eform(e^coef)" }
   }

   if ("`1'"=="") {
     if ("$S_E_cmd"~="glm") { error 301 }
     di in gr "$S_1"
     di in gr "$S_2" 
     regress, noheader level(`level') `eform'
     exit
   }
   
   qui {
       tempvar mu eta w z dev oldev chi2 wt m touse pv r 
       local y "`1'"
       mac shift
       local m=1
       
* OFFSET/EXPOSURE OPTION
       tempvar LNEX
       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 "`exp'"~="" {
              gen float `wt' `exp' if `touse'
              replace `wt'=. if `wt'<=0
              replace `touse'=0 if `wt'==.
              if ("`weight'"=="aweight") {
                  summ `wt'
                  replace `wt' = `wt' / _result(3)
              }
        }
        else gen float `wt' = `touse' if `touse'

        if "`offset'"~="" {
             replace `touse'=0 if `offset'==.
             replace `wt'=. if `offset'==.
             local poffset "+`offset'"
             local moffset "-`offset'"
             local offopt "offset(`offset')"
        }

* GROUPED DATA SET OPTION
        if "`group'"~="" {
              tempvar ym
              local m "`1'"
              mac shift
              gen `ym'=`y'/`m' if `touse'
        }

*INITIALIZATION
       qui regress `y' `*'  `if' `in'
       mac list
       local nobs   = _result(1)
       if `ltolera'<0  { local ltolera 0.0001 }  /* DEFAULT TOLERANCE  */
       if `iterate'<=0 { local iterate 50 }      /* DEFAULT ITERATIONS */
       egen `mu'    = mean(`y')                  /* INITIALIZE MU  */
       replace `mu' = (`y' +`mu')/(`m'+1)          
       
       local i    = 1
       gen `w'    = 0 if `touse' 
       gen `z'    = 0 if `touse'  
       gen `dev'  = 0 if `touse'
       gen `oldev'= 1 if `touse'  
       local ddev = 1  

* SET LINK
       if "`fam'"=="bin" & "`link'"~="pow"  { 
            replace `mu'=(`y'+.5)/(`m'+1)
            if "`link'"=="l" { gen `eta' = log(`mu'/(1-`mu')) if `touse' }
            if "`link'"=="p" { gen `eta' = log(`mu'/(1-`mu')) if `touse' }
            if "`link'"=="c" { gen `eta' = log(-log(1-`mu')) if `touse' }
       }
       if "`link'"~="pow" {
            if "`fam'"=="gau"  { gen `eta' = `mu' if `touse' }
            if "`fam'"=="poi"  { gen `eta' = log(`mu') if `touse' }
            if "`fam'"=="gam"  { gen `eta' = 1/`mu' if `touse' }
            if "`fam'"=="ivg"  { gen `eta' = 1/`mu'^2 if `touse' }
            if "`fam'"=="nb"   { gen `eta' = log(`mu'/(`mu'+`k')) if `touse' }       
*            if "`fam'"=="nb"   { gen `eta' = - log(1+ 1/(`mu'*`k')) if `touse' }       
       }     
       if "`link'"=="pow" {  
            if `pow'==0 { gen `eta' =  log(`mu') if `touse' }
            else { gen `eta' = `mu'^`pow' if `touse' }
       }

* LOCAL SCORING ALGORITHM: IRLS   
       while ((abs(`ddev') > `ltolera'& `i'<`iterate')| (`i'== 0))   {      
           
           * CALL FOR VARIANCE FUNCTION
           _glmwgt `fam' "`link'" `eta' `mu' `w' `pow' `k' `touse' `disp' 
           
           * CALCULATE DUMMY RESPONSE
           if ("`fam'"=="gam"  & "`link'"~="pow")  {
              replace  `z' = `eta' - (`y'-`mu')/`w' `moffset' if `touse' 
           }
           else if ("`fam'"=="ivg" & "`link'"~="pow") {
              replace `z' = `eta'-(2*(`y'-`mu'))/`w' `moffset' if `touse'
           }  
           else if ("`fam'"=="bin" & "`group'"~="") {
              replace `z' = `eta' + (`ym'-`mu')/`w' `moffset' if `touse' 
           }
           else  { replace `z' = `eta' + (`y'-`mu')/`w' `moffset' if `touse'}
           
           * ADJUST WEIGHTING
           if ("`link'"=="pow" & "`fam'"=="gau") { replace `w'=`w'^2 }
           if (("`link'"=="pow"|"`link'"=="p"|"`link'"=="c")&"`fam'"=="bin"){ 
              replace `w'=`w'^2/(`mu'*(1-`mu')) 
           }
           if ("`link'"=="pow") { 
             if "`fam'"=="poi" { replace `w'=`w'^2/`mu' }
             if "`fam'"=="gam" { replace `w'=`w'^2/`mu'^2 }
             if "`fam'"=="ivg" { replace `w'=`w'^2/`mu'^3 }
             if "`fam'"=="nb" { replace `w'=`w'^2/(`mu'+`k'*`mu'^2)}
           }
           if "`group'"~="" { replace `w' = `w'*`m' }

           * WEIGHTED REGRESSION
           regress `z' `*' [iw=`w'* `wt'], mse1 dof(100000) dep(`y')
           
           * OBTAIN LINEAR PREDICTION
           cap drop `eta'
           predict `eta' if `touse' 
           replace `eta' = `eta' `poffset'
           
           * CALL TO INVERSE LINK
           mac def S_E_fam "`fam'"
           mac def S_E_link "`link'"
           _glmilnk `eta' `mu' `pow' `k'
           
           * DEVIANCE
           replace `oldev' = `dev' 
           if "`fam'"=="gau" {
                replace `dev' = (`y'-`mu')^2
                replace `dev' = sum(`dev'*`wt')
                replace `dev' = `dev'[_N]/`disp'
           }
           if "`fam'"=="bin" & "`group'" == "" {
                replace `dev'= cond(`y',log(`mu'),log(1-`mu'))
                replace `dev'= sum(`dev'*`wt')
                replace `dev'=-2*`dev'[_N]/`disp'
           }
           if "`fam'"=="bin" & "`group'" ~= "" {
                #delimit ;
                replace `dev' = (`m'-`y')*log((`m'-`y')/(`m'-(`m'*`mu'))) 
                                if (`y'==0);
                replace `dev' = `y'*log(`y'/(`m'*`mu')) if (`y'==`m');
                replace `dev' = `y'*log(`y'/(`m'*`mu')) +
                                 (`m'-`y')*log((`m'-`y')/(`m'-(`m'*`mu'))) 
                                 if (`y'~=0 & `y'~=`m');
                #delimit cr
                replace `dev' = sum(`dev'*`wt')
                replace `dev' = 2*`dev'[_N]/`disp'
           }
           if "`fam'" == "poi" { 
                replace `dev' = (`y'*log(`y'/`mu') - (`y'-`mu')) if `y'>0
                replace `dev' = `mu' if `y'==0
                replace `dev' = sum(`dev'*`wt')
                replace `dev' = 2*`dev'[_N]/`disp'
           }
           if "`fam'" == "gam" { replace `dev' = -log(`y'/`mu') + (`y'-`mu')/`mu' 
                replace `dev' = sum(`dev'*`wt')
                replace `dev' = 2*`dev'[_N]/`disp'
           }
           if "`fam'" == "ivg" { replace `dev' = (`y'-`mu')^2 / (`mu'^2*`y')  
                replace `dev' = sum(`dev'*`wt')
                replace `dev' = `dev'[_N]/`disp'
           }
           if "`fam'" == "nb"  {
                #delimit ;
                replace `dev' = `y'*log(`y'/`mu')-(1+`k'*`y')/`k' * 
                                 log((1+`k'*`y')/(1+`k'*`mu')) if `y'>0;
                replace `dev' = log(1+`k'*`mu')/`k' if `y'==0;
                #delimit cr
                replace `dev' = sum(`dev'*`wt')
                replace `dev' = 2*`dev'[_N]/`disp'
           }      
           local ddev    = `dev' - `oldev'
           noi di in gr "Iter `i' : Dev = " in ye %9.4f `dev'
           local i = `i'+1
        }
 

* CALCULATE CHI2 AND DISPERSIONS

if "`fam'"=="gau" { gen `chi2' = (`y'-`mu')^2 * `disp' if `touse'}
if "`fam'"=="bin" { gen `chi2' = (`y'-`mu')^2/(`mu'*(1-`mu')*`m'*`disp') if `touse' }
if "`fam'"=="poi" { gen `chi2' = (`y'-`mu')^2/(`mu'*`disp') if `touse' }
if "`fam'"=="gam" { gen `chi2' = (`y'-`mu')^2/(`mu'^2*`disp') if `touse'}
if "`fam'"=="ivg" { gen `chi2' = (`y'-`mu')^2/(`mu'^3*`disp') if `touse' }
if "`fam'"=="nb"  { gen `chi2' = (`y'-`mu')^2/((`mu'+`k'*`mu'*`mu')*`disp') if `touse' }
replace `chi2' = sum(`chi2'*`wt')
replace `chi2' = `chi2'[_N]
local df = `nobs'-_result(3)-1
local dispc = `chi2'/`df'
local dispd = `dev'[_N]/`df'
}


* DISPLAY
di in gr _col(57) "No obs.    = "  in ye %9.0g `nobs'
di in gr "Chi2       = " in ye %9.0g `chi2' in gr _col(57) "Deviance   = "  in ye %9.0g `dev'
di in gr "Prob>chi2  = " in ye %9.0g chiprob(`df',`chi2') in gr _col(57) "Prob>chi2  = "  in ye %9.0g chiprob(`df',`dev')
di in gr "Dispersion = " in ye %9.0g `dispc' in gr _col(57) "Dispersion = "  in ye %9.0g `dispd' _n

mac def S_E_fam "`fam'" 
mac def S_E_link "`link'"
_glmfl `pow'
di in gr "$S_1"
di in gr "$S_2"


if ("`scale'"=="d") {qui regress `z' `*' [iw=`w'*`wt'/`dispd'], mse1 dof(100000) dep(`y') } 
if ("`scale'"=="c") {qui regress `z' `*' [iw=`w'*`wt'/`dispc'], mse1 dof(100000) dep(`y') } 
noi regress, noheader level(`level') `eform'

if ("`scale'"=="d") { di in bl "(standard errors scaled by square root of deviance dispersion)" }
if ("`scale'"=="c") { di in bl "(standard errors scaled by square root of chi2 dispersion)" }
if (`disp'~=1) { di in bl "(Quasi-likelihood model)" }
local delta = 1
if ("`scale'"=="d") { local delta = `dispd' } 
if ("`scale'"=="c") { local delta = `dispc' }
  local dev = `dev'
  mac def S_E_cmd  "glm"
  mac def S_E_fam  "`fam'"
  mac def S_E_grp  "`group'"
  mac def S_E_m    "`m'"
  mac def S_E_depv "`y'"
  mac def S_E_link "`link'"
  mac def S_E_pow  "`pow'"
  mac def S_E_po   "`poffset'"
  mac def S_E_ef   "`eform'"
  mac def S_E_k    "`k'"  
  mac def S_E_disp "`disp'"
  mac def S_E_del  "`delta'"
  mac def S_E_dev  "`dev'"  
  mac def S_E_cdp  "`dispc'"  
 
end