/* mvtubes.sas */

title 'Multivariate tests on Little Fungus Tube data';

data mould;
     infile '/home/u1407221/441s24/data/LittleTubes2.data.txt' firstobs=2; 
     input line id mcg length10 sclr10 weight rate;
     label mcg      = 'Mycelial Compatibility Group'
           length10 = 'Length of fungus in tube at day 10'
           sclr10   = 'Number of sclerotia at day 10'
           weight   = 'Sclerotial weight at end of experiment' 
           rate     = 'Linear growth rate';
     /* Make indicator dummy variables using arrays */
     /* v1-v6 will have these values for every case. */
     v1=198; v2=205; v3=213; v4=221; v5=223; v6=225;
     array dummy{6} mcg198 mcg205 mcg213 mcg221 mcg223 mcg225;
     array value{6} v1-v6;
     do j = 1 to 6;
        if mcg = . then dummy{j} = .;
        else if mcg = value{j} then dummy{j}=1;
        else dummy{j}=0;
     end; 
     /* Effect coding */
     array fungus{5} fungus1-fungus5; 
     do i = 1 to 5;
        if mcg=225 then fungus{i} = -1;
        else fungus{i}=dummy{i};
     end;
/*
proc freq;
     title2 'Check dummy variable creation';
     tables (mcg198--mcg225 fungus1-fungus5)*mcg / norow nocol nopercent missing;
*/

proc means;
     title2 'Potential response variables broken down by MCG';
     class mcg;
     var length10 sclr10 weight rate;

/* Using just length10 and weight to limit the output */

proc glm;
     title2 'Manova with proc glm';
     class mcg;
     model length10 weight = mcg;
     /* Test equality of mcgs excluding 198: a COLLECTION of contrasts */
     contrast 'AllBut198'      mcg   0   1  -1   0   0   0,
                               mcg   0   0   1  -1   0   0,
                               mcg   0   0   0   1  -1   0,
                               mcg   0   0   0   0   1  -1;
     manova H = _all_; /* Test all main effects and interactions */
     means mcg;
     /* lsmeans output is strictly univariate */
  
proc reg plots=none;
     title2 'Manova with proc reg and effect coding;';
     model length10 weight = fungus1-fungus5;
     overall: mtest fungus1=fungus2=fungus3=fungus4=fungus5=0;

proc reg plots=none;
     title2 'Manova with proc reg and cell means coding';
     model length10 weight = mcg198 mcg205 mcg213 mcg221 mcg223 mcg225 / noint;
     Overall:     test mcg198=mcg205=mcg213=mcg221=mcg223=mcg225; 
     mOverall:   mtest mcg198=mcg205=mcg213=mcg221=mcg223=mcg225; 
     AllBut198:   test mcg205=mcg213=mcg221=mcg223=mcg225;
     mAllBut198: mtest mcg205=mcg213=mcg221=mcg223=mcg225;
run;
quit;