3

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?

neither-nor
  • 145

1 Answers1

1

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.

Community
  • 1
juod
  • 2,240
  • 13
  • 20
  • 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? – neither-nor Mar 14 '17 at 21:58
  • 1
    @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 Mar 14 '17 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 Mar 14 '17 at 22:17
  • Would Bender's method work if I want to assume an exponential distribution for survival time, in which case hazard is constant and the inverse function is not defined? – neither-nor Mar 17 '17 at 16:26
  • 1
    @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
  • Thanks! One follow-up question: To find, for example, fitted Weibull baseline hazards parameters, do I simply run flexsurvreg(Surv(lifespan, censor) ~ size + resource, data = pop.surv, dist = "weibull") and then use the est values for shape and scale? – neither-nor Mar 17 '17 at 18:11
  • 1
    @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