options nodate nonumber;
/*import the data file*/
proc import datafile="d:/barley/nd.csv" out=d replace dbms=csv;
run;

/*look at the first 37 rows of the data*/
proc print data=d (obs=37);
run;

/* Prepare data for mixed model analysis */
data ndat;  set d;
run;

data ndat; set ndat;
  gene= _n_;
  rep=1;
  geno="ML6";
  time=0;
  intensity=V1;
  output;
  gene= _n_;
  rep=1;
  geno="ML6";
  time=8;
  intensity=V2;
  output;
  gene= _n_;
  rep=1;
  geno="ML6";
  time=16;
  intensity=V3;
  output;
  gene= _n_;
  rep=1;
  geno="ML6";
  time=20;
  intensity=V4;
  output;
  gene= _n_;
  rep=1;
  geno="ML6";
  time=24;
  intensity=V5;
  output;
  gene= _n_;
  rep=1;
  geno="ML6";
  time=32;
  intensity=V6;
  output;
  gene= _n_;
  rep=1;
  geno="ML13";
  time=0;
  intensity=V7;
  output;
  gene= _n_;
  rep=1;
  geno="ML13";
  time=8;
  intensity=V8;
  output;
  gene= _n_;
  rep=1;
  geno="ML13";
  time=16;
  intensity=V9;
  output;
  gene= _n_;
  rep=1;
  geno="ML13";
  time=20;
  intensity=V10;
  output;
  gene= _n_;
  rep=1;
  geno="ML13";
  time=24;
  intensity=V11;
  output;
  gene= _n_;
  rep=1;
  geno="ML13";
  time=32;
  intensity=V12;
  output;

/*rep2*/
  gene= _n_;
  rep=2;
  geno="ML6";
  time=0;
  intensity=V13;
  output;
  gene= _n_;
  rep=2;
  geno="ML6";
  time=8;
  intensity=V14;
  output;
  gene= _n_;
  rep=2;
  geno="ML6";
  time=16;
  intensity=V15;
  output;
  gene= _n_;
  rep=2;
  geno="ML6";
  time=20;
  intensity=V16;
  output;
  gene= _n_;
  rep=2;
  geno="ML6";
  time=24;
  intensity=V17;
  output;
  gene= _n_;
  rep=2;
  geno="ML6";
  time=32;
  intensity=V18;
  output;
  gene= _n_;
  rep=2;
  geno="ML13";
  time=0;
  intensity=V19;
  output;
  gene= _n_;
  rep=2;
  geno="ML13";
  time=8;
  intensity=V20;
  output;
  gene= _n_;
  rep=2;
  geno="ML13";
  time=16;
  intensity=V21;
  output;
  gene= _n_;
  rep=2;
  geno="ML13";
  time=20;
  intensity=V22;
  output;
  gene= _n_;
  rep=2;
  geno="ML13";
  time=24;
  intensity=V23;
  output;
  gene= _n_;
  rep=2;
  geno="ML13";
  time=32;
  intensity=V24;
  output; 

/*rep3*/
  gene= _n_; 
  rep=3;
  geno="ML6";
  time=0;
  intensity=V25;
  output;
  gene= _n_;
  rep=3;
  geno="ML6";
  time=8;
  intensity=V26;
  output;
  gene= _n_;
  rep=3;
  geno="ML6";
  time=16;
  intensity=V27;
  output;
  gene= _n_;
  rep=3;
  geno="ML6";
  time=20;
  intensity=V28;
  output;
  gene= _n_;
  rep=3;
  geno="ML6";
  time=24;
  intensity=V29;
  output;
  gene= _n_;
  rep=3;
  geno="ML6";
  time=32;
  intensity=V30;
  output;
  gene= _n_;
  rep=3;
  geno="ML13";
  time=0;
  intensity=V31;
  output;
  gene= _n_;
  rep=3;
  geno="ML13";
  time=8;
  intensity=V32;
  output;
  gene= _n_;
  rep=3;
  geno="ML13";
  time=16;
  intensity=V33;
  output;
  gene= _n_;
  rep=3;
  geno="ML13";
  time=20;
  intensity=V34;
  output;
  gene= _n_;
  rep=3;
  geno="ML13";
  time=24;
  intensity=V35;
  output;
  gene= _n_;
  rep=3;
  geno="ML13";
  time=32;
  intensity=V36;
  output;
  drop V1-V36;
run; 

proc print data=ndat (obs=36);
run;

/*the data for the first gene*/
data temp;  set ndat;
  if gene=1;
run;

ods trace on; /* turns on the writing of the trace record */
run;

/*fitting mixed linear model for the first gene*/
proc mixed data=temp covtest IC;
  class gene rep geno time intensity; /*classification variables*/
  model intensity=geno time geno*time/DDFM=KENWARDROGER;
  random rep rep*geno rep*time;
  lsmeans geno*time;
run;

ods trace off; /* turns off the writing of the trace record */
run;

ods noresults;
run;

ods listing close; /*do not print in the output window*/
run;

options nomlogic nonotes; /*Using SAS System Options to Suppress Log Output*/

proc mixed data=ndat covtest IC;
  class gene rep geno time intensity;
  by gene;
  ods output InfoCrit=IC;
  ods output Tests3=results;
  ods output CovParms=varcomps;
  ods output lsmeans=means;
  model intensity=geno time geno*time/DDFM=KENWARDROGER outp=resids;
  random rep rep*geno rep*time;
  lsmeans geno*time;
run;

ods listing;  /*to print in the output window*/
run;

proc print data=IC (obs=10);
run;

proc print data=varcomps (obs=10);
run;

proc print data=results (obs=10);
run;

proc print data=means (obs=10);
run;

/*export estimates of the variance components*/
proc export data=varcomps outfile="D:\barley\vars.txt" dbms=TAB replace;
run;

/*export information criteria*/
proc export data=IC outfile="D:\barley\IC.txt" dbms=TAB replace;
run;

/*export p-values for the test of fixed effects, genotype, time, and genotype-time interaction */
data results; set results;
  pval=1*probf;
  drop probf;
run;

proc export data=results outfile="D:\barley\results.txt" dbms=TAB replace;
run;

/*export lsmeans*/
proc export data=means outfile="D:\barley\lsmeans.txt" dbms=TAB replace;
run;

/*export residuals*/
proc export data=resids outfile="D:\barley\resids.txt" dbms=TAB replace;
run;