data INITIAL; 
 infile "C:\anti.dat";
 input anti1 anti2 anti3 anti4 gen homecog id;
run;

**************************************************************************************************
** Macro variables input section                                                                **
**************************************************************************************************;

%let idvarname = id;
%let varlist = anti1 anti2 anti3 anti4;
%let timepts = 0 1 2 3;
%let inclmiss = n;
%let level = both;
%let regtype = lin;
%let numplot = ;
%let hist = y;
%let outds = ;

**************************************************************************************************
** END of user input section                                                                    **
**************************************************************************************************;

**************************************************************************************************
**** PROGRAM NAME:  OLStraj.sas                                                                 **
**** AUTHOR:        Madeline Carrig                                                             **
**** LAST MODIFIED: October 3, 2002                                                             **
**** FOR FULL DOCUMENTATION:  Visit http://www.unc.edu/~curran/OLStraj.htm                      **
**************************************************************************************************;
dm 'out;clear;log;clear;pgm;recall;';
options nocenter linesize=100 nodate nonumber nomprint nosymbolgen;

%MACRO OLStraj;

  %global _varnum;
  %let _varnum = 0;
  %do %until(%scan(&varlist,%eval(&_varnum+1),' ') eq );
    %let _varnum = %eval(&_varnum+1);
  %end;

  %let _timenum = 0;
  %do %until(%scan(&timepts,%eval(&_timenum+1),' ') eq );
    %let _timenum = %eval(&_timenum+1);
  %end;

  %if &_varnum ne &_timenum %then %do;
    %put +++ERROR: NUMBER OF TIME POINTS DOES NOT EQUAL NUMBER OF REPEATED MEASURES;
    %put +++RERUN PROGRAM AFTER CORRECTING PARAMETER ASSIGNMENTS;
    %goto exit;
  %end;

  %let inclmiss = %upcase(%left(&inclmiss));
  %let level = %upcase(%left(&level));
  %let regtype = %upcase(%left(&regtype));
  %let hist = %upcase(%left(&hist));

  data initial;
    set initial;
    array varlist {&_varnum} &varlist;

    %if &numplot ne %then if _n_ le &numplot;;

    numpres = 0;
    do i = 1 to &_varnum;
      if varlist{i} ne . then numpres = numpres + 1;
    end;

    %if &inclmiss = N %then if numpres = &_varnum then output;;
  run;

  proc sort data = initial;
    by &idvarname;
  run;

  data initial;
    keep &idvarname _newid time score %if &regtype ne LIN %then timesq;;

    set initial end = last;
    retain _newid;

    if _n_ = 1 then _newid = 1;
    else _newid = _newid + 1;

    array varlist {&_varnum} &varlist;

    do i = 1 to &_varnum;
      time = input(scan("&timepts",i,' '),32.);
      %if &regtype ne LIN %then timesq = time**2;;
      score = varlist{i};
      output;
    end;

    %if &level ne IND %then if last then call symput('_numobs',_newid);;
  run;

  proc reg data = initial outest = ahatbhat rsquare;
    by &idvarname;
    id _newid;
    model score = time
      %if &regtype ne LIN %then timesq;
      / noprint;
    output out = olsdat p = y_est;
  run;

  data ahatbhat (keep = &idvarname intercept linear %if &regtype ne LIN %then quad; rsquared);
    set ahatbhat;
    rename time  = linear
           %if &regtype ne LIN %then timesq = quad;
           _rsq_ = rsquared;
  run;

  %if &level ne IND %then %do;
    %let _numobs = %left(&_numobs);

    %do i = 1 %to &_numobs;
      data new (keep = time score&i pred&i);
        set olsdat;
        if _newid = &i;

        rename score = score&i y_est = pred&i;
      run;

      proc sort data = new;
        by time;
      run;

      data all;
        %if &i = 1 %then %do;
          set new;
        %end;
        %else %do;
          merge all
                new;
          by time;
        %end;
      run;
    %end;

    goptions reset=all ftext=swiss interpol=join;

    symbol1;
    symbol2;
    symbol3;

    proc gplot data = all;
      title "Simple Join";
      %if &numplot ne %then %do;
        title2 "First &numplot Observations";
      %end;
      %else %do;
        title2 "All Observations";
      %end;
      plot (score1-score&_numobs) * time / overlay vaxis = axis1 haxis = axis2;
      label time = 'Time';
      axis1 label = ("Score");
      axis2 label = ("Time");
    run;

    %if &regtype ne QUAD %then %do;
      goptions interpol = rl;

      proc gplot data = all;
        title "Estimated OLS Linear Trajectories";
        %if &numplot ne %then %do;
          title2 "First &numplot Observations";
        %end;
        %else %do;
          title2 "All Observations";
        %end;
        plot (pred1-pred&_numobs) * time / overlay vaxis = axis1 haxis = axis2;
        label time = 'Time';
        axis1 label = ("Score");
        axis2 label = ("Time");
      run;
    %end;

    %if &regtype ne LIN %then %do;
      goptions interpol = rq;

      proc gplot data = all;
        title "Estimated OLS Quadratic Trajectories";
        %if &numplot ne %then %do;
          title2 "First &numplot Observations";
        %end;
        %else %do;
          title2 "All Observations";
        %end;
        plot (pred1-pred&_numobs) * time / overlay vaxis = axis1 haxis = axis2;
        label time = 'Time';
        axis1 label = ("Score");
        axis2 label = ("Time");
      run;
    %end;

  %end;

  %if &level ne GRP %then %do;
    options nobyline;
    goptions reset=all ftext=swiss;

    symbol1 color = red
            value = dot
            height = 1;

    symbol2 color = blue width=2 interpol = rl;

    symbol3 color = blue width=2
            %if &regtype ne QUAD %then line=2;
            interpol = rq;

    legend1 frame label=none
            %if &regtype = LIN %then %do;
              value=(tick=1 'Observed Data' tick=2 'OLS Linear Trajectory');
            %end;
            %else %if &regtype = QUAD %then %do;
              value=(tick=1 'Observed Data' tick=2 'OLS Quadratic Trajectory');
            %end;
            %else %do;
              value=(tick=1 'Observed Data' tick=2 'OLS Linear Trajectory'
                     tick=3 'OLS Quadratic Trajectory');
            %end;

    proc gplot data = olsdat nocache uniform;
      by &idvarname;
      title "%upcase(%left(&idvarname)) = #byval1";
      title2;
      plot score * time = 1
           %if &regtype ne QUAD %then score * time = 2;
           %if &regtype ne LIN %then score * time = 3;
           / overlay vaxis  = axis1 haxis = axis2 legend = legend1;
      label time = 'Time';
      axis1 label = ("Score");
      axis2 label = ("Time");
    run;
  %end;

  quit;

  %if &outds ne %then %do;
    data &outds;
      set ahatbhat;
    run;
  %end;

  %if &hist = Y %then %do;
    data ahatbhat;
      keep &idvarname est_grp ols_est;
      length est_grp $ 9;
      set ahatbhat;

      est_grp = 'Intercept';
      ols_est = intercept;
      output;

      est_grp = 'Linear';
      ols_est = linear;
      output;

      %if &regtype ne LIN %then %do;
        est_grp = 'Quadratic';
        ols_est = quad;
        output;
      %end;
    run;

    proc sort data = ahatbhat;
      by est_grp;
    run;

    proc univariate data = ahatbhat noprint;
      by est_grp;
      var ols_est;
      output out = minmax min = minval max = maxval;
    run;

    goptions reset=all ftext=swiss;
    symbol v=dot c=salmon;
    axis1 minor=none color=black label=(angle=90 rotate=0);

    %macro rept(grp,lowpt,hipt,fact);
      proc univariate data = ahatbhat noprint;
        where est_grp = "&grp";
        title "Histogram of OLS Estimates:  Parameter = &grp";
        %if &numplot ne %then %do;
          title2 "First &numplot Observations";
        %end;
        %else %do;
          title2 "All Observations";
        %end;
        histogram ols_est / cfill=salmon cbarline=black vscale=count
                            midpoints = &lowpt to &hipt by &fact;
        label ols_est = "Value";
      run;
    %mend rept;

    data _null_;
      set minmax;
      lowpt = round((minval + (maxval-minval)/36),.01);
      hipt = round((maxval - (maxval-minval)/36),.01);
      fact = round(((maxval-minval)/18),.01);
      call execute('%rept('||est_grp||','||lowpt||','||hipt||','||fact||')');
    run;

    proc boxplot data = ahatbhat;
      plot ols_est * est_grp /
                   boxstyle = schematicid
                   nohlabel
                   cboxes   = dagr
                   cboxfill = ywh
                   cframe   = vligb
                   idcolor  = salmon
                   vaxis    = axis1;
      id &idvarname;
      label ols_est = 'OLS Estimate';
      title 'Box Plots of OLS Estimates';
      %if &numplot ne %then %do;
        title2 "First &numplot Observations";
      %end;
      %else %do;
        title2 "All Observations";
      %end;
      title3 '(Any Extreme Observations are Identified by Value of ID Variable)';
    run;
  %end;

%EXIT:  %MEND OLStraj;

%olstraj;