***********************************************************************
* This macro computes agreement statistics on CCC(rc), precision(r),
* accuracy(ca), TDI(k), CP, and their one sided confidence limits.
* They are computed when target values are random or fixed, and when the
* error is constant or proportional.  It also produce the agreement plot
* around the identity line.  

* Written by Lawrence Lin and verified by Min Yang.
* Updated by Lawrence Lin on 6/14/2007
***********************************************************************;

* Parameters are:;
*********************************************************************
* Dataset: name of your data;
* Y: reading 1, verticle axis 
* V_label: label for reading 1, verticle axis 
* X: reading 2, horizontal axis
* H_label: label for reading 2, horizontal axis
* Min: minimum of the plotting range.  
* Max: maximum of the plotting range.  
* By: increment of the plotting range for constant error,
  or the log scale increment in between min and max for proportional
  error case
* CCC_a: CCC allowance
* CP_a: coverage probability (CP) allowance
* TDI_a: TDI allowance, must be in % when Error="prop" is specified.
* Error = "const" for constant error structure.  Here, TDI is the absolute
          difference.
        = "prop" for proportional error structure.  Here, TDI is the
          % change.  Log transformation to data will be applied.
* Target: "random" or "fixed".  
* Alpha: 100(1-alpha)% one-tailed confidence limit  
* Dec: significant digits after the decimal point printed for TDI when the
  "const" error is specified (eg. 2).
***********************************************************************;

* This macro creats a .rtf file containing a graph and a table.
* Need to specify the location of this rtf file.
* Need 'ODS rtf close' statement after each macro run.
* You can use "insert picture" in MS-WORD to import the created .rtf file.
***********************************************************************;

%macro agreement(dataset=, y=, x=, V_label=, H_label=, min=, max=, by=, 
       CCC_a=, CP_a=, TDI_a=, error=, target=, dec=, alpha=, );

**** Set ods escape character to insert simple formatting text into ODS output         ****; 
ods escapechar='^';
%***** Create symbol for turning on subscript. ******;
%LET SUBS=%str(^\sub);

data c; 
 set &dataset; 
 dummy=1;
 y=&y; x=&x; 
 if y=. or x=. then delete;
 if &error="prop" then do;
   y=log(y); x=log(x); 
 end; 
 d=y-x; 
proc means mean var uss n noprint; var y x d; 
output out=t mean=m1 m2 yyy var=v1 v2 yy uss=xx xxx e2 n=xxxx xxxxx n;
run;

data t; 
 set t; 
 call symputx('obs', n);
 dummy=1;
 mu_d=m1-m2; d2=mu_d**2;
 e2=e2/n; sd2=(e2-d2)*n/(n-3); s12=v1*(n-1)/n; s22=v2*(n-1)/n; 
 U=mu_d/sqrt(SQRT(s12*S22)); V=sqrt(S12/S22); Ca=2/(V+1/V+U**2); 
 rc=1-e2/(d2+s12+s22); 
 r=(rc/ca); 
 d2_s2=(d2/sd2); b1=r*v; b0=m1-b1*m2; 
 se2=s12*(1-r**2)*n/(n-3); e2=e2*n/(n-1);
 kp=&TDI_a;  if &error="prop" then kp=log(&TDI_a/100+1);

* For CP when the target value is fixed;
%if &target="fixed" %then %do;
 data cc; merge c t; by dummy;
 di=b0+(b1-1)*x; d2i=di**2;
 cpi=probchi(kp**2/se2,1,d2i/se2); 
 se=sqrt(se2); kpm=(kp+di)/se; kmm=(kp-di)/se;
 C0=pdf('normal',-kpm)-pdf('normal',kmm); C1=C0*x;
 C2=-kpm*pdf('normal',-kpm)-kmm*pdf('normal',kmm);
 proc means noprint mean; 
  var cpi C0 C1 C2; 
  output out=cpf mean=cp C0 C1 C2;
  run;

 data t; merge t cpf;
  s_p=C0**2+(C0*m2-C1)**2/s22+C2**2/2;
  s_T=sqrt(s_p/((n-3)*cp**2*(1-cp)**2));
  T=log(cp/(1-cp));
  T_l=T-probit(1-&alpha)*s_T;
  cp_l=exp(T_l)/(1+exp(T_l));
 run;
