# Joint distribution of log hazard ratio estimates for two outcomes [closed]

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).

library(survival)
set.seed(144)
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:

library(MVN)
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,])  • I may not follow all the code, but won't censoring have an effect on correlation? Jan 25, 2020 at 16:06 • @ToddD there is no censoring in the example provided -- we observe both outcomes for all individuals. Jan 25, 2020 at 23:52 • 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 Jul 16, 2023 at 2:54 • 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. Mar 27 at 17:24