/******************************************************
 Analysis using Proc Glimmix
*******************************************************/
data veneer;
  infile "C:\temp\veneer.dat" firstobs=2 dlm="09"X;
  input patient tooth age base_gcf cda time gcf;
run;
proc sort data = veneer;
  by patient tooth time;
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;

/*************************************************************
 Model 7.1: loaded mean structure. 
            Residuals have default 
            covariance matrix (Ri) 
            structure = sigma_squared I.
**************************************************************/
title "Model 7.1";
proc glimmix data = veneer noclprint ;
   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 g;
   random intercept / subject = tooth(patient) g;
   covtest "Hypothesis 7.1" general 0 0 0 1 0 / est;
run;


/***************************************************************
  Model the Residual Covariance Structure (Ri)
  Create new variable: CATTIME = Time 
  because need a categorical  time variable 
  for the R-side Random statement.
****************************************************************/
data veneer;
   set veneer;
   cattime = time;
run;

/*No warnings in Log.

Covariance Parameter Estimates 
Cov Parm Subject Estimate Standard
Error 
UN(1,1) patient 546.62 279.35 
UN(2,1) patient -148.65 74.3885 
UN(2,2) patient 44.6431 21.1000 
Intercept tooth(patient) 7.7465 16.5274 
UN(1,1) tooth(patient) 101.55 18.8107 
UN(2,1) tooth(patient) 39.1689 .      <-------No S.E.
UN(2,2) tooth(patient) 76.1174 15.3007 
*/
title "Model 7.2A: Unstructured R Matrix";
proc glimmix data = veneer noclprint;
   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 g;
   random intercept / subject = tooth(patient) g;
   random cattime / subject = tooth(patient) type = un residual ;
run;

/*NOTE: At least one element of the gradient is greater than 1e-3.

Covariance Parameter Estimates 
Cov Parm Subject Estimate Standard
Error 
UN(1,1) patient 555.41 279.77 
UN(2,1) patient -149.77 74.5492 
UN(2,2) patient 44.7250 21.1510 
Intercept tooth(patient) 40.8338 16.6689 
CS tooth(patient) 6.1272 .        <-------No S.E.
Residual   49.6887 10.9217 
*/
title "Model 7.2B: Compound Symmetric R Matrix";
proc glimmix data = veneer noclprint;
   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;
   random intercept / subject = tooth(patient) solution;
   random cattime / subject = tooth(patient) type = cs residual;
run;

/*Test Hypothesis 7.2: Constant Residual Variance*/
/*All covariance params and S.E. are computed.

Covariance Parameter Estimates 
Cov Parm Subject Estimate Standard
Error 
UN(1,1) patient 546.61 279.34 
UN(2,1) patient -148.64 74.3846 
UN(2,2) patient 44.6418 21.0988 
Intercept tooth(patient) 46.9150 16.5273 
Var(1) tooth(patient) 62.3778 18.8108 
Var(2) tooth(patient) 36.9485 15.3007 
*/
title "Model 7.2C: Heterogeneous Variance, Toeph(1), R Matrix";
proc glimmix data = veneer noclprint;
   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 
	v = 1 vcorr = 1;
   random intercept / subject = tooth(patient) ;
   random cattime / subject = tooth(patient) type = toeph(1) residual;
   covtest "hypothesis 7.2" general 0 0 0 0 1 -1 / est;
run;

/****************************************************************
  Reduce Model by Removing non-significant fixed effects
*****************************************************************/
 
title "Model 7.1: ML Estimation";
proc glimmix data = veneer noclprint noreml;
   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 
	v = 1 vcorr = 1 g;
   random intercept / subject = tooth(patient) g;
run;

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

title "P-Value for Hypothesis 7.3: Remove 3 FE parameters";
data _null_;
  p_value = 1 - probchi(1.84,3);
  format p_value 10.8;
  put "Hypothesis 7.3:  "  p_value= ;
run;

/*******************************************************************
  Final Model 7.3
********************************************************************/
title "Model 7.3: REML Estimation";
ods exclude solutionr;
ods output solutionr=eblupdat;
proc glimmix data = veneer noclprint
    plots =(residualpanel(blup) studentpanel(blup) boxplot(student blup));
 class patient cattime tooth;
 model gcf = time base_gcf cda age / solution;
 random intercept time / subject = patient type = un  
	                     v=1 vcorr=1 solution cl g;
 random intercept / subject = tooth(patient) solution cl g;
 output out = pdat 
   pred(blup)=p_condit resid(blup)=r_condit student(blup)=s_condit
   pred(noblup)=p_marg resid(noblup)=r_marg  student(noblup)=s_marg; 
run;

/***************************************************************************
  Check distribution of EBLUPs
****************************************************************************/ 
title "Eblups for Random Patient Intercepts";
proc univariate data=eblupdat(where =(effect="Intercept"
                     and substr(subject,1,7)="patient")) normal;
   var estimate;
   histogram ;
   qqplot / normal(mu=est sigma=est);
run;

title "Eblups for Random Time Effects";
proc univariate data=eblupdat(where =(effect="time"
                     and substr(subject,1,7)="patient")) normal;
   var estimate;
   histogram ;
   qqplot / normal(mu=est sigma=est);
run;

title "Eblups for Random Tooth Effects";
proc univariate data=eblupdat(where =(effect="Intercept" 
      and substr(subject,1,14)="tooth(patient)")) normal;
   var estimate;
   histogram ;
   qqplot / normal(mu=est sigma=est);
run;

/***************************************************************************
  Model 7.3 Alternative Approach
****************************************************************************/ 
title "Model 7.3 Alternate syntax";
proc glimmix data = veneer noclprint;
   class patient tooth cattime;
   model gcf = time base_gcf cda age / solution;
   random intercept time / subject = patient type = un 
	                       v = 1 vcorr = 1 solution ;
   random cattime / subject = tooth(patient) type=cs residual g;
run;