# 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?

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.
• @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 Mar 17 '17 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 Mar 17 '17 at 19:14