/*
Peter H. Van Ness
January 12, 2009
Purpose:  	This program provides SAS code for estimating an intraclass correlation coefficient using PROC MIXED
			and then calculating a 90% confidence interval using PROC NLMIXED. This
			program and corresponding data serve for demonstrational purposes only. 

   
   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
   
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
   
   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/



*Intraclass Correlation Analysis;

*Twenty subjects are tested twice for a response whose range runs from 1 to 6.
 The data is fabricated but might be the sort of data generated by a small
 test-retest reliability study;

data icc;
      input subject time response @@;
	  datalines; 
    1 1 0   1 2 1   2 1 0   2 2 2 
    3 1 1   3 2 3   4 1 2   4 2 2 
    5 1 2   5 2 3   6 1 2   6 2 2 
    7 1 5   7 2 6   8 1 6   8 2 6 
    9 1 5   9 2 5  10 1 5  10 2 5 
   11 1 2  11 2 3  12 1 1  12 2 1 
   13 1 0  13 2 4  14 1 0  14 2 0 
   15 1 1  15 2 1  16 1 2  16 2 2 
   17 1 4  17 2 4  18 1 2  18 2 2 
   19 1 3  19 2 5  20 1 6  20 2 5 
   ;
run;

proc print data=icc; run;

proc freq data=icc; tables subject time response; run;

*The off diagonal value from the correlation matrix is the intraclass correlation coefficient.
 It is equal to the estimate for the covariance parameter for the intercept 2.7994 divided by the
 the sum of this value and the estimate for the covariance parameter of the residual 0.8250:
 2.7994/(2.7994+0.8250) = 0.7724;


proc sort data=icc; by subject; run;

proc mixed data=icc method=ml cl covtest;
	title	'Random-Effects ANOVA MIXED Model';
	model response =  / solution ddfm=bw;
	random intercept / subject = subject type=vc v vcorr;
run; 


*While confidence intervals are provided for the off diagonal variance estimate they are not provided for the
 off diagonal value from the correlation matrix. NLMIXED can provide the appropriate confidence interval;


*The intercept value from the MIXED model is used as the initial intercept value for the NLMIXED model;
*The ESTIMATE statement calculates the intraclass correlation coefficient and provides a confidence interval
 according to the significance level indicated in the 'alpha' option;
*Here an alpha level of 0.10 is used to generate a 90% confidence interval
because one is really only interested in the lower bound being greater than 
some acceptable level acording to a Type I error level of 0.05;


proc nlmixed data=icc /*corr cov*/;
	title	'Random-Effects ANOVA NLMIXED Model';
	parms beta0 = 2.875 g11 = 0 to 3 by 0.5 s2 = 1;
	eta = beta0 + b1;
	model response ~ normal(eta, s2);
	random b1 ~ normal(0,g11) subject=subject;
	estimate 'ICC' g11/(g11 + s2) alpha=0.10;
run;

*Note that g11 is equal to the Intercept covariance parameter estimate from the MIXED model: 
 both represent the between-groups variance.
 Also, s2 is equal to the Residual covariance parameter estimate from the MIXED model:
 both represent the within-group variance;