*! version 1.1.1 2-21-93  J. Hilbe  STB-12 sg16.1
* RESIDUAL DIAGNOSTICS FOR GLM  
    
cap program drop _glmresd    
program define _glmresd
    version 3.0
    tempvar y mu m class stdp hat STDpr STDdev w
    
qui  {    

    gen `y'    = ___y 
    gen `mu'   = ___mu
    gen `class'= ___c
    gen `w'    = ___w
    if (`class'==2) {gen `m' = ___m}
    drop ___*

/*  CLASS 
    1 : family = bin, bernoulli
    2 : family = bin, binomial
    3 : family = poi
    4 : family = gam
    5 : family = invg
    6 : family = gau
*/

    
    if (`class' == 1) {                         /* bernoulli */
        tempvar stdp hat STDpr STDdev STDbeta
        gen _Presid     = (`y'-`mu')/sqrt(`mu'*(1-`mu'))
        gen _Dresid     = `y'*log(`y'/`mu') if `y'==1
        replace _Dresid = (1+`y')*log((1-`y')/(1-`mu')) if `y'==0
        replace _Dresid = sqrt(2*_Dresid)
        replace _Dresid = -_Dresid if (`y'-`mu')<0
        predict `stdp', stdp
        gen `hat'   =`stdp'^2*`mu'*(1-`mu')   
        gen `STDpr' = (_Presid/sqrt(1-`hat'))^2 
        gen `STDdev'= (_Dresid/sqrt(1-`hat'))^2  
        gen `STDbeta' = (_Presid^2*`hat')/(1-`hat')^2
        gen _Lresid = sqrt(`hat'*`STDpr'+(1-`hat')*`STDdev')
        replace _Lresid = -_Lresid if (`y'-`mu')<0
        gen _Aresid = (2.0533902*ibeta(2/3,2/3,`y'))-(2.0533902*ibeta(2/3,2/3,`mu'))
        replace _Aresid = _Aresid/(`mu'*(1-`mu'))^(1/6) 
        gen _hat = `hat'        
        gen _Dpr = `STDpr'    
        gen _Ddev = `STDdev'
        gen _Dbeta = `STDbeta'
        noi di in gr "Variables created: " in ye "  _Presid, _Dresid, _Lresid, _Aresid, _hat"
        noi di in ye "                     _Dpr, _Ddev, _Dbeta"
        lab var _Presid "Pearson residual"
        lab var _Dresid "Deviance residual"
        lab var _Lresid "Likelihood residual"
        lab var _Aresid "Anscombe residual"
        lab var _hat "Hat"
        lab var _Dpr    "Delta Pearson residual"
        lab var _Ddev   "Delta deviance residual"
        lab var _Dbeta  "Delta beta residual"
    }
    if (`class' == 2) {                              /* binomial  */
        tempvar yhat pp hat stdp STDpr STDdev STDbeta
        gen `yhat' = `mu'*`m'
        gen `pp' = `yhat'/`m'
        predict `stdp',stdp
        gen `hat'=`stdp'^2 * `yhat'*(1-`mu') 
        gen _Presid = (`y'-`yhat')/sqrt(`yhat'*(1-`mu'))
        gen _Dresid = `y'*log(`y'/`yhat') + (`m'-`y')*log((`m'-`y')/(`m'-`yhat'))
        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 `STDbeta' = (_Presid^2*`hat')/(1-`hat')^2  
        gen `STDpr'=_Presid/sqrt(1-`hat')
        gen `STDdev'=_Dresid/sqrt(1-`hat')
        gen _Lresid = sqrt(`hat'*`STDpr'^2+(1-`hat')*`STDdev'^2)
        replace _Lresid = -_Lresid if (`y'-`mu')<0
        noi di in gr "Variable created: " in ye " _Presid _Dresid  _Aresid _Lresid" 
        noi di in ye "                   _hat _Dpr _Ddev _Dbeta"
        gen _Dpr = `STDpr'
        gen _Ddev = `STDdev'
        gen _Dbeta = `STDbeta'
        gen _hat = `hat'
        lab var _hat "Hat"
        lab var _Presid "Pearson residual"
        lab var _Dresid "Deviance residual"
        lab var _Aresid "Anscombe residual"
        lab var _Lresid "Likelihood residual"
        lab var _Dpr "Delta Pearson residual"
        lab var _Ddev   "Delta deviance residual"
        lab var _Dbeta  "Delta beta residual"
    }
    if (`class' == 3)   {                           /* poisson  */
        gen _Presid = (`y'-`mu')/sqrt(`mu') * sqrt(`w')
        gen _Aresid = (3/2)*(`y'^(3/2) - `mu'^(3/2)) / `mu'^(1/6) * sqrt(`w')
        gen _Dresid = (`y'-`mu')*sqrt(2*(`y'*log(`y'/`mu')-`y'+`mu')) * sqrt(`w')
        replace _Dresid = -_Dresid if (`y'-`mu')<0
        noi di in gr "Variables created: " in ye " _Presid, _Aresid, _Dresid"
        lab var _Presid "Pearson residual"
        lab var _Dresid "Deviance residual"
        lab var _Aresid "Anscombe residual"
    }
    if (`class' == 4)   {                           /* gamma  */
        gen _Presid = (`y'-`mu')/`mu' * sqrt(`w')           
        gen _Aresid = (3*(`y'^(1/3) - `mu'^(1/3))/`mu'^(1/3)) * sqrt(`w')
        gen _Dresid = (`y'-`mu')*sqrt(2*(-log(`y'/`mu')+(`y'-`mu')/`mu')) * sqrt(`w')
        replace _Dresid = -_Dresid if (`y'-`mu')<0
        noi di in gr "Variables created: " in ye " _Presid, _Aresid, _Dresid"
        lab var _Presid "Pearson residual"
        lab var _Dresid "Deviance residual"
        lab var _Aresid "Anscombe residual"
    }        
    if (`class' == 5)   {                    /* inverse gaussian  */
        gen _Presid = (`y'-`mu')/`mu'^(3/2) * sqrt(`w')           
        gen _Aresid = (log(`y')-log(`mu'))/`mu'^(1/2) * sqrt(`w')
        gen _Dresid = (`y'-`mu')*sqrt((`y'-`mu')^2/(`y'*`mu'^2)) * sqrt(`w')
        replace _Dresid = -_Dresid if (`y'-`mu')<0
        noi di in gr "Variables created: " in ye " _Presid, _Aresid, _Dresid"
        lab var _Presid "Pearson residual"
        lab var _Dresid "Deviance residual"
        lab var _Aresid "Anscombe residual"
    }        
    if (`class' == 6)  {                     /* gaussian */
        gen _Presid = `y'-`mu'
        noi di in gr "Variable created :  " in ye " _Presid"
        lab var _Presid "Pearson residual"
    }

}    
end