/********************* uvtubes.sas ***********************/
title 'Fungus Tube data in univariate format'; 

data mould1;
     infile '/folders/myfolders/tubes.data.txt' firstobs=2  missover; 
      /* missover makes data omitted at the end of a line missing. */
     input line mcg rep day  AML AMS AMld   PML PMS PMld   AMslp PMslp SWeight; 
     label mcg     = 'Fungus Type'
           rep     = 'Replication (1-4)'
           AML     = 'AM Length'    AMS = 'AM Number of Sclerotia'  AMld = 'AM Leading Edge'
           PML     = 'PM Length'    PMS = 'PM Number of Sclerotia'  PMld = 'PM Leading Edge'
           AMslp   = 'AM Least-squares Slope'
           PMslp   = 'PM Least-squares Slope'
           SWeight = 'Slerotial Weight';
     Tube = 1000*rep + mcg; /* 4-digit tube identifier */
     Length = (AML+PML)/2;
     Sclerotia = (AMS+PMS)/2; 
     /* Zero sclerotia were recorded as missing. */
     if Sclerotia=. then Sclerotia=0;
     Slope = (AMslp+PMslp)/2;
     if Tube ne 1213; /* Getting rid of the outlier, based on science. */
     
/* Fungus grew past the end of some tubes after day 10 so it was really a 10-day
   experiment, but day 14 has slopes and weight of sclerotia. These variables will
   still be available in mould1, with one non-missing value for each tube. */

data mould2;
     set mould1;
     if day < 11;

proc print;
     var Tube mcg rep day Length Sclerotia;

proc freq;
     title2 'Sample sizes';
     tables day*mcg / norow nocol nopercent missing;

proc tabulate;
     title2 'Mean Length';
     class day mcg;
     var Length;
     table (day all),(mcg all) * (mean*Length);

proc glm;
     title2 'Proc glm just for the mean plot';
     class day mcg;
     model Length = day|mcg;

proc mixed cl;
     title2 'Unstructured Covariance Matrix';
     class day mcg;
     model Length = day|mcg;
     repeated / type=un subject=tube;
     lsmeans mcg / pdiff adjust=bon;
/* Test linear growth starting with day 2. H0 is
                          mu1  mu2  mu3  mu4  mu5  mu6  mu7  mu8  mu9  mu10
   mu2-mu3 = mu3-mu4  <=>  0    1   -2    1    0    0    0    0    0    0
   mu3-mu4 = mu4-mu5  <=>  0    0    1   -2    1    0    0    0    0    0
   mu4-mu5 = mu5-mu6  <=>  0    0    0    1   -2    1    0    0    0    0
   mu5-mu6 = mu6-mu7  <=>  0    0    0    0    1   -2    1    0    0    0
   mu6-mu7 = mu7-mu8  <=>  0    0    0    0    0    1   -2    1    0    0
   mu7-mu8 = mu8-mu9  <=>  0    0    0    0    0    0    1   -2    1    0
   mu8-mu9 = mu9-mu10 <=>  0    0    0    0    0    0    0    1   -2    1
*/
     contrast 'LinearAfter2' 
                       day 0    1   -2    1    0    0    0    0    0    0,
                       day 0    0    1   -2    1    0    0    0    0    0,
                       day 0    0    0    1   -2    1    0    0    0    0,
                       day 0    0    0    0    1   -2    1    0    0    0,
                       day 0    0    0    0    0    1   -2    1    0    0,
                       day 0    0    0    0    0    0    1   -2    1    0,
                       day 0    0    0    0    0    0    0    1   -2    1;

/* Finally, test whether an ar(1) is good enough */

proc mixed method=ml;
     title2 'Full model for Sigma is unstructured';
     class day mcg;
     model Length = day|mcg;
     repeated / type=un subject=tube;

proc mixed method=ml;
     title2 'Reduced model for Sigma is AR(1)';
     class day mcg;
     model Length = day|mcg;
     repeated / type=ar(1) subject=tube;
     
proc iml;
     title2 'Difference between chi-squared values is chi-squared';
     Chisquare = 382.99-219.86; df=54-1; pvalue = 1-probchi(Chisquare,df);
     print Chisquare df pvalue;