/***********************************************************************************************/
/*                                                                                             */
/* aic_macro_SAS9                                                                              */
/*                                                                                             */
/*IMPORTANT!!! Use these macros with SAS version 9 and greater                                 */
/*                                                                                             */
/* SAS macro for computing aic and aicc from PROC GENMOD output.  Results from multiple        */
/* runs of GENMOD are accumulated in a dataset named by the user.  To use the macro, first     */
/* run PROC GENMOD and include the following statements after the model statement:             */
/*                                                                                             */
/*     ODS output modelfit=modelfit;                                                           */
/*     ODS output modelinfo=modelinfo;                                                         */
/*     ODS output ParameterEstimates=ParameterEstimates;                                       */
/*                                                                                             */
/* These statements invoke the Output Delivery System to create SAS datasets containing        */
/* GENMOD output.  Next, invoke the macro:                                                     */
/*                                                                                             */
/*     %AICC(dsn=chat_aic,estdsn=estimate,model=sex age date**2 avgprecip,model_dimension=1);  */
/*                                                                                             */
/* The macro has 4 parameters.  Values are assigned by the user. Parameters are as follows:    */
/*     DSN names a dataset in which AIC results are accumulated                                */
/*     ESTDSN names a dataset in which parameter estimates from the model are stored           */
/*     MODEL is a label used to identify the model                                             */
/*     MODEL_DIMENSION takes on an integer value that identifies the highest order interaction */
/*        term (involving only CLASS variables) that is specified in ANY model in the candidate*/
/*        set, not just the current model.  For example, suppose the candidate set of models   */
/*        includes two CLASS variables (variables that are names in a CLASS statement), say    */
/*        YEAR and SEX. If the candidate set includes models that contain YEAR*SEX interaction */
/*        terms, then MODEL_DIMENSION=2.  If the candidate set includes terms for SEX and AGE  */
/*        that involve only main effects or interactions with non-CLASS variables, then        */
/*        MODEL_DIMENSION=1.                                                                   */
/*                                                                                             */
/* IMPORTANT: This macro calls a 1-observation dataset named "N_eff" that contains a variable  */
/*            named "N_eff."  The value of that variable is the effective sample size for      */
/*            computing AICc.  This data set must be created by the user prior to calling      */
/*            the aicc macro.  Sample code for creating the data set appears below:            */
/*                                                                                             */
/*            The data set Nests contains one observation for each observation interval.       */
/*            The variable t is the length of the interval.                                    */
/*            Surv=1 if the clutch survived the interval, Surv=0 if it did not.                */
/*                                                                                             */
/*       Data N_eff(keep=n_eff);                                                               */
/*         Set Nests end=lastobs;                                                              */
/*         if surv=0 then n_eff+1;                                                             */
/*         else if surv=1 then n_eff+t;                                                        */
/*         if lastobs;                                                                         */
/*       run;                                                                                  */
/*                                                                                             */
/* NOTE:  Missing values for one or more explantory variables can result in different subsets  */
/*        of the data being used for different models.  Before running analyses (and before    */
/*        creating the N_eff data set) I strongly recommend that any observation with missing  */
/*        values for any explanatory variable in the canadidate set be deleted.                */
/*                                                                                             */
/***********************************************************************************************/

