* Reference:
* Passing H, Bablok W (1983) A new biometrical procedure for testing the equality
* of measurements from two different analytical methods. Application of linear
* regression procedures for method comparison studies in Clinical Chemistry,
* Part I. J. Clin. Chem. Clin. Biochem., 21:709-720.

*************************************************
*        PASSING-BABLOK REGRESSION              *
*************************************************
* Runs with SPSS 14/15 (both graphs might give  *
* problems with SPSS 16/17)                     *
*************************************************
* (C) Marta García-Granero 04/2009              *
* send questions to: mgarciagranero@gmail.com   *
* Downloaded from: http://gjyp.nl/marta/        *
* Feel free to use or modify this code, but     *
* acknowledge the author and the web site       *
*************************************************
* The asymptotic 2-tailed p-value for linearity *
* test is estimated using the first three terms *
* of the Smirnov (1948) formula (see SPSS help  *
* for details)                                  *
*************************************************
* The accuracy has been checked with MedCalc    *
* and Method Validator                          *
*************************************************
* WARNING: variable names lenght should be only *
* 8 characters at most (MATRIX limitation)      *
*************************************************

PRESERVE.
SET RESULTS=NONE MESSAGES=NONE ERRORS=NONE.
HOST COMMAND=['IF not exist C:\Temp MD C:\Temp'].
RESTORE.

DEFINE PASSINGBABLOK(x=!TOKENS(1) / y=!TOKENS(1)).
DATASET NAME OriginalData.
DATASET COPY WorkingData   WINDOW=HIDDEN.
DATASET DECLARE Parameters WINDOW=HIDDEN.
DATASET ACTIVATE WorkingData.
FREQUENCIES
  VARIABLES=!x !y
  /FORMAT=NOTABLE
  /STATISTICS=STDDEV SEMEAN MINIMUM MAXIMUM MEAN MEDIAN.
TEMPORARY.
SET MXLOOPS=100000.
MATRIX.
PRINT /TITLE="**** PASSING-BABLOK REGRESSION ****".
GET pair
 /VARIABLES= !x !y
 /NAMES=vname
 /MISSING OMIT.
PRINT {T(vname)}
 /FORMAT='A8'
 /RLABELS='X=','Y='
 /TITLE='Variables compared'.
