Two linear predictors interact significantly (see below). How can I visualize this interaction in a plot?
> data(pbc)
> d <- pbc
> rm(pbc, pbcseq)
> d$status <- ifelse(d$status != 0, 1, 0)
>
> dd = datadist(d)
> options(datadist='dd')
> (m <- cph(Surv(time, status) ~ bili * alk.phos, data=d))
Cox Proportional Hazards Model
cph(formula = Surv(time, status) ~ bili * alk.phos, data = d)
Frequencies of Missing Values Due to Each Variable
Surv(time, status) bili alk.phos
0 0 106
Model Tests Discrimination
Indexes
Obs 312 LR chi2 94.76 R2 0.264
Events 144 d.f. 3 Dxy 0.565
Center 0.676 Pr(> chi2) 0.0000 g 0.641
Score chi2 193.72 gr 1.898
Pr(> chi2) 0.0000
Coef S.E. Wald Z Pr(>|Z|)
bili 0.2280 0.0300 7.59 <0.0001
alk.phos 0.0001 0.0000 1.83 0.0667
bili * alk.phos 0.0000 0.0000 -2.86 0.0043
One way I can think of is to dichotomize one predictor and plot the high values with the low values as two line plots in one figure. However, I cannot reproduce the example in the rms package under plot.Predict()
using the example data above.
> d$alk.phos.high <- ifelse(d$alk.phos > 1259, 1, 0)
> (m <- cph(Surv(time, status) ~ bili * alk.phos.high, data=d))
Cox Proportional Hazards Model
cph(formula = Surv(time, status) ~ bili * alk.phos.high, data = d)
Frequencies of Missing Values Due to Each Variable
Surv(time, status) bili alk.phos.high
0 0 106
Model Tests Discrimination
Indexes
Obs 312 LR chi2 97.95 R2 0.272
Events 144 d.f. 3 Dxy 0.540
Center 0.81 Pr(> chi2) 0.0000 g 0.727
Score chi2 194.00 gr 2.069
Pr(> chi2) 0.0000
Coef S.E. Wald Z Pr(>|Z|)
bili 0.2374 0.0277 8.57 <0.0001
alk.phos.high 0.5667 0.2214 2.56 0.0105
bili * alk.phos.high -0.1139 0.0309 -3.69 0.0002
UPDATE #1
My trying to figure out how to plot both groups of a dichotomized predictor in a single figure, I figured out how to plot lines for several values of one interacting predictor in a plot of the other predictor against the hazard ratio.
I thinks this kind of plot shows the effect of the interaction term in an easy to understand way (which is especially important for a physician g).
- Does this kind of plot have a special name? How would you call this kind of plot?
- How can I interpret this interaction? Would it be correct to say that the prognostic impact of bilirubin increases with lower values for alkaline phosphatase?
library(rms)
data(pbc)
d <- pbc
rm(pbc, pbcseq)
d$status <- ifelse(d$status != 0, 1, 0)
dd = datadist(d)
options(datadist='dd')
m1 <- cph(Surv(time, status) ~ bili * alk.phos, data=d)
p1 <- Predict(m1, bili, alk.phos=c(850, 1250, 2000), conf.int=FALSE, fun=exp)
plot(p1, ylab="Hazard Ratio")
m2 <- cph(Surv(time, status) ~ bili + alk.phos, data=d)
p2 <- Predict(m2, bili, alk.phos=c(850, 1250, 2000), conf.int=FALSE, fun=exp)
plot(p2, ylab="Hazard Ratio")
first figure: model m2
without interaction
second figure: model m1
with interaction