%Macro aicc(dsn=,estdsn=,model=,model_dimension=);

	/* Obtain dfg (number of observations or intervals used in the analysis minus */
	/* the number of parameters estimated), loglike (maximized value of the       */
	/* Log-likelihood function), and the model deviance.                          */

	data temp;
	   retain dfg deviance loglike;
	   set modelfit;
	   if _n_=1 then dfg=df;
	   if criterion='Log Likelihood' then loglike=value;
	   if criterion='Deviance' then deviance=value;
	   if criterion='Log Likelihood' then output;
	   keep dfg loglike deviance;
	run;
	 
	/* obsused is the number of observations (intervals) used in the analysis. */
	/* ness is the effective sample size used in AICc claculations             */
	/* model is a label that is assigned by the user.                          */

	data temp2;
	   retain N_eff;
	   if _n_=1 then set N_eff;
	   set nobs;
	   length model $ 60;
	   model="&model";
	   if label='Number of Observations Used';
	   Obsused=n;
	   keep obsused model n_eff;
	run;
	
	/* param2 is the parameter name for continuous covariates, and the parameter */
	/* name concatenated with level for variables named in a class statement.    */
	/* estimate is the estimated value of the parameter.                         */
	/* StdErr is the estimated standard error for the parameter named by param2. */

	data temp4(keep=param2 estimate) temp4a(keep=param2 StdErr)
	     temp4b(keep=level1-level&model_dimension df parameter estimate StdErr);
	    length Param2 $24 parameter $16 level1-level&model_dimension $16;
		retain level1-level&model_dimension ' ';
	    Set ParameterEstimates;
	    if parameter='Scale' then delete;
	    if parameter='Intercept' then level1=' ';
	    if level1= '.' then param2=trim(parameter); 
	    else param2=trim(parameter)||trim(level1); 
	run;

	data temp6;
	   retain model obsused n_eff;
	   if _n_=1 then set temp2;
	   set temp4b;
	run;
	    
	/* Compute AIC (Burnham and Anderson 2002:p.61) and  */
	/* AICC (Burnham and Anderson 2002:p.66.             */

	data temp3;
	   merge temp temp2;
	   n=n_eff;
	   k=obsused-dfg;
	   aic = -2*loglike + 2*K;
	   aicc = aic +(2*K*(K+1))/(n-k-1);
	run;
	  
	/* Save the values of AIC and AICC for this particular model to a dataset that    */
	/* that will hold the results from all models.  It is important that the dataset  */
	/* named in the base= statement be an empty data set when the first model is run. */
	  
	Proc append base=&dsn data=temp3 force;
	run;

	*proc print data=&dsn; run;

	/* Save the parameter estimates for this particular model */

	Proc append base=&estdsn data = temp6 force;
	run;

	*proc print data=&estdsn; run;

%Mend aicc;


/***********************************************************************************************/
/*                                                                                             */
/* delta_aic_macro.sas                                                                         */
/*                                                                                             */
/* SAS macro for computing values of delta AICC and Akaike weights.  To use the macro, first   */
/* use the AICC macro to generate a dataset containing AIC results from                        */
/* multiple models.  Then invoke the deltaaic macro as follows:                  :             */
/*                                                                                             */
/*     %DeltaAICDataIn=chat_aic,DataOut=chat_delta_aicc,Vari=aicc);                            */
/*                                                                                             */
/* The macro has 3 parameters.  Values are assigned by the user. Parameters are as follows:    */
/*     DATAIN names a dataset created with the AICC macro                                      */
/*     DATAOUT names a dataset in which delta AICC values and Akaike weights will be stored    */
/*     VARI=AIC request an analysis on AIC values, whereas VARI=AICC requests an analaysis on  */
/*        AICC values.                                                                         */
/*                                                                                             */
/***********************************************************************************************/

%Macro DeltaAic(DataIn=,DataOut=,Vari=);

	/** obtain minimum aicc value **/
	proc univariate noprint data=&DataIn;
	 var &vari;
	 output out=min min=minvari;
	run;

	/** add minaicc to aic2 data set and compute delta aicc (delta_aicc)**/
	data aic2;
	 retain minvari;
	 set &DataIn;
	 if _n_=1 then set min;
	 delta_&vari = &vari - minvari;
	 wnum=exp((-0.5 * delta_&vari));
	run;

	proc sort; by delta_&vari;
	run;

	proc univariate noprint data=aic2;
	  Var wnum;
	  output out=wdenominator sum=wden;
	run;

	data &DataOut;
	  Retain wden;
	  if _n_ = 1 then set wdenominator;
	  set aic2;
	  w_Akaike= wnum/wden;
	run;

	proc print data=&DataOut;
	  var model loglike deviance n k aic aicc delta_&vari w_Akaike;
	  title "Delta values and Akaike weights based on &vari";
	run;

%Mend DeltaAic;