COMPUTE NData=NROW(pair).
COMPUTE X=pair(:,1).
COMPUTE Y=pair(:,2).
* Max number of slopes *.
COMPUTE N=NData*(NData-1)/2.
* Compute all valid Sij *.
COMPUTE SData=MAKE(N,1,0).
COMPUTE count=0.
LOOP i=1 TO NData-1.
. LOOP j=i+1 TO NData.
* Skip 0/0 *.
.  DO IF X(i) NE X(j).
.   COMPUTE slope=(Y(i)-Y(j))/(X(i)-X(j)).
* Record it only if Sij=-1 *.
.   DO IF slope NE -1.
.    COMPUTE count=count+1.
.    COMPUTE SData(count)=slope.
.   END IF.
* Add values if X(i)-X(j)EQ 0 but Y(i)-Y(j) NE 0 (+/- infinite, or almost) *.
.  ELSE IF (X(i) EQ X(j)) AND(Y(i) NE Y(j)).
.   COMPUTE count=count+1.
.   COMPUTE SData(count)=(Y(i)-Y(j))*9.9E50.
.  END IF.
. END LOOP.
END LOOP.
* Trim to real number of slopes *.
COMPUTE N=count.
COMPUTE Sij=SData(1:N).
* Finding offset (K=Nr of Sij<-1) *.
COMPUTE K=0.
LOOP i=1 TO N.
. DO IF Sij(i) LT -1.
.  COMPUTE K=K+1.
. END IF.
END LOOP.
* Sort Sij to compute shifted median (from offset) *.
COMPUTE SortedS=Sij.
COMPUTE SortedS(GRADE(Sij))=Sij.
* Compute b *.
COMPUTE even=(TRUNC(N/2) EQ N/2).
* Formula for N odd *.
DO IF even EQ 0. 
. COMPUTE b=SortedS(K+(N+1)/2).
ELSE.
* Formula for N even *.
COMPUTE b1=SortedS(K+N/2).
COMPUTE b2=SortedS(1+K+N/2).
COMPUTE b=(b1+b2)/2.
END IF.
* Confidence interval for b *.
COMPUTE Zalpha2=1.96.
COMPUTE Calpha=Zalpha2*SQRT((NData*(NData-1)*(2*NData+5))/18).
COMPUTE M1=RND((N-Calpha)/2).
COMPUTE M2=N-M1+1.
COMPUTE LowB=SortedS(M1+K).
COMPUTE UppB=SortedS(M2+K).
* Compute Aij *.
COMPUTE Aij=Y-b&*X.
COMPUTE LowAij=Y-UppB&*X.
COMPUTE UppAij=Y-LowB&*X.
* Sort Aij to compute median *.
COMPUTE SortedA=Aij.
COMPUTE SortedA(GRADE(Aij))=Aij.
* Sort LowAij to compute Lower CL for a *.
COMPUTE SortedAL=LowAij.
COMPUTE SortedAL(GRADE(LowAij))=LowAij.
* Sort UppAij to compute Upper CL for a *.
COMPUTE SortedAU=UppAij.
COMPUTE SortedAU(GRADE(UppAij))=UppAij.
* Compute a & CI *.
COMPUTE even=(TRUNC(NData/2) EQ NData/2).
* Formula for NData odd *.
DO IF even EQ 0. 
. COMPUTE a=SortedA((NData+1)/2).
. COMPUTE LowA=SortedAL((NData+1)/2).
. COMPUTE UppA=SortedAU((NData+1)/2).
ELSE.
* Formula for N even *.
COMPUTE a1=SortedA(NData/2).
COMPUTE a2=SortedA(1+NData/2).
COMPUTE a=(a1+a2)/2.
COMPUTE Lowa1=SortedAL(NData/2).
COMPUTE Lowa2=SortedAL(1+NData/2).
COMPUTE LowA=(Lowa1+Lowa2)/2.
COMPUTE Uppa1=SortedAU(NData/2).
COMPUTE Uppa2=SortedAU(1+NData/2).
COMPUTE UppA=(Uppa1+Uppa2)/2.
END IF.
PRINT {a,LowA,UppA;b,LowB,UppB}
 /FORMAT='F8.4'
 /RLABELS='A=','B='
 /CLABEL='Value','Low95%CL','Upp95%CL'
 /TITLE='Intercept and slope'.
