I have a Cox proportional hazards model in R (see made-up example below) that models the effect of some variable, say weight. From this model, I'd like to extrapolate what a change in weight from say 90 to 60 would mean to survival, taking into account the fact that for such a change occurring at say age 40, certain amount of risk has already accumulated (and assuming weight change is instantaneous).
I've attached some code which involves
- fitting the Cox model (using age as the time scale);
- extracting the predicted cumulative survival $S(t)$ using survfit for weight=90 and 60;
- getting the cumulative hazard $H(t) = -\log(S(t))$;
- getting the "instantaneous" hazard $h(t)$ via differencing $H(t)$ (plus small fudge factor to avoid zero hazard), which seems to do the job but probably a bit hacky;
- adding a constant to the $\log(h(t))$ for all timepoints after the change, equivalent to the $\beta$ coefficient from the Cox regression times the difference in weights (90-60=30);
- get the new survival functions $S^\prime(t)$ as $\exp(-{\rm cumsum}(\exp(\log(h^\prime(t)))))$.
This procedure produces reasonable results (plotted as $1 - S(t)$), but is it correct or am I just lucky?
library(survival)
set.seed(1)
rm(list=ls())
# Simulate some semi-realistic data
n <- 1e3
age <- round(runif(n, 1, 60))
weight <- round(rnorm(n, 70, 10))
height <- round(runif(n, 1.3, 1.9), 2)
sex <- sample(c("M", "F"), length(age), replace=TRUE, prob=c(0.7, 0.3))
d.time <- ceiling(rexp(n, weight / 1e4))
cens <- round(runif(n, 1, 60))
death <- d.time <= cens
d.time <- pmin(d.time, cens)
d <- data.frame(age=age, weight=weight, height=height, difftime=d.time,
time=d.time + age, sex=sex, death=death)
s <- coxph(Surv(age, time, death) ~ height + weight, data=d)
d.new <- data.frame(weight=c(60, 90), height=1.7)
sf <- survfit(s, d.new)
# The cumulative hazard function H is -log(S(t)) where S(t) is the survivor function
# (aka cumulative survival)
S <- sf$surv[,2]
# Assume we start off with high weight
H <- -log(S)
# The hazard is the derivative (here, finite difference) of the cumulative hazard H
# But the hazard can't be zero exactly as when we take log hazard, won't make sense
h <- diff(c(0, H)) + 1e-6
# We introduce a changepoint in the hazard, but must make sure that the
# hazard does not become negative - this is naturally achieved because the
# Cox model is linear in the log-hazard. This means that the final survivor
# function will always be monotonically decreasing for any value of delta in
# (-Inf, +Inf); delta > 0 increases hazard, delta < 0 decreases hazard
delta <- coef(s)["weight"] * (d.new$weight[1] - d.new$weight[2])
logh <- log(h)
age <- 40
logh[sf$time > age] <- logh[sf$time > age] + delta
h <- exp(logh)
# Get the new cumulative hazard and new survivor functions
H <- cumsum(h)
S <- exp(-H)
# Compare original survivor function with modified one
plot(sf, lwd=5, col=1:2, conf.int=FALSE, mark=NA, fun="event",
xlab="Age", ylab="Cumulative risk")
lines(c(0, sf$time), 1 - c(1, S), type="s", col=3, lwd=5)
abline(v=age, lty=2)
legend(x="topleft", legend=c("Weight=60", "Weight=90", "Weight decreased 90 to 60"),
col=1:3, lwd=5)