libname focus_mi "C:\Users\tneilands\Box Sync\My Documents\CAPS\Methods Core\Presentations\Missing Data 2015\Part 2\Example 2"; 

ods rtf style = journal file = "C:\Users\tneilands\Box Sync\My Documents\CAPS\Methods Core\Presentations\Missing Data 2015\Part 2\Example 2\Duncan_MI.rtf" ; 

** --> Generate multiple imputations <-- ** ; 

proc sort data = focus_mi.wide_demo out = sorted ;
     by randstat ;
run ; 
 
ods graphics on ;
proc mi data = sorted nimpute = 25 out = outmi seed = 998129000 ;
     var ln_sumsymarv1 ln_sumsymarv2 ln_sumsymarv3 avgbother0 avgbother3 avgbother6 ;
     fcs plots = trace
        nbiter = 20
        regpmm (/details);
     by randstat ;
run ;
ods graphics off ; 
 
** --> Tranpose data to long format required by PROC MIXED <-- ** ; 
 
data dmv_mi ;
     set outmi ;
     array mvar(2,3) ln_sumsymarv1 ln_sumsymarv2 ln_sumsymarv3 avgbother0 avgbother3 avgbother6 ;
     do measure = 1 to 2 ;
          do time = 1 to 3 ;
               arv_sym_both = mvar(measure, time) ;
               output ;
          end ;
     end ;
     keep arv_sym_both subject randstat measure time _imputation_ ; run ; 
 
proc sort data = dmv_mi out = sorted4mixed ;
     by _imputation_ ;
run ; 
 
** --> Perform doubly-mv analyses on the long data retaining cases with missing data <-- ** ; 

ods exlude all ; 
proc mixed data = sorted4mixed namelen=100;
     title "Focus: Doubly-Multivariate Analysis: ARV symptoms and bother" ;
     title2 "All cases" ;
     class measure randstat time subject ;
     model arv_sym_both = measure|randstat|time / solution covb ;
     repeated measure time / subject = subject type = un@un ;
	 slice randstat*measure*time / sliceby = measure*time diff noF ;  
	 lsmeans randstat*measure*time ; 
     by _Imputation_;
	 ods output SolutionF=mixparms CovB=mixcovb ; 
	 ods output slicediffs=sliced ; 
	 ods output lsmeans=lsms ; 
run ; 
ods exclude none ; 

** --> Combine datasets and obtain inference for the randstat*time joint effect <-- ** ; 

** --> First drop rows in the output PARMS and MIXCOVB data sets b/c MIANALYZE won't work with overparameterized design matrices <-- ** ; 

data mixparms;
   set mixparms;
   if (StdErr ~= .);
run;
 
data mixcovb;
   set mixcovb;
   if (measure  ~= 2);
   if (randstat ~= 1);
   if (time     ~= 3);
run;
 
** --> Obtain the inference for the randstat*time joint effect <-- ** ; 

proc mianalyze parms(classvar=full)=mixparms covb(effectvar=rowcol)=mixcovb mult; 
    title "Focus: Doubly-Multivariate Analysis: ARV symptoms and bother" ;
    title2 "All cases" ;
	title3 "Test of interaction of MBSR group and time"; 
	class randstat time ; 
	modeleffects randstat*time ; 
run;  

** --> Combine data sets and obtain inferences for group comparisons within time and measure <-- ** ; 

** --> First re-arrange the slices output into a form PROC MIANALYZE can use <-- ** ; 

data sliced2;
	set sliced;
	comparison=randstat||' vs '||left(_randstat); 
run;
 
proc sort data=sliced2;
	by slice comparison _imputation_;
run;
 
proc mianalyze data=sliced2;
    title "Focus: Doubly-Multivariate Analysis: ARV symptoms and bother" ;
    title2 "All cases" ;
	title3 "Tests of MBSR group differences within each time point"; 
	by slice comparison;
	modeleffects estimate;
	stderr StdErr;
run; 

** --> Obtain MI-based estimates of LSMeans <-- ** ; 

proc sort data=lsms ;
	by  measure randstat time _imputation_;
run;

proc mianalyze data=lsms;
    title "Focus: Doubly-Multivariate Analysis: ARV symptoms and bother" ;
    title2 "All cases" ;
	title3 "Estimated means for each measure, group, and time point"; 
	by measure randstat time ;
	modeleffects estimate;
	stderr StdErr;
run; 

ods rtf close ;
ods graphics on ;