%MACRO twoway;
****************************************************************************;
* Hays, R. D. Wang, E., & Sonksen, M.  (1995, September).
*  General Reliability and Intraclass Correlation Program (GRIP).
*  Proceedings of the 3rd Annual Western Users of SAS Conference,
*  220-223.
****************************************************************************;
* Modified by Sally Carson to add REL70;
* Modified by Karen Spritzer to add REL80, REL90, and REL95 and to run 
* under SAS 8.0;
*
* Modified by Karen Spritzer on 2/18/04 to put "upcase" on both sides of the 
* comparison between "source" and "target" in the 2-way portion of code;  
* SAS 8 is case-sensitive with respect to variable names and that was mixing 
* up the procedure;
****************************************************************************;
proc glm data=&indata outstat=stats noprint;
  class &targetv &repeatv;
  model &dv = &targetv &repeatv ;
  run;


proc sort data=stats;
  by _name_ _SOURCE_;
  run;

data allrel;
  retain bdf bms jdf edf wms jms ems k;
  set stats;
  by _name_;
  if _type_='SS1' then delete;
  if _source_='ERROR' then do;
     ems=ss/df;
     edf=df;
  end;
  if upcase(_source_)="%upcase(&targetv)" then do; * modified by KS on 2-18-04;  
     bms=ss/df;
     bdf=df;
  end;
  if upcase(_source_)="%upcase(&repeatv)" then do;  * modified by KS on 2-18-04 ;
     jms=ss/df;
     jdf=df;
     k=df+1;
  end;
  if last._name_ then do;
    wms=((ems*edf)+(jms*jdf))/(edf+jdf);
    n=bdf+1;
m=(n*(k-1))/(n*(k-1)-2);
theta=(bms-(m*wms))/(k*m*wms);
rii=theta/(1+theta);
rtt=(k*theta)/(1+(k*theta));


fixed=(bms-ems)/(bms+((k-1)*ems));
fixedk=(bms-ems)/bms;
biased=(bms-wms)/(bms+(k-1)*wms);
random=(bms-ems)/((bms)+((k-1)*ems)+((k*(jms-ems))/n));
rk=(bms-ems);
randk=rk/(bms+((jms-ems)/n));
k=(bms-wms)/bms;
    output;
  end;
run;



data _null_;
file print header=hea1 ps=64 notitles;
set  allrel;
kkm=jdf+edf;
kk3=edf+bdf+jdf;
*+++++++++1+++++++++2+++++++++3+++++++++4+++++++++5+++++++++6+++++++++7+++++++++8+++++++++9++++++++10;
put
@20                'Degrees of     mean      Label  for'/
@5  'Source          freedom' @35 'square' @45 'mean Square'/
@5 '----------------------------------------------------------';


put @5 'Ratees  (N-1)'  @24 bdf  3.   @34 bms 7.2  @49 'BMS'/
    @5 'Within'         @24 kkm 3.   @34 wms 7.2  @49 'WMS'/
    @7 'Raters  (K-1)'  @25 jdf  3.   @34 jms 7.2  @49 'JMS'/
    @7 'Raters x Ratees ' @25 edf 3. @34 ems 7.2  @49 'EMS'/
@5 '----------------------------------------------------------'/
@5 'Total'            @24 kk3    3. /;
return;


hea1:
   do;
   put @5 "&t1 "/;
   end;
       return;
run;



data _null_;
file print header=hea1 ps=64 notitles;
set  allrel;
*+++++++++1+++++++++2+++++++++3+++++++++4+++++++++5+++++++++6+++++++++7+++++++++8+++++++++9++++++++10;
put
@5 'Model            Reliability            Intraclass Correlation'/
@5 '---------------------------------------------------------------';
put @5 'One way'/
    @7 'Biased'       @24 k 5.3   @53 biased 5.3 /
    @7 'Unbiased'     @24  rtt 5.3   @53 rii 5.3 //
    @5 'Two-way'/
    @7 'Fixed effects'        @24  fixedk 5.3   @53 fixed 5.3 /
    @7 'Random effects'       @24 randk 5.3   @53 random 5.3 /
@5 '---------------------------------------------------------------'/;
return;


hea1:
   do;
   put @5 " &t2 "//;
   end;
       return;
run;
%MEND twoway;
**********************************************;


**********************************************;
%MACRO oneway;
**********************************************;
title1 "&t1";
title2 "&t2";

proc anova data=&indata outstat=est1 noprint;
 class &targetv;
 model &dv= &targetv;
run;


data est;
 set est1;
retain;
if _type_='ERROR' then wms=ss/df;
if _type_='ERROR' then n=df;
if _type_='ERROR' then errdf=df;
if _type_='ANOVA' then bms=ss/df;
if _type_='ANOVA' then betdf=df;
if _type_='ANOVA' then nrated=df+1;
if _type_='ANOVA' then nn=n+df+1;
if _type_='ANOVA' then k=nn/nrated;
OUTPUT;
data est;
 set est;
m=(n*(k-1))/(n*(k-1)-2);
theta=(bms-(m*wms))/(k*m*wms);
rii=theta/(1+theta);
rtt=(k*theta)/(1+(k*theta));
frtt=(bms-wms)/bms;
frii=frtt/(k * (1 + (frtt * (1/k - 1))));


/* added by sally carson ----> */
*bn4rel70=(frtt-(frii*frtt))/(frii*(1-frtt));
 bn4rel70=(.70 -(frii*.70 ))/(frii*(1-.70 ));
*un4rel70=(rtt-(rii*rtt))/(rii*(1-rtt));
 un4rel70=(.70-(rii*.70))/(rii*(1-.70));

/* added by karen spritzer ----> */
 bn4rel80=(.80 -(frii*.80 ))/(frii*(1-.80 ));
 un4rel80=(.80-(rii*.80))/(rii*(1-.80));
 
 bn4rel90=(.90 -(frii*.90 ))/(frii*(1-.90 ));
 un4rel90=(.90-(rii*.90))/(rii*(1-.90));

 bn4rel95=(.95 -(frii*.95 ))/(frii*(1-.95 ));
 un4rel95=(.95-(rii*.95))/(rii*(1-.95));

OUTPUT;
run;
DATA EST;
 SET EST;
if k ne .;
RUN;


data _null_;
file print ps=64 /*notitles*/;
set  est;
k2=k-1;
n2=n-1;
k3=n2+n;
kk=betdf+errdf;
*+++++++++1+++++++++2+++++++++3+++++++++4+++++++++5+++++++++6+++++++++7+++++++++8+++++++++9++++++++10;
put ///// @5 'GRIP' //
@20                'Degrees of     mean      Label  for'/
@5  'Source          freedom' @35 'square' @45 'mean Square'/
@5 '----------------------------------------------------------';


put @5 'Between    '  @24 betdf  5.   @34 bms 7.2  @49 'BMS'/
    @5 'Within'       @24 errdf  5.   @34 wms 7.2  @49 'WMS'/
@5 '----------------------------------------------------------'/
@5 'Total'            @24 kk     5. /;
PUT /////;
*return;
/*
data _null_;
file print ps=64 /*notitles**;
set  est;                    */
*+++++++++1+++++++++2+++++++++3+++++++++4+++++++++5+++++++++6+++++++++7+++++++++8+++++++++9++++++++10;
put
@5 'Model            Reliability   ICC    N4rel.70  N4rel.80  N4rel.90    N4rel.95'/
@5 '-------------------------------------------------------------------------------------------------';
put @5 'One way'/
    @7 'Biased'       @24 frtt 5.3   @33 frii 5.3    @43 bN4rel70 7.2   @53 bN4rel80 7.2 @63 bN4rel90 7.2 @76 bN4rel95 7.2/
    @7 'Unbiased'     @24 rtt  5.3   @33  rii 5.3    @43 uN4rel70 7.2   @53 uN4rel80 7.2 @63 bN4rel90 7.2 @76 bN4rel95 7.2//
@5 '-------------------------------------------------------------------------------------------------'/;
return;
RUN;
%MEND oneway;
**********************************************;


%MACRO grip(indata=,targetv=,repeatv=,dv=,nrepeatv=,type=,t1=,t2=);
 %IF %EVAL(&type) %then %DO;
  %twoway;
 %end;
%ELSE
 %DO;
  %oneway;
 %end;
%MEND grip;



******************************************************;
** Program and macro call starts here:
******************************************************;

data all;
input plant $ 1-2 expcon rater height0 height1;
cards;
A1 1 1 120 121
A1 1 2 118 120
A2 2 1  84  85
A2 2 2  96  88
B1 2 1 107 108 
B1 2 2 105 104
B2 1 1  94 100
B2 1 2  97 104
C1 2 1  85  88
C1 2 2  91  96
C2 1 1  79  86
C2 1 2  78  92
D1 1 1  70  76
D1 1 2  72  80
D2 2 1  54  56
D2 2 2  56  60
E1 1 1  85 101
E1 1 2  97 108
E2 2 1  90  84
E2 2 2  92  96
;
run;
******************************************************;
%GRIP(indata=all,targetv=plant,repeatv=rater,dv=height0,
      type=1,t1=test of GRIP macro,t2=);
******************************************************;