# 21  Case Study in Cox Regression

## 21.1 Choosing the Number of Parameters and Fitting the Model

• Clinical trial of estrogen for prostate cancer
• Response is time to death, all causes
• Base analysis on Cox proportional hazards model
• $$S(t | X)$$ = probability of surviving at least to time $$t$$ given set of predictor values $$X$$
• $$S(t | X) = S_{0}(t)^{\exp(X\beta)}$$
• Censor time to death at time of last follow-up for patients still alive at end of study (treat survival time for pt.
censored at 24m as 24m+)
• Use simple, partial approaches to data reduction
• Use transcan for single imputation
• Again combine last 2 categories for ekg,pf
• See if we can use a full additive model (4 knots for continuous $$X$$)

A

Predictor Name d.f. Original Levels
Dose of estrogen rx 3 placebo, 0.2, 1.0, 5.0 mg estrogen
Age in years age 3
Weight index: wt(kg)-ht(cm)+200 wt 3
Performance rating pf 2 normal, in bed <50% of time, in bed >50%, in bed always
History of cardiovascular disease hx 1 present/absent
Systolic blood pressure/10 sbp 3
Diastolic blood pressure/10 dbp 3
Electrocardiogram code ekg 5 normal, benign, rhythm disturb., block, strain, old myocardial infarction, new MI
Serum hemoglobin (g/100ml) hg 3
Tumor size (cm$$^2$$) sz 3
Stage/histologic grade combination sg 3
Serum prostatic acid phosphatase ap 3
Bone metastasis bm 1 present/absent
• Total of 36 candidate d.f.
• Impute missings and estimate shrinkage B
Code
require(rms)
options(prType='html')    # for print, summary, anova, validate
getHdata(prostate)
levels(prostate$ekg)[levels(prostate$ekg) %in%
c('old MI','recent MI')] <- 'MI'
# combines last 2 levels and uses a new name, MI

prostate$pf.coded <- as.integer(prostate$pf)
# save original pf, re-code to 1-4
levels(prostate$pf) <- c(levels(prostate$pf)[1:3],
levels(prostate\$pf)[3])
# combine last 2 levels

w <- transcan(~ sz + sg + ap + sbp + dbp + age +
wt + hg + ekg + pf + bm + hx,
imputed=TRUE, data=prostate, pl=FALSE, pr=FALSE)

attach(prostate)
sz  <- impute(w, sz, data=prostate)
sg  <- impute(w, sg, data=prostate)
age <- impute(w, age,data=prostate)
wt  <- impute(w, wt, data=prostate)
ekg <- impute(w, ekg,data=prostate)

units(dtime) <- 'Month'
S <- Surv(dtime, status!='alive')

f <- cph(S ~ rx + rcs(age,4) + rcs(wt,4) + pf + hx +
rcs(sbp,4) + rcs(dbp,4) + ekg + rcs(hg,4) +
rcs(sg,4) + rcs(sz,4) + rcs(log(ap),4) + bm)
print(f, coefs=FALSE)

Cox Proportional Hazards Model

cph(formula = S ~ rx + rcs(age, 4) + rcs(wt, 4) + pf + hx + rcs(sbp,
4) + rcs(dbp, 4) + ekg + rcs(hg, 4) + rcs(sg, 4) + rcs(sz,
4) + rcs(log(ap), 4) + bm)

