*! Version 1.0.5  9-30-92  J. Hilbe  STB-11 sg16
* Generalized Linear Models 

program define glm
  version 3.0
  local options "LEvel(integer $S_level) EForm"
  local varlist "req ex"
  local options "`options' LTolerance(real -1) ITerate(integer 0) Fam(string) Link(string) Group Stats Resid Offset(string) EXposure(string)"
  local in "opt"
  local if "opt"
  local weight "fweight"
  parse "`*'"
  parse "`varlist'",parse(" ")
  qui   {
       tempvar mu eta w z dev oldev m wt touse pv
       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'==.
        }
        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')"
        }
  
       if "`group'"~="" { 
              tempvar ym
              local m "`1'"
              mac shift
              gen `ym'=`y'/`m' if `touse' 
        }
        qui regress `y' `*' if `touse'


        count  if `touse'
        local nobs = _result(1)
        if `ltolera'<0  { local ltolera 0.0001 }
        if `iterate'<=0 { local iterate 50 }
        
        if ("`fam'"=="gau" | "`fam'"=="") {  
              qui reg `y' `*', dep(`y')  
              tempvar eta mu lrchi2 lrdev
              predict `mu' if `touse'
              gen `lrchi2'=(`y'-`mu')^2 if `touse'
              replace `lrchi2'=sum(`lrchi2')
              local lrchi2=`lrchi2'[_N]
              local lrdev=`lrchi2'
              gen `eta'=`mu' if `touse'
       }
        