* Testing for linearity *.
COMPUTE Lsmall=0.
COMPUTE Lbig=0.
LOOP i=1 TO NData.
. DO IF Y(i) GT a+b*X(i).
.  COMPUTE Lsmall=Lsmall+1.
. ELSE IF Y(i) LT a+b*X(i).
.  COMPUTE Lbig=Lbig+1.
. END IF.
END LOOP.
COMPUTE Ri=MAKE(NData,1,0).
LOOP i=1 TO NData.
. DO IF Y(i) GT a+b*X(i).
.  COMPUTE Ri(i)=SQRT(Lbig/Lsmall).
. ELSE IF Y(i) LT a+b*X(i).
.  COMPUTE Ri(i)=-SQRT(Lsmall/Lbig).
. END IF.
END LOOP.
COMPUTE Di=(Y-a+X/b)/SQRT(1+(1/b**2)).
* Sort Ri by ascending Di *.
COMPUTE SortedRi=Ri.
COMPUTE SortedRi(grade(Di))=Ri.
COMPUTE cusum=MAKE(NData,1,0).
LOOP i=1 TO NData.
. COMPUTE cusum(i)=CSUM(SortedRi(1:i)).
END LOOP.
* KS critical values for alpha=0.10, 0.05, 0.025, 0.01 & 0.005 (one tailed) *.
DO IF NData LE 40. /* KS critical values for n LE 40 *.
. COMPUTE KSValues={0.9000,0.9500,0.9750,0.9900,0.9950;
                    0.6838,0.7764,0.8419,0.9000,0.9293;
                    0.5648,0.6360,0.7076,0.7846,0.8290;
                    0.4927,0.5652,0.6239,0.6889,0.7342;
                    0.4470,0.5094,0.5633,0.6272,0.6685;
                    0.4104,0.4680,0.5193,0.5774,0.6166;
                    0.3815,0.4361,0.4834,0.5384,0.5758;
                    0.3583,0.4096,0.4543,0.5065,0.5418;
                    0.3391,0.3875,0.4300,0.4796,0.5133;
                    0.3226,0.3687,0.4092,0.4566,0.4889;
                    0.3083,0.3524,0.3912,0.4367,0.4677;
                    0.2958,0.3382,0.3754,0.4192,0.4490;
                    0.2847,0.3255,0.3614,0.4036,0.4325;
                    0.2748,0.3142,0.3489,0.3897,0.4176;
                    0.2659,0.3040,0.3376,0.3771,0.4042;
                    0.2578,0.2947,0.3273,0.3657,0.3920;
                    0.2504,0.2863,0.3180,0.3553,0.3809;
                    0.2436,0.2785,0.3094,0.3457,0.3706;
                    0.2373,0.2714,0.3014,0.3369,0.3612;
                    0.2316,0.2647,0.2941,0.3287,0.3524;
                    0.2262,0.2586,0.2872,0.3210,0.3443;
                    0.2212,0.2528,0.2809,0.3139,0.3367;
                    0.2165,0.2475,0.2749,0.3073,0.3295;
                    0.2120,0.2424,0.2693,0.3010,0.3229;
                    0.2079,0.2377,0.2640,0.2952,0.3166;
                    0.2040,0.2332,0.2591,0.2896,0.3106;
                    0.2003,0.2290,0.2544,0.2844,0.3050;
                    0.1968,0.2250,0.2499,0.2794,0.2997;
                    0.1935,0.2212,0.2457,0.2747,0.2947;
                    0.1903,0.2176,0.2417,0.2702,0.2899;
                    0.1873,0.2141,0.2379,0.2660,0.2853;
                    0.1844,0.2108,0.2342,0.2619,0.2809;
                    0.1817,0.2077,0.2308,0.2580,0.2768;
                    0.1791,0.2047,0.2274,0.2543,0.2728;
                    0.1766,0.2018,0.2242,0.2507,0.2690;
                    0.1742,0.1991,0.2212,0.2473,0.2653;
                    0.1719,0.1965,0.2183,0.2440,0.2618;
                    0.1697,0.1939,0.2154,0.2409,0.2584;
                    0.1675,0.1915,0.2127,0.2379,0.2552;
                    0.1655,0.1891,0.2101,0.2349,0.2521}.
. COMPUTE KSVal=SQRT(Ndata)*KSValues(Ndata,:).
ELSE. /* KS critical values for n GT 40 *.
. COMPUTE KSVal={1.07,1.22,1.36,1.52,1.63}.
END IF.
COMPUTE CritVal=KSVal*SQRT(Lbig+1).
COMPUTE LinStat=CMAX(ABS(cusum)).
PRINT {LinStat,CritVal}
 /FORMAT='F8.3'
 /CLABEL='Stat(*)','(0.10)','(0.05)','(0.025)','(0.01)','(0.005)'
 /TITLE='Linearity: Max(Abs(cusum)) & Critical values for different 1-tail alpha levels'.
DO IF Linstat LE CritVal(1).
. PRINT /TITLE='(*) P>0.10'.
ELSE IF Linstat LE CritVal(2).
. PRINT /TITLE='(*) 0.05<P<0.10'.
ELSE IF Linstat LE CritVal(3).
. PRINT /TITLE='(*) 0.025<P<0.05'.
ELSE IF Linstat LE CritVal(4).
. PRINT /TITLE='(*) 0.01<P<0.025'.
ELSE IF Linstat LE CritVal(5).
. PRINT /TITLE='(*) 0.005<P<0.01'.
ELSE IF Linstat GT CritVal(5).
. PRINT /TITLE='(*) P<0.005'.
END IF.
* KS exact p value using 3 terms of Smirnov formula (1948) *.
COMPUTE ZValue=LinStat/SQRT(Lbig+1).
DO IF ZValue LT 0.27.
. COMPUTE PValue2=1.
ELSE IF (ZValue GE 0.27) AND (ZValue LT 1).
. COMPUTE Q=EXP(-1.233701*(ZValue**(-2))).
. COMPUTE PValue2=1-(2.5006628/ZValue)*(Q+Q**9+Q**25).
ELSE IF (ZValue GE 1) AND (ZValue LT 3.1).
. COMPUTE Q=EXP(-2*(ZValue**2)).
. COMPUTE PValue2=2*(Q-Q**4+Q**9-Q**16).
ELSE.
. COMPUTE PValue2=0.
END IF.
PRINT {ZValue,PValue2,PValue2/2}
 /FORMAT='F8.3'
 /CLABELS='Z Stat.','2-tail p','1-tail p'
 /TITLE='Linearity: Z statistic and asymptotic significance [Smirnov (1948)]'.
* Preparing & exporting data needed for the graphs *.
COMPUTE Residual=Y-(a+b*X).
COMPUTE vname={vname,'Residual'}.
COMPUTE Export={X,Y,Residual}.
COMPUTE Param={a,b,LowA,LowB,UppA,UppB}.
SAVE Export /OUTFILE=WorkingData /NAMES=vname.
SAVE Param  /OUTFILE=Parameters.
END MATRIX.
DATASET ACTIVATE Parameters.
PRESERVE.
SET LOCALE=ENGLISH.
STRING #var1 TO #var6 (A12).
DO REPEAT A=#var1 TO #var6 /B=col1 TO col6.
- COMPUTE A = STRING(B,F8.2).
END REPEAT.
!LET !a=   !UNQUOTE(#var1).
!LET !b=   !UNQUOTE(#var2).
!LET !alow=!UNQUOTE(#var3).
!LET !blow=!UNQUOTE(#var4).
!LET !aupp=!UNQUOTE(#var5).
!LET !bupp=!UNQUOTE(#var6).
* Write chart templates *.
WRITE OUTFILE 'C:\Temp\Residual.sgt'
 /'<?xml version="1.0" encoding="UTF-8" standalone="no"?>'
 /'<template SPSS-Version="1.4" date="2009-04-16" description="" selectPath="12 13 900 "'
 /' xmlns="http://xml.spss.com/spss/visualization"'
 /'xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"'
 /'xsi:schemaLocation="http://xml.spss.com/spss/visualization'
 /'http://xml.spss.com/spss/visualization/vizml-template-3.0.xsd">'
 /'	<setAxisStyle categorical="false" role="x">'
 /'		<style color="#000000" stroke-width="0.5pt" visible="true"/>'
 /'	</setAxisStyle>'
 /'	<setAxisStyle categorical="false" role="y">'
 /'		<style color="#000000" stroke-width="0.5pt" visible="true"/>'
 /'	</setAxisStyle>'
 /'	<addReferenceLine styleOnly="false" y="0.0" ycategorical="false">'
 /'		<style stroke-dasharray="3px,2px" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addFrame count="1" styleOnly="true" type="visualization">'
 /'		<style color="#ffffff" color2="transparent" number="0" visible="true"/>'
 /'		<style font-family="SansSerif" font-size="8pt" number="1" pattern="0" stroke-linecap="butt" text-fit="true" visible="true"/>'
 /'	</addFrame>'
 /'	<addFrame count="1" styleOnly="true" type="graph">'
 /'		<style color="transparent" color2="transparent" visible="true"/>'
 /'		<style color="#ffffff" color2="#000000" number="1" visible="true"/>'
 /'	</addFrame>'
 /'	<addFrame count="1" styleOnly="true" type="title">'
 /'		<style color="transparent" color2="transparent" visible="true"/>'
 /'	</addFrame>'
 /' <setStyle subtype="simple" type="scatter">'
 /'  <style color="#000000" color2="#000000" size="4px" symbol="circle" visible="true"/>'
 /' </setStyle>'
 /'</template>'.
WRITE OUTFILE 'C:\Temp\Regression.sgt'
 /'<?xml version="1.0" encoding="UTF-8" standalone="no"?>'
 /'<template SPSS-Version="1.4" date="2009-04-16" description="" selectPath="2 12 900 901 902 903 "'
 /'xmlns="http://xml.spss.com/spss/visualization"'
 /'xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"'
 /'xsi:schemaLocation="http://xml.spss.com/spss/visualization'
 /'http://xml.spss.com/spss/visualization/vizml-template-3.0.xsd">'
 /'	<addFrame count="1" type="visualization">'
 /'		<location bottom="533px" left="0px" right="519px" top="0px"/>'
 /'		<style color="#ffffff" color2="transparent" number="0" visible="true"/>'
 /'		<style font-family="SansSerif" font-size="8pt" number="1" pattern="0" stroke-linecap="butt" text-fit="true" visible="true"/>'
 /'	</addFrame>'
 /'	<addReferenceLine numberPoints="200" styleOnly="false" y="'!bupp' * x + '!aupp'" ycategorical="false">'
 /'		<style color="#0000ff" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addReferenceLine numberPoints="200" styleOnly="false" y="'!b' * x +  '!a'" ycategorical="false">'
 /'		<style color="#000000" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addReferenceLine numberPoints="200" styleOnly="false" y="'!blow' * x + '!alow'" ycategorical="false">'
 /'		<style color="#0000ff" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addReferenceLine numberPoints="200" styleOnly="false" y="1 * x + 0" ycategorical="false">'
 /'		<style color="#ff0000" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addFrame count="1" styleOnly="true" type="graph">'
 /'		<style color="transparent" color2="transparent" visible="true"/>'
 /'		<style color="#ffffff" color2="#000000" number="1" visible="true"/>'
 /'	</addFrame>'
 /'	<addFrame count="1" styleOnly="true" type="title">'
 /'		<style color="transparent" color2="transparent" visible="true"/>'
 /'	</addFrame>'
 /' <setStyle subtype="simple" type="scatter">'
 /'  <style color="#000000" color2="#000000" size="4px" symbol="circle" visible="true"/>'
 /' </setStyle>'
 /'</template>'.
DATASET ACTIVATE WorkingData.
DATASET CLOSE Parameters.
GRAPH /TITLE='Regression plot'
 /SCATTERPLOT(BIVAR)=!x WITH !y
 /TEMPLATE='C:\Temp\Regression.sgt'.
ECHO 'Red:   Equality line (y=x)'.
ECHO 'Black: Regression line (y=a+b*x)'.
ECHO 'Blue:  Confidence limits'.
EXE.
GRAPH /TITLE='Residual plot'
 /SCATTERPLOT(BIVAR)=!x WITH Residual
 /TEMPLATE='C:\Temp\Residual.sgt'.
RESTORE.
DATASET ACTIVATE OriginalData WINDOW=ASIS.
DATASET CLOSE WorkingData.
!ENDDEFINE.

* Sample dataset #1 (Wright meter 1st & 2nd measures from Bland&Altman (1986) Lancet paper) *.
DATA  LIST LIST/wright1 wright2(2 F8).
BEGIN DATA
494 490
395 397
516 512
434 401
476 470
557 611
413 415
442 431
650 638
433 429
417 420
656 633
267 275
478 492
178 165
423 372
427 421
END DATA.
VAR LEVEL ALL(SCALE).

* Macro call (1st example tested with Medcalc) *.
PASSINGBABLOK x=wright2 y=wright1.
PASSINGBABLOK x=wright1 y=wright2.

* Sample dataset #2 (Fleiss JL. The design and analysis of clinical experiments,
  New York: John Wiley and Sons, 1986, p 20, table 1.7) *.

DATA LIST LIST/examin1 examin2 (2 F8).
BEGIN DATA
 8  7
13 11
 0  0
 3  6
13 13
19 23
 0  0
 2  0
18 20
 5  3
END DATA.
VAR LEVEL ALL(SCALE).

* Macro call (also tested with Medcalc) *.
PASSINGBABLOK x=examin1 y=examin2.