Model Tests Discrimination
Indexes
Obs 502 LR χ2 136.22 R2 0.238
Events 354 d.f. 36 R236,502 0.181
Center -2.9933 Pr(>χ2) 0.0000 R236,354 0.247
Score χ2 143.62 Dxy -0.333
Pr(>χ2) 0.0000
• Global LR $$\chi^2$$ is 135 and very significant $$\rightarrow$$ modeling warranted C
• AIC on $$\chi^2$$ scale = $$136.2 - 2 \times 36 = 64.2$$
• Rough shrinkage: 0.74 ($$\frac{136.2 - 36}{136.2}$$)
• Informal data reduction (increase for ap) D
Data reduction strategy
Variables Reductions d.f. Saved
wt Assume variable not important enough for 4 knots; use 3 1
pf Assume linearity 1
hx,ekg Make new 0,1,2 variable and assume linearity: 2=hx and ekg not normal or benign, 1=either, 0=none 5
sbp,dbp Combine into mean arterial bp and use 3 knots: map=$$\frac{2}{3}$$ dbp $$+ \frac{1}{3}$$ sbp 4
sg Use 3 knots 1
sz Use 3 knots 1
ap Look at shape of effect of ap in detail, and take log before expanding as spline to achieve stability: add 1 knot -1
Code
heart <- hx + ekg %nin% c('normal','benign')
label(heart) <- 'Heart Disease Code'
map   <- (2*dbp + sbp)/3
label(map) <- 'Mean Arterial Pressure/10'

f <- cph(S ~ rx + rcs(age,4) + rcs(wt,3) + pf.coded +
heart + rcs(map,3) + rcs(hg,4) +
rcs(sg,3) + rcs(sz,3) + rcs(log(ap),5) + bm,
x=TRUE, y=TRUE, surv=TRUE, time.inc=5*12)
print(f, coefs=FALSE)

Cox Proportional Hazards Model

cph(formula = S ~ rx + rcs(age, 4) + rcs(wt, 3) + pf.coded +
heart + rcs(map, 3) + rcs(hg, 4) + rcs(sg, 3) + rcs(sz, 3) +
rcs(log(ap), 5) + bm, x = TRUE, y = TRUE, surv = TRUE, time.inc = 5 *
12)

Model Tests Discrimination
Indexes
Obs 502 LR χ2 118.37 R2 0.210
Events 354 d.f. 24 R224,502 0.171
Center -2.4307 Pr(>χ2) 0.0000 R224,354 0.234
Score χ2 125.58 Dxy -0.321
Pr(>χ2) 0.0000
Code
# x, y for anova LR, predict, validate, calibrate;
# surv, time.inc for calibrate
anova(f, test='LR')
Likelihood Ratio Statistics for S
χ2 d.f. P
rx 8.46 3 0.0373
age 11.98 3 0.0074
Nonlinear 8.31 2 0.0157
wt 7.80 2 0.0202
Nonlinear 2.45 1 0.1176
pf.coded 3.49 1 0.0616
heart 23.82 1 <0.0001
map 0.05 2 0.9777
Nonlinear 0.04 1 0.8339
hg 11.44 3 0.0095
Nonlinear 7.53 2 0.0231
sg 1.65 2 0.4386
Nonlinear 0.05 1 0.8295
sz 11.81 2 0.0027
Nonlinear 0.06 1 0.7998
ap 6.22 4 0.1836
Nonlinear 5.82 3 0.1208
bm 0.03 1 0.8674
TOTAL NONLINEAR 22.73 11 0.0193
TOTAL 118.37 24 <0.0001
• Savings of 12 d.f. E
• AIC=70, shrinkage 0.80

## 21.2 Checking Proportional Hazards

• This is our tentative model
• Examine distributional assumptions using scaled Schoenfeld residuals
• Complication arising from predictors using multiple d.f.
• Transform to 1 d.f. empirically using $$X\hat{\beta}$$
• cox.zph does this automatically
• Following analysis approx. since internal coefficients estimated
Code
z <- predict(f, type='terms')
# required x=T above to store design matrix
f.short <- cph(S ~ z, x=TRUE, y=TRUE)
# store raw x, y so can get residuals
• Fit f.short has same LR $$\chi^2$$ of 118 as the fit f, F but with falsely low d.f.
• All $$\beta=1$$
Code
require(survival)   # or use survival::cox.zph(...)
phtest <- cox.zph(f, transform='identity')
phtest
                   chisq df    p
