I want to set up a Cox proportional hazards model with two correlated covariates $x$ and $y$. The two covariates are defined as this in R:
x <- rexp(n)
y <- 2x + rnorm(n)
Does the hazard function look like this:
$h(t|x) = h_0(t)\times e^{\beta_1x+\beta_2y}$ or $h(t|x) = h_0(t)\times e^{\beta_1x+\beta_2y+\beta_1\beta_2xy}$?
This depends on how you want to model survival: with x and y independently associated with outcome (although they are correlated) or if you want their product (interaction) also to be associated with outcome. In other words: do you want the effect of x on outcome to depend on the value of y (and vice-versa) or not?
For the simple no-interaction situation, your first formula is correct.
An interaction would most generally be modeled as
$$h(t|x) = h_0(t)\times e^{\beta_1x+\beta_2y+\beta_3 xy}.$$
If you are simulating data and want to define $\beta_3=\beta_1 \beta_2$ exactly, your second formula would be OK. If you are fitting a model with an interaction to data, then the above form is what is assumed: $\beta_3$ is modeled separately from either of $\beta_1$ or $\beta_2$.
