/*************************************************************************************** PRDS Reliability Analysis Tools Program: BlandAltmanPlot.sas Title: Bland & Altman Plot ((Bland & Altman, 1986 and 1999) Description: The Bland & Altman plot is made in which the differences between the two techniques are plotted against the averages of the two techniques. Input: inds -- input dataset; X, Y -- paired measurements in the input data Author: NANCY CHENG Date: 7/27/2004 Modified Date: 8/1/2006 /***************************************************************************************/ options nocenter linesize=80; %macro BlandAltmanPlot(inds=, x=, y=); * compute difference, and mean of two method measurements; data ba; set &inds; dif=x-y; avg=(x+y)/2; type=1; *for merging; label dif='Method Difference' avg='Method Mean'; run; * get mean and std of the difference, note that Bland Altman analysis use std dev not std err for the mean; proc means data=ba mean std noprint; var dif avg; output out=mout (drop=_type_ _freq_) mean(dif)=mdif std(dif)=stddif max(dif avg)=maxdif maxavg min(dif)=mindif n=nmax; run; data mout; set mout; type=1; ulmdif=mdif+1.96*stddif; * mean+1.96SD; llmdif=mdif-1.96*stddif; * mean-1.96SD; maxmdif=mdif+3*stddif; * for display purpose; minmdif=mdif-3*stddif; vaxisMx=max(maxmdif, maxdif); vaxisMn=min(minmdif, mindif); call symput('ulmdif', put(ulmdif, 8.3)); call symput('llmdif', put(llmdif, 8.3)); call symput('vaxisMx', put(vaxisMx, 8.3)); call symput('vaxisMn', put(vaxisMn, 8.3)); label mdif='mean difference' stddif='std of difference' maxdif='max difference' mindif='min difference' maxavg='max average' vaxisMx='Max value of VAxis' vaxisMn='Min value of VAxis'; run; proc sort data=ba;by type;run; proc sort data=mout; by type;run; * merge to get mdif, uldif, lldif...into the dataset ; data ba; merge ba mout; by type; run; * add a new variable, haxisv, for setting the range of method mean to get a better plot variable ref0 for drawing reference 0 line; data ba; set ba; if _n_=1 then haxisv=0; if _n_=nmax then haxisv=maxavg; ref0=0; run; data temp; ulmdif=&ulmdif; llmdif=&llmdif; vaxisMx=&vaxisMx; vaxisMn=&vaxisMn; * get the range for vertical axis to show the bias ; if ulmdif <= 0.0 then vaxisMx=-vaxisMn; *need to change max value in vaxis; if llmdif >= 0.0 then vaxisMn=-vaxisMx; *need to change min value in vaxis; inc0=(vaxisMx - &vaxisMn)/10; *get initial increment; if 0.001<= inc0 <0.01 then rdunit=0.001; else if 0.01<= inc0 <0.1 then rdunit=0.01; else if 0.1<= inc0 <1.0 then rdunit=0.1; else if 1.0<= inc0 <10.0 then rdunit=1; else if 10.0<= inc0 <100.0 then rdunit=10; else if 100.0<= inc0 <1000.0 then rdunit=100; else if 1000.0 <= inc0 <10000.0 then rdunit=1000; *adjust the max and min value in vaxis; newVaxisMx=ceil(vaxisMx*0.1/rdunit) /(0.1/rdunit); newVaxisMn=floor(vaxisMn*0.1/rdunit) /(0.1/rdunit); incpunit=inc0/rdunit; *get final increment; if 1<=incpunit<2 then inc=rdunit; else if 2<=incpunit<3 then inc=2*rdunit; else if 3<=incpunit<4 then inc=3*rdunit; else if 4<=incpunit<6 then inc=5*rdunit; else if 6<=incpunit<8 then inc=8*rdunit; else inc=10*rdunit; call symput('newVaxisMx', put(newVaxisMx, 8.3)); call symput('newVaxisMn', put(newVaxisMn, 8.3)); call symput('inc', put(inc, 8.3)); run; %global vtaxis; %let vtaxis=%str(label=(angle=90) order=(&newVaxisMn to &newVaxisMx by &inc)); * vertical axis definition; *set grahics options; goptions ftext=swissb ftitle=swissb htitle=2 htext=1 ; symbol1 cv=orange v=dot i=none height=0; symbol2 v=none i=join line=2 ci=blue height=3; *l=1, dash line; symbol3 v=none i=join line=1 ci=green height=3; *l=2, solid line; symbol4 v=none i=join line=1 ci=red height=3; *l=2, solid line; axis1 %str(&vtaxis); legend1 label=none value=(tick=1 'Method Difference vs Method Mean' tick=2 'Mean+1.96SD' tick=3 'Mean Difference' tick=4 'Mean-1.96SD' tick=5 'Method Difference =0') position=(bottom center outside); proc gplot data=ba; title "Bland & Altman Plot"; plot dif*avg=1 ulmdif*haxisv=2 mdif*haxisv=3 llmdif*haxisv=2 ref0*haxisv=4 /overlay vaxis=axis1 hminor = 0 vminor = 0 legend=legend1; footnote1 " "; footnote2 " "; run; quit; %mend BlandAltmanPlot; data test1; input x y; datalines; 0.49 0.39 0.83 0.69 0.71 0.55 0.38 0.05 0.57 0.31 0.68 0.51 0.69 0.62 0.5 0.6 0.7 0.55 0.62 0.65 0.54 0.46 0.58 0.54 0.32 0.41 0.33 0.3 0.44 0.42 0.47 0.45 0.45 0.47 0.51 0.48 0.56 0.37 0.43 0.39 0.53 0.49 0.75 0.8 0.65 0.7 0.73 0.7 0.42 0.41 0.49 0.42 0.59 0.55 0.58 0.6 0.49 0.43 ; run; %BlandAltmanPlot(inds=test1, x=x, y=y); data test2; input x y; datalines; 0.49 0.39 0.83 0.69 0.71 0.55 0.38 0.05 0.57 0.31 0.68 0.51 ; run; %BlandAltmanPlot(inds=test2, x=x, y=y); data test3; input x y; datalines; 0.020 0.025 0.030 0.028 0.040 0.042 0.050 0.055 0.060 0.065 0.070 0.068 0.080 0.086 ; run; %BlandAltmanPlot(inds=test3, x=x, y=y); data test4; input x y; datalines; 2.0 2.5 3.0 2.8 4.0 4.2 5.0 5.5 6.0 6.5 7.0 6.8 8.0 8.6 ; run; %BlandAltmanPlot(inds=test4, x=x, y=y); *runit=0.1; data test5; input x y; datalines; 20 25 30 28 40 42 50 55 60 65 70 68 80 86 ; run; %BlandAltmanPlot(inds=test5, x=x, y=y); *runit=1; data test6; input x y; datalines; 200 250 300 280 400 420 500 550 600 650 700 680 800 860 ; run; %BlandAltmanPlot(inds=test6, x=x, y=y); *runit=10; data test7; input x y; datalines; 2000 2500 3000 2800 4000 4200 5000 5500 6000 6500 7000 6800 8000 8600 ; run; %BlandAltmanPlot(inds=test7, x=x, y=y);*runit=100;