%macro varmod(version,

  /* Identify input data set, response variable and data set to score */
       data=_last_,     /* Specify the data set to be used.  If omitted,
                           the last-created data set is used.         */
       response=,       /* REQUIRED.  The response variable for the mean
                           model.                                     */
       score=,          /* Specify an optional data set of observations
                           for which predicted means and variances will
                           be computed.  These observations do not
                           affect the model fit.  The response variable,
                           if present, is copied to the variable
                           _OBSRESP.  All explanatory variable values
                           must be nonmissing for predicted values to be
                           computed.  The OUTSCORE= option should also
                           be specified when SCORE= is specified.     */

  /* Specify mean and variance models */
       mclass=,         /* Specify all CLASS variables (if any) used in
                           the mean model.                            */
       mmodel=,         /* Specify all terms in the mean model using
                           PROC GLM MODEL statement syntax.           */
       vclass=,         /* Specify all CLASS variables (if any) used in
                           the variance model.                        */
       vmodel=,         /* Specify all terms in the variance model using
                           PROC GLM MODEL statement syntax.           */
       mnoint=,         /* If this is non-blank, fit a no-intercept mean
                           model.                                     */

  /* Output data sets */
       outmdiag=,       /* Specify name of output data set of predicted
                           values of the mean and diagnostics.        */
       outvdiag=,       /* Specify name of output data set of predicted
                           values of the variance and diagnostics.    */
       outpred=,        /* Specify name of output data set containing
                           the DATA= data set plus a variable _PREDM
                           giving the predicted mean values and a
                           variable _PREDV giving the predicted
                           variances.                                 */
       outscore=,       /* Specify name of output data set containing
                           the SCORE= data set plus a variable _PREDM
                           giving the predicted mean values and a
                           variable _PREDV giving the predicted
                           variances.  SCORE= should be specified when
                           OUTSCORE= is specified.                    */
       outmparm=,       /* Specify name of output data set containing
                           parameter estimates of the mean model.     */
       outvparm=,       /* Specify name of output data set containing
                           parameter estimates of the variance model. */
       outlr=,          /* Specify name of output data set containing
                           the likelihood ratio test for equal
                           variances.                                 */


  /* Control printing */
       noprint=,        /* If non-blank, suppress all printed output.
                           Useful when you just want to create output
                           data sets.                                 */
       noeqvar=,        /* If non-blank, suppress printing equal
                           variances mean model results.              */
       noiter=,         /* If non-blank, suppress printing iteration
                           history.                                   */
       nolr=,           /* If non-blank, suppress printing likelihood
                           ratio test for equal variances.            */
       noparms=,        /* If non-blank, suppress printing parameters
                           of the mean and variance models.           */
       mgof=,           /* If non-blank, prints goodness of fit table
                           for final mean model.                      */
       vgof=,           /* If non-blank, prints goodness of fit table
                           for final variance model.                  */

  /* Control maximum likelihood process */
       maxiter=20,      /* Maximum number of iterations allowed.  An
                           iteration is one fitting of both models.   */
       converge=1e-3    /* Convergence is declared when the maximum
                           change in loglikelihood is less than this. */
       );

%if &version ne %then %put VARMOD macro Version 1.7;

%let opts = %sysfunc(getoption(notes))
            _last_=%sysfunc(getoption(_last_));
%if &data=_last_ %then %let data=&syslast;
options nonotes;

/* Verify that RESPONSE= option is specified. */
%if &response= %then %do;
    %put ERROR: Response variable must be specified in the RESPONSE= argument.;
    %goto exit;
%end;

%if &noprint ne %then %do;
    %let noeqvar=no; %let noiter=no; %let nolr=no; %let noparms=no;
%end;

/* Add a variable to identify the SCORE= observations when concatenated
   to the DATA= observations.  _OBSRESP contains the observed responses,
   or is missing.  The response variable is set to missing to ensure
   these observations do not affect the model fit.
*/
%if &score ne %then %do;
data _score;
  set &score;
  _scorobs=1;
  _obsresp=&response;
  &response=.;
  run;
%end;

/* For 0th iteration (equal variances model fit):
   - move data to temporary file (_work)
   - set weights to 1 (_wgt)
   - initialize -2LL (_sum) and iteration (_iter) variables
*/
data _work;
  set &data
      %if &score ne %then _score;
      ;
  _wgt = 1;
  _sum = 0;
  run;

