Libname quirk 'c:\sasdata\';
proc iml;
   seed = 42109;
   u = j(80,3);
   age = j(80,1);
     do i = 1 to 80;
       do j = 1 to 3;
   	   u[i,j] = rannor(seed);
   	end;
     age[i,1] = 20*ranuni(seed)+65;
     end;
   cov = {13 3 0, 3 4 0, 0 0 0.10};
   covcho = root(cov);
   uvec = u * covcho;
   rand = uvec||age;
   vars= {u0 u1 u2 age};
   create quirk.gsafirst from rand [colname=vars];
   append from rand;
run;

data quirk.gsatwo;
  set quirk.gsafirst;
  id=_n_;
  age=round(age);
  beta0 = 115+u0;
  beta1 = -0.25+u1;
  beta2 = -0.08+u2;
  wave=1;
  output;
  wave=2;
  output;
  wave=3;
  output;
  wave=4;
  output;

run;

proc iml;
  seed = 32291;
  e = j(320,1);
  do i=1 to 320;
  e[i,1] = rannor(seed);
  end;
  vars={eit};
  create quirk.gsaerr from e [colname=vars];
  append from e;

run;
  

data quirk.gsadata;
  merge quirk.gsatwo quirk.gsaerr;
  age=age+wave-1;
  y = beta0 + (beta1*(age-75)) + (beta2*((age-75)**2)) + (4*eit);
  y=round(y);
  agesq=age*age;
  if age le 75 then old=0;
      else old=1;

run;

proc print data=quirk.gsadata;
 var y age wave beta0 beta1 beta2 eit u0 u1 u2;

run;



data timeone;
  set quirk.gsadata;
  if wave=1;



proc ttest data=timeone;
  class old;
  var y u0 u1 u2;

proc corr data=timeone;
  var u0 u1 u2;


proc glm data=timeone;
  class old;
  model y=old/nouni;
  estimate 'group difference' old 1 -1;



run;


data quirk.gsamult;
  set quirk.gsadata;
  by id;
  if wave=1 then y1=y;
  if wave=2 then y2=y;
  if wave=3 then y3=y;
  if wave=4 then y4=y;
  retain y1 y2 y3 y4;
  if last.id=1 then output;




proc corr;
  var y1 y2 y3 y4;

run;



run;

data quirk.gsamult3;
  set quirk.gsamult;
  age=age-3;
  lin=y3-y1;
  keep id old age y1 y2 y3 lin;




run;

/* everything up to here is simply creating the data, and doing some checks to see whether the general pattern of results is as desired */


/*actually begin here with my part of workshop  */


proc print data=quirk.gsamult3;
  var age y1 y2 y3 old lin;


proc corr data=quirk.gsamult3;
  var y1 y2 y3;

data quirk.gsauni;
  set quirk.gsamult3;
  wave=1;
  y=y1;
  output;
  wave=2;
  age=age+1;
  y=y2;
  output;
  wave=3;
  age=age+1;
  y=y3;
  output;

data quirk.gsauniv;
  set quirk.gsauni;
  wavezero=wave-1;
  wavcmean=wave-2;
  agec75=age-75;
  agec75sq=agec75*agec75;

proc print data=quirk.gsauniv;
  var id y wave wavezero wavcmean age agec75 agec75sq;

run;

proc reg data=quirk.gsauniv;
  model y = age;

proc reg data=quirk.gsauniv;
  model y = wave;

proc glm data=quirk.gsamult3;
  model y1 y2 y3=/nouni;
  repeated time 3 (0 1 2) polynomial / summary;


proc mixed data=quirk.gsauniv;
  model y=age/s;

proc mixed data=quirk.gsauniv;
  model y=wave/s;


proc mixed data=quirk.gsauniv;
  model y=wave/s;
  random int/sub=id;

proc mixed data=quirk.gsauniv;
  model y=wave/s;
  random int wave/sub=id type=un gcorr;
  

proc mixed data=quirk.gsauniv;
  model y=wavcmean;
  random int wavcmean/sub=id type=un gcorr;

proc mixed data=quirk.gsauniv;
  model y=agec75;
  random int agec75/sub=id type=un gcorr;

proc corr data=quirk.gsamult3;
  var y1 lin;

proc mixed data=quirk.gsauniv;
  class old;
  model y=age old old*age;
  random int/sub=id;

proc mixed data=quirk.gsauniv;
  class old;
  model y=age old old*age/s;
  random int age/sub=id type=un gcorr;
  


proc glm data=quirk.gsamult3;
  class old;
  model y1 y2 y3=old/nouni;
  repeated time 3 (0 1 2) polynomial / summary;
  estimate 'group difference' old 1 -1;


proc ttest data=quirk.gsamult3;
  class old;
  var lin;



proc print data=quirk.gsauniv;
  var id y wavcmean old;

proc mixed data=quirk.gsauniv;
  class old;
  model y=wavcmean old old*wavcmean;
  random int wavcmean/sub=id type=un gcorr;

proc mixed data=quirk.gsauniv;
  class old;
  model y=wavcmean old old*wavcmean/s;
  repeated/sub=id type=un rcorr;


data quirk.gsadata2;
  set quirk.gsadata;
  wavezero=wave-1;
  wavcmean=wave-2;
  agec75=age-75;
  agec75sq=agec75*agec75;



proc mixed data=quirk.gsadata2;
  model y=agec75 agec75sq;
  random int agec75/sub=id type=un gcorr;


proc mixed data=quirk.gsadata2;
  model y=age agec75sq;
  random int agec75 agec75sq/sub=id type=un gcorr;



proc mixed data=quirk.gsadata2;
  model y=age agesq;
  random int age/sub=id type=un gcorr;


proc mixed data=quirk.gsadata2;
  model y=age agesq;
  random int age agesq/sub=id type=un gcorr;

run;



goptions reset=symbol;



symbol1 color=black interpol=join
        width=1 value=none
        height=0 line=1;

symbol2 color=black interpol=join
        width=1 value=none
        height=0 line=1;

symbol3 color=black interpol=join
        width=1 value=none
        height=0 line=1;

symbol4 color=black interpol=join
        width=1 value=none
        height=0 line=1;

symbol5 color=black interpol=join
        width=1 value=none
        height=0 line=1;

symbol6 color=black interpol=join
        width=1 value=none
        height=0 line=1;

symbol7 color=black interpol=join
        width=1 value=none
        height=0 line=1;

symbol8 color=black interpol=join
        width=1 value=none
        height=0 line=1;

symbol9 color=black interpol=join
        width=1 value=none
        height=0 line=1;

symbol10 color=black interpol=join
        width=1 value=none
        height=0 line=1;


data tempmix;
  set quirk.gsauniv;
  if id le 10;


proc gplot data=tempmix;
  plot y*wave=id / nolegend haxis=1 to 3 hminor=0;
run;




quit;