/* calculate model selection statistics for specified models */
/* to choose appropriate correlation structure */

options formdlim = '-';

/* read the data file */
data asparagus;
  infile 'asparagus.txt' firstobs = 2;
  input trt year block yield;

%macro model(type);
/* the code to fit a proc mixed model goes here */
/* it needs to specify everything except the */
/*   type of correlation structure */
/* do need to explicitly indicate the data set in the */
/*  proc mixed statement */
/* you probably needs to change anything in CAPS */
/* all output is suppressed, so no need for any bells and whistles */

proc mixed data = ASPARAGUS;
  class BLOCK YEAR TRT;
  model YIELD =  BLOCK TRT YEAR YEAR*TRT;
  repeated YEAR /type = &type subject = BLOCK*TRT;
  
  ods output fitstatistics = fit;
  ods listing exclude all;
  
run;
       
data fit2;
  set fit;
  model =  symget('type');
  criterion = substr(descr,1,4);
   
  drop descr;
   
proc datasets;
  append base = allfit data = fit2 force;

run;
%mend;

/* if you want to consider models with an additional */
/*   random effect, add the appropriate proc mixed code in */
/*   this macro */

%macro modelrandom(type);
/* the code to fit a proc mixed model goes here */
/* it needs to specify everything except the */
/*   type of correlation structure */
/* do need to explicitly indicate the data set in the */
/*  proc mixed statement */
/* you probably need to change anything in CAPS */
/* all output is suppressed, so no need for any bells and whistles */

proc mixed data = ASPARAGUS;
  class BLOCK YEAR TRT;
  model YIELD =  BLOCK TRT YEAR YEAR*TRT;
  repeated YEAR /type = &type subject = BLOCK*TRT;
  random BLOCK*TRT;
  
  ods output fitstatistics = fit;
  ods listing exclude all;
  
run;
       
data fit2;
  set fit;
  model = trim(symget('type'))||'+RE';
  criterion = substr(descr,1,4);
  drop descr;
  
  
proc datasets;
  append base = allfit data = fit2 force;

run;
%mend;

%macro modelindep;

/* the code to fit a proc mixed model goes here */
/* it needs to specify everything except the */
/*   type of correlation structure */
/* do need to explicitly indicate the data set in the */
/*  proc mixed statement */
/* all output is suppressed, so no need for any bells and whistles */
/* you probably needs to change anything in CAPS */

proc mixed data = ASPARAGUS;
  class BLOCK YEAR TRT;
  model YIELD =  BLOCK TRT YEAR YEAR*TRT;
              
  ods output fitstatistics = fit;
  ods listing exclude all;

data fit2;
  set fit;
  model = 'Indep';
  criterion = substr(descr,1,4);
  drop descr;
    
proc datasets;
  append base = allfit data = fit2 force;

run;
%mend;

proc datasets;
  delete allfit;
  run;

/* then add one line per desired correlation model */
/*  current code does not fit models with additional */
/*  random term */
/* unfortunately, can't do the independence model */
/*  at least without adding more code*/

%model(cs);
%model(ar(1));
%model(arh(1));
%model(ante(1));
%model(un);


/* if want models with additional random effects */
/*  invoke them now */
%modelrandom(ar(1));
%modelrandom(ante(1));

/* if want independence model, run that too */
%modelindep;


/* here's the final reorganization and printing of results */
/* you shouldn't need to change anything below here */

proc sort data = allfit;
  by model criterion;

/* make a nice table */
proc transpose data = allfit out = tfit;
  by model;
  id criterion;
  var value;
run;

data tfit2;
  set tfit;
  Nparam = (AIC - _2_R)/2;
              
/* and print it */     
proc print data = tfit2;
  ods listing select all;
  id model;
  var Nparam aic aicc bic ;
  
  title 'Fit criteria for different correlation models';
  run;