data _tmp;
  _sum = 0; _iter=0;
  run;

%let conv = 0;
%let iter = 0;


/********************************************************************
   Iteratively fit mean and variance models until convergence.
   0th iteration fits the equal variances model (intercept-only).
   Convergence declared if change in -2*logLikelihood < &converge in
   fewer than &maxiter iterations.
********************************************************************/
%do %while( &conv = 0 );

%if &iter > &maxiter - 1 %then %do;
  %put WARNING: Convergence not attained in &maxiter iterations.  Terminating.;
  %goto exit;
%end;

/* Store previous iteration's -2LL for comparison with this iteration */
data _old;
  set _tmp;
  _devold = _sum;
  keep _devold;
  run;

/*************************** MEAN MODEL *****************************
  set NOSCALE so that scale is not estimated
  select OBSTATS to get residuals
  SCWGT selects dispersion parameter weights
*/
%put VARMOD Iteration &iter: Fitting mean model;
ods exclude all;
%if &iter = 0 and &noeqvar eq %then %str(ods exclude obstats;);
title "Equal variances model";
proc genmod data=_work;
  scwgt _wgt;
  ods output modelfit=_mnfit;
  ods output obstats=_mndiag;
  ods output parameterestimates=_meanest;
  %if &mclass ne %then %str(class &mclass;);
  model &response = %quote(&mmodel)
     / noscale obstats
       %if &mnoint ne %then noint;
     ;
  run;
%if &iter = 0 and &noeqvar eq %then %str(ods exclude all;);

/* Obtain squared residuals
   Add 1e-12 to the squared residuals to keep them off the zero
   boundary.
*/
data _work;
  drop _sum resraw;
  set _mndiag (keep=resraw);
  set _work;
  _rsquare = resraw*resraw + 1e-12;
  run;

/*********************** VARIANCE MODEL *****************************
  set scoring=100 to use Fisher scoring
  set scale = .5 for 1 df in gamma distribution
  first time thru, fit intercept-only model to get -2LL for the equal
    variances model
*/
%put VARMOD Iteration &iter: Fitting variance model;
proc genmod data=_work;
  ods output modelfit=_varfit;
  ods output obstats=_vardiag;
  ods output parameterestimates=_varests;
  %if &vclass ne %then %str(class &vclass;);
  model _rsquare =   %if &iter>0 %then %quote(&vmodel);
        / dist=gamma link=log obstats noscale scale=.5 scoring=100;
  run;

/* Get DF of equal variances model */
%if &iter=0 %then %do;
  data _null_; set _varfit;
    if _n_=1 then call symput("eqdf",df);
    stop;
    run;
%end;

/* Compute weights for mean model
   Compute _sum = -2*log(likelihood)
*/
data _work;
  drop xbeta pred std hesswgt upper lower resraw reschi resdev;
  set _work;
  set _vardiag;
  _wgt = 1./pred;
  _sum + (_rsquare/pred + log(pred));
  /* Aitkin adds log(2*pi) to the likelihood contribution above.  It
     isn't necessary since it isn't a function of the parms, but it can
     be added in by inserting the following before the last parenthesis
     above:          + log(2*arcos(-1))
  */
  run;

/* Update -2LL */
data _tmp;
  set _work nobs = _last;
  keep _sum _iter;
  if _n_ = _last;
  _iter = &iter ;
  run;

/* Update iteration history */
data _iter;
  set %if &iter>0 %then _iter;
      _tmp;
  run;

/* Check for convergence */
data _NULL_;
  set _tmp; set _old;
  put "VARMOD Iteration &iter: -2*LogLikelihood = " _sum;
  if ( abs(_sum - _devold) <= &converge ) then do;
     _conv = 1;
     put "VARMOD: Convergence attained.";
  end;
  else _conv = 0;
  call symput( 'dev', left(put(_sum,12.4)) );
  call symput( 'conv', left(put(_conv,3.)) );
  run;

%let iter=%eval( &iter+1 );

%end;
/******************* End of model-fitting loop *********************/

ods select all;

/*************** Print iteration history ******************/
%if &noiter= %then %do;
title 'Iteration History';
proc print data=_iter label noobs;
  var _iter _sum;
  format _sum best15.;
  label _iter="Iteration"
        _sum="-2*LogLikelihood";
  run;
%end;

