%macro autocor(data=_last_, y=, time=, id=);
/*********************************************************************
This Macro is designed to create a scatterplot matrix of the residuals 
for logitudinal data. In the output window we get the autocorrelation 
matrix and the autocorrelation function under the stationarity assumption.
The first graph window gives the autocorrelation function plot and the 
second graph window gives the scatterplot matrix of residuals. Since 
we make the scatterplot matrix by SAS/Insight software, another data
window will be opened, just disregard it.
There are four parameters:

DATA     is the name of the data set containing logitudinal data.
         The default is the most recently created data set.

Y        is the response variable.

TIME     is the time when the response was measured, later was convert
         to indicator of time point. 

ID       is a person identifier.

Example of usage:

%autocor(data=b,y=dist,time= age, id=id);

Author:  Hongfei Guo, hfguo@jhsph.edu
         Department of Biostatistics 
         Johns Hopkins University  
*************************************************************************/
options LS=80  NODATE ; 
/*READ THE DATA AND RENAME THE VARIABLES*/
data _inset_;
 set &data;
 y=&y; t=&time; id=&id;
 keep y t id;
/*CREATE LAG TIME VARIABLE: NUM */
proc sort data=_inset_; by id t;
data _inset2_; set _inset_; by id;
where y ne . and t ne .;
	if first.id then num=1;
  	retain num;
 	if first.id ne 1 then
        do;num=num+1;end;
con=1;
run;
/*MACRO TO CREATE VARIABLES LIST*/
%macro names(name= ,number= );
   %do n=1 %to &number; &name&n  %end;
%mend names;
/*CALCULATE THE RESIDUALS OF THE REGRESSION*/
proc iml;
use _inset2_;
read all var {y id num}; 
read all var {con t} into X;
res=y-X*inv(X`*X)*X`*y;
/*SAVE TO A NEW DATASET*/
create _inset3_  var {id num res};
append;
close _inset2_;
close _inset3_;
/*FIND THE MAXIMUM NUMBER OF TIME POINT*/
numax=max(num);
create _nm_ var{numax};
append;
close _nm_;
quit;
/*DEFINE THE MAXIMUM NUMBER OF TIME MACRO*/
data _null_; set _nm_;
   call symput('nm',numax);
/*RESHAPE THE RESIDUALS SUCH THAT EACH SUBJECT HAS ONE OBS*/
proc sort data=_inset3_; by id num;
data _inset4_ ; set _inset3_; by id num;
	array tt{*} %names(name=time,number=&nm);
	retain %names(name=time,number=&nm);
	if first.id then 
        do i=1 to &nm; tt{i}= .; end;
	tt{num}= res;		
	if last.id then output;
keep  id %names(name=time,number=&nm);
/*CORRELATION OF THE RESIDUALS*/
proc corr data=_inset4_ noprint outp=_corr_;
	var %names(name=time,number=&nm) ;
data _corrp_; set _corr_;
where _type_ eq 'CORR';
run;
proc iml;
use _corrp_;
read all var{%names(name=time,number=&nm)} into Autocorrelation_matrix;
names = {%names(name=time,number=&nm)};
print Autocorrelation_matrix[rowname=names colname=names];
/*REARRANGE THE RESIDUALS AND CALCULATE THE AUTOCORRELATION FUNCTION*/
use _inset4_;
read all var{%names(name=time,number=&nm)} into X;
mm=nrow(X);
lag=shape(X`,1);
nn=ncol(lag);
lag=lag`;
acf=repeat(.,&nm-1,2);
do i=1 to &nm-1;
	y = repeat(.,nn,1);
	y[(i*mm+1):nn]=lag[1:(nn-i*mm)];
	create _t1_ var{lag y};
	append;
	delete all where (lag =. | y =.);
	read all into acfd;
	acf[i,1]=i;
	acf[i,2]=corr(acfd)[1,2];
	close _t1_;
end;
print acf;
close _inset4_;
/*AUTOCORRELATION FUNCTION PLOT*/
create _acfplot_ from acf[colname={'Lag' 'acf'}];
append from acf;
close _acfplot_;
quit;
goptions reset = all;
axis1 order = (0 to &nm) minor = none;
axis2 order = (0 to 1 by .2) minor=none label=(a=90);
symbol1 i = none v = circle h = .5;
proc gplot data=_acfplot_;
	plot acf*lag/haxis=axis1 vaxis=axis2;
	label acf= "Autocorrelation function";
run;
goptions reset = all;
quit;
/*SCATTER MATRIX PLOTS*/
proc insight data=_inset4_ NOMENU  NOBUTTON NOCONFIRM;
	scatter %names(name=time,number=&nm)*%names(name=time,number=&nm);
run;
quit;
%mend autocor;