I can calculate the motion of heavenly bodies but not the madness of people. -Isaac Newton

Estimate Quasi-Binomial Model with GENMOD Procedure in SAS

In my previous post (https://statcompute.wordpress.com/2015/11/01/quasi-binomial-model-in-sas/), I’ve shown why there is an interest in estimating Quasi-Binomial models for financial practitioners and how to estimate with GLIMMIX procedure in SAS.

In the demonstration below, I will show an example how to estimate a Quasi-Binomial model with GENMOD procedure by specifying VARIANCE and DEVIANCE. While the CPU time for model estimation is a lot faster with GENMOD than with GLIMMIX, additional steps are necessary to ensure the correct statistical inference.

```ods listing close;
ods output modelfit = fit;
ods output parameterestimates = parm1;
proc genmod data = kyphosis;
model y = age number start / link = logit noscale;
variance v = _mean_ * (1 - _mean_);
deviance d = (-2) * log((_mean_ ** _resp_) * ((1 - _mean_) ** (1 - _resp_)));
run;
ods listing;

proc sql noprint;
select
distinct valuedf format = 12.8, df format = 8.0 into :disp, :df
from
fit
where
index(criterion, "Pearson") > 0;

select
value format = 12.4 into :ll
from
fit
where
criterion = "Log Likelihood";

select
sum(df) into :k
from
parm1;
quit;

%let aic = %sysevalf((-&ll + &k) * 2);
%let bic = %sysevalf(-&ll * 2 + &k * %sysfunc(log(&df + &k)));

data parm2 (keep = parameter estimate stderr2 df t_value p_value);
set parm1;
where df > 0;

stderr2 = stderr * (&scale ** 0.5);
df = &df;
t_value = estimate / stderr2;
p_value = (1 - probt(abs(t_value), &df)) * 2;
run;

title;
proc report data = parm2 spacing = 1 headline nowindows split = "*";
column(" Parameter Estimate of Quasi-Binomial Model "
parameter estimate stderr2 t_value df p_value);
compute before _page_;
line @5 "Fit Statistics";
line " ";
line @3 "Quasi Log Likelihood: %sysfunc(round(&ll, 0.01))";
line @3 "Quasi-AIC (smaller is better): %sysfunc(round(&aic, 0.01))";
line @3 "Quasi-BIC (smaller is better): %sysfunc(round(&bic, 0.01))";
line @3 "(Dispersion Parameter for Quasi-Binomial is %sysfunc(round(&disp, 0.0001)))";
line " ";
endcomp;
define parameter / "Parmeter" width = 10 order order = data;
define estimate  / "Estimate" width = 10 format = 10.4;
define stderr2   / "Std Err"  width = 10 format = 10.4;
define t_value   / "T Value"  width = 10 format = 10.2;
define df        / "DF"       width = 5  format = 4.0;
define p_value   / "Pr > |t|" width = 10 format = 10.4;
run;
quit;

/*
Fit Statistics

Quasi Log Likelihood: -30.69
Quasi-AIC (smaller is better): 69.38
Quasi-BIC (smaller is better): 78.96
(Dispersion Parameter for Quasi-Binomial is 0.9132)

Parameter Estimate of Quasi-Binomial Model
Parmeter     Estimate    Std Err    T Value    DF   Pr > |t|
------------------------------------------------------------
Intercept     -2.0369     1.3853      -1.47    77     0.1455
Age            0.0109     0.0062       1.77    77     0.0800
Number         0.4106     0.2149       1.91    77     0.0598
Start         -0.2065     0.0647      -3.19    77     0.0020
*/
```