* Survival Analysis Example.

* Posted to SPSSX-L on 2004/05/13 by Marta Garcia-Granero.

* Dataset slightly modified (some leukocytes data changed)
  from Selvin S (1996) "Statistical analysis of epidemiological data"
  Oxford University Press *
* Survival times of 33 patients with acute mieloid leukhaemia *.

DATA LIST FREE /time(F8.1) status auer_r leuko (3 F8.0).
BEGIN DATA
65.0    1 0 2300 156.0 1 0 750 100.0 1 0 4300
134.0   1 0 2600 16.0   1 0 16000 108.0 1 0 10500
121.0   1 0 10000 4.0   1 0 17000 69.0 1 0 5400
143.0   1 0 7000 56.0   1 0 9400 26.0   1 0 32000
22.0    1 0 35000 1.0   1 0 100000 1.0 1 0 130000
5.0     1 0 52000 65.0  1 0 100000 65.0 1 1 3000
17.0    1 1 4000 7.0    1 1 1500 16.0   1 1 9000
22.0    1 1 5300 3.0    1 1 10000 4.0   1 1 19000
2.0     1 1 27000 3.0   1 1 28000 8.0   1 1 31000
4.0     1 1 26000 3.0   1 1 21000 30.0 1 1 79000
4.0     1 1 100000 43.0 1 1 100000 56.0 1 1 4400
END DATA.

VAR LABELS time "Follow-up time (weeks)"
 /status "Censoring status"
 /auer_r "Auer Rods"
 /leuko "Leucocytes (counts/mmі)".

VALUE LABELS status 0 "Alive" 1 "Dead"
 /auer_r 0 "Presents" 1 "Absents".

* Kaplan-Meier analysis *.

KM time  BY auer_r  /STATUS=status(1)
  /PRINT TABLE MEAN
  /PLOT SURVIVAL
  /TEST LOGRANK
  /COMPARE OVERALL POOLED .

RANK
  VARIABLES=leuko (A) /NTILES (3) /PRINT=NO
  /TIES=CONDENSE .

VAR LABELS nleuko "Leukocytes in 3 cat.".
VALUE LABELS nleuko 1 "Lower 3rd (<7.000 mmі)"
                    2 "Middle 3rd (7.000-27.000 mmі)"
                    3 "Upper 3rd (>27.000 mmі)".

KM time  BY nleuko  /STATUS=status(1)
  /PRINT TABLE MEAN
  /PLOT SURVIVAL
  /TEST LOGRANK
  /COMPARE OVERALL POOLED .

* Pair-wise comparisons *.

KM  time  BY nleuko /STATUS=status(1)
  /PRINT NONE
  /TEST LOGRANK
  /COMPARE PAIRWISE POOLED .

* Stratified analysis *.

KM  time  BY  auer_r /STRATA=nleuko /STATUS=status(1)
  /PRINT TABLE MEAN
  /PLOT SURVIVAL
  /TEST LOGRANK
  /COMPARE OVERALL POOLED .

KM  time  BY  nleuko/STRATA=auer_r  /STATUS=status(1)
  /PRINT TABLE MEAN
  /PLOT SURVIVAL
  /TEST LOGRANK
  /COMPARE OVERALL POOLED .

* Comments: from the graphs, interaction is suspected *.

* Cox univariate *.

COXREG
  time /STATUS=status(1)
  /CONTRAST (auer_r)=Indicator(1)
  /METHOD=ENTER auer_r
  /PRINT=CI(95)
  /CRITERIA=PIN(.05) POUT(.10) ITERATE(20) .

COXREG
  time /STATUS=status(1)
  /METHOD=ENTER leuko
  /PRINT=CI(95)
  /CRITERIA=PIN(.05) POUT(.10) ITERATE(20) .

* Comment: the increase in HR associated to 1 leukocyte count change is small,
  it's better to work with 1000 counts increase *.

COMPUTE thouleuk = leuko / 1000 .
EXECUTE .

COXREG
  time /STATUS=status(1)
  /METHOD=ENTER thouleuk
  /PRINT=CI(95)
  /CRITERIA=PIN(.05) POUT(.10) ITERATE(20) .

* Cox main effects, BSTEP(LR) with PIN(1)&POUT(1) is used to get LR adjusted tests
  (better - more sensitive - than Wald tests )*.

COXREG
  time  /STATUS=status(1)
  /CONTRAST (auer_r)=Indicator(1)
  /METHOD=BSTEP(LR) auer_r thouleuk
  /PRINT=CI(95)
  /CRITERIA=PIN(1) POUT(1) ITERATE(20) .

* Model with interaction *.

COXREG
  time  /STATUS=status(1)
  /CONTRAST (auer_r)=Indicator(1)
  /METHOD=BSTEP(LR) auer_r thouleuk auer_r*thouleuk
  /PRINT=CI(95)
  /CRITERIA=PIN(1) POUT(1) ITERATE(20) .

* Graphical presentation of the interaction: *.

COMPUTE hr=EXP(1.924*auer_r + 0.022*thouleuk -0.022*auer_r*thouleuk).
* (coefficients read from the output ) *.
EXEC.
GRAPH /SCATTERPLOT(BIVAR)=thouleuk WITH hr BY auer_r.

* Edit the graph to add fit lines to both subgroups *.

* Extras: verification of the proportional hazards assumption *.

* (a) LML (Log Minus Log) graph (qualitative predictors or categorized quantitative) *.

* If we try to get the LML graph this way: *.

COXREG
  time /STATUS=status(1)
  /PATTERN BY auer_r
  /CONTRAST (auer_r)=Indicator(1)
  /METHOD=ENTER auer_r
  /PLOT LML.

* We'll ALWAYS get an image of absolute proportionality because the ESTIMATED LML are plotted,
  not the OBSERVED.

* The correct way is: *.
COXREG
  time  /STATUS=status(1)  /STRATA=auer_r
  /PLOT LML.

* This will give the observed LML for both groups, not the predicted by the regression model *.

* (b) For quantitative predictors, the way of testing the proportionality assumption is Cox model
  with time-dep covariate. Check it in 2 ways: using time or LN(time) for the T_COV_ expression: *.

TIME PROGRAM.
COMPUTE T_COV_ = T_*thouleuk .
COXREG
  time  /STATUS=status(1)
  /METHOD=BSTEP(LR) T_COV_ thouleuk
  /PRINT=CI(95)
  /CRITERIA=PIN(1) POUT(1) ITERATE(20) .

* LR test p-value is=0.781 (ns), the proportionality assumption is not violated for leukocytes *.