/*  Tukey test for additivity */

/******************* MACRO BEGINS HERE ****************************/
%MACRO tukeyaddit(dataset=, factor1=, factor2=, response=);
ods listing close;
proc glm data=&dataset;
  class &factor1 &factor2;
  model &response = &factor1 &factor2;
  ods output overallanova=overall modelanova=model;
run;
quit;
ods listing;
ods output close;
data _null_;
  set overall;
  if source='Corrected Total' then call symput('overall', ss);
run;
data _null_;
  set model ;
  if hypothesistype=1 and source="&factor2" then call symput('ssa', ss);
  if hypothesistype=1 and source="&factor1" then call symput('ssb', ss);
  if hypothesistype=1 and source="&factor2" then call symput('dfa', df);
  if hypothesistype=1 and source="&factor1" then call symput('dfb', df);
run;
%put here is &overall &ssa &ssb &dfa &dfb; /* the statement will appear in the log file so you can check the calculations */
proc sql;
  create table temp1 as
  select &response, &factor1, &factor2 , mean(&response) as yj
  from &dataset
  group by &factor1;
quit;
proc sql;
  create table temp2 as
  select *, mean(&response) as yi
  from temp1
  group by &factor2;
quit;
proc sql noprint; 
   select mean(&response) into :meanp from temp1;
quit;
%put here is  &meanp; /* check value in log file */
proc sql noprint;
  select sum( (yi - &meanp)*(yj - &meanp)*&response ) into :total from temp2;
quit;
%put here is &total; /*check value in log file*/
data final;
  msa = &ssa/(&dfb+1);
  msb = &ssb/(&dfa+1);
  ssab = (&total*&total) / ( msa*msb );
  ssrem = &overall - &ssa - &ssb - ssab;
  f = ssab/( ssrem/((&dfa+1)*(&dfb+1) - (&dfa+1) - (&dfb+1)) );
  p_value = 1- cdf('F',f,  1, (&dfa+1)*(&dfb+1) - (&dfa+1) - (&dfb+1) );
run;
proc print data=final;
run;
%MEND tukeyaddit;
/******************** MACRO ENDS HERE *****************************/

/* 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 =  ; /* input number of observations per cell */
r =  ; /* input number of levels of the random factor */
MSTR =  ; /* input MSTR from PROC GLM */
MSE =   ; /* input MSE from PROC GLM  */

/* 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;

************************************************************************