* INITIALIZATION        
        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
        
        if "`fam'"=="bin" {                 
              gen `mu' = (`y'+.5)/(`m'+1) if `touse'
              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 "`fam'" == "poi"  { 
              egen `mu' = mean(`y') if `touse'
              if "`link'"=="id" { replace `z' = `y' }
              else { gen `eta' = log(`mu') if `touse'}
        }
        if "`fam'" == "gam"  { 
              egen `mu' = mean(`y') if `touse'
              if "`link'"=="log" {
                 gen `eta' = log(`mu') if `touse' }  
              else if "`link'"=="id" {  
                 replace `mu' = `y' 
                 replace `z' = `y' } 
              else { gen `eta' = 1/`mu' if `touse'}  
        }
        if "`fam'" == "invg" { 
              egen `mu' = mean(`y') if `touse'
              if "`link'"=="id" { 
                 replace `mu' = `y'  
                 replace `z' = `y' }
              else if "`link'"=="log" { 
                 gen `eta' = log(`mu') if `touse' }
              else { gen `eta' = 1/`mu'^2  if `touse'}
        }

* LOCAL SCORING ALGORITHM (IRLS)
if ("`fam'"~= "gau" & "`fam'"~="") {
        while ((abs(`ddev') >`ltolera' & `i'<`iterate')| (`i'== 0))   {
             
             *  WEIGHTING  (calculate initial w)
             if "`fam'" == "bin" { 
                if "`link'" == "l" { replace `w' = `mu'*(1-`mu') if `touse'}
                if "`link'" == "p" { replace `w' = 1/(sqrt(2*_pi)*exp(-`eta'^2/2)) if `touse'}
                if "`link'" == "c" { replace `w' = exp(`eta')*exp(- exp(`eta')) if `touse'}
             }
             if "`fam'" == "poi"  { 
                if "`link'"=="id" { replace `w' = `wt'/`mu' }
                else { replace `w' = `mu' if `touse' } 
             }
             if "`fam'" == "gam"  { 
                 if "`link'"=="log" { replace `w' = `mu' if `touse'} 
                 else if "`link'"=="id"  { replace `w' = `wt'/`mu'^2 if `touse' }
                 else  { replace `w' = `mu'^2 if `touse'}  
             } 
             if "`fam'" == "invg" { 
                 if "`link'"=="log" { replace `w' =  `mu' if `touse' }
                 else if "`link'"=="id" {replace `w' = `wt'/`mu'^3  if `touse'}
                 else { replace `w' = `mu'^3 if `touse'} 
             }
             
             * ARTIFICIAL RESPONSE (calculate z)
             if "`fam'" == "bin" & "`group'"=="" { 
                 replace `z' = `eta' + (`y' -`mu')/`w' `moffset' if `touse' 
             }
             if "`fam'"== "bin" & "`group'" ~= ""  {  
                 replace `z' = `eta' + (`ym'-`mu')/`w' `moffset' if `touse'
             }
             if "`fam'" == "poi" & "`link'"=="" { 
                 replace `z' = `eta' + (`y' -`mu')/`w' `moffset' if `touse' 
             }
             if "`fam'" == "gam" { 
                if "`link'"=="log" { replace `z' = `eta' + (`y' -`mu')/`w' `moffset' if `touse' } 
                else if "`link'"== ""   { replace `z' = `eta' - (`y' -`mu')/`w' `moffset' if `touse' }
             }
             if "`fam'" == "invg" & "`link'"~="id" { 
                if "`link'"=="log" { replace `z' = `eta' + (2*(`y'-`mu'))/`w' `moffset' if `touse' }
                else if "`link'"=="" { replace `z'=`eta'- (2*(`y'-`mu'))/`w' `moffset' if `touse' } 
             }  
            
             * ADJUSTED WEIGHTING FOR NONCANONICAL LINKS AND GROUPING
             if "`link'" == "p" { replace `w' = `w'^2/(`mu'*(1-`mu')) }
             if "`link'" == "c" { replace `w' = `w'^2/(`mu'*(1-`mu')) }
             if "`link'" == "log" { replace `w' = 1 }
             if "`group'"~= ""  { replace `w' = `w'*`m' } 
             
             * REGRESS
             regress `z' `*' [iw=`w'*`wt'], mse1 dof(100000) dep(`y')

             * INVERSE LINK
             cap drop `eta'
             predict `eta' if `touse' 
             replace `eta' = `eta' `poffset' 
             if "`fam'" == "bin"  { 
                    if "`link'" == "l"  {                
                         replace `mu' = 1/(1+exp(-`eta'))  
                    }
                    if "`link'" == "p"  {                
                         replace `mu' = normprob(`eta') 
                    }
                    if "`link'" == "c"  {                
                         replace `mu' = 1-exp(-exp(`eta')) 
                    }
             }
             if "`fam'" == "poi"  { 
                 if "`link'"=="id" { replace `mu' = `eta'}
                 else { replace `mu' = exp(`eta') }
             }
             if "`fam'" == "gam"  { 
                 if "`link'" == "log" { replace `mu' = exp(`eta') }
                 else if "`link'" == "id"  { replace `mu' = `eta'}
                 else { replace `mu' = 1/`eta' }
             }
             if "`fam'" == "invg" { 
                 if "`link'" == "log" { replace `mu' = exp(`eta') }
                 else if "`link'" == "id" { replace `mu' = `eta'}
                 else { replace `mu' = 1/sqrt(`eta') } 
             }
             
             * ITERATIONS IN DEVIANCE 
             replace `oldev' = `dev' 
             if "`fam'" == "bin" & "`group'" == ""  { 
                  replace `dev' = cond(`y',log(`mu'), log(1-`mu'))
             }     
             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
             }
             if "`fam'" == "poi" { 
                  replace `dev' = (`y'*log(`y'/`mu') - (`y'-`mu'))
             }
             if "`fam'" == "gam" {
                  replace `dev' = -log(`y'/`mu') + (`y'-`mu')/`mu'
             }
             if "`fam'" == "invg" { 
                  replace `dev' = (`y'-`mu')^2 / (`mu'^2*`y')
             }
             
             replace `dev' = sum(`dev'*`wt') 
             if ("`fam'"=="bin" & "`group'"=="") {
                     replace `dev' = -2*`dev'[_N] 
             }
             else if ("`fam'"=="invg") {
                     replace `dev' = `dev'[_N] 
             }
             else { replace `dev' = 2*`dev'[_N] }
             local ddev = `dev'-`oldev'
             noi di in gr "Iter `i' : Dev = " in ye %9.4f `dev'
             local i = `i' + 1

        }
* CHI SQUARE CALCULATION
   tempvar chi2
     if "`fam'"=="bin"  { gen `chi2' = (`y'-`mu')^2/(`m'*`mu'*(1-`mu')) if `touse' }
     if "`fam'"=="poi"  { gen `chi2' = (`y'-`mu')^2/`mu' if `touse' }
     if "`fam'"=="gam"  { gen `chi2' = (`y'-`mu')^2/`mu'^2  if `touse'}
     if "`fam'"=="invg" { gen `chi2' = (`y'-`mu')^2/`mu'^3 if `touse'}
     replace `chi2' = sum(`wt'*`chi2')
     replace `chi2' = `chi2'[_N]
   }
}
local df = `nobs'-_result(3)-1    

if ("`fam'"=="gau" |"`fam'"=="") {
   local   dev  =`lrdev'
   local   chi2 =`lrchi2'
}
di in gr _col(57) "No obs.    = " in ye %9.0g `nobs' 
di in gr _col(57) "Deviance   = " in ye %9.0g `dev'
di in gr _col(57) "Prob>chi2  = " in ye %9.0g chiprob(`df',`dev')
di in gr _col(57) "Dispersion = " in ye %9.0g `chi2'/`df'

if ("`fam'"=="gau" | "`fam'"=="")  {di in gr "Gaussian distribution"} 
if ("`fam'"=="bin" & "`group'"=="") {di in gr "Bernoulli distribution"} 
if ("`fam'"=="bin" & "`group'"~="") {di in gr "Binomial distribution"}
if "`fam'"=="poi" {di in gr "Poisson distribution"}
if "`fam'"=="gam" {di in gr "Gamma distribution"}
if "`fam'"=="invg" {di in gr "Inverse Gaussian distribution"}
if "`link'"=="l" {di in gr "Logit link"}
if "`link'"=="p" {di in gr "Probit link"}
if "`link'"=="c" {di in gr "Cloglog link"}
if "`link'"=="log" {di in gr "Log link"}
if "`link'"=="id" {di in gr "Identity link"}
if "`link'"=="" {di in gr "Canonical link"}

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

*  GAMMA AND INVERSE GAUSSIAN
if ("`fam'"=="gam"|"`fam'"=="invg") {       
di in gr _dup(78) "-"
di in gr _col(1) "`y'" _col(10) "|" _col(17) "Coef." _col(25) "Std.Err." _col(41) "t" _col(47) "P>|t|" _col(59) "[`level'% Conf. Interval]"
di in gr _dup(78) "-"
while "`1'"~="" {
local j=1
local se`j'=sqrt(`chi2'/`df')*_se[`1']
local t`j' = _b[`1']/`se`j''
local tpb`j' = tprob(`df',`t`j'')
local lci`j' = _b[`1']-(invnorm(`level'/100+(1-`level'/100)/2))*(`se`j'')
local uci`j' = _b[`1']+(invnorm(`level'/100+(1-`level'/100)/2))*(`se`j'')
di in gr _col(1) "`1'" _col(10) "|" in ye   _col(13) %9.0g _b[`1']  _col(24) %9.0g `se`j'' _col(35) %9.3f `t`j''  _col(47) %4.3f `tpb`j'' _col(58) %9.0g `lci`j'' _col(70) %9.0g `uci`j''
local j=`j'+1
mac shift
}
local secon = sqrt(`chi2'/`df')*_se[_cons]
local tcon = _b[_cons]/`secon'
local tpbcon = tprob(`df',`tcon') 
local lcicon = _b[_cons]-(invnorm(`level'/100+(1-`level'/100)/2))*(`secon')
local ucicon = _b[_cons]+(invnorm(`level'/100+(1-`level'/100)/2))*(`secon')
di in gr _col(1) "_cons" _col(10) "|" in ye _col(13) %9.0g _b[_cons] _col(24) %9.0g `secon' _col(35) %9.3f `tcon' _col(47) %4.3f `tpbcon' _col(58) %9.0g `lcicon' _col(70) %9.0g `ucicon'
di in gr _dup(78) "-"
di in gr "Standard errors adjusted by sqrt(dispersion)"
}
else { noi regress, noheader level(`level') `eform' }


* CALCULATION OF LINEAR PREDICTOR (ETA) AND PREDICTED (MU)
if "`stats'"~="" {
    qui   { 
        gen _eta = `eta' if `touse'
        gen _mu  = `mu' if `touse'
           if "`group'"~="" {
               replace _mu = _mu*`m'
           }
    }
    di in gr "Variables created:  "  in ye " _eta, _mu"
    lab var _eta "Linear predictor"
    lab var _mu "Predicted"
}

if "`resid'"~="" {
   qui {
        if ("`fam'"=="bin" & "`group'"=="")  { gen ___c = 1 if `touse'}
        else if ("`fam'"=="poi")  { gen ___c = 3 if `touse' }
        else if ("`fam'"=="gam")  { gen ___c = 4 if `touse'}
        else if ("`fam'"=="invg") { gen ___c = 5 if `touse' }
        else if ("`fam'"=="bin" & "`group'"~= "") {
              gen ___m = `m' if `touse'
              gen ___c = 2 if `touse'
        }
        gen ___y  = `y' if `touse'
        gen ___mu = `mu' if `touse'
   }
   _glmresd
}


end