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:

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?

• @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