*****************************************************************
*   KALBFLEISCH-PRENTICE 95% CI INTERVALS FOR SURVIVAL CURVES   *
*****************************************************************
* Kalbfleisch JD and Prentice RL(1980). The Statistical         *
* Analysis of Failure Time Data. New York: John Wiley & Sons    *
*****************************************************************
* (c) Marta Garcia-Granero (june 2005 - original code) modified *
* (turnt into macro) september 2008.  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     *
         ***********************************************

* Unlike Greenwood's, the Kalbfleisch–Prentice approach has the attractive
  property that the upper and lower bounds are always between 0 and 1.

* MACRO definition *.
DEFINE KPSURVIVAL(!POS=!TOKENS(1)/!POS=!CHAREND('(')/!POS=!CHAREND(')')).
* Creating the chart template *.
PRESERVE.
SET RESULTS=NONE MESSAGES=NONE ERRORS=NONE.
HOST COMMAND=['IF not exist C:\Temp MD C:\Temp'].
HOST COMMAND=['DEL C:\Temp\*.sgt'].
RESTORE.
DO IF $casenum EQ 1.
WRITE OUTFILE='C:\Temp\KPTemplate.sgt'
/'<?xml version="1.0" encoding="UTF-8" standalone="no"?>'
/'<template SPSS-Version="2.2" selectPath="12 47 "'
/' 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">'
/'	<setAxisMargin categorical="false" lowerMargin="0%" role="x" upperMargin="5%"/>'
/'	<setAxisMargin categorical="false" lowerMargin="0%" role="y" upperMargin="5%"/>'
/'	<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="transparent" number="1" visible="true"/>'
/'	</addFrame>'
/'	<addFrame count="1" styleOnly="true" type="title">'
/'		<style color="transparent" color2="transparent" visible="true"/>'
/'	</addFrame>'
/'</template>'.
END IF.
EXECUTE.
DATASET NAME OriginalData.
* Copy data to leave original unmodified *.
DATASET COPY WorkingData1 WINDOW=MINIMIZED.
DATASET DECLARE AuxData WINDOW=HIDDEN.
DATASET ACTIVATE WorkingData1.
COMPUTE status=(!2 EQ !3).
RENAME VARIABLES (!1 = time).
* Get rid of every other variable *.
MATCH FILES FILE=*
 /KEEP=time status.
* Save sorted file to match CI limits later *.
SORT CASES BY time(A).
DATASET COPY WorkingData2  WINDOW=MINIMIZED.
* Saving KM survival estimates *.
PRESERVE.
SET ERRORS=NONE RESULTS=NONE.
KM time /STATUS=status(1)
  /PRINT NONE
  /PLOT SURVIVAL
  /SAVE SURVIVAL  .
* Total sample size *.
RANK VARIABLES=time(A) /N INTO totaln/PRINT=NO.
* Preparing data to compute pj and Std Error of Log(survival) *.
AGGREGATE /OUTFILE = *  MODE=REPLACE
  /BREAK = time
  /aj = SUM(status)
  /sj = FIRST(SUR_1)
  /totaln = FIRST(totaln)
  /N = n.
* Auxiliary file #1 *.
MATRIX.
COMPUTE time=0.
COMPUTE aj=0.
COMPUTE sj=1.
COMPUTE totaln=0.
COMPUTE n=0.
COMPUTE vnames={'time','aj','sj','totaln','n'}.
COMPUTE newdata={time,aj,sj,totaln,n}.
SAVE newdata /OUTFILE='AuxData'/NAMES=vnames.
END MATRIX.
* Adding time 0 data (from auxiliary file)*.
ADD FILES /FILE=*
 /FILE='AuxData'
 /IN=source.
IF source EQ 1 totaln=LAG(totaln).
SORT CASES BY time(A).
* Computing auxiliary statistics *.
COMPUTE rj=totaln.
DO IF $casenum GT 1.
- COMPUTE n = n+LAG(n) .
- COMPUTE rj= totaln-LAG(n).
END IF.
EXECUTE ./* Don't eliminate it (necessary) *.
* Getting rid of censored data *.
SELECT IF NOT MISSING(sj).
* More statistics and reordering dataset *.
COMPUTE j= $casenum-1.
COMPUTE pj = (rj-aj)/rj .
MATCH FILES FILE = *
 /KEEP = j time aj rj pj sj .
* 95% CI limits *.
COMPUTE sum = ((1-pj)/(pj*rj)).
DO IF $casenum GT 1.
- COMPUTE sum=sum+LAG(sum).
- COMPUTE sdlog = SQRT(sum/((LN(sj))**2)).
* 95% Limits using Log(-Log(Sj)) & back transform
* Replace 1.959964 by 2.575829 if 99%CI are wanted *.
- COMPUTE lower = EXP(-EXP(LN(-LN(sj))+1.959964*sdlog)).
- COMPUTE upper = EXP(-EXP(LN(-LN(sj))-1.959964*sdlog)).
END IF.
RESTORE.
* Report *.
FORMAT j TO rj (F8.0) sj pj lower upper (F8.3).
SUMMARIZE
 /TABLES=time sj lower upper
 /FORMAT=LIST NOCASENUM NOTOTAL
 /TITLE='Kaplan-Meier Estimates and Kalbfleisch-Prentice 95% CI'
 /CELLS=NONE.

* Adding survival & CI limits to original dataset *.
MATCH FILES FILE=*
 /KEEP = time sj lower upper.
MATCH FILES
 /FILE='WorkingData2'
 /TABLE=*
 /BY time.
* Adding time=0 survival data (auxiliary file #2) *.
PRESERVE.
SET ERRORS=NONE RESULTS=NONE.
MATRIX.
COMPUTE time=0.
COMPUTE sj=1.
COMPUTE lower=1.
COMPUTE upper=1.
COMPUTE vnames={'time','sj','lower','upper'}.
COMPUTE newdata={time,sj,lower,upper}.
SAVE newdata /OUTFILE='AuxData'/NAMES=vnames.
END MATRIX.
ADD FILES /FILE=*
 /FILE='AuxData'.
SORT CASES BY time(a).
* Filling data for censored cases after the first event *.
DO IF MISSING(sj).
- COMPUTE sj=LAG(sj).
- COMPUTE lower=LAG(lower).
- COMPUTE upper=LAG(upper).
END IF.
RESTORE.
* Survival plot (fully formatted, thanks to ViAnn Beadle) *.
FORMAT sj lower upper (F8.1).
!LET !gdataset=!CONCAT('"','graphdataset','"').
!LET !time=!CONCAT('"','time','"').
!LET !sj=!CONCAT('"','sj','"').
!LET !lower=!CONCAT('"','lower','"').
!LET !upper=!CONCAT('"','upper','"').
!LET !ylabel=!CONCAT('"','Cumulative Survival','"').
!LET !xlabel=!CONCAT('"','Follow-up Time','"').
!LET !title=!CONCAT('"','KM estimate with Kalbfleish-Prentice 95%CI','"').
!LET !px=!CONCAT('"','2px','"').
GGRAPH
  /GRAPHDATASET NAME=!gdataset VARIABLES=time sj lower upper
  MISSING=LISTWISE REPORTMISSING=NO
  /GRAPHSPEC SOURCE=INLINE TEMPLATE="C:\Temp\KPTemplate.sgt".
BEGIN GPL
  SOURCE: s=userSource(id(!gdataset))
  DATA: time=col(source(s), name(!time))
  DATA: sj=col(source(s), name(!sj))
  DATA: lower=col(source(s), name(!lower))
  DATA: upper=col(source(s), name(!upper))
  SCALE: linear(dim(2), include(0)) 
  GUIDE: axis(dim(1), label(!xlabel))
  GUIDE: axis(dim(2), label(!ylabel), delta(0.1))
  GUIDE: legend(aesthetic(aesthetic.shape.interior), null())
  GUIDE: text.title(label(!title))
  ELEMENT: line(position(smooth.step.left(time*lower)), missing.wings(), 
           shape.interior(shape.dash), color(color.black))
  ELEMENT: line(position(smooth.step.left(time*upper)), missing.wings(), 
           shape.interior(shape.dash), color(color.black))
  ELEMENT: line(position(smooth.step.left(time*sj)), missing.wings(), 
           shape.interior(shape.solid), size(size.!px))
END GPL.

* Tidying up *.
DATASET ACTIVATE OriginalData WINDOW=ASIS.
DATASET CLOSE AuxData.
DATASET CLOSE WorkingData1.
DATASET CLOSE WorkingData2.
!ENDDEFINE.


*****************************************************************
*                     EXAMPLE DATASET                           *
* Table 9.1 Biostatistical Methods in Epidemiology (pp 177-178) *
* Stephen C Newman (2001); JOHN WILEY&SONS                      *
*****************************************************************.
DATA LIST FREE/FuTime Event stage reclevel (4 F8.0).
BEGIN DATA
50 0 1 1 51 1 1 1 51 0 1 1 53 0 1 1 53 0 1 1 54 0 1 1
54 0 1 1 55 0 1 1 56 1 1 1 56 0 1 1 57 0 1 1 60 0 1 1
10 1 1 2 34 1 1 2 34 0 1 2 47 1 1 2 47 1 1 2 49 0 1 2
49 0 1 2 50 0 1 2 50 0 1 2 50 0 1 2 50 0 1 2 50 0 1 2
50 0 1 2 50 0 1 2 51 0 1 2 51 0 1 2 51 0 1 2 51 0 1 2
51 0 1 2 51 0 1 2 52 0 1 2 52 0 1 2 52 0 1 2 52 0 1 2
52 0 1 2 53 0 1 2 53 0 1 2 53 0 1 2 53 0 1 2 53 0 1 2
53 0 1 2 54 0 1 2 54 0 1 2 54 0 1 2 54 0 1 2 54 0 1 2
55 0 1 2 55 0 1 2 56 0 1 2 56 0 1 2 57 0 1 2 57 0 1 2
57 0 1 2 57 0 1 2 57 0 1 2 58 0 1 2 58 0 1 2 58 0 1 2
58 0 1 2 58 0 1 2 59 0 1 2 59 0 1 2 59 0 1 2 59 0 1 2
60 0 1 2 60 0 1 2 60 0 1 2  4 0 2 1  9 1 2 1 13 1 2 1
21 1 2 1 29 1 2 1 29 1 2 1 40 1 2 1 46 1 2 1 49 0 2 1
49 0 2 1 52 0 2 1 52 0 2 1 53 0 2 1 54 0 2 1 55 0 2 1
55 0 2 1 56 0 2 1 57 1 2 1 57 0 2 1 58 0 2 1 58 0 2 1
59 0 2 1 60 0 2 1 11 1 2 2 16 1 2 2 21 1 2 2 23 1 2 2
23 1 2 2 24 1 2 2 33 1 2 2 33 1 2 2 36 1 2 2 36 1 2 2
36 0 2 2 37 1 2 2 45 1 2 2 46 1 2 2 49 0 2 2 49 0 2 2
50 0 2 2 50 0 2 2 50 0 2 2 50 0 2 2 50 0 2 2 50 0 2 2
51 0 2 2 51 0 2 2 51 0 2 2 51 0 2 2 52 0 2 2 52 0 2 2
52 0 2 2 52 0 2 2 52 0 2 2 53 0 2 2 53 0 2 2 53 0 2 2
53 0 2 2 53 0 2 2 54 0 2 2 54 0 2 2 54 0 2 2 54 0 2 2
55 0 2 2 55 0 2 2 55 0 2 2 55 0 2 2 56 0 2 2 56 0 2 2
56 0 2 2 56 0 2 2 56 0 2 2 56 0 2 2 57 0 2 2 57 0 2 2
57 0 2 2 57 0 2 2 58 1 2 2 58 1 2 2 58 0 2 2 58 0 2 2
58 0 2 2 58 0 2 2 58 0 2 2 58 0 2 2 58 0 2 2 58 0 2 2
59 0 2 2 59 0 2 2 59 0 2 2 59 0 2 2 59 0 2 2 60 0 2 2
60 0 2 2 60 0 2 2 60 0 2 2 60 0 2 2 60 0 2 2  9 1 3 1
12 1 3 1 14 1 3 1 15 1 3 1 15 0 3 1 17 1 3 1 21 1 3 1
22 1 3 1 23 1 3 1 23 1 3 1 31 1 3 1 34 1 3 1 35 1 3 1
53 0 3 1 60 0 3 1  7 0 3 2  9 1 3 2 17 1 3 2 21 0 3 2
22 1 3 2 22 1 3 2 34 1 3 2 34 1 3 2 41 1 3 2 49 0 3 2
52 0 3 2 55 1 3 2 56 0 3 2 58 0 3 2 58 0 3 2 59 0 3 2
59 0 3 2
END DATA.
VAR LABEL stage' Tumour Stage'/ reclevel 'Receptor level'
/FuTime 'Survival time (months)' /Event 'Censoring Status'.
VALUE label stage 1'I' 2'II' 3'III'
 /reclevel 1'Low' 2'High' /Event 0'Censored' 1'Event'.

KPSURVIVAL FuTime Event(1).