I have a study where individuals are randomized to a treatment or control group, and there are two time-to-event outcomes $S_i$ and $T_i$ measured for each individual $i$. The two outcomes have some arbitrary joint distribution in the control group and in the treatment group, with correlations $\rho_C$ and $\rho_T$.

To measure the treatment effect, I plan to estimate the log hazard ratio (e.g. with the Cox proportional hazards model) for each outcome, yielding two log hazard ratio estimates $\hat\mu_S$ and $\hat\mu_T$. What is the joint distribution of $\hat\mu_S$ and $\hat\mu_T$?

I ask because in a simple numerical example (see below) I simulate a scenario with $\rho_C=\rho_T=0.7$ and no censoring, yet the log hazard ratio estimates have simulated correlation of only 0.484, far from the 0.7 I was expecting. As expected, the estimates are approximately bivariate normal, centered at the true log hazard ratios for the two outcomes.

Numerical example (simulation): I simulate 10,000 cohorts with 20,000 individuals each (10,000 in the control group and 10,000 in the treatment group). I simulate the two outcomes for each arm using a bivariate exponential distribution, controlling the correlation to be 0.7 in each arm. In the control group the outcomes have rates 5 and 1, while in the treatment group they have rates 1.8 and 0.05 (so the log hazard ratios for the two outcomes are -1.02 and -3.00).

bivariate.exp <- function(n, rate1, rate2, rho) {
  joint <- rate1*(rho-1)/(1-rate2/rate1-rate1/rate2+rho)
  e <- rexp(n) ; c <- rexp(n, joint)
  data.frame(o1=pmin(e/(rate1-joint), c), o2=pmin(e/(rate2-joint), c))
dat <- replicate(10000, {
  n <- 10000
  ctl <- bivariate.exp(n, 5, 1, 0.7)
  exp <- bivariate.exp(n, 1.8, 0.05, 0.7)
  cph.o1 <- unname(coxph(Surv(c(ctl$o1, exp$o1), rep(1, 2*n)) ~ 
    rep(c("c","e"), c(n, n)))$coefficients)
  cph.o2 <- unname(coxph(Surv(c(ctl$o2, exp$o2), rep(1, 2*n)) ~ 
    rep(c("c","e"), c(n, n)))$coefficients)
  c(cph.o1, cph.o2)

The correlation of the estimates is only 0.484:

cor(dat[1,], dat[2,])
# [1] 0.4836657

The estimates are approximately bivariate normal, centered at the true log hazard ratios:

mvn(t(dat)[1:5000,], mvnTest="hz")
# $multivariateNormality
#            Test        HZ   p value MVN
# 1 Henze-Zirkler 0.9636181 0.2502943 YES
# ...
plot(dat[1,], dat[2,])

enter image description here

  • $\begingroup$ I may not follow all the code, but won't censoring have an effect on correlation? $\endgroup$
    – Todd D
    Jan 25, 2020 at 16:06
  • $\begingroup$ @ToddD there is no censoring in the example provided -- we observe both outcomes for all individuals. $\endgroup$
    – josliber
    Jan 25, 2020 at 23:52
  • $\begingroup$ For those dealing with a problem similar to this one, please, search around competing risks models. Check the following references: academic.oup.com/ndt/article/28/11/2670/1823847 journals.sagepub.com/doi/pdf/10.1177/1536867X1301300209 $\endgroup$
    – Jacob
    Jul 16, 2023 at 2:54
  • $\begingroup$ How are you estimating "log hazard ratios" from a Cox model? Are $\mu_s$ and $\mu_t$ in fact constant (log) hazard rates for processes $S$ and $T$? If so this is an exponential model and you can solve for the joint covariance fairly easily. But I have some doubt you are using notation and modeling concepts precisely here. $\endgroup$
    – AdamO
    Mar 27 at 17:24