* References:

* 1) Cornbleet, PJ, Gochman N. Incorrect least-square regression coefficients in
*    method-comparison analysis. Clin Chem, 1979; 25: 432-8
* 2) Linnet K. Estimation of the Linear Relationship Between the Measurements of Two
*    methods with Proportional Errors. Stat Med, 1990; 9: 1463-73.
* 3) Linnet K. Evaluation of Regression Procedures for Methods Comparison Studies.
*    Clin Chem, 1993; 39(3): 424-32.
* 4) Linnet K. Reference Manual CBStat 5.1 (1997-2004); p 23-32.

*************************************************
*   DEMING REGRESSION (UNWEIGHTED & WEIGHTED)   *
*************************************************
* All macros tested with SPSS 15 & PASW 17.0.2  *
*************************************************
* (C) Marta García-Granero 04/30/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 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.

*****************************************************
***************** MACRO DEFINITIONS *****************
*****************************************************

* A) MACRO FOR UNWEIGHTED DEMING REGRESSION (SINGLE MEASUREMENTS, CONSTANT ERRORS) *.

DEFINE DEMINGREG1 (x=!TOKENS(1) / y=!TOKENS(1) /lambda=!DEFAULT(1) !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.
SELECT IF (NOT MISSING(!x)) AND (NOT MISSING(!y)).
RANK VARIABLES=!x(A) /N INTO N1@ /PRINT=NO .
DO IF $casenum EQ 1.
. COMPUTE TValue = IDF.T(0.975,N1@-2) .
END IF.
TEMPORARY.
SET MXLOOPS=10000.
MATRIX.
PRINT /TITLE='****** UNWEIGHTED DEMING REGRESSION (WITH JACKKNIFE ESTIMATION OF SE&95%CI) ******'.
GET T95 /VAR=TValue /MISSING=OMIT.
* Compute user specified Lambda (variance ratio) *.
COMPUTE L=!lambda.
GET pair /VAR=!x !y /NAMES=vname /MISSING=OMIT.
PRINT {T(vname)}
 /FORMAT='A8'
 /RLABELS='Method X','Method Y'
 /TITLE='Variables compared'.
PRINT L
 /FORMAT='F8.3'
 /TITLE='Lambda (fixed by user)'.
* 1st variable is X, 2nd is Y *.
COMPUTE X=pair(:,1).
COMPUTE Y=pair(:,2).
COMPUTE n=NROW(pair).
COMPUTE MeanX=CSUM(X)/n.
COMPUTE MeanY=CSUM(Y)/n.
COMPUTE u=CSSQ(X)-((CSUM(X))**2)/n.
COMPUTE q=CSSQ(Y)-((CSUM(Y))**2)/n.
COMPUTE p=CSUM(X&*Y)-CSUM(X)*CSUM(Y)/n.
COMPUTE R=p/SQRT(u*q).
COMPUTE b=((L*q-u)+SQRT((u-L*q)**2+4*L*p**2))/(2*L*p).
COMPUTE a=MeanY.
COMPUTE a0=a-b*MeanX.
* Compute empty matrix to store JK statistics (a0i&bi) *.
COMPUTE jack=MAKE(n,2,0).
COMPUTE nj=n-1.
* Cycle thru all values *.
LOOP i=1 TO n.
. DO IF (i EQ 1).                          /* JK sample for 1st position *.
.  COMPUTE sample=pair(2:n,:).
. ELSE IF (i GT 1) AND (i LT n). /* JK samples for 2nd to last-1 postion *.
.  COMPUTE sample={pair(1:(i-1),:);
                   pair((i+1):n,:)}.
. ELSE IF (i EQ n).                       /* JK sample for last position *.
.  COMPUTE sample=pair(1:(n-1),:).
. END IF.
. COMPUTE Xi=sample(:,1).
. COMPUTE Yi=sample(:,2).
. COMPUTE MeanXi=CSUM(Xi)/nj.
. COMPUTE MeanYi=CSUM(Yi)/nj.
. COMPUTE ui=CSSQ(Xi)-((CSUM(Xi))**2)/nj.
. COMPUTE qi=CSSQ(Yi)-((CSUM(Yi))**2)/nj.
. COMPUTE pi=CSUM(Xi&*Yi)-CSUM(Xi)*CSUM(Yi)/nj.
. COMPUTE ai=MeanYi.
. COMPUTE bi=((L*qi-ui)+SQRT((ui-L*qi)**2+4*L*pi**2))/(2*L*pi).
. COMPUTE jack(i,1)=n*(a-b*MeanX)-(n-1)*(ai-bi*MeanXi).
. COMPUTE jack(i,2)=n*b-(n-1)*bi.
END LOOP.
COMPUTE MeanJack=CSUM(jack)/n.
COMPUTE VarJack=(CSSQ(jack)-n*(MeanJack&**2))/(n-1).
COMPUTE SEJack=SQRT(VarJack/n).
COMPUTE Low={a0,b}-T95*SEJack.
COMPUTE Upp={a0,b}+T95*SEJack.
COMPUTE TValue=({a0,b}-{0,1})/SEJack.
COMPUTE PValue=2*(1-TCDF(ABS(TValue),n-2)).
PRINT {T({a0,b}),T(SEJack),T(Low),T(Upp),T(TValue),T(PValue)}
 /FORMAT='F8.5'
 /RLABEL='A=','B='
 /CLABEL='Coeff.','SE(coef)','Low95%CL','Upp95%CL','T Value','2-tail p'
 /TITLE='Regression line & Jackknife estimators of SE & 95CI'.
PRINT R
 /FORMAT='F8.3'
 /TITLE='Product-moment coefficient of correlation'.
* Preparing data to export *.
* 1) Standardized residuals *.
COMPUTE d=Y-(a0+b*X).
COMPUTE Xesti=X+(L*b*d)/(1+L*(b**2)).
COMPUTE Yesti=Y-(d/(1+L*(b**2))).
COMPUTE Sign=(Y-Yesti)/ABS(Y-Yesti).
COMPUTE Ri=Sign&*SQRT((X-Xesti)&**2+L*(Y-Yesti)&**2).
COMPUTE SDr=SQRT(CSSQ(Ri)/(n-2)).
COMPUTE Residual=Ri/SDr.
COMPUTE Mean=(Xesti+L*Yesti)/(1+L).
COMPUTE vname={'Method1','Method2','Mean','Residual'}.
COMPUTE Export={X,Y,Mean,Residual}.
SAVE Export /OUTFILE=WorkingData /NAMES=vname.
* 2) Regression line *.
COMPUTE Param={a0,b}.
SAVE Param  /OUTFILE=Parameters.
END MATRIX.
DATASET ACTIVATE Parameters.
PRESERVE.
SET LOCALE=ENGLISH.
STRING #var1 #var2 (A12).
DO REPEAT A=#var1 #var2 /B=col1 col2.
- COMPUTE A = STRING(B,F8.2).
END REPEAT.
!LET !a=   !UNQUOTE(#var1).
!LET !b=   !UNQUOTE(#var2).
* 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">'
 /'	<setAxisInterval lowerMapping="exact" max="5.0" min="-5.0" role="y" upperMapping="exact"/>'
 /'	<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>'
 /'	<addReferenceLine styleOnly="false" y="4.0" ycategorical="false">'
 /'		<style stroke-dasharray="3px,2px" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addReferenceLine styleOnly="false" y="-4.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="'!b' * x +  '!a'" ycategorical="false">'
 /'		<style color="#000000" 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.
VAR LABEL Method1'Method X'/Method2 'Method Y'.
GRAPH /TITLE='Regression plot'
 /SCATTERPLOT(BIVAR)=Method1 WITH Method2
 /TEMPLATE='C:\Temp\Regression.sgt'.
ECHO 'Red:   Equality line (y=x)'.
ECHO 'Black: Regression line (y=a+b*x)'.
EXE.
VAR LABEL Residual 'Standardized Residual' /Mean 'Average'.
GRAPH /TITLE='Residual plot'
 /SCATTERPLOT(BIVAR)=Mean WITH Residual
 /TEMPLATE='C:\Temp\Residual.sgt'.
RESTORE.
COMPUTE ID=$casenum.
FORMAT ID(F8).
COMPUTE AbsRes=ABS(Residual).
SORT CASES BY AbsRes(D).
ECHO '5 biggest residuals (absolute values)'.
LIST VARIABLES=ID Residual /CASES=FROM 1 TO 5.
ECHO 'Values over 4 are outliers'.
ECHO '-------------------------------------------'.
ECHO 'Linearity: RUNS test'.
SORT CASES BY Mean (A) .
NPAR TESTS /RUNS(MEDIAN)=Residual.
EXAMINE
  VARIABLES=Residual
  /PLOT BOXPLOT NPPLOT.
DATASET ACTIVATE OriginalData WINDOW=ASIS.
DATASET CLOSE WorkingData.
!ENDDEFINE.

* B) MACRO FOR WEIGHTED DEMING REGRESSION (SINGLE MEASUREMENTS, PROPORTIONAL ERRORS) *.

DEFINE DEMINGREG2 (x=!TOKENS(1) / y=!TOKENS(1) /lambda=!DEFAULT(1) !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.
COUNT nmiss = !x !y (SYSMIS) .
SELECT IF (nmiss=0). 
RANK VARIABLES=!x(A) /N INTO N1@ /PRINT=NO .
DO IF $casenum EQ 1.
. COMPUTE TValue = IDF.T(0.975,N1@-2) .
END IF.
TEMPORARY.
SET MXLOOPS=10000.
MATRIX.
PRINT /TITLE='****** WEIGHTED DEMING REGRESSION (WITH JACKKNIFE ESTIMATION OF SE&95%CI) ******'.
PRINT/TITLE='Single measurements & proportional measurement errors'.
GET T95 /VAR=TValue /MISSING=OMIT.
* Maximum number of iterations for estimation of weighted coefficients (usually only 3-4 are needed) *.
COMPUTE MxIter=25.
* Compute user specified Lambda (variance ratio) *.
COMPUTE L=!lambda.
GET pair /VAR=!x !y/NAMES=vname /MISSING=OMIT.
PRINT {T(vname)}
 /FORMAT='A8'
 /RLABELS='Method X','Method Y'
 /TITLE='Variables compared'.
PRINT L
 /FORMAT='F8.3'
 /TITLE='Lambda (fixed by user)'.
* 1st variable is X, 2nd is Y *.
COMPUTE X=pair(:,1).
COMPUTE Y=pair(:,2).
COMPUTE n=NROW(pair).
COMPUTE MeanX=CSUM(X)/n.
COMPUTE MeanY=CSUM(Y)/n.
COMPUTE u=CSSQ(X)-((CSUM(X))**2)/n.
COMPUTE q=CSSQ(Y)-((CSUM(Y))**2)/n.
COMPUTE p=CSUM(X&*Y)-CSUM(X)*CSUM(Y)/n.
COMPUTE b=((L*q-u)+SQRT((u-L*q)**2+4*L*p**2))/(2*L*p).
COMPUTE a=MeanY.
COMPUTE a0=a-b*MeanX.
* Weighted regression: begin iteration process *.
LOOP jk=1 TO MxIter.
. COMPUTE B0=b.
. COMPUTE di=Y-(a0+b*X).
. COMPUTE Xhati=X+(L*b*di)/(1+L*(b**2)).
. COMPUTE Yhati=Y-(di/(1+L*(b**2))).
. COMPUTE Whati=((Xhati+L*Yhati)/(1+L))&**(-2).
. COMPUTE AverXiYi=(Xhati+Yhati)/2.
. COMPUTE MeanXw=CSUM(Whati&*X)/CSUM(Whati).
. COMPUTE MeanYw=CSUM(Whati&*Y)/CSUM(Whati).
. COMPUTE uw=CSUM(Whati&*((X-MeanXw)&**2)).
. COMPUTE qw=CSUM(Whati&*((Y-MeanYw)&**2)).
. COMPUTE pw=CSUM(Whati&*(X-MeanXw)&*(Y-MeanYw)).
. COMPUTE b=((L*qw-uw)+SQRT((uw-L*qw)**2+4*L*pw**2))/(2*L*pw).
. COMPUTE a=MeanYw.
. COMPUTE a0=a-b*MeanXw.
. COMPUTE diff=ABS(B0-b).
END LOOP IF (diff LT 1.0E-5).
COMPUTE Rw=pw/SQRT(uw*qw).
* Jackknife estimation of standard errors *.
* Compute empty matrix to store JK statistics (a0i&bi) *.
COMPUTE jack=MAKE(n,2,0).
COMPUTE nj=n-1.
* Cycle thru all values *.
LOOP i=1 TO n.
. DO IF (i EQ 1).                          /* JK sample for 1st position *.
.  COMPUTE sample=pair(2:n,:).
. ELSE IF (i GT 1) AND (i LT n). /* JK samples for 2nd to last-1 postion *.
.  COMPUTE sample={pair(1:(i-1),:);
                   pair((i+1):n,:)}.
. ELSE IF (i EQ n).                       /* JK sample for last position *.
.  COMPUTE sample=pair(1:(n-1),:).
. END IF.
. COMPUTE Xi=sample(:,1).
. COMPUTE Yi=sample(:,2).
. COMPUTE MeanXi=CSUM(Xi)/nj.
. COMPUTE MeanYi=CSUM(Yi)/nj.
. COMPUTE ui=CSSQ(Xi)-((CSUM(Xi))**2)/nj.
. COMPUTE qi=CSSQ(Yi)-((CSUM(Yi))**2)/nj.
. COMPUTE pi=CSUM(Xi&*Yi)-CSUM(Xi)*CSUM(Yi)/nj.
. COMPUTE bi=((L*qi-ui)+SQRT((ui-L*qi)**2+4*L*pi**2))/(2*L*pi).
. COMPUTE a0i=MeanYi-bi*Meanxi.
* Begin iteration process for each jackknife pseudovariate *.
. LOOP jk=1 TO MxIter.
.  COMPUTE B0=bi.
.  COMPUTE di=Yi-(a0i+bi*Xi).
.  COMPUTE Xhati=Xi+(L*bi*di)/(1+L*(bi**2)).
.  COMPUTE Yhati=Yi-(di/(1+L*(bi**2))).
.  COMPUTE Whati=((Xhati+L*Yhati)/(1+L))&**(-2).
.  COMPUTE AverXiYi=(Xhati+Yhati)/2.
.  COMPUTE MeanXwi=CSUM(Whati&*Xi)/CSUM(Whati).
.  COMPUTE MeanYwi=CSUM(Whati&*Yi)/CSUM(Whati).
.  COMPUTE uwi=CSUM(Whati&*((Xi-MeanXwi)&**2)).
.  COMPUTE qwi=CSUM(Whati&*((Yi-MeanYwi)&**2)).
.  COMPUTE pwi=CSUM(Whati&*(Xi-MeanXwi)&*(Yi-MeanYwi)).
.  COMPUTE bi=((L*qwi-uwi)+SQRT((uwi-L*qwi)**2+4*L*pwi**2))/(2*L*pwi).
.  COMPUTE a0i=MeanYwi-bi*MeanXwi.
.  COMPUTE diff=ABS(B0-bi).
. END LOOP IF (diff LT 1.0E-5).
. COMPUTE jack(i,1)=n*a0-(n-1)*a0i.
. COMPUTE jack(i,2)=n*b-(n-1)*bi.
END LOOP.
COMPUTE MeanJack=CSUM(jack)/n.
COMPUTE VarJack=(CSSQ(jack)-n*(MeanJack&**2))/(n-1).
COMPUTE SEJack=SQRT(VarJack/n).
COMPUTE Low={a0,b}-T95*SEJack.
COMPUTE Upp={a0,b}+T95*SEJack.
COMPUTE TValue=({a0,b}-{0,1})/SEJack.
COMPUTE PValue=2*(1-TCDF(ABS(TValue),n-2)).
PRINT {T({a0,b}),T(SEJack),T(Low),T(Upp),T(TValue),T(PValue)}
 /FORMAT='F8.5'
 /RLABEL='A=','B='
 /CLABEL='Coeff.','SE(coef)','Low95%CL','Upp95%CL','T Value','2-tail p'
 /TITLE='Weighted regression line & Jackknife estimators of SE & 95CI'.
PRINT Rw
 /FORMAT='F8.3'
 /TITLE='Weighted product-moment coefficient of correlation'.
* Preparing data to export *.
* 1) Standardized residuals *.
COMPUTE d=Y-(a0+b*X).
COMPUTE Xesti=X+(L*b*d)/(1+L*(b**2)).
COMPUTE Yesti=Y-(d/(1+L*(b**2))).
COMPUTE Wi=((Xesti+L*Yesti)/(1+L))&**(-2).
COMPUTE Sign=(Y-Yesti)/ABS(Y-Yesti).
COMPUTE Ri=Sign&*SQRT(Wi&*(X-Xesti)&**2+L*Wi&*(Y-Yesti)&**2).
COMPUTE SDr=SQRT(CSSQ(Ri)/(n-2)).
COMPUTE Residual=Ri/SDr.
COMPUTE Mean=(Xesti+L*Yesti)/(1+L).
COMPUTE vname={'Method1','Method2','Mean','Residual'}.
COMPUTE Export={X,Y,Mean,Residual}.
SAVE Export /OUTFILE=WorkingData /NAMES=vname.
* 2) Regression line *.
COMPUTE Param={a0,b}.
SAVE Param  /OUTFILE=Parameters.
END MATRIX.
DATASET ACTIVATE Parameters.
PRESERVE.
SET LOCALE=ENGLISH.
STRING #var1 #var2 (A12).
DO REPEAT A=#var1 #var2 /B=col1 col2.
- COMPUTE A = STRING(B,F8.2).
END REPEAT.
!LET !a=   !UNQUOTE(#var1).
!LET !b=   !UNQUOTE(#var2).
* 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">'
 /'	<setAxisInterval lowerMapping="exact" max="5.0" min="-5.0" role="y" upperMapping="exact"/>'
 /'	<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>'
 /'	<addReferenceLine styleOnly="false" y="4.0" ycategorical="false">'
 /'		<style stroke-dasharray="3px,2px" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addReferenceLine styleOnly="false" y="-4.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="'!b' * x +  '!a'" ycategorical="false">'
 /'		<style color="#000000" 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.
VAR LABEL Method1'Method X'/Method2 'Method Y'.
GRAPH /TITLE='Regression plot'
 /SCATTERPLOT(BIVAR)=Method1 WITH Method2
 /TEMPLATE='C:\Temp\Regression.sgt'.
ECHO 'Red:   Equality line (y=x)'.
ECHO 'Black: Regression line (y=a+b*x)'.
EXE.
VAR LABEL Residual 'Weighted Standardized Residual' /Mean 'Weighted Average'.
GRAPH /TITLE='Residual plot'
 /SCATTERPLOT(BIVAR)=Mean WITH Residual
 /TEMPLATE='C:\Temp\Residual.sgt'.
RESTORE.
COMPUTE ID=$casenum.
FORMAT ID(F8).
COMPUTE AbsRes=ABS(Residual).
SORT CASES BY AbsRes(D).
ECHO '5 biggest residuals (absolute values)'.
LIST VARIABLES=ID Residual /CASES=FROM 1 TO 5.
ECHO 'Values over 4 are outliers'.
ECHO '-------------------------------------------'.
ECHO 'Linearity: RUNS test'.
SORT CASES BY Mean (A) .
NPAR TESTS /RUNS(MEDIAN)=Residual.
EXAMINE
  VARIABLES=Residual
  /PLOT BOXPLOT NPPLOT.
DATASET ACTIVATE OriginalData WINDOW=ASIS.
DATASET CLOSE WorkingData.
!ENDDEFINE.

* C) MACRO FOR UNWEIGHTED DEMING REGRESSION (PAIRS OF MEASUREMENTS WITH CONSTANT ERRORS) *.

DEFINE DEMINGREG3 (x=!TOKENS(2) / y=!TOKENS(2)).
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.
COUNT nmiss = !x !y (SYSMIS) .
SELECT IF (nmiss=0). 
RANK VARIABLES=!x(A) /N INTO N1@ /PRINT=NO .
DO IF $casenum EQ 1.
. COMPUTE TValue = IDF.T(0.975,N1@-2) .
END IF.
TEMPORARY.
SET MXLOOPS=10000.
MATRIX.
PRINT /TITLE='*** UNWEIGHTED DEMING REGRESSION (WITH JACKKNIFE ESTIMATION OF SE&95%CI) ***'.
PRINT /TITLE='Pairs of measurements (direct estimation of Lambda)'.
PRINT /TITLE='Constant error of measurement is assumed for both methods'.
GET T95 /VAR=TValue /MISSING=OMIT.
GET pair1 /VAR=!x /NAMES=vname1.
GET pair2 /VAR=!y /NAMES=vname2.
PRINT {vname1;vname2}
 /FORMAT='A8'
 /RLABELS='Method X','Method Y'
 /CLABELS='--1st---','--2nd---'
 /TITLE='Variables compared'.
* 1st pair of variables form X, 2nd pair form Y *.
COMPUTE X=RSUM(pair1)/2.
COMPUTE Y=RSUM(pair2)/2.
COMPUTE n=NROW(pair1).
COMPUTE MeanX=CSUM(X)/n.
COMPUTE MeanY=CSUM(Y)/n.
COMPUTE DiffX=pair1(:,1)-pair1(:,2).
COMPUTE DiffY=pair2(:,1)-pair2(:,2).
COMPUTE Varx=CSUM(DiffX&**2)/(2*n).
COMPUTE Vary=CSUM(DiffY&**2)/(2*n).
COMPUTE Lhat=VarX/VarY.
PRINT {MeanX,VarX,100*SQRT(VarX)/MeanX;MeanY,VarY,100*SQRT(VarY)/MeanY}
 /FORMAT='F8.3'
 /RLABELS='Method X','Method Y'
 /CLABELS='Mean','SDa²','CV(%)'
 /TITLE='Statistics for both methods'.
PRINT Lhat
 /FORMAT='F8.3'
 /TITLE='Variance ratio (Lambda)'.
COMPUTE u=CSSQ(X)-((CSUM(X))**2)/n.
COMPUTE q=CSSQ(Y)-((CSUM(Y))**2)/n.
COMPUTE p=CSUM(X&*Y)-CSUM(X)*CSUM(Y)/n.
COMPUTE R=p/SQRT(u*q).
COMPUTE b=((Lhat*q-u)+SQRT((u-Lhat*q)**2+4*Lhat*p**2))/(2*Lhat*p).
COMPUTE a=MeanY.
COMPUTE a0=a-b*MeanX.
* Compute empty matrix to store JK statistics (a0i&bi) *.
COMPUTE AllSamp={X,Y,DiffX,DiffY}.
COMPUTE jack=MAKE(n,2,0).
COMPUTE nj=n-1.
* Cycle thru all values *.
LOOP i=1 TO n.
. DO IF (i EQ 1).                          /* JK sample for 1st position *.
.  COMPUTE sample=AllSamp(2:n,:).
. ELSE IF (i GT 1) AND (i LT n). /* JK samples for 2nd to last-1 postion *.
.  COMPUTE sample={AllSamp(1:(i-1),:);
                   AllSamp((i+1):n,:)}.
. ELSE IF (i EQ n).                       /* JK sample for last position *.
.  COMPUTE sample=AllSamp(1:(n-1),:).
. END IF.
. COMPUTE Xi=sample(:,1).
. COMPUTE Yi=sample(:,2).
. COMPUTE DiffXi=sample(:,3).
. COMPUTE DiffYi=sample(:,4).
. COMPUTE Varxi=CSUM(DiffXi&**2)/(2*n).
. COMPUTE Varyi=CSUM(DiffYi&**2)/(2*n).
. COMPUTE Lhati=VarXi/VarYi.
. COMPUTE MeanXi=CSUM(Xi)/nj.
. COMPUTE MeanYi=CSUM(Yi)/nj.
. COMPUTE ui=CSSQ(Xi)-((CSUM(Xi))**2)/nj.
. COMPUTE qi=CSSQ(Yi)-((CSUM(Yi))**2)/nj.
. COMPUTE pi=CSUM(Xi&*Yi)-CSUM(Xi)*CSUM(Yi)/nj.
. COMPUTE ai=MeanYi.
. COMPUTE bi=((Lhati*qi-ui)+SQRT((ui-Lhati*qi)**2+4*Lhati*pi**2))/(2*Lhati*pi).
. COMPUTE jack(i,1)=n*(a-b*MeanX)-(n-1)*(ai-bi*MeanXi).
. COMPUTE jack(i,2)=n*b-(n-1)*bi.
END LOOP.
COMPUTE MeanJack=CSUM(jack)/n.
COMPUTE VarJack=(CSSQ(jack)-n*(MeanJack&**2))/(n-1).
COMPUTE SEJack=SQRT(VarJack/n).
COMPUTE Low={a0,b}-T95*SEJack.
COMPUTE Upp={a0,b}+T95*SEJack.
COMPUTE TValue=({a0,b}-{0,1})/SEJack.
COMPUTE PValue=2*(1-TCDF(ABS(TValue),n-2)).
PRINT {T({a0,b}),T(SEJack),T(Low),T(Upp),T(TValue),T(PValue)}
 /FORMAT='F8.4'
 /RLABEL='A=','B='
 /CLABEL='Coeff.','SE(coef)','Low95%CL','Upp95%CL','T Value','2-tail p'
 /TITLE='Regression line & Jackknife estimators of SE & 95CI'.
PRINT R
 /FORMAT='F8.3'
 /TITLE='Product-moment coefficient of correlation'.
* Preparing data to export *.
* 1) Standardized residuals *.
COMPUTE d=Y-(a0+b*X).
COMPUTE Xesti=X+Lhat*b*d/(1+Lhat*(b**2)).
COMPUTE Yesti=Y-d/(1+Lhat*(b**2)).
COMPUTE Sign=(Y-Yesti)/ABS(Y-Yesti).
COMPUTE Ri=Sign&*SQRT((X-Xesti)&**2+Lhat*(Y-Yesti)&**2).
COMPUTE SDr=SQRT(CSSQ(Ri)/(n-2)).
COMPUTE Residual=Ri/SDr.
COMPUTE Mean=(Xesti+Lhat*Yesti)/(1+Lhat).
COMPUTE vname={'Method1','Method2','Mean','Residual'}.
COMPUTE Export={X,Y,Mean,Residual}.
SAVE Export /OUTFILE=WorkingData /NAMES=vname.
* 2) Regression line *.
COMPUTE Param={a0,b}.
SAVE Param  /OUTFILE=Parameters.
END MATRIX.
DATASET ACTIVATE Parameters.
PRESERVE.
SET LOCALE=ENGLISH.
STRING #var1 #var2 (A12).
DO REPEAT A=#var1 #var2 /B=col1 col2.
- COMPUTE A = STRING(B,F8.2).
END REPEAT.
!LET !a=   !UNQUOTE(#var1).
!LET !b=   !UNQUOTE(#var2).
* 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">'
 /'	<setAxisInterval lowerMapping="exact" max="5.0" min="-5.0" role="y" upperMapping="exact"/>'
 /'	<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>'
 /'	<addReferenceLine styleOnly="false" y="4.0" ycategorical="false">'
 /'		<style stroke-dasharray="3px,2px" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addReferenceLine styleOnly="false" y="-4.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="'!b' * x +  '!a'" ycategorical="false">'
 /'		<style color="#000000" 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.
VAR LABEL Method1'Method X'/Method2 'Method Y'.
GRAPH /TITLE='Regression plot'
 /SCATTERPLOT(BIVAR)=Method1 WITH Method2
 /TEMPLATE='C:\Temp\Regression.sgt'.
ECHO 'Red:   Equality line (y=x)'.
ECHO 'Black: Regression line (y=a+b*x)'.
EXE.
VAR LABEL Residual 'Weighted Standardized Residual' /Mean 'Weighted Average'.
GRAPH /TITLE='Residual plot'
 /SCATTERPLOT(BIVAR)=Mean WITH Residual
 /TEMPLATE='C:\Temp\Residual.sgt'.
RESTORE.
COMPUTE ID=$casenum.
FORMAT ID(F8).
COMPUTE AbsRes=ABS(Residual).
SORT CASES BY AbsRes(D).
ECHO '5 biggest residuals (absolute values)'.
LIST VARIABLES=ID Residual /CASES=FROM 1 TO 5.
ECHO 'Values over 4 are outliers'.
ECHO '-------------------------------------------'.
ECHO 'Linearity: RUNS test'.
SORT CASES BY Mean (A) .
NPAR TESTS /RUNS(MEDIAN)=Residual.
EXAMINE
  VARIABLES=Residual
  /PLOT BOXPLOT NPPLOT.
DATASET ACTIVATE OriginalData WINDOW=ASIS.
DATASET CLOSE WorkingData.
!ENDDEFINE.

* D) MACRO FOR WEIGHTED DEMING REGRESSION (PAIRS OF MEASUREMENTS WITH PROPORTIONAL ERRORS) *.

DEFINE DEMINGREG4 (x=!TOKENS(2) / y=!TOKENS(2)).
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.
COUNT nmiss = !x !y (SYSMIS) .
SELECT IF (nmiss=0). 
RANK VARIABLES=!x(A) /N INTO N1@ /PRINT=NO .
DO IF $casenum EQ 1.
. COMPUTE TValue = IDF.T(0.975,N1@-2) .
END IF.
TEMPORARY.
SET MXLOOPS=10000.
MATRIX.
PRINT /TITLE='****** WEIGHTED DEMING REGRESSION (WITH JACKKNIFE ESTIMATION OF SE&95%CI) ******'.
PRINT /TITLE='Pairs of measurements (direct estimation of Lambda)'.
PRINT /TITLE='Measurement error proportional to squared mean value'.
* Maximum number of iterations for estimation of weighted coefficients (usually only 3-4 are needed) *.
COMPUTE MxIter=25.
* Get all working data *.
GET T95 /VAR=TValue /MISSING=OMIT.
GET pair1 /VAR=!x /NAMES=vname1.
GET pair2 /VAR=!y /NAMES=vname2.
PRINT {vname1;vname2}
 /FORMAT='A8'
 /RLABELS='Method X','Method Y'
 /CLABELS='--1st---','--2nd---'
 /TITLE='Variables compared'.
* 1st pair of variables form X, 2nd pair form Y *.
COMPUTE X=RSUM(pair1)/2.
COMPUTE Y=RSUM(pair2)/2.
COMPUTE AverXY=(X+Y)/2.
COMPUTE n=NROW(pair1).
COMPUTE DiffX=pair1(:,1)-pair1(:,2).
COMPUTE DiffY=pair2(:,1)-pair2(:,2).
COMPUTE VarX=CSUM(DiffX&**2)/(2*n).
COMPUTE VarY=CSUM(DiffY&**2)/(2*n).
COMPUTE Fhatx2=CSUM((DiffX&/AverXY)&**2)/(2*n).
COMPUTE Fhaty2=CSUM((DiffY&/AverXY)&**2)/(2*n).
COMPUTE Lhat=Fhatx2/Fhaty2.
COMPUTE MeanX=CSUM(X)/n.
COMPUTE MeanY=CSUM(Y)/n.
COMPUTE u=CSSQ(X)-((CSUM(X))**2)/n.
COMPUTE q=CSSQ(Y)-((CSUM(Y))**2)/n.
COMPUTE p=CSUM(X&*Y)-CSUM(X)*CSUM(Y)/n.
PRINT {MeanX,VarX;MeanY,VarY}
 /FORMAT='F8.3'
 /RLABELS='Method X','Method Y'
 /CLABELS='Mean','SDa²'
 /TITLE='Statistics for both methods'.
* Unweighted estimation (starting point) *.
COMPUTE b=((Lhat*q-u)+SQRT((u-Lhat*q)**2+4*Lhat*p**2))/(2*Lhat*p).
COMPUTE a=MeanY.
COMPUTE a0=a-b*MeanX.
PRINT {100*SQRT(Fhatx2),100*SQRT(Fhaty2),Lhat,a0,b}
 /FORMAT='F8.3'
 /CLABEL='Fx(%)','Fy(%)','Lambda','A','B'
 /TITLE='Initial (unweighted) estimates'.
* Begin iteration process *.
LOOP jk=1 TO MxIter.
. COMPUTE B0=b.
. COMPUTE di=Y-(a0+b*X).
. COMPUTE Xhati=X+Lhat*b*di/(1+Lhat*(b**2)).
. COMPUTE Yhati=Y-di/(1+Lhat*(b**2)).
. COMPUTE Whati=((Xhati+Lhat*Yhati)/(1+Lhat))&**(-2).
. COMPUTE AverXiYi=(Xhati+Yhati)/2.
. COMPUTE Fhatx2=CSUM((DiffX&/AverXiYi)&**2)/(2*n).
. COMPUTE Fhaty2=CSUM((DiffY&/AverXiYi)&**2)/(2*n).
. COMPUTE Lhat=Fhatx2/Fhaty2.
. COMPUTE MeanXw=CSUM(Whati&*X)/CSUM(Whati).
. COMPUTE MeanYw=CSUM(Whati&*Y)/CSUM(Whati).
. COMPUTE uw=CSUM(Whati&*((X-MeanXw)&**2)).
. COMPUTE qw=CSUM(Whati&*((Y-MeanYw)&**2)).
. COMPUTE pw=CSUM(Whati&*(X-MeanXw)&*(Y-MeanYw)).
. COMPUTE b=((Lhat*qw-uw)+SQRT((uw-Lhat*qw)**2+4*Lhat*pw**2))/(2*Lhat*pw).
. COMPUTE a=MeanYw.
. COMPUTE a0=a-b*MeanXw.
. COMPUTE diff=ABS(B0-b).
END LOOP IF (diff LT 1.0E-5).
COMPUTE Rw=pw/SQRT(uw*qw).
* Report *.
PRINT {100*SQRT(Fhatx2),100*SQRT(Fhaty2),Lhat,a0,b}
 /FORMAT='F8.3'
 /CLABEL='Fx(%)','Fy(%)','Lambda','A','B'
 /TITLE='Weighted estimates (*)'.
PRINT jk
 /FORMAT='F8.0'
 /TITLE='(*) Estimation converged (diff. LT 1E-5) at iteration:'.
* Jackknife estimation of standard errors *.
* Compute empty matrix to store JK statistics (a0i&bi) *.
COMPUTE AllSamp={X,Y,DiffX,DiffY}.
COMPUTE jack=MAKE(n,2,0).
COMPUTE nj=n-1.
* Cycle thru all values *.
LOOP i=1 TO n.
. DO IF (i EQ 1).                          /* JK sample for 1st position *.
.  COMPUTE sample=AllSamp(2:n,:).
. ELSE IF (i GT 1) AND (i LT n). /* JK samples for 2nd to last-1 postion *.
.  COMPUTE sample={AllSamp(1:(i-1),:);
                   AllSamp((i+1):n,:)}.
. ELSE IF (i EQ n).                       /* JK sample for last position *.
.  COMPUTE sample=AllSamp(1:(n-1),:).
. END IF.
. COMPUTE Xi=sample(:,1).
. COMPUTE Yi=sample(:,2).
* Unweightedjackknife estimator *.
. COMPUTE AverXiYi=(Xi+Yi)/2.
. COMPUTE DiffXi=sample(:,3).
. COMPUTE DiffYi=sample(:,4).
. COMPUTE Fhatix2=CSUM((DiffXi&/AverXiYi)&**2)/(2*nj).
. COMPUTE Fhatiy2=CSUM((DiffYi&/AverXiYi)&**2)/(2*nj).
. COMPUTE Lhati=Fhatix2/Fhatiy2.
. COMPUTE MeanXi=CSUM(Xi)/nj.
. COMPUTE MeanYi=CSUM(Yi)/nj.
. COMPUTE ui=CSSQ(Xi)-((CSUM(Xi))**2)/nj.
. COMPUTE qi=CSSQ(Yi)-((CSUM(Yi))**2)/nj.
. COMPUTE pi=CSUM(Xi&*Yi)-CSUM(Xi)*CSUM(Yi)/nj.
. COMPUTE bi=((Lhati*qi-ui)+SQRT((ui-Lhati*qi)**2+4*Lhati*pi**2))/(2*Lhati*pi).
. COMPUTE a0i=MeanYi-bi*MeanXi.
* Begin iteration process for each jackknife weighted pseudovariate *.
. LOOP jk=1 TO MxIter.
.  COMPUTE B0i=bi.
.  COMPUTE di=Yi-(a0i+bi*Xi).
.  COMPUTE Xhati=Xi+Lhati*bi*di/(1+Lhati*(bi**2)).
.  COMPUTE Yhati=Yi-di/(1+Lhati*(bi**2)).
.  COMPUTE Whati=((Xhati+Lhati*Yhati)/(1+Lhati))&**(-2).
.  COMPUTE AverXiYi=(Xhati+Yhati)/2.
.  COMPUTE Fhatix2=CSUM((DiffXi&/AverXiYi)&**2)/(2*nj).
.  COMPUTE Fhatiy2=CSUM((DiffYi&/AverXiYi)&**2)/(2*nj).
.  COMPUTE Lhati=Fhatix2/Fhatiy2.
.  COMPUTE MeanXwi=CSUM(Whati&*Xi)/CSUM(Whati).
.  COMPUTE MeanYwi=CSUM(Whati&*Yi)/CSUM(Whati).
.  COMPUTE uwi=CSUM(Whati&*((Xi-MeanXwi)&**2)).
.  COMPUTE qwi=CSUM(Whati&*((Yi-MeanYwi)&**2)).
.  COMPUTE pwi=CSUM(Whati&*(Xi-MeanXwi)&*(Yi-MeanYwi)).
.  COMPUTE bi=((Lhati*qwi-uwi)+SQRT((uwi-Lhati*qwi)**2+4*Lhati*pwi**2))/(2*Lhati*pwi).
.  COMPUTE a0i=MeanYwi-bi*MeanXwi.
.  COMPUTE diffi=ABS(B0i-bi).
. END LOOP IF (diffi LT 1.0E-5).
. COMPUTE jack(i,1)=n*a0-(n-1)*a0i.
. COMPUTE jack(i,2)=n*b-(n-1)*bi.
END LOOP.
COMPUTE MeanJack=CSUM(jack)/n.
COMPUTE VarJack=(CSSQ(jack)-n*(MeanJack&**2))/(n-1).
COMPUTE SEJack=SQRT(VarJack/n).
COMPUTE Low={a0,b}-T95*SEJack.
COMPUTE Upp={a0,b}+T95*SEJack.
COMPUTE TValue=({a0,b}-{0,1})/SEJack.
COMPUTE PValue=2*(1-TCDF(ABS(TValue),n-2)).
PRINT {T({a0,b}),T(SEJack),T(Low),T(Upp),T(TValue),T(PValue)}
 /FORMAT='F8.4'
 /RLABEL='A=','B='
 /CLABEL='Coeff.','SE(coef)','Low95%CL','Upp95%CL','T Value','2-tail p'
 /TITLE='Regression line & Jackknife estimators of SE & 95CI'.
PRINT Rw
 /FORMAT='F8.3'
 /TITLE='Weighted product-moment coefficient of correlation'.
* Preparing data to export *.
* 1) Standardized residuals *.
COMPUTE d=Y-(a0+b*X).
COMPUTE Xesti=X+Lhat*b*d/(1+Lhat*(b**2)).
COMPUTE Yesti=Y-d/(1+Lhat*(b**2)).
COMPUTE Wi=((Xesti+Lhat*Yesti)/(1+Lhat))&**(-2).
COMPUTE Sign=(Y-Yesti)/ABS(Y-Yesti).
COMPUTE Ri=Sign&*SQRT(Wi&*(X-Xesti)&**2+Lhat*Wi&*(Y-Yesti)&**2).
COMPUTE SDr=SQRT(CSSQ(Ri)/(n-2)).
COMPUTE Residual=Ri/SDr.
COMPUTE Mean=(Xesti+Lhat*Yesti)/(1+Lhat).
COMPUTE vname={'Method1','Method2','Mean','Residual'}.
COMPUTE Export={X,Y,Mean,Residual}.
SAVE Export /OUTFILE=WorkingData /NAMES=vname.
* 2) Regression line *.
COMPUTE Param={a0,b}.
SAVE Param  /OUTFILE=Parameters.
END MATRIX.
DATASET ACTIVATE Parameters.
PRESERVE.
SET LOCALE=ENGLISH.
STRING #var1 #var2 (A12).
DO REPEAT A=#var1 #var2 /B=col1 col2.
- COMPUTE A = STRING(B,F8.2).
END REPEAT.
!LET !a=   !UNQUOTE(#var1).
!LET !b=   !UNQUOTE(#var2).
* 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">'
 /'	<setAxisInterval lowerMapping="exact" max="5.0" min="-5.0" role="y" upperMapping="exact"/>'
 /'	<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>'
 /'	<addReferenceLine styleOnly="false" y="4.0" ycategorical="false">'
 /'		<style stroke-dasharray="3px,2px" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addReferenceLine styleOnly="false" y="-4.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="'!b' * x +  '!a'" ycategorical="false">'
 /'		<style color="#000000" 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.
VAR LABEL Method1 'Method X'/Method2 'Method Y'.
GRAPH /TITLE='Regression plot'
 /SCATTERPLOT(BIVAR)=Method1 WITH Method2
 /TEMPLATE='C:\Temp\Regression.sgt'.
ECHO 'Red:   Equality line (y=x)'.
ECHO 'Black: Regression line (y=a+b*x)'.
EXE.
VAR LABEL Residual 'Weighted Standardized Residual' /Mean 'Weighted Average'.
GRAPH /TITLE='Residual plot'
 /SCATTERPLOT(BIVAR)=Mean WITH Residual
 /TEMPLATE='C:\Temp\Residual.sgt'.
RESTORE.
COMPUTE ID=$casenum.
FORMAT ID(F8).
COMPUTE AbsRes=ABS(Residual).
SORT CASES BY AbsRes(D).
ECHO '5 biggest residuals (absolute values)'.
LIST VARIABLES=ID Residual /CASES=FROM 1 TO 5.
ECHO 'Values over 4 are outliers'.
ECHO '-------------------------------------------'.
ECHO 'Linearity: RUNS test'.
SORT CASES BY Mean (A) .
NPAR TESTS /RUNS(MEDIAN)=Residual.
EXAMINE
  VARIABLES=Residual
  /PLOT BOXPLOT NPPLOT.
DATASET ACTIVATE OriginalData WINDOW=ASIS.
DATASET CLOSE WorkingData.
!ENDDEFINE.

* E) MACRO FOR WEIGHTED DEMING REGRESSION (SINGLE MEASUREMENTS, SQUARE ROOT WEIGHTING) *.

DEFINE DEMINGREG5 (x=!TOKENS(1) / y=!TOKENS(1) /lambda=!DEFAULT(1) !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.
COUNT nmiss = !x !y (SYSMIS) .
SELECT IF (nmiss=0). 
RANK VARIABLES=!x(A) /N INTO N1@ /PRINT=NO .
DO IF $casenum EQ 1.
. COMPUTE TValue = IDF.T(0.975,N1@-2) .
END IF.
TEMPORARY.
SET MXLOOPS=10000.
MATRIX.
PRINT /TITLE='****** WEIGHTED DEMING REGRESSION (WITH JACKKNIFE ESTIMATION OF SE&95%CI) ******'.
PRINT/TITLE='Single measurements & non constant measurement errors'.
PRINT/TITLE='Square root weighting (e.g. used when counting isotopes)'.
GET T95 /VAR=TValue /MISSING=OMIT.
* Maximum number of iterations for estimation of weighted coefficients (usually only 3-4 are needed) *.
COMPUTE MxIter=25.
* Compute user specified Lambda (variance ratio) *.
COMPUTE L=!lambda.
GET pair /VAR=!x !y/NAMES=vname /MISSING=OMIT.
PRINT {T(vname)}
 /FORMAT='A8'
 /RLABELS='Method X','Method Y'
 /TITLE='Variables compared'.
PRINT L
 /FORMAT='F8.3'
 /TITLE='Lambda (fixed by user)'.
* 1st variable is X, 2nd is Y *.
COMPUTE X=pair(:,1).
COMPUTE Y=pair(:,2).
COMPUTE n=NROW(pair).
COMPUTE MeanX=CSUM(X)/n.
COMPUTE MeanY=CSUM(Y)/n.
COMPUTE u=CSSQ(X)-((CSUM(X))**2)/n.
COMPUTE q=CSSQ(Y)-((CSUM(Y))**2)/n.
COMPUTE p=CSUM(X&*Y)-CSUM(X)*CSUM(Y)/n.
COMPUTE b=((L*q-u)+SQRT((u-L*q)**2+4*L*p**2))/(2*L*p).
COMPUTE a=MeanY.
COMPUTE a0=a-b*MeanX.
* Weighted regression: begin iteration process *.
LOOP jk=1 TO MxIter.
. COMPUTE B0=b.
. COMPUTE di=Y-(a0+b*X).
. COMPUTE Xhati=X+(L*b*di)/(1+L*(b**2)).
. COMPUTE Yhati=Y-(di/(1+L*(b**2))).
. COMPUTE Whati=((Xhati+L*Yhati)/(1+L))&**(-1).
. COMPUTE AverXiYi=(Xhati+Yhati)/2.
. COMPUTE MeanXw=CSUM(Whati&*X)/CSUM(Whati).
. COMPUTE MeanYw=CSUM(Whati&*Y)/CSUM(Whati).
. COMPUTE uw=CSUM(Whati&*((X-MeanXw)&**2)).
. COMPUTE qw=CSUM(Whati&*((Y-MeanYw)&**2)).
. COMPUTE pw=CSUM(Whati&*(X-MeanXw)&*(Y-MeanYw)).
. COMPUTE b=((L*qw-uw)+SQRT((uw-L*qw)**2+4*L*pw**2))/(2*L*pw).
. COMPUTE a=MeanYw.
. COMPUTE a0=a-b*MeanXw.
. COMPUTE diff=ABS(B0-b).
END LOOP IF (diff LT 1.0E-5).
COMPUTE Rw=pw/SQRT(uw*qw).
* Jackknife estimation of standard errors *.
* Compute empty matrix to store JK statistics (a0i&bi) *.
COMPUTE jack=MAKE(n,2,0).
COMPUTE nj=n-1.
* Cycle thru all values *.
LOOP i=1 TO n.
. DO IF (i EQ 1).                          /* JK sample for 1st position *.
.  COMPUTE sample=pair(2:n,:).
. ELSE IF (i GT 1) AND (i LT n). /* JK samples for 2nd to last-1 postion *.
.  COMPUTE sample={pair(1:(i-1),:);
                   pair((i+1):n,:)}.
. ELSE IF (i EQ n).                       /* JK sample for last position *.
.  COMPUTE sample=pair(1:(n-1),:).
. END IF.
. COMPUTE Xi=sample(:,1).
. COMPUTE Yi=sample(:,2).
. COMPUTE MeanXi=CSUM(Xi)/nj.
. COMPUTE MeanYi=CSUM(Yi)/nj.
. COMPUTE ui=CSSQ(Xi)-((CSUM(Xi))**2)/nj.
. COMPUTE qi=CSSQ(Yi)-((CSUM(Yi))**2)/nj.
. COMPUTE pi=CSUM(Xi&*Yi)-CSUM(Xi)*CSUM(Yi)/nj.
. COMPUTE bi=((L*qi-ui)+SQRT((ui-L*qi)**2+4*L*pi**2))/(2*L*pi).
. COMPUTE a0i=MeanYi-bi*Meanxi.
* Begin iteration process for each jackknife pseudovariate *.
. LOOP jk=1 TO MxIter.
.  COMPUTE B0=bi.
.  COMPUTE di=Yi-(a0i+bi*Xi).
.  COMPUTE Xhati=Xi+(L*bi*di)/(1+L*(bi**2)).
.  COMPUTE Yhati=Yi-(di/(1+L*(bi**2))).
.  COMPUTE Whati=((Xhati+L*Yhati)/(1+L))&**(-1).
.  COMPUTE AverXiYi=(Xhati+Yhati)/2.
.  COMPUTE MeanXwi=CSUM(Whati&*Xi)/CSUM(Whati).
.  COMPUTE MeanYwi=CSUM(Whati&*Yi)/CSUM(Whati).
.  COMPUTE uwi=CSUM(Whati&*((Xi-MeanXwi)&**2)).
.  COMPUTE qwi=CSUM(Whati&*((Yi-MeanYwi)&**2)).
.  COMPUTE pwi=CSUM(Whati&*(Xi-MeanXwi)&*(Yi-MeanYwi)).
.  COMPUTE bi=((L*qwi-uwi)+SQRT((uwi-L*qwi)**2+4*L*pwi**2))/(2*L*pwi).
.  COMPUTE a0i=MeanYwi-bi*MeanXwi.
.  COMPUTE diff=ABS(B0-bi).
. END LOOP IF (diff LT 1.0E-5).
. COMPUTE jack(i,1)=n*a0-(n-1)*a0i.
. COMPUTE jack(i,2)=n*b-(n-1)*bi.
END LOOP.
COMPUTE MeanJack=CSUM(jack)/n.
COMPUTE VarJack=(CSSQ(jack)-n*(MeanJack&**2))/(n-1).
COMPUTE SEJack=SQRT(VarJack/n).
COMPUTE Low={a0,b}-T95*SEJack.
COMPUTE Upp={a0,b}+T95*SEJack.
COMPUTE TValue=({a0,b}-{0,1})/SEJack.
COMPUTE PValue=2*(1-TCDF(ABS(TValue),n-2)).
PRINT {T({a0,b}),T(SEJack),T(Low),T(Upp),T(TValue),T(PValue)}
 /FORMAT='F8.5'
 /RLABEL='A=','B='
 /CLABEL='Coeff.','SE(coef)','Low95%CL','Upp95%CL','T Value','2-tail p'
 /TITLE='Weighted regression line & Jackknife estimators of SE & 95CI'.
PRINT Rw
 /FORMAT='F8.3'
 /TITLE='Weighted product-moment coefficient of correlation'.
* Preparing data to export *.
* 1) Standardized residuals *.
COMPUTE d=Y-(a0+b*X).
COMPUTE Xesti=X+(L*b*d)/(1+L*(b**2)).
COMPUTE Yesti=Y-(d/(1+L*(b**2))).
COMPUTE Wi=((Xesti+L*Yesti)/(1+L))&**(-2).
COMPUTE Sign=(Y-Yesti)/ABS(Y-Yesti).
COMPUTE Ri=Sign&*SQRT(Wi&*(X-Xesti)&**2+L*Wi&*(Y-Yesti)&**2).
COMPUTE SDr=SQRT(CSSQ(Ri)/(n-2)).
COMPUTE Residual=Ri/SDr.
COMPUTE Mean=(Xesti+L*Yesti)/(1+L).
COMPUTE vname={'Method1','Method2','Mean','Residual'}.
COMPUTE Export={X,Y,Mean,Residual}.
SAVE Export /OUTFILE=WorkingData /NAMES=vname.
* 2) Regression line *.
COMPUTE Param={a0,b}.
SAVE Param  /OUTFILE=Parameters.
END MATRIX.
DATASET ACTIVATE Parameters.
PRESERVE.
SET LOCALE=ENGLISH.
STRING #var1 #var2 (A12).
DO REPEAT A=#var1 #var2 /B=col1 col2.
- COMPUTE A = STRING(B,F8.2).
END REPEAT.
!LET !a=   !UNQUOTE(#var1).
!LET !b=   !UNQUOTE(#var2).
* 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">'
 /'	<setAxisInterval lowerMapping="exact" max="5.0" min="-5.0" role="y" upperMapping="exact"/>'
 /'	<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>'
 /'	<addReferenceLine styleOnly="false" y="4.0" ycategorical="false">'
 /'		<style stroke-dasharray="3px,2px" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addReferenceLine styleOnly="false" y="-4.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="'!b' * x +  '!a'" ycategorical="false">'
 /'		<style color="#000000" 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.
VAR LABEL Method1 'Method X'/Method2 'Method Y'.
GRAPH /TITLE='Regression plot'
 /SCATTERPLOT(BIVAR)=Method1 WITH Method2
 /TEMPLATE='C:\Temp\Regression.sgt'.
ECHO 'Red:   Equality line (y=x)'.
ECHO 'Black: Regression line (y=a+b*x)'.
EXE.
VAR LABEL Residual 'Weighted Standardized Residual' /Mean 'Weighted Average'.
GRAPH /TITLE='Residual plot'
 /SCATTERPLOT(BIVAR)=Mean WITH Residual
 /TEMPLATE='C:\Temp\Residual.sgt'.
RESTORE.
COMPUTE ID=$casenum.
FORMAT ID(F8).
COMPUTE AbsRes=ABS(Residual).
SORT CASES BY AbsRes(D).
ECHO '5 biggest residuals (absolute values)'.
LIST VARIABLES=ID Residual /CASES=FROM 1 TO 5.
ECHO 'Values over 4 are outliers'.
ECHO '-------------------------------------------'.
ECHO 'Linearity: RUNS test'.
SORT CASES BY Mean (A) .
NPAR TESTS /RUNS(MEDIAN)=Residual.
EXAMINE
  VARIABLES=Residual
  /PLOT BOXPLOT NPPLOT.
DATASET ACTIVATE OriginalData WINDOW=ASIS.
DATASET CLOSE WorkingData.
!ENDDEFINE.

* F) MACRO FOR WEIGHTED DEMING REGRESSION (PAIRS OF MEASUREMENTS, SQUARE ROOT WEIGHTING) *.

DEFINE DEMINGREG6 (x=!TOKENS(2) / y=!TOKENS(2)).
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.
COUNT nmiss = !x !y (SYSMIS) .
SELECT IF (nmiss=0). 
RANK VARIABLES=!x(A) /N INTO N1@ /PRINT=NO .
DO IF $casenum EQ 1.
. COMPUTE TValue = IDF.T(0.975,N1@-2) .
END IF.
TEMPORARY.
SET MXLOOPS=10000.
MATRIX.
PRINT /TITLE='****** WEIGHTED DEMING REGRESSION (WITH JACKKNIFE ESTIMATION OF SE&95%CI) ******'.
PRINT /TITLE='Pairs of measurements (direct estimation of Lambda)'.
PRINT/TITLE='Non constant measurement errors & square root weighting'.
PRINT /TITLE='(e.g. used when counting isotopes)'.
* Maximum number of iterations for estimation of weighted coefficients (usually only 3-4 are needed) *.
COMPUTE MxIter=25.
* Get all working data *.
GET T95 /VAR=TValue /MISSING=OMIT.
GET pair1 /VAR=!x /NAMES=vname1.
GET pair2 /VAR=!y /NAMES=vname2.
PRINT {vname1;vname2}
 /FORMAT='A8'
 /RLABELS='Method X','Method Y'
 /CLABELS='--1st---','--2nd---'
 /TITLE='Variables compared'.
* 1st pair of variables form X, 2nd pair form Y *.
COMPUTE X=RSUM(pair1)/2.
COMPUTE Y=RSUM(pair2)/2.
COMPUTE AverXY=(X+Y)/2.
COMPUTE n=NROW(pair1).
COMPUTE DiffX=pair1(:,1)-pair1(:,2).
COMPUTE DiffY=pair2(:,1)-pair2(:,2).
COMPUTE VarX=CSUM(DiffX&**2)/(2*n).
COMPUTE VarY=CSUM(DiffY&**2)/(2*n).
COMPUTE Fhatx2=CSUM((DiffX&/AverXY)&**2)/(2*n).
COMPUTE Fhaty2=CSUM((DiffY&/AverXY)&**2)/(2*n).
COMPUTE Lhat=Fhatx2/Fhaty2.
COMPUTE MeanX=CSUM(X)/n.
COMPUTE MeanY=CSUM(Y)/n.
COMPUTE u=CSSQ(X)-((CSUM(X))**2)/n.
COMPUTE q=CSSQ(Y)-((CSUM(Y))**2)/n.
COMPUTE p=CSUM(X&*Y)-CSUM(X)*CSUM(Y)/n.
PRINT {MeanX,VarX;MeanY,VarY}
 /FORMAT='F8.3'
 /RLABELS='Method X','Method Y'
 /CLABELS='Mean','SDa²'
 /TITLE='Statistics for both methods'.
* Unweighted estimation (starting point) *.
COMPUTE b=((Lhat*q-u)+SQRT((u-Lhat*q)**2+4*Lhat*p**2))/(2*Lhat*p).
COMPUTE a=MeanY.
COMPUTE a0=a-b*MeanX.
PRINT {100*SQRT(Fhatx2),100*SQRT(Fhaty2),Lhat,a0,b}
 /FORMAT='F8.3'
 /CLABEL='Fx(%)','Fy(%)','Lambda','A','B'
 /TITLE='Initial (unweighted) estimates'.
* Begin iteration process *.
LOOP jk=1 TO MxIter.
. COMPUTE B0=b.
. COMPUTE di=Y-(a0+b*X).
. COMPUTE Xhati=X+Lhat*b*di/(1+Lhat*(b**2)).
. COMPUTE Yhati=Y-di/(1+Lhat*(b**2)).
. COMPUTE Whati=((Xhati+Lhat*Yhati)/(1+Lhat))&**(-1).
. COMPUTE AverXiYi=(Xhati+Yhati)/2.
. COMPUTE Fhatx2=CSUM((DiffX&/AverXiYi)&**2)/(2*n).
. COMPUTE Fhaty2=CSUM((DiffY&/AverXiYi)&**2)/(2*n).
. COMPUTE Lhat=Fhatx2/Fhaty2.
. COMPUTE MeanXw=CSUM(Whati&*X)/CSUM(Whati).
. COMPUTE MeanYw=CSUM(Whati&*Y)/CSUM(Whati).
. COMPUTE uw=CSUM(Whati&*((X-MeanXw)&**2)).
. COMPUTE qw=CSUM(Whati&*((Y-MeanYw)&**2)).
. COMPUTE pw=CSUM(Whati&*(X-MeanXw)&*(Y-MeanYw)).
. COMPUTE b=((Lhat*qw-uw)+SQRT((uw-Lhat*qw)**2+4*Lhat*pw**2))/(2*Lhat*pw).
. COMPUTE a=MeanYw.
. COMPUTE a0=a-b*MeanXw.
. COMPUTE diff=ABS(B0-b).
END LOOP IF (diff LT 1.0E-5).
COMPUTE Rw=pw/SQRT(uw*qw).
* Report *.
PRINT {100*SQRT(Fhatx2),100*SQRT(Fhaty2),Lhat,a0,b}
 /FORMAT='F8.3'
 /CLABEL='Fx(%)','Fy(%)','Lambda','A','B'
 /TITLE='Weighted estimates (*)'.
PRINT jk
 /FORMAT='F8.0'
 /TITLE='(*) Estimation converged (diff. LT 1E-5) at iteration:'.
* Jackknife estimation of standard errors *.
* Compute empty matrix to store JK statistics (a0i&bi) *.
COMPUTE AllSamp={X,Y,DiffX,DiffY}.
COMPUTE jack=MAKE(n,2,0).
COMPUTE nj=n-1.
* Cycle thru all values *.
LOOP i=1 TO n.
. DO IF (i EQ 1).                          /* JK sample for 1st position *.
.  COMPUTE sample=AllSamp(2:n,:).
. ELSE IF (i GT 1) AND (i LT n). /* JK samples for 2nd to last-1 postion *.
.  COMPUTE sample={AllSamp(1:(i-1),:);
                   AllSamp((i+1):n,:)}.
. ELSE IF (i EQ n).                       /* JK sample for last position *.
.  COMPUTE sample=AllSamp(1:(n-1),:).
. END IF.
. COMPUTE Xi=sample(:,1).
. COMPUTE Yi=sample(:,2).
* Unweightedjackknife estimator *.
. COMPUTE AverXiYi=(Xi+Yi)/2.
. COMPUTE DiffXi=sample(:,3).
. COMPUTE DiffYi=sample(:,4).
. COMPUTE Fhatix2=CSUM((DiffXi&/AverXiYi)&**2)/(2*nj).
. COMPUTE Fhatiy2=CSUM((DiffYi&/AverXiYi)&**2)/(2*nj).
. COMPUTE Lhati=Fhatix2/Fhatiy2.
. COMPUTE MeanXi=CSUM(Xi)/nj.
. COMPUTE MeanYi=CSUM(Yi)/nj.
. COMPUTE ui=CSSQ(Xi)-((CSUM(Xi))**2)/nj.
. COMPUTE qi=CSSQ(Yi)-((CSUM(Yi))**2)/nj.
. COMPUTE pi=CSUM(Xi&*Yi)-CSUM(Xi)*CSUM(Yi)/nj.
. COMPUTE bi=((Lhati*qi-ui)+SQRT((ui-Lhati*qi)**2+4*Lhati*pi**2))/(2*Lhati*pi).
. COMPUTE a0i=MeanYi-bi*MeanXi.
* Begin iteration process for each jackknife weighted pseudovariate *.
. LOOP jk=1 TO MxIter.
.  COMPUTE B0i=bi.
.  COMPUTE di=Yi-(a0i+bi*Xi).
.  COMPUTE Xhati=Xi+Lhati*bi*di/(1+Lhati*(bi**2)).
.  COMPUTE Yhati=Yi-di/(1+Lhati*(bi**2)).
.  COMPUTE Whati=((Xhati+Lhati*Yhati)/(1+Lhati))&**(-1).
.  COMPUTE AverXiYi=(Xhati+Yhati)/2.
.  COMPUTE Fhatix2=CSUM((DiffXi&/AverXiYi)&**2)/(2*nj).
.  COMPUTE Fhatiy2=CSUM((DiffYi&/AverXiYi)&**2)/(2*nj).
.  COMPUTE Lhati=Fhatix2/Fhatiy2.
.  COMPUTE MeanXwi=CSUM(Whati&*Xi)/CSUM(Whati).
.  COMPUTE MeanYwi=CSUM(Whati&*Yi)/CSUM(Whati).
.  COMPUTE uwi=CSUM(Whati&*((Xi-MeanXwi)&**2)).
.  COMPUTE qwi=CSUM(Whati&*((Yi-MeanYwi)&**2)).
.  COMPUTE pwi=CSUM(Whati&*(Xi-MeanXwi)&*(Yi-MeanYwi)).
.  COMPUTE bi=((Lhati*qwi-uwi)+SQRT((uwi-Lhati*qwi)**2+4*Lhati*pwi**2))/(2*Lhati*pwi).
.  COMPUTE a0i=MeanYwi-bi*MeanXwi.
.  COMPUTE diffi=ABS(B0i-bi).
. END LOOP IF (diffi LT 1.0E-5).
. COMPUTE jack(i,1)=n*a0-(n-1)*a0i.
. COMPUTE jack(i,2)=n*b-(n-1)*bi.
END LOOP.
COMPUTE MeanJack=CSUM(jack)/n.
COMPUTE VarJack=(CSSQ(jack)-n*(MeanJack&**2))/(n-1).
COMPUTE SEJack=SQRT(VarJack/n).
COMPUTE Low={a0,b}-T95*SEJack.
COMPUTE Upp={a0,b}+T95*SEJack.
COMPUTE TValue=({a0,b}-{0,1})/SEJack.
COMPUTE PValue=2*(1-TCDF(ABS(TValue),n-2)).
PRINT {T({a0,b}),T(SEJack),T(Low),T(Upp),T(TValue),T(PValue)}
 /FORMAT='F8.4'
 /RLABEL='A=','B='
 /CLABEL='Coeff.','SE(coef)','Low95%CL','Upp95%CL','T Value','2-tail p'
 /TITLE='Regression line & Jackknife estimators of SE & 95CI'.
PRINT Rw
 /FORMAT='F8.3'
 /TITLE='Weighted product-moment coefficient of correlation'.
* Preparing data to export *.
* 1) Standardized residuals *.
COMPUTE d=Y-(a0+b*X).
COMPUTE Xesti=X+Lhat*b*d/(1+Lhat*(b**2)).
COMPUTE Yesti=Y-d/(1+Lhat*(b**2)).
COMPUTE Wi=((Xesti+Lhat*Yesti)/(1+Lhat))&**(-2).
COMPUTE Sign=(Y-Yesti)/ABS(Y-Yesti).
COMPUTE Ri=Sign&*SQRT(Wi&*(X-Xesti)&**2+Lhat*Wi&*(Y-Yesti)&**2).
COMPUTE SDr=SQRT(CSSQ(Ri)/(n-2)).
COMPUTE Residual=Ri/SDr.
COMPUTE Mean=(Xesti+Lhat*Yesti)/(1+Lhat).
COMPUTE vname={'Method1','Method2','Mean','Residual'}.
COMPUTE Export={X,Y,Mean,Residual}.
SAVE Export /OUTFILE=WorkingData /NAMES=vname.
* 2) Regression line *.
COMPUTE Param={a0,b}.
SAVE Param  /OUTFILE=Parameters.
END MATRIX.
DATASET ACTIVATE Parameters.
PRESERVE.
SET LOCALE=ENGLISH.
STRING #var1 #var2 (A12).
DO REPEAT A=#var1 #var2 /B=col1 col2.
- COMPUTE A = STRING(B,F8.2).
END REPEAT.
!LET !a=   !UNQUOTE(#var1).
!LET !b=   !UNQUOTE(#var2).
* 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">'
 /'	<setAxisInterval lowerMapping="exact" max="5.0" min="-5.0" role="y" upperMapping="exact"/>'
 /'	<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>'
 /'	<addReferenceLine styleOnly="false" y="4.0" ycategorical="false">'
 /'		<style stroke-dasharray="3px,2px" visible="true"/>'
 /'	</addReferenceLine>'
 /'	<addReferenceLine styleOnly="false" y="-4.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="'!b' * x +  '!a'" ycategorical="false">'
 /'		<style color="#000000" 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.
VAR LABEL Method1 'Method X'/Method2 'Method Y'.
GRAPH /TITLE='Regression plot'
 /SCATTERPLOT(BIVAR)=Method1 WITH Method2
 /TEMPLATE='C:\Temp\Regression.sgt'.
ECHO 'Red:   Equality line (y=x)'.
ECHO 'Black: Regression line (y=a+b*x)'.
EXE.
VAR LABEL Residual 'Weighted Standardized Residual' /Mean 'Weighted Average'.
GRAPH /TITLE='Residual plot'
 /SCATTERPLOT(BIVAR)=Mean WITH Residual
 /TEMPLATE='C:\Temp\Residual.sgt'.
RESTORE.
COMPUTE ID=$casenum.
FORMAT ID(F8).
COMPUTE AbsRes=ABS(Residual).
SORT CASES BY AbsRes(D).
ECHO '5 biggest residuals (absolute values)'.
LIST VARIABLES=ID Residual /CASES=FROM 1 TO 5.
ECHO 'Values over 4 are outliers'.
ECHO '-------------------------------------------'.
ECHO 'Linearity: RUNS test'.
SORT CASES BY Mean (A) .
NPAR TESTS /RUNS(MEDIAN)=Residual.
EXAMINE
  VARIABLES=Residual
  /PLOT BOXPLOT NPPLOT.
DATASET ACTIVATE OriginalData WINDOW=ASIS.
DATASET CLOSE WorkingData.
!ENDDEFINE.

DISPLAY MACROS.

***********************************************************
***************** EXAMPLES OF MACRO CALLS *****************
***********************************************************

* Sample dataset (from Bland&Altman Lancet(1986) 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 calls (some examples were tested with Medcalc & Method Validator) *

***************
* FIRST MACRO *
***************

* 2 arguments (x & y) with the names of the methods; there is an optional argument
 (Lambda=error variances ratio), with a default value of 1 *.

DEMINGREG1 x=wright1 y=mini1.

* If we have external information that tells us that the Wright Meter has better precision,
  (for instance: it has around 60% the measurement error of the Mini Meter), we can use that
  as an estimate of Lambda *.

DEMINGREG1 x=wright1 y=mini1 lambda=0.6.

* Here we can safely assume a variance ratio = 1 (two consecutive measurements of same method) *.

DEMINGREG1 x=wright1 y=wright2.

* Using averaged values and external estimation for Lambda (it would be better to use the 3rd macro!) *.

COMPUTE wright=(wright1+wright2)/2.
COMPUTE mini=(mini1+mini2)/2.
DEMINGREG1 x=wright y=mini lambda=0.6.

****************
* SECOND MACRO *
****************

* 2 arguments (x & y) with the names of the methods; there is an optional argument
 (Lambda=error variances ratio), with a default value of 1 *.

DEMINGREG2 x=wright1 y=mini1.

* If we have external information that tells us that the Wright Meter has better precision,
  (for instance: it has around 60% the measurement error of the Mini Meter), we can use that
  as an estimate of Lambda *.

DEMINGREG2 x=wright1 y=mini1 lambda=0.6.

***************
* THIRD MACRO *
***************

* Pairs of measurements for both methods are needed; since Lambda is estimated from the repeated data,
  it is not supplied as an extra argument *.

DEMINGREG3 x=wright1 wright2 y=mini1 mini2 .

***************
*  4TH MACRO  *
***************

* Pairs of measurements for both methods are needed; since Lambda is estimated from the repeated data,
  it is not supplied as an extra argument *.

DEMINGREG4 x=wright1 wright2 y=mini1 mini2 .

***************
*  5TH MACRO  *
***************

* Same examples as the second macro. The only difference is the weighting used (square-root instead of proportional) *.

DEMINGREG5 x=wright1 y=mini1 lambda=0.6.

***************
*  6TH MACRO  *
***************

* Same examples as the 4th macro. The only difference is the weighting used (square-root instead of proportional) *.

DEMINGREG6 x=wright1 wright2 y=mini1 mini2 .