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



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.