PROC IMPORT OUT= WORK.a 
            DATAFILE= "D:\e\hpp\VEROEFF\13\AR1 and confounding\p rep sim\CycDesigN.csv" 
            DBMS=CSV REPLACE;
     GETNAMES=YES;
     DATAROW=2; 
RUN;

/*Simulate AR(1) for p-rep design; can also use the same dataset to simulate RCBD with 200 complete blocks of
 nine treatments. Just some recoding needed and randomization within blocks.*/

%let seed=2373;   
%let var_b=1;
%let var_e=1;

data b;
se_b=sqrt(&var_b);
se_e=sqrt(&var_e);
do loc_=1 to 5;
  do block_=1 to 40;
    b=se_b*normal(&seed);
	do plot_=1 to 9;
	  e=se_e*normal(&seed);
	  y=b+e;
	  y=round(10000*y)/10000;
	  output;
	end;
  end;
end;
run;

data c;
merge a b;
by loc_ block_ plot_;
if block_<21 then rep_=1; else rep_=2;
keep loc_ rep_ block_ plot_ treatment_ y;
run;

/*randomization-based analysis*/

proc mixed data=c lognote;
class loc_ rep_ plot_ treatment_;
model y=loc_*rep_ treatment_;
random int/sub=loc_*rep_*block_;
run;

/*randomization-based analysis plus AR(1) instead of nugget*/

proc mixed data=c lognote;
class loc_ rep_ plot_ treatment_;
model y=loc_*rep_ treatment_;
random int/sub=loc_*rep_*block_;
repeated plot_/sub=loc_*rep_*block_ type=sp(pow) (plot_);
parms (1)(.1 to .9 by .1)(1)/lowerb=0,0,0;
run;

/*randomization-based analysis plus AR(1) with nugget*/
proc mixed data=c lognote asycorr;
class loc_ rep_ plot_ treatment_;
model y=loc_*rep_ treatment_;
random int/sub=loc_*rep_*block_;
repeated plot_/sub=loc_*rep_*block_ type=sp(pow) (plot_) local;
parms (1)(.5)(.1 to .9 by .1 .999)(.5)/lowerb=0,0,0,0;
run;

proc mixed data=c lognote asycorr;
class loc_ rep_ plot_ treatment_;
model y=loc_*rep_ treatment_;
random int/sub=loc_*rep_*block_;
repeated plot_/sub=loc_*rep_*block_ type=sp(exp) (plot_) local;
parms (1)(.5)(.1 1 10)(.5)/lowerb=0,0,0,0;
run;

PROC EXPORT DATA= WORK.C 
            OUTFILE= "D:\e\hpp\VEROEFF\13\AR1 and confounding\p rep sim\C2273.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="Table1"; 
     NEWFILE=YES;
RUN;

/*LV model*/
data lv;
input parm row col1-col9;
datalines;
1 1 8 7 6 5 4 3 2 1 0
1 2 7 8 7 6 5 4 3 2 1
1 3 6 7 8 7 6 5 4 3 2
1 4 5 6 7 8 7 6 5 4 3
1 5 4 5 6 7 8 7 6 5 4
1 6 3 4 5 6 7 8 7 6 5
1 7 2 3 4 5 6 7 8 7 6
1 8 1 2 3 4 5 6 7 8 7
1 9 0 1 2 3 4 5 6 7 8
;
proc mixed data=c lognote asycorr;
class loc_ rep_ plot_ treatment_;
model y=loc_*rep_ treatment_;
random int/sub=loc_*rep_*block_;
repeated plot_/sub=loc_*rep_*block_ type=lin(1) ldata=lv local;
parms (1)(.1 to .9 by .1)(.5)/lowerb=0,0,0,0;
run;