# Simulating Time-Dependent Covariate and Applying coxph Function to Estimate Association Parameters

Let $$T$$ denote a failure time random variable drawn from a proportional hazards model with a time-invariant covariate $$X\sim U(2, 7)$$ and a time-varying covariate $$Z(t) = kt$$ for $$k > 0$$. Assuming the baseline hazard function is Exponentially distributed with parameter $$\lambda = 1$$, I can generate $$200$$ observations of right-censored failure time data with covariates from this model using the R code:

library(survival)
# Set proportional hazards model parameters
gamma=c(0.15,0.05)
k = 5
# Generate failure/censoring times and covariate data
failure=c(); censor=c()
obs=c(); ind=c()
time_invariant=c(); time_variant=c()

## Survival Data Simulation
for (i in 1:200){
temp.Z=runif(1, 2, 7)
temp.U= runif(1,0,1)

temp.T=(1/(gamma[2]*k))*log(1+(-(log(temp.U)*gamma[2]*k)/(exp(gamma[1]*temp.Z))))

temp.C = rexp(1, 0.5)
temp.obs = min(temp.T,temp.C)

## Append empty list
failure = append(failure, temp.T); censor = append(censor,temp.C)
obs = append(obs,temp.obs); ind = append(ind, 1*(temp.T==temp.obs))
time_invariant = append(time_invariant,temp.Z)
time_variant = append(time_variant, k*temp.obs)
}


After generating the simulated data, I applied the coxph function treating the variable time_invariant as a time-invariant covariate and I estimate the gamma[1] parameter correctly.

mod = coxph(Surv(obs, ind) ~ time_invariant)
summary(mod)


I want to include the time-varying covariate time_variant into the coxph function to obtain an estimate of gamma[2]. How do I code the variable $$Z(t)$$ so that I can include it as a time-dependent covariate in the coxph function?

• Are you familiar with the vignette on time-dependence in survival models? See Section 5 on "Predictable time-dependent covariates," which seems to be what you need for a time-varying covariate like your Z(t).
– EdM
Commented Oct 7, 2023 at 12:06
• I have tried using the tt code but the estimates do not match with the simulation parameters. My question is how to setup the code to include a time-varying covariate of the form $Z(t) = kt$. The code I used was mod = coxph(Surv(obs, ind)~time_invariant + tt(time_variant), tt = function(x, t, ...) {time_variant <- x*t}) summary(mod) The estimates that were determined were: 4.442e-02 (time_invariant) and -1.887e+01 (time_variant) Commented Oct 7, 2023 at 13:50
• Have you tried expanding the data set yourself to use the counting process format Surv(startTime, stopTime, event) instead of relying on the tt() function? As discussed in the vignette, that also allows for more extensive post-processing of the model fit.
– EdM
Commented Oct 7, 2023 at 15:25
• I have tried to include longitudinal measurements for the time-dependent covariate and included them in the "start", "stop", "event" format but the coxph functions runs into numerical errors. I have essentially the same problem as link but the question of how to include the $Z(t)$ function into the coxph function was not resolved. Commented Oct 7, 2023 at 17:00
• If $Z(t)$ is deterministic, i.e. does not change for different observations, then it's just part of the baseline hazard. Same as a constant covariate is just part of the intercept in a OLS Commented Oct 13, 2023 at 12:48

First of, the baseline hazard is not really distributed, but an arbitrary function from time to $$[0, \infty)$$. A Cox-Model with one time invariant covariable $$x_1$$ and a time variant $$x_2$$ can be setup with a hazard function like this: $$\lambda(t, x_1, x_2) = \lambda_0(t) \cdot \exp(\gamma_1x_1 + \gamma_2f(t, x_2))$$

$$\lambda_0(t)$$ being the baseline hazard. From here you could find the inverse CDF of the failure times for the case specific values of $$x_{1/2}$$ and then use that to sample failure times efficiently or just use your CPU and simulate with tiny time-steps $$\Delta$$ and with probability of failure $$\Delta\lambda(t, x_1, x_2)$$. I have made a few choices in my model to demonstrate the algorithm better, the big one being $$f(t, x_2) = \sin(2\pi t)\cdot x_2$$

library(parallel)
library(survival)

gamma <- c(0.25, 0.5) # i want them a bit larger
base_hazard <- function(t){1}
loc_lambda <- function(t, x1, x2){
base_hazard(t) * exp(gamma[1]*x1 + gamma[2] * sin(t*pi*2) * x2)
}

step_size <- 1/10^5 # very small 0.01 would be perfectly sufficient
sim_case <- function(x1, x2){
failure <- FALSE
t <- 0
while(!failure){
t <- t + step_size
lambda_t <- loc_lambda(t, x1, x2)*step_size
failure <- sample(c(FALSE, TRUE), size = 1, prob = c(1-lambda_t, lambda_t))
}
return(t)
}
sim_case(2, 3)

# setup data, n = 200 is too small for checking the algorithm
n <- 2000

dat <- data.frame(x1 = runif(n), x2 = runif(n))

# simulate with all but 2 of your CPU cores. takes about 1 minute on my machine(AMD Ryzen 5 3600)
start_time <- Sys.time()
print(start_time)
failure_times <- mcmapply(1:n, mc.cores = detectCores() - 2, FUN = function(i){
sim_case(dat[i, "x1"], dat[i, "x2"])
})
print(Sys.time() - start_time)

censor_times <- rexp(n, 2)
dat$$failure <- failure_times < censor_times dat$$time <- pmin(failure_times, censor_times)
summary(dat)
vfit3 <- coxph(Surv(time, failure) ~ x1 + tt(x2),
data=dat,
tt = function(x, t, ...) x * sin(t*pi*2))

summary(vfit3)


We recover the parameter values of $0.25, 0.5$ reasonably well:

Call:
coxph(formula = Surv(time, failure) ~ x1 + tt(x2), data = dat,
tt = function(x, t, ...) x * sin(t * pi * 2))

n= 2000, number of events= 779

coef exp(coef) se(coef)     z Pr(>|z|)
x1     0.1576    1.1707   0.1237 1.274  0.20271
tt(x2) 0.4785    1.6137   0.1808 2.646  0.00814 **



Now your question seems to imply, that you don't actually have a time variant covariate, but instead a feature that increases the hazard equally for all cases. This would be part of the baseline hazard and Cox-models explicitly try to avoid modeling the baseline hazard. I would point you towards Discrete-Time Survival, where the baseline hazard is model and can be tested: https://www.rensvandeschoot.com/tutorials/discrete-time-survival/