/**********************************************************************
 Analysis using Proc Mixed
***********************************************************************/
data veneer;
  infile "C:\temp\veneer.dat" firstobs=2 dlm="09"X;
  input patient tooth age base_gcf cda time gcf;
run;

/*Plots*/
title "Panels ordered by Age";
proc sgpanel data=veneer;
  panelby age / columns = 4;
  series y=gcf x=time/ group=tooth;
run;

title "Panels ordered by Patient Number";
proc sgpanel data=veneer;
  panelby patient / columns = 4;
  series y=gcf x=time / group=tooth;
run;

/*******************************************************************
  Models
********************************************************************/
title "Model 7.1";
proc mixed data = veneer noclprint covtest;
   class patient tooth;
   model gcf = time base_gcf cda age 
               time*base_gcf  time*cda time*age / solution;
   random intercept time / subject = patient type = un solution;
   random intercept / subject = tooth(patient) g;
run;

title "Model 7.1A";
proc mixed data = veneer noclprint covtest;
   class patient tooth;
   model gcf = time base_gcf cda age 
                     time*base_gcf  time*cda time*age / solution;
   random intercept time / subject = patient type = un solution;
   *random intercept / subject = tooth(patient) g;
run;

title "P-Value for Hypothesis 7.1: Model 7.1A vs 7.1";
data _null_;
  lrtstat = 858.3 - 847.1;
  df = 1;
  p_value = 0.5*(1-probchi(lrtstat,df));
  format p_value 10.8;
  put "Hypothesis 7.1:  "  lrtstat=   p_value= ;
run;


/***************************************************************
  Model the Residual Covariance Structure (Rij)
  Create new variable: CATTIME = Time 
  because need a categorical  time variable 
  for the "repeated" statement.
****************************************************************/
data veneer;
   set veneer;
   cattime = time;
run;

/*This model has an aliasing problem*/
title "Model 7.2A";
proc mixed data = veneer noclprint covtest;
   class patient tooth cattime;
   model gcf = time base_gcf cda age 
               time*base_gcf  time*cda time*age / solution;
   random intercept time / subject = patient type = un solution;
   random intercept / subject = tooth(patient) g;
   repeated cattime / subject = tooth(patient) type = un r rcorr;
run;

/*This model has an aliasing problem*/
title "Model 7.2B";
proc mixed data = veneer noclprint covtest;
   class patient tooth cattime;
   model gcf = time base_gcf cda age 
               time*base_gcf  time*cda time*age / solution;
   random intercept time / subject = patient type = un solution;
   random intercept / subject = tooth(patient) g;
   repeated cattime / subject = tooth(patient) type = cs r rcorr;
run;

/*This model does not have an aliasing problem*/
title "Model 7.2C";
proc mixed data = veneer noclprint covtest;
   class patient tooth cattime;
   model gcf = time base_gcf cda age 
               time*base_gcf  time*cda time*age / solution;
   random intercept time / subject = patient type = un solution;
   random intercept / subject = tooth(patient) g;
   repeated cattime / subject = tooth(patient) type = toeph(1) r rcorr;
run;

title "P-Value for Hypothesis 7.2: Model 7.1 vs 7.2C";
data _null_;
  lrtstat = 847.1 - 846.2;
  df = 1;
  p_value =(1-probchi(lrtstat,df));
  format p_value 10.8;
  put "Hypothesis 7.2:  "  df=  p_value= ;
run;


/****************************************************************
  Reduce Model by Removing non-significant fixed effects
  Use method = ML to estimate the models.
*****************************************************************/ 
title "Model 7.1: ML Estimation";
proc mixed data = veneer noclprint covtest method = ml;
   class patient tooth;
   model gcf = time base_gcf cda age 
               time*base_gcf  time*cda time*age / solution;
   random intercept time / subject = patient type = un solution;
   random intercept / subject = tooth(patient) g;
run;

title "Model 7.3: ML Estimation";
proc mixed data = veneer noclprint covtest method = ml;
   class patient tooth;
   model gcf = time base_gcf cda age / solution;
   random intercept time / subject = patient type = un solution;
   random intercept / subject = tooth(patient) g;
run;

title "P-Value for Hypothesis 7.3: Remove 3 Interactions with Time";
data _null_;
  lrtstat = 845.5 - 843.7;
  df = 3;
  p_value = 1 - probchi(lrtstat,df);
  format p_value 10.8;
  put "Hypothesis 7.3:  " lrtstat=  df=  p_value= ;
run;

/*******************************************************************
  Model 7.3
********************************************************************/
proc sort data = veneer;
  by patient tooth time;
run;

title "Model 7.3: REML Estimation";
proc mixed data = veneer noclprint covtest 
         plots = (residualpanel(blup) studentpanel(blup));
   class patient tooth;
   model gcf = time base_gcf cda age / solution outpred = resids residual;
   random intercept time / subject = patient type = un solution v = 1 vcorr = 1;
   random intercept / subject = tooth(patient) g;
   ods output solutionR = eblupdat; 
run;

/*Different Sorting order to make Vi Matrix more readable*/
proc sort data = veneer;
  by patient time tooth;
run;

ods exclude solutionR;
title "Model 7.3: REML Estimation";
proc mixed data = veneer noclprint covtest 
         plots = (residualpanel(blup) studentpanel(blup));
   class patient tooth;
   model gcf = time base_gcf cda age / solution outpred = resids residual;
   random intercept time / subject = patient type = un solution v = 1 vcorr = 1;
   random intercept / subject = tooth(patient) g;
   ods output solutionR = eblupdat; 
run;

proc print data=eblupdat(obs=10);
run;
/*Note: The eblupdat data set uses ._ as the missing value indicator*/
title "Eblups for Random Intercepts for Patients";
proc univariate data=eblupdat(where =(effect="Intercept"
                     and patient not=._ and tooth = ._)) normal;
   var estimate;
   histogram ;
   qqplot / normal(mu=est sigma=est);
run;

title "Eblups for Random Time Slope for Patients";
proc univariate data=eblupdat(where =(effect="time")) normal;
   var estimate;
   histogram ;
   qqplot / normal(mu=est sigma=est);
run;

title "Eblups for Random Intercept for Teeth Nested Within Patients";
proc univariate data=eblupdat(where =(effect="Intercept"
                     and patient not=._ and tooth not= ._)) normal;
   var estimate;
   histogram ;
   qqplot / normal(mu=est sigma=est);
run;

/*Model 7.3A,same results as Model 7.3 but specified differently*/
title "Alternative Model 7.3A";
proc mixed data = veneer noclprint covtest 
         plots = (residualpanel(blup) studentpanel(blup));
   class patient tooth cattime;
   model gcf = time base_gcf cda age / solution outpred = resids residual;
   random intercept time / subject = patient type = un solution v = 1 vcorr = 1;
   repeated cattime / subject = tooth(patient) type=cs;
run;