* Quasilikelihood Chapter
drop _all

use accident, clear
xi: glm y I.site I.variety, family(binom) link(logit) nolog 

xi: glm y I.site I.variety, family(binom) link(logit) scale(x2) nolog 

predict double eta, eta
predict double pearson, pearson
predict double mu
gen double logvar = log(mu*mu*(1-mu)*(1-mu))
label var pearson "Pearson residuals"
label var eta     "Linear predictor"
label var logvar  "Log(variance)"
twoway scatter pearson eta, yline(0)
twoway scatter pearson logvar, yline(0) 

twoway scatter pearson eta, yline(0)
graph export quasi1.eps, replace

twoway scatter pearson logvar, yline(0) 
graph export quasi2.eps, replace


capture drop eta
capture drop mu
capture drop pearson
capture drop logvar
                                                                                                                             
capture program drop binsq
program define binsq
        args todo eta mu return touse
                                                                                                                             
        if `todo' == -1 {
                if "$SGLM_m" != "1" {                           /* Title */
                        global SGLM_vt "Quasi"
                        global SGLM_vf "u^2*(1-u/$SGLM_m)^2"
                }
                else {
                        global SGLM_vt "Quasi"
                        global SGLM_vf "u^2*(1-u)^2"
                }
                                                                                                                             
                local y     "$SGLM_y"
                local m     "$SGLM_m"
                local touse "`eta'"
                                                                                                                             
                capture assert `m'>0 if `touse'         /* sic, > */
                if _rc {
                        di in red `"`m' has nonpositive values"'
                        exit 499
                }
                capture assert `m'==int(`m') if `touse'
                if _rc {
                        di in red `"`m' has noninteger values"'
                        exit 499
                }
                capture assert `y'>=0 if `touse'
                if _rc {
                        di in red `"dependent variable `y' has negative values"'
                        exit 499
                }
                capture assert `y'<=`m' if `touse'
                if _rc {
                        noi di in red `"`y' > `m' in some cases"'
                        exit 499
                }
                                                                                                                             
                global SGLM_mu "glim_mu 0 $SGLM_m"
                exit
        }
        if `todo' == 0 {                        /* Initialize eta */
                if "$SGLM_L" == "glim_l01" {                    /* Identity */
                        gen double `eta' = `mu'/$SGLM_m
                }
                else if "$SGLM_L" == "glim_l02" {               /* Logit */
                        gen double `eta' = ln(`mu'/($SGLM_m-`mu'))
                }
                else if "$SGLM_L" == "glim_l03" {               /* Log */
                        gen double `eta' = ln(`mu'/$SGLM_m)
                }
                else if "$SGLM_L" == "glim_l05" {               /* Log compl */
                        gen double `eta' = ln(1-`mu'/$SGLM_m)
                }
                else if "$SGLM_L" == "glim_l06" {               /* Loglog */
                        gen double `eta' = -ln(-ln(`mu'/$SGLM_m))
                }
                else if "$SGLM_L" == "glim_l07" {               /* Cloglog */
                        gen double `eta' = ln(-ln(1-`mu'/$SGLM_m))
                }
                 else if "$SGLM_L" == "glim_l08" {               /* Probit */
                        gen double `eta' = invnormal(`mu'/$SGLM_m)
                }
                else if "$SGLM_L" == "glim_l09" {               /* Reciprocal */
                        gen double `eta' = $SGLM_m/`mu'
                }
                else if "$SGLM_L" == "glim_l10" {               /* Power(-2) */
                        gen double `eta' = ($SGLM_m/`mu')^2
                }
                else if "$SGLM_L" == "glim_l11" {               /* Power(a) */
                        gen double `eta' = (`mu'/$SGLM_m)^$SGLM_a
                }
                else if "$SGLM_L" == "glim_l12" {               /* OPower(a) */
                        gen double `eta' = /*
                                */ ((`mu'/($SGLM_m-`mu'))^$SGLM_a-1) / $SGLM_a
                }
                else {
                        gen double `eta' = $SGLM_m*($SGLM_y+.5)/($SGLM_m+1)
                }
                exit
        }
        if `todo' == 1 {                        /* V(mu) */
                gen double `return' =  `mu'*`mu'*(1-`mu'/$SGLM_m)*(1-`mu'/$SGLM_m)
                exit
        }
        if `todo' == 2 {                        /* (d V)/(d mu) */
                gen double `return' =  2*`mu' - 6*`mu'*`mu'/$SGLM_m + /*
                        */ 4*`mu'*`mu'*`mu'/($SGLM_m*$SGLM_m)
                exit
        }
        if `todo' == 3 {                        /* deviance */
                noi di in red "deviance calculation not supported"
                gen double `return' = 0
                exit
        }
        if `todo' == 4 {                        /* Anscombe */
                noi di in red "Anscombe residuals not supported"
                exit 198
         }
        if `todo' == 5 {                        /* ln-likelihood */
                local y "$SGLM_y"
                local m "$SGLM_m"
                if "`y'" == "" {
                        local y "`e(depvar)'"
                }
                if "`m'" == "" {
                        local m "`e(m)'"
                }
                gen double `return' = (2*`y'-1)*ln(`mu'/(`m'-`mu')) - `y'/`mu' - /*
                        */ (`m'-`y')/(`m'-`mu')
                exit
        }
        if `todo' == 6 {
                noi di in red "adjusted deviance residuals not supported"
                exit 198
        }
        noi di in red "Unknown call to glim variance function"
        error 198
end

xi: glm y I.site I.variety, family(binsq) link(logit) nolog 

predict double eta, eta
predict double pearson, pearson
predict double mu
gen double logvar = log(mu*mu*(1-mu)*(1-mu))
label var eta     "Linear predictor"
label var pearson "Pearson residuals"
label var logvar  "Log(variance)"
twoway scatter pearson eta, yline(0) 
graph export quasi3.eps, replace

twoway scatter pearson logvar, yline(0) 
graph export quasi4.eps, replace

* The gamfit.exe executable no longer works

drop _all   
mata: mata clear                                                                                                                        
set mem 64m
use fast2, clear
gen ante = anterior
capture program drop gam

glm death anterior hcabg kk2 kk3 age2-age4 cases, family(bin) irls nolog 

if 0 {
gam death anterior hcabg kk2 kk3 age2-age4 cases, family(bin) df(cases:4)
}

use medpar, clear
glm los hmo white age, family(poisson) irls nolog