Libname quirk 'k:\sboker\temp\GSA2001\';

proc iml;
   seed = 33000;
   er = j(80,3);
   age = j(80,1);
   do i = 1 to 80;
      do j = 1 to 3;
         er[i,j] = rannor(seed);
      end;
      age[i,1] = 20*ranuni(seed)+65;
   end;
   rand = er||age;
   vars= {er1 er2 er3 age};
   create quirk.exp1 from rand [colname=vars];
   append from rand;
run;

data quirk.exp1b;
  set quirk.exp1;
  id=_n_;
  age=round(age);
  beta0 = 50 + (10 * er1);
  beta1 = .5 + er2;
  trial=1;
  output;
  trial=2;
  output;
run;


data quirk.exp1data;
  set quirk.exp1b;
  agec=age-75;
  y = beta0 + (beta1*trial) + (4*er3);
  y=round(y);
  agecat=0;
  if age > 75 then agecat=1;
run;

proc print data=quirk.exp1data;
 var id trial y age agec beta0 beta1 agecat;

run;


title1 '***********************************************';
title2 ' Univariate statistics for EXP1DATA            ';
title3 '***********************************************';

PROC MEANS DATA=quirk.exp1data N MEAN STD MIN MAX;
RUN;


PROC SORT DATA=quirk.exp1data;
    BY age;
RUN;

PROC MEANS DATA=quirk.exp1data N MEAN STD MIN MAX;
    BY agecat;
RUN;

title1 '***********************************************';
title2 ' PROC GLM for EXP1DATA                         ';
title3 '***********************************************';

PROC GLM DATA=quirk.exp1data;
     MODEL y = trial;
RUN;

title1 '***********************************************';
title2 ' PROC MIXED for EXP1DATA                       ';
title3 '***********************************************';

PROC MIXED DATA=quirk.exp1data;
  MODEL y=trial / S;
  RANDOM int / SUB=id TYPE=un gcorr;
RUN;

title1 '***********************************************';
title2 ' PROC GLM for EXP1DATA  (with interaction)     ';
title3 '***********************************************';


PROC GLM DATA=quirk.exp1data;
     MODEL y = trial agec trial*agec;
RUN;

title1 '***********************************************';
title2 ' PROC MIXED for EXP1DATA  (with interaction)   ';
title3 '***********************************************';

PROC MIXED DATA=quirk.exp1data;
  MODEL y=trial agec trial*agec / S;
  RANDOM int agec / SUB=id TYPE=un gcorr;
RUN;