* Dataset and formulae from:
* 1) http://www-users.york.ac.uk/~mb55/meas/ba.htm
* 2) Bland JM, Altman DG. Measuring agreement in method comparison studies.
*    Statistical Methods in Medical Research 1999; 8: 135-160

* READ CAREFULLY BOTH PAPERS BEFORE USING THE MACRO *.

*******************************************************
*                  B A L O A R E P                    *
*******************************************************
* The MACRO needs that 'C:\Temp\' folder exists       *
* Runs with SPSS 14 or newer                          *
*******************************************************
* (C) Marta García-Granero 03/2009                    *
* send questions to: mgarciagranero@gmail.com         *
*******************************************************.

DEFINE BALOAREP(!POSITIONAL=!ENCLOSE('(',')') /!POSITIONAL=!ENCLOSE('(',')')/norm= !DEFAULT(H) !TOKENS(1)).
DATASET NAME OriginalData.
DATASET COPY WorkingData.
DATASET DECLARE Results WINDOW=HIDDEN.
DATASET ACTIVATE WorkingData.
COUNT nmiss = !1 !2 (SYSMIS) .
SELECT IF (nmiss=0). 
* Compute critical value for Student's t distribution *.
RANK VARIABLES=!1(A) /N INTO N1@ /PRINT=NO .
DO IF $casenum EQ 1.
- COMPUTE TValue = IDF.T(0.975,N1@-1) .
END IF.
MATRIX.
PRINT /TITLE='BLAND&ALTMAN LOA WITH REPEATED MEASURES'.
GET T95   /VAR=TValue /MISSING=OMIT.
GET data1 /VARIABLES= !1 /NAMES=vname1.
GET data2 /VARIABLES= !2 /NAMES=vname2.
COMPUTE n=NROW(data1).
* Statistics for 1st method *.
COMPUTE m1=NCOL(data1).
COMPUTE mean1=RSUM(data1)/m1.
COMPUTE Iroofs=MSSQ(data1).
COMPUTE Aroofs=MSUM((RSUM(data1))&**2)/m1. 
COMPUTE Sw2.1=(Iroofs-Aroofs)/(n*(m1-1)).
* Statistics for 2st method *.
COMPUTE m2=NCOL(data2).
COMPUTE mean2=RSUM(data2)/m2.
COMPUTE Iroofs=MSSQ(data2).
COMPUTE Aroofs=MSUM((RSUM(data2))&**2)/m2. 
COMPUTE Sw2.2=(Iroofs-Aroofs)/(n*(m2-1)).
* BALOA analysis with average values *.
COMPUTE data={mean1,mean2}.
COMPUTE diff=data(:,1)-data(:,2).
COMPUTE averdiff=MSUM(diff)/n.
COMPUTE SD2diff=(CSSQ(diff)-n*averdiff&**2)/(n-1).
COMPUTE CorrSD2=SD2Diff+((1-(1/m1))*Sw2.1)+((1-(1/m2))*Sw2.2).
COMPUTE sddiff=SQRT(CorrSD2).
COMPUTE mean=RSUM(data)/2.
COMPUTE gmean=MSUM(mean)/n.
COMPUTE tvalue=averdiff/(sddiff/SQRT(n)).
COMPUTE pvalue=2*(1-TCDF(ABS(tvalue),n-1)).
COMPUTE lower=averdiff-1.96*sddiff.
COMPUTE upper=averdiff+1.96*sddiff.
COMPUTE minrange=CMIN(diff)-0.1*(CMAX(diff)-CMIN(diff)).
DO IF lower LT minrange.
- COMPUTE minrange=lower-0.1*(CMAX(diff)-CMIN(diff)).
END IF.
COMPUTE maxrange=CMAX(diff)+0.1*(CMAX(diff)-CMIN(diff)).
DO IF upper GT maxrange.
- COMPUTE maxrange=upper+0.1*(CMAX(diff)-CMIN(diff)).
END IF.
COMPUTE cr=100*sddiff/gmean.
* Confidence interval for bias *.
COMPUTE low95=averdiff-T95*sddiff/SQRT(n).
COMPUTE upp95=averdiff+T95*sddiff/SQRT(n).
* Confidence interval for agreement limits *.
COMPUTE SELOA=SQRT(((1/n)+(1.96**2)/(2*(n-1)))*sddiff**2).
COMPUTE low95low=lower-T95*SELOA.
COMPUTE upp95low=lower+T95*SELOA.
COMPUTE low95upp=upper-T95*SELOA.
COMPUTE upp95upp=upper+T95*SELOA.
PRINT {averdiff,sddiff,gmean}
 /FORMAT='F8.3'
 /CLABELS='MeanDiff','SD(Diff)','G-Mean'
 /TITLE='Statistics (*)'.
PRINT {SQRT(Sw2.1);SQRT(Sw2.2);SQRT(SD2diff)}
 /FORMAT='F8.3'
 /RLABELS='Sw(1)','Sw(2)','SD(Dif)'
 /TITLE='Within subjects SD & observed SD for differences'.