/***********************************************************************************************/
/*                                                                                             */
/* model_avg_macro.sas                                                                         */
/*                                                                                             */
/* SAS macro for computing model-averaged parameter estimates and unconditional standard       */
/* errors.  To use the macro, first use the AICC macro to generate a dataset                   */
/* containing AIC results from multiple models. Then invoke the DELTAAIC                       */
/* to compute Akaike weights for each model. Next, invoke the MODELAVG macro as follows:       */
/*                                                                                             */
/*     %MODELAVG(Estimate=Params, Akaike=Aweights, DataOut=Model_Averaged,.model_dimension=1)  */
/*                                                                                             */
/* The macro has 4 parameters.  Values are assigned by the user. Parameters are as follows:    */
/*     ESTIMATE names a dataset containing parameter estimates, and is created with the AICC   */
/*          macro                                                                              */
/*     AKAIKE names a dataset containing Akaike weights and is created with the DELTAAIC macro */
/*     DATAOUT names a dataset to hold the results of model-averaging                          */
/*     MODEL_DIMENSION is an integer that gives the dimension of the highest-order model (see  */
/*         the AICC macro.for a detailed description                                           */
/*                                                                                             */
/***********************************************************************************************/

%Macro Modelavg(estimate=, Akaike=, DataOut=, model_dimension=);

	proc sort data=&Estimate nodupkey out=allcombos(keep=parameter level1-level&model_dimension);
	  by parameter level1-level&model_dimension;
	  run; 

	 proc sort data=&estimate;
	   by model;
	 run;

	 data estimator;
	   retain modlnum 0;
	   set &Estimate;
	   by model;
	   if first.model then modlnum+1;
	 run;

	 proc sort data=estimator;
           by modlnum parameter level1-level&model_dimension;
	 run;

	 proc univariate noprint data=estimator;
	   var modlnum;
	   output out=nmodels max=nmodels;
	   run;
	 
	 data allcombos2;
	   retain nmodels;
	   if _n_=1 then set nmodels;
	   set allcombos;
	   do modlnum=1 to nmodels;
	     output;
	   end;
	 run;

	proc sort data=allcombos2;
	  by modlnum parameter level1-level&model_dimension;
	run;

	data estimate2;
	   merge allcombos2(in=inall) estimator(in=insome);
	   by modlnum parameter level1-level&model_dimension;
	   if inall and not insome then do;
	     estimate=0;  StdErr=0;
	   end;
	   drop model obsused;
	run;

	proc sort data=&Akaike out=w1(keep=model w_akaike);
	  by model;
	run;

	data w2;
	  retain modlnum 0;
	  set w1;
	  by model;
	  if first.model then modlnum+1;
	run;

	data estimate3;
	  merge w2 estimate2;
	  by modlnum;
	run;

	/* Compute model-averaged parameter estimates (estimate_model_avgd)    */
	/* according to equation 4.1 of Burnham and Anderson (2002, page 150). */

	proc summary missing nway data=estimate3;
	 class parameter level1-level&model_dimension;
	 var estimate;
	 weight  w_Akaike;
	 output out=model_avg(drop=_freq_ _type_) mean=estimate_model_avgd;
	run;

	proc sort data=estimate3;
	  by parameter level1-level&model_dimension;
	run;

	data estimate4;
	  merge model_avg estimate3;
	  by parameter level1-level&model_dimension;
	  var_piece = w_Akaike * sqrt(stdErr**2 + (estimate-estimate_model_avgd)**2);
	    /* above is partial formula for unconditional variance, equation 4.9 of */
	    /* Burnham and Anderson 2002: page 162.                                 */
	run;
	
	proc summary nway noprint data=estimate4;
	  by parameter level1-level&model_dimension;
	  var var_piece;
	  output out=OutData(drop=_type_ rename=(_freq_=no_models)) sum=UnconditionalSE;
	  id estimate_model_avgd;
	run;
	
	data &Dataout;
	  set outdata;
	  *Unconditional_Variance=UnconditionalSE**2;
	  l95cl = estimate_model_avgd - 2*unconditionalSE;
	  u95cl = estimate_model_avgd + 2*unconditionalSE;
	run;

	proc print data=&DataOut; 
	  title 'Model-averaged parameter estimates along with unconditional standard errors and confidence intervals';
    run;

%Mend Modelavg;