Dear ,
with the help of Rick Wicklin's book on simulating in SAS, I managed to simulate 1 dataset for a longitudinal analysis with three timepoints, 2 treatment groups and 5 subjects in each treatment group.
Now, I want to simulate 1000 datasets. I found something like "do SampleID = 1 to &NumSamples, then I could run the proc mixed with a by variable, but I don't succeed in completing this.
This is the code for 1 dataset. Any suggestions where I could find information?
data Interactions; WtGain = 0; /* the value of y does not matter */
do Week = 1 to 3;
do Treatment = 1 to 2;
do Indiv = 1 to 5;
output;
end;
end;
end;
run;
proc print data=Interactions;
run;
proc sort data=Interactions;
by Treatment Indiv Week;
run;
proc print data=Interactions;
run;
proc glmmod data = Interactions noprint outparm = Parm outdesign = Design( drop = WtGain); /* DROP y */
class Week Treatment;
model WtGain = Week | Treatment;
run;
proc iml;
call randseed( 1);
use Design;
read all var _NUM_ into X;
close Design; /* Intcpt |--Week--| Treatment |-----Interactions-----| */
beta = {1, -1, -0.5, 0, -0.5, 0, -0.25, -0.25, 0, 0, 0, 0};
k = 3; /* 3 weeks thus 3 timepoints */
g=2; /*2 treatments*/
s1 = 5; /* 5 individuals in each group*/
s = g*s1; /* 10 individuals in total*/
B = { 1 0.5 0.25 ,
0.5 1 0.5 ,
0.25 0.5 1 } ;
print B;
R = B; /* first block */
do i = 2 to s; /* create block-diagonal matrix */
R = block( R, b); /* with s blocks */
end;
print R;
Week = T(repeat(1:k, 1,s)); /* column: 1,2,3,1,2,3,… */
print Week ;
Indiv = colvec(repeat(T(1:s),1,k)); /* 1,1,1,2,2,2,3,3,3,… */
print Indiv ;
Treatment = colvec(repeat(T(1:g),1,s1*k));
print Treatment ;
call randseed( 1234);
create Mix var {"WtGain" "Week" "Treatment" "Indiv"}; /* name the variables */
/* random effects */
zero = j( 1, k*s, 0); /* the zero vector */
eps = RandNormal( 1, zero, R); /* eps ~ MVN( 0, R) */
print eps ;
fixed = X* beta ;
print fixed;
WtGain = fixed + T(eps) ; /* fixed effects, correlated errors */
append;
close;
quit;
* analyze the data;
proc mixed data = Mix CL;
class Indiv Week Treatment;
model WtGain = Week | Treatment/ s;
repeated / type = ar(1) subject = Indiv R;
*ods select Dimensions R CovParms;
run;size