PRINT {vname1;vname2}
 /FORMAT='A8'
 /RLABELS='Method1:','Method2:'
 /TITLE='(*) Methods compared:'.
PRINT {averdiff,SDDiff/SQRT(n),low95,upp95,tvalue,pvalue}
 /FORMAT='F8.3'
 /CLABEL='MeanDiff','Se(Diff)','Low95%CL','Upp5%CL','t-value','2-tail p'
 /TITLE='Paired samples t-test for the significance of bias'.
PRINT{lower,low95low,upp95low;upper,low95upp,upp95upp}
 /FORMAT='F8.3'
 /RLABELS='Lower','Upper'
 /CLABELS='LOA','Low95%CL','Upp95%CL'
 /TITLE='Limits of agreement (with 95%CI)'.
PRINT cr
 /FORMAT='F8.1'
 /RLABEL='RSDR(%)='
 /TITLE='Reproducibility relative standard deviation'.
SAVE {averdiff,lower,upper,minrange,maxrange} /OUTFILE=Results.
COMPUTE namevec={'Mean','Diff'}.
SAVE {mean,diff} /OUTFILE=WorkingData /NAMES=namevec.
END MATRIX.
PRESERVE.
SET LOCALE=ENGLISH.
DATASET ACTIVATE Results.
STRING #var1 TO #var5 (A12).
DO REPEAT A=#var1 TO #var5 /B=col1 TO col5.
- COMPUTE A = STRING(B,E11.3).
END REPEAT.
!LET !reflin1=!UNQUOTE(#var1).
!LET !reflin2=!UNQUOTE(#var2).
!LET !reflin3=!UNQUOTE(#var3).
!LET !ymin=   !UNQUOTE(#var4).
!LET !ymax=   !UNQUOTE(#var5).
* Write Graph to syntax file *.
WRITE OUTFILE 'c:\Temp\LOAGraph.sps'
 /"IGRAPH"
 /"  /VIEWNAME='Bland-Altman LOA plot'"
 /"  /X1 = VAR(Mean) TYPE = SCALE"
 /"  /Y = VAR(Diff)  TYPE = SCALE (MIN="!ymin" MAX="!ymax")"
 /"  /COORDINATE = VERTICAL"
 /"  /X1LENGTH=4.0"
 /"  /YLENGTH= 3.0"
 /"  /REFLINE Diff "!reflin1
 /"  /REFLINE Diff "!reflin2
 /"  /REFLINE Diff "!reflin3
 /"  /TITLE='                Bland-Altman plot (with LOA)'"
 /"  /CAPTION='Methods compared: " !QUOTE(!1)" & " !QUOTE(!2) "'"
 /"  /SCATTER COINCIDENT = NONE.".
DATASET ACTIVATE WorkingData.
DATASET CLOSE Results.
VARIABLE LEVEL Mean Diff (SCALE).
INCLUDE 'C:\Temp\LOAGraph.sps'. 
RESTORE.
* Normality of diferences *.
!IF (!UPCASE(!norm)='H') !THEN
FREQUENCIES VARIABLES=Diff
  /STATISTICS=SKEWNESS SESKEW KURTOSIS SEKURT
  /HISTOGRAM  NORMAL.
!ELSE
EXAMINE VARIABLES=Diff
  /PLOT BOXPLOT NPPLOT.
!IFEND.
* Correlation between ABS(Diff)&Mean (looking for error proportional to the mean) *.
COMPUTE Diff=ABS(Diff).
VAR LABEL Diff'Abs. differences'.
NONPAR CORR
  /VARIABLES=Mean Diff
  /PRINT=SPEARMAN TWOTAIL NOSIG.
* Get rid of copy of dataset *.
DATASET ACTIVATE OriginalData.
DATASET CLOSE WorkingData.
!ENDDEFINE.

* Sample dataset (from 1st paper) *.
DATA  LIST LIST/subject wright1 wright2 mini1 mini2 (5 F8).
BEGIN DATA
 1 494 490 512 525
 2 395 397 430 415
 3 516 512 520 508
 4 434 401 428 444
 5 476 470 500 500
 6 557 611 600 625
 7 413 415 364 460
 8 442 431 380 390
 9 650 638 658 642
10 433 429 445 432
11 417 420 432 420
12 656 633 626 605
13 267 275 260 227
14 478 492 477 467
15 178 165 259 268
16 423 372 350 370
17 427 421 451 443
END DATA.
VAR LEVEL ALL(SCALE).

* Macro call: 2 arguments are mandatory (the set of repetitions of
  2 methods to be compared, enclosed betwenn parenthesis)
  and one (norm, with 2 possible values, H and SW) is optional:
  if sample size is small, use SW to check for normality (Shapiro-Wilk&boxplot)
  otherwise use H (histogram with normal curve plotted) .

BALOAREP (wright1 wright2) (mini1 mini2) norm=SW.