/* SAS Examples of Analysis of Random Effects Model.     */
/* We analyze the Apex Enterprises data from Chapter 25. */

DATA apex;
  INPUT rating officer candidate;
cards;
  76  1  1
  65  1  2
  85  1  3
  74  1  4
  59  2  1
  75  2  2
  81  2  3
  67  2  4
  49  3  1
  63  3  2
  61  3  3
  46  3  4
  74  4  1
  71  4  2
  85  4  3
  89  4  4
  66  5  1
  84  5  2
  80  5  3
  79  5  4
;
run;

PROC GLM data = apex;
  CLASS officer;
  MODEL rating = officer;
  RANDOM officer;
run;
/* We use the RANDOM statement to tell SAS that "officer" has random levels. */

/* We see there is significant variation in rating among the population */
/* of officers.  (F* = 5.39, P-value = 0.0068.)                         */

/*********************************************************************************/          

/* 90% CI for mu-dot and 90% CI for sigma^2 */

PROC MIXED data = apex method = reml ASYCOV CL COVTEST ALPHA = 0.10;
/* The ASYCOV CL COVTEST ALPHA = 0.10 part        */
/* of the above line gives the 90% CI for sigma^2 */
CLASS officer;
MODEL rating = / CL ALPHA = 0.10; 
/* The above line gives the 90% CI for mu-dot        */
/* (for the mean rating of all personnel officers).  */
RANDOM officer;
run;

/* The 90% CI for sigma^2 is (43.98, 151.39).  (Find this on the SAS output.) */

/* The 90% CI for mu-dot is (61.98, 80.92).  (Find this on the SAS output.) */

/*********************************************************************************/

/* CI for intraclass correlation coefficient and CI for sigma_mu^2 */

DATA ciinfo;
ALPHA = 0.10; /*This gives a 90% CI, can be changed depending on needed confidence level*/
n = 4; /* input number of observations per cell */
r = 5; /* input number of levels of the random factor */
MSTR = 394.925; /* input MSTR from PROC GLM output above */
MSE = 73.28333; /* input MSE from PROC GLM output above */

/* lots of calculations */
fl=finv(1-alpha/2,r-1,r*(n-1)); fu=finv(alpha/2,r-1,r*(n-1));
l=1/n*((mstr/mse)*(1/fl)-1); u=1/n*((mstr/mse)*(1/fu)-1);
LCL_ICC=l/(1+l); UCL_ICC=u/(1+u);
c1=1/n; c2=-1/n; lhat=(mstr-mse)/n;
f1=finv(1-alpha/2,r-1,9999); f2=finv(1-alpha/2,r*(n-1),9999); f3=finv(1-alpha/2,9999,r-1);
f4=finv(1-alpha/2,9999,r*(n-1)); f5=finv(1-alpha/2,r-1,r*(n-1)); f6=finv(1-alpha/2,r*(n-1),r-1);
g1=1-(1/f1); g2= 1-(1/f2); g3=((f5-1)**2-(g1*f5)**2-(f4-1)**2)/f5;
g4=f6*(((f6-1)/f6)**2-(((f3-1)/f6)**2)-g2**2);
hl=(((g1*c1*mstr)**2)+(((f4-1)*c2*mse)**2)-g3*c1*c2*mstr*mse)**.5;
hu=(((f3-1)*c1*mstr)**2+((g2*c2*mse)**2)-(g4*c1*c2*mstr*mse))**.5;
LCL_sigma_mu_sq=lhat-hl; 
UCL_sigma_mu_sq=lhat+hu; run;

/* Printing out the Confidence Limits for the CIs            */
/* for intraclass correlation coefficient and for sigma_mu^2 */
PROC PRINT data=ciinfo;
VAR LCL_ICC UCL_ICC LCL_sigma_mu_sq UCL_sigma_mu_sq;
run;


/* A 90% CI for the intraclass correlation coefficient is (0.160, 0.884).  */
/* Compare to book's result, at top of page 1041. */
/* A 90% CI for sigma_mu^2 is (20.22, 536.32). */
/* Compare to book's result, at bottom of page 1046. */

/*********************************************************************************/