/***************** Compute/print/output LR test ********************/
title 'Likelihood Ratio Test of Equal Variances';
data _lreqv;
  keep lrchi lrdf lrp;
  label lrchi="LR ChiSquare"
        lrdf="DF"
        lrp="Pr>Chi";
  format lrp 6.4;
  set _iter end=_eof;
  retain eq2ll;
  if _n_=1 then do;
     eq2ll=_sum;
     set _varfit;
     retain lrdf;
     lrdf=&eqdf-df;
  end;
  if _eof then do;
     lrchi=eq2ll-_sum;
     if lrdf>0 then lrp=1-probchi(lrchi,lrdf);
     output;
  end;
  run;

%if &nolr= %then %do;
proc print data=_lreqv label noobs;
  var lrchi lrdf lrp;
  format lrp pvalue.;
  run;
%end;

/*************** Print mean model parameter estimates **************/
title 'Final Mean Model';
%if &mgof ne %then %do;
proc print data=_mnfit noobs label;
  where criterio ne "Log Likelihood";
  title2 "Criteria for Assessing Goodness of Fit";
  run;
%end;
%if &noparms= %then %do;
proc print data=_meanest noobs label;
  title2 "Analysis of Parameter Estimates";
  run;
%end;

/************** Print variance model parameter estimates *************/
title 'Final Variance Model';
%if &vgof ne %then %do;
proc print data=_varfit noobs label;
  where criterio ne "Log Likelihood";
  title2 "Criteria for Assessing Goodness of Fit";
  run;
%end;
/* Compute variance ratios */
data _varests;
  set _varests;
  if parameter ne "SCALE" and df ne 0 then
     varratio=exp(estimate);
  label varratio="Variance Ratio";
  run;
%if &noparms= %then %do;
proc print data=_varests noobs label;
  title2 "Analysis of Parameter Estimates";
  run;
%end;

/********************** Output data sets *********************/
options notes;
/* OUTMDIAG= data set -- mean model diagnostics */
%if &outmdiag ne %then %do;
  %if &score= %then %do;
  data &outmdiag;
    set _mndiag;
    run;
  %end;
  %else %do;
  data &outmdiag;
    merge _work (keep=_scorobs)
          _mndiag;
    drop _scorobs;
    if _scorobs ne 1;
    run;
  %end;
%end;

/* OUTVDIAG= data set -- variance model diagnostics */
%if &outvdiag ne %then %do;
  %if &score= %then %do;
  data &outvdiag;
    set _vardiag;
    run;
  %end;
  %else %do;
  data &outvdiag;
    merge _work (keep=_scorobs)
          _vardiag;
    drop _scorobs;
    if _scorobs ne 1;
    run;
  %end;
%end;

/* OUTPRED= data set -- predicted mean and variance for DATA= */
%if &outpred ne %then %do;
  data &outpred;
    merge _work
          _mndiag  (keep=pred rename=(pred=_predm))
          _vardiag (keep=pred rename=(pred=_predv));
    %if &score ne %then %str(if _scorobs ne 1;);
    drop _wgt _sum _rsquare
      %if &score ne %then _scorobs _obsresp;
      ;
    run;
%end;

/* OUTSCORE= data set -- predicted mean and variance for SCORE= */
%if &outscore ne %then %do;
  %if &score= %then %do;
    %put ERROR: No SCORE= data set specified.  OUTSCORE= is ignored.;
    %goto clean;
  %end;
  data &outscore;
    merge _work
          _mndiag  (keep=pred rename=(pred=_predm))
          _vardiag (keep=pred rename=(pred=_predv));
    drop _wgt _sum _rsquare _scorobs &response;
    if _scorobs=1;
    run;
%end;

/* OUTMPARM= data set -- mean model parameter estimates */
%if &outmparm ne %then %do;
  data &outmparm;
    set _meanest;
    run;
%end;

/* OUTVPARM= data set -- variance model parameter estimates */
%if &outvparm ne %then %do;
  data &outvparm;
    set _varests;
    run;
%end;

/* OUTLR= data set -- Likelihood ratio test of equal variances */
%if &outlr ne %then %do;
  data &outlr;
    set _lreqv;
    run;
%end;
options nonotes;

/******************** Clean up temporary data sets ******************/
%clean:
proc datasets nolist;
  delete _tmp _work _old _iter _mnfit _mndiag _meanest
         _varfit _vardiag _varests _lreqv
         %if &score ne %then _score;
         ;
  quit;

title;
%exit:
options &opts;
ods select all;
%mend varmod;