%end;

data tt; 
 set t; 
 label rc='CCC' rc_l='Lower CL of CCC'
       r='Precision' r_l='Lower CL of Precision'
       ca='Accuracy' ca_l='Lower CL of Accuracy'
       e2='MSD' e2_u='Upper CL of MSD'
       k='TDI' k_u='Upper CL of TDI'
       cp='CP' cp_l='Lower CL of CP' d2_s2='Relative Bias Square';

 * For CCC, accuracy and precision;
 L_ca=log(ca/(1-ca));
 z_rc=LOG((1+RC)/(1-RC))/2; 
 z_r=LOG((1+R)/(1-R))/2; 

 * For MSD & approx TDI & CP;
 W=LOG(e2);
 k=probit(1-(1-&CP_a)/2)*sqrt(e2);
 cp_a=probchi(kp**2/e2,1); 

 * Standard errors of CCC, accuracy, precision, and MSD;
 if &target="random" then do;
  SZ_rc=((1-R**2)*RC**2/(1-RC**2)/R**2 + 2*RC**3*(1-RC)*U**2/R/(1-RC**2)**2
     - RC**4*U**4/R**2/(1-RC**2)**2/2)/(n-2);
  SZ_r=1/(n-3);
  SL_ca=((ca**2*u**2*(v+1/v-2*r) + ca**2*(v**2+1/v**2+2*r**2)/2
     + (1+r**2)*(ca*u**2-1))/((n-2)*(1-ca)**2));
  s_W=2*(1-d2**2/e2**2)/(n-2);
 end;
 if &target="fixed" then do;
  SZ_rc=(1-r**2)*rc**2/((n-2)*(1-rc**2)**2*r**2)*(v*u**2*rc**2+(1-rc*r*v)**2
     + v**2*rc**2*(1-r**2)/2); 
  SZ_r=(1-r**2/2)/(n-3);
  SL_ca=(u**2*v*ca**2*(1-r**2) + (1-v*ca)**2*(1-r**4)/2)/((n-2)*(1-ca)**2);
  s_W=2*(1-(d2+s22*(1-b1)**2)**2/e2**2)/(n-2);
* Add case when there is only one X level. 7/18/03;
  if s22=0 then s_W=2*(1-d2**2/e2**2)/(n-2);
 end;

 SZ_rc=SQRT(SZ_rc); SZ_r=SQRT(SZ_r); SL_ca=SQRT(SL_ca); S_W=sqrt(S_W);
 z_rc_l=z_rc-probit(1-&alpha)*sz_rc; 
 z_r_l=z_r-probit(1-&alpha)*sz_r;
 l_ca_l=L_ca-probit(1-&alpha)*sl_ca;
 rc_l=tanh(z_rc_l); r_l=tanh(z_r_l); ca_l=exp(l_ca_l)/(1+exp(l_ca_l));

 W_u=W+probit(1-&alpha)*s_W;
 e2_u=exp(W_u);
 k_u=probit(1-(1-&CP_a)/2)*sqrt(e2_u);
 cp_a_l=probchi(kp**2/e2_u,1); 
 if &error="prop" then do; k=100*(exp(k)-1); k_u=100*(exp(k_u)-1); end;

* For CP when the target value is random;;
 if &target="random" then do;
  cp=probchi(kp**2/sd2,1,d2_s2); 
  sd=sqrt(sd2); kpm=(kp+mu_d)/sd; kmm=(kp-mu_d)/sd;
  s_T=((pdf('normal',-kpm)-pdf('normal',kmm))**2
     + (kmm*pdf('normal',kmm)+kpm*pdf('normal',-kpm))**2/2)
     /((n-3)*cp**2*(1-cp)**2); s_T=sqrt(s_T);
  T=log(cp/(1-cp)); T_l=T-probit(1-&alpha)*s_T;
  cp_l=exp(T_l)/(1+exp(T_l));
 end;
 conf=(1-&alpha)*100;
 call symputx('conf', conf);
 keep rc rc_l r r_l ca ca_l e2 e2_u k k_u cp cp_l d2_s2;
*proc print data=tt;
run;

data tb;
 set tt;
 S='Estimate            ';
 CCC=rc; Precision=r; Accuracy=ca; TDI=k; CP=cp; RBS=d2_s2; output;
 S="&conf.% Conf. Limit";
 CCC=rc_l; Precision=r_l; Accuracy=ca_l; TDI=k_u; CP=cp_l; RBS=.; output;
 S='Allowance';
 CCC=&CCC_a; Precision=.; Accuracy=.; TDI=&TDI_a; CP=&CP_a; RBS=.; output;
keep s ccc Precision Accuracy TDI CP RBS;
run;

proc report data=tb headline headskip split='|' nowindows;
  column s CCC Precision Accuracy TDI CP RBS;
 define s/ 'Statistics' format=$20. left;
 define CCC/ 'CCC' format=10.4 right;
 define Precision / 'Precision' format=10.4 right;
 define Accuracy / 'Accuracy' format=10.4 right;
  %if &error="const" %then %do;
   define TDI / "TDI &subs &CP_a" format=10.&dec right;
   define CP / "CP &subs &TDI_a" format=10.4 right;
  %end;
  %else %do;
   define TDI / "TDI% &subs &CP_a" format=10.2 right;
   define CP / "CP &subs &TDI_a.%" format=10.4 right;
  %end;
  define RBS / "RBS*" format=10.2 right;
  compute after _page_;
  line "n=&obs";
  line "*: The relative bias squared (RBS) must be less than 1 or 8 for CP_a of 0.9 or 0.8,";
  line "respectively, in order for the approximated TDI values to be valid.";
  endcomp;
 run;
 quit;

* For graph;
goptions reset=all ftext=swiss htext=2 device=cgmof97L; 
AXIS1 label=(a=90 "&V_label") offset=(1 cm) order=&min to &max by &by ;
AXIS2 label=("&H_label") offset=(1 cm) order=&min to &max by &by ;
AXIS3 label=(a=90 "&V_label") offset=(1 cm) logbase=2 order=&min &by &max ;
AXIS4 label=("&H_label") offset=(1 cm) logbase=2 order=&min &by &max ;

data p; 
 set c;
%if &error="const" %then %do; %let na=1; %let nb=2; %end;
%if &error="prop" %then %do; %let na=3; %let nb=4; %end;
 y0=.; if _n_=1 then y0=&min; if _n_=2 then y0=&max; 
 if &error="prop" then do; y=exp(y); x=exp(x); end; 
PROC GPLOT;
PLOT y*x y0*y0/
overlay FR NOLEGEND VAXIS=AXIS&na HAXIS=AXIS&nb; 
SYMBOL1 h=2 c=black V=circle I=NONE;
SYMBOL2 W=1 c=black V=' ' h= I=join l=1;
*title ' ';
*footnote ' ';
run;
quit;

%mend agreement;

* The macro calling examples for the JASA paper are provided below:.;
/*
* Data must be arranged as one row per sample with paired measurements
* per row;
 
* Calling macro for Example 1;
ods rtf file='x:\xx\xxx\AgreementEX1.rtf' 
    style=styles.TablesRTF NOTOC_DATA;
agreement(dataset=a, y=HemoCue, x=Sigma, V_label=HemoCue, H_label=Sigma,   
          min=0, max=2000, by=250,
          CCC_a=0.9775, CP_a=0.9, TDI_a=150, 
          error="const", target="random", dec=1, alpha=0.05); 
ods rtf close;

* Calling macro for Example 2;
ods rtf file='x:\xx\xxx\AgreementEX2.rtf' 
    style=styles.TablesRTF NOTOC_DATA;
%agreement(dataset=b,y=Y1, x=x, V_label=1:5, H_label=Actual FVIII %, min=1, max=256,
            by=2 4 8 16 32 64 128, 
            CCC_a=0.9775, CP_a=0.8, TDI_a=50,
            error="prop", target="fixed", alpha=0.05); 
%agreement(dataset=b,y=Y2, x=x, V_label=1:10, H_label=Actual FVIII %, min=1, max=256,
            by=2 4 8 16 32 64 128, 
            CCC_a=0.9775, CP_a=0.8, TDI_a=50,
            error="prop", target="fixed", alpha=0.05); 
ods rtf close;
*/