# Generating random survival values using Cox model estimates

Background

I'm fresh to survival analysis and I'm using R's survival and coxme libraries to evaluate the effects of two covariates -- population size and resource level -- on the lifespan (in weeks) of local populations.

I scaled down population size by 100 and resource measure by 10. From the subsequent censored data frame:

> head(pop.surv)
location lifespan censor size resource
1       13        2      1 3.10      0.0
2       13        1      1 0.68      0.0
3       26        2      1 2.02      0.0
4       26        2      1 2.04      0.0
5       30        3      1 5.23      0.1
6       13        1      1 5.22      0.0


I ran a mixed-effect cox-proportional hazard model:

res <- coxme(Surv(lifespan, censor) ~ size + resource + (1|location), data=pop.surv)


Based on the result,

> summary(res)
Cox mixed-effects model fit by maximum likelihood
Data: pop.surv
events, n = 1940, 1940
Iterations= 23 165
NULL Integrated    Fitted
Log-likelihood -12751.36  -12318.14 -12288.69

Chisq    df p    AIC    BIC
Integrated loglik 866.45  3.00 0 860.45 843.74
Penalized loglik 925.35 21.51 0 882.33 762.50

Model:  Surv(lifespan, censor) ~ size + resource + (1 | location)
Fixed coefficients
coef exp(coef)    se(coef)      z       p
size     -0.01693793 0.9832047 0.003612058  -4.69 2.7e-06
resource -0.15943564 0.8526248 0.007610163 -20.95 0.0e+00

Random effects
Group    Variable  Std Dev   Variance
location Intercept 0.3320527 0.1102590


I interpret that, holding the other covariate constant, an additional 100 members in a population reduces the weekly hazard of extinction by a factor of 0.9832 on average -- that is, by 1.68 percent. Similarly, each 10 unit increase in resource level reduces the hazard by a factor of 0.8526, or 14.74 percent.

Question

Based on this knowledge, I now want to write a predictive function survfunc(s,r) that takes the arguments of population size s and resource level r, then outputs a survival distribution with a covariate-dependent hazard rate and randomly samples a lifespan value from it. How would I do that?

Cox regression is an attractive choice for survival analysis because it focuses on comparing covariate levels without assuming much about the overall survival distribution. However, simulating the survival times is not so easy - you will need to specify the baseline hazard rate. Without it, your model does not know anything about the scale and shape of the distribution that should result, only that it should become "steeper" in the high-risk groups.

Typically, baseline hazards are chosen so that the resulting survival time distributions would be exponential, Weibull, or Gompertz. My advice would be to try fitting these functions, using, for example, flexsurv package. Once you choose the function, you can take the parameter values obtained from the fit, and then use the inverse hazard function to generate the survival times. Bender et al. present a nice overview of this method and include formulas for the typical distributions, and this answer will be helpful as well.

• Thanks! I didn't think the process would be this complicated. So far, I have fitted the entire lifespan data to a negative binomial due to the data's high number of zeros. Would the fitted parameters from that be appropriate? Commented Mar 14, 2017 at 21:58
• @neither-nor You can apply the inverse hazard method as long as you have the inverse hazard function for the underlying distribution. But I don't have an intuition on how complicated it would be to do that for a neg binomial. Worst case scenario - you can always write a simulation to loop through each day and "roll a die" for each survivor on that day.
– juod
Commented Mar 14, 2017 at 22:15
• By the way, where do the zeros come from? Is that a very short lifespan, or is it something else - maybe censoring could be more appropriate then?
– juod
Commented Mar 14, 2017 at 22:17
• @neither-nor exponential distribution works fine, the inverse is $t/\lambda$, or $t e^{-\beta x} / \lambda$, with $t = -log(Unif)$. (I should have written that you need the inverse of cumulative hazard function $H(t)$, not $h(t)$.)
– juod
Commented Mar 17, 2017 at 17:05
• @neither-nor yes, that would be my approach. For Weibull in particular additional care is needed because it has several common parametrizations - if your simulations suddenly look weird, you might need to fiddle with the parameters, e.g. use $1/\lambda$ instead of $\lambda$.
– juod
Commented Mar 17, 2017 at 19:14