rx              4.07e+00  3 0.25
rcs(age, 4)     4.27e+00  3 0.23
rcs(wt, 3)      2.22e-01  2 0.89
pf.coded        5.34e-02  1 0.82
heart           4.95e-01  1 0.48
rcs(map, 3)     3.20e+00  2 0.20
rcs(hg, 4)      5.26e+00  3 0.15
rcs(sg, 3)      1.01e+00  2 0.60
rcs(sz, 3)      3.07e-01  2 0.86
rcs(log(ap), 5) 3.59e+00  4 0.47
bm              2.11e-06  1 1.00
GLOBAL          2.30e+01 24 0.52
Code
plot(phtest[1])  # plot only the first variable
• None of the effects significantly change over time
• Global test of PH $$P=0.52$$

## 21.3 Testing Interactions

• Will ignore non-PH for dose even though it makes sense
• More accurate predictions could be obtained using stratification or time dep. cov.
• Test all interactions with dose
Reduce to 1 d.f. as before
Code
z.dose <- z[,"rx"]  # same as saying z[,1] - get first column
z.other <- z[,-1]   # all but the first column of z
f.ia <- cph(S ~ z.dose * z.other, x=TRUE, y=TRUE)
anova(f.ia, test='LR')
Likelihood Ratio Statistics for S
χ2 d.f. P
z.dose (Factor+Higher Order Factors) 20.65 11 0.0371
All Interactions 11.89 10 0.2926
z.other (Factor+Higher Order Factors) 121.19 20 <0.0001
All Interactions 11.89 10 0.2926
z.dose × z.other (Factor+Higher Order Factors) 11.89 10 0.2926
TOTAL 130.26 21 <0.0001

## 21.4 Describing Predictor Effects

G

• Plot relationship between each predictor and $$\log \lambda$$
Code
ggplot(Predict(f), sepdiscrete='vertical', nlevels=4,
vnames='names')

## 21.5 Validating the Model

H

• Validate for $$D_{xy}$$ and slope shrinkage
Code
set.seed(1)  # so can reproduce results
v <- validate(f, B=300)
v
Index Original
Sample
Training
Sample
Test
Sample
Optimism Corrected
Index
Successful
Resamples
Dxy -0.3208 -0.3494 -0.2953 -0.0541 -0.2667 300
R2 0.2101 0.2481 0.1756 0.0724 0.1377 300
Slope 1 1 0.7863 0.2137 0.7863 300
D 0.0292 0.0354 0.0239 0.0116 0.0176 300
U -5e-04 -5e-04 0.0024 -0.0029 0.0024 300
Q 0.0297 0.0359 0.0215 0.0144 0.0153 300
g 0.7174 0.7999 0.629 0.1708 0.5466 300
• Shrinkage surprisingly close to heuristic estimate of 0.79
• Now validate 5-year survival probability estimates
Code
cal <- calibrate(f, B=300, u=5*12, maxdim=3)
Using Cox survival estimates at 60 Months
Code
plot(cal)

## 21.6 Presenting the Model

• Display hazard ratios, overriding default for ap
Code
spar(top=1)
plot(summary(f, ap=c(1,20)), log=TRUE, main='')
• Draw nomogram, with predictions stated 4 ways
Code
spar(ps=8)
surv  <- Survival(f)
surv3 <- function(x) surv(3*12,lp=x)
surv5 <- function(x) surv(5*12,lp=x)
quan  <- Quantile(f)
med   <- function(x) quan(lp=x)/12
ss    <- c(.05,.1,.2,.3,.4,.5,.6,.7,.8,.9,.95)

nom <- nomogram(f, ap=c(.1,.5,1,2,3,4,5,10,20,30,40),
fun=list(surv3, surv5, med),
funlabel=c('3-year Survival','5-year Survival',
'Median Survival Time (years)'),
fun.at=list(ss, ss, c(.5,1:6)))
plot(nom, xfrac=.65, lmgp=.35)