Survival analyses for complicated models

Question

I've just began looking into survival analyses, and I'm having some difficulty interpreting the results.

I have a model linear regression model where time ~ x1 * x2 * x3... x9, this model seems to fit well (adjusted R${^2}$ of ~0.9) however it was suggested that survival analyses may be a better direction.

I have tried using the R tool for this, but I can't interpret the data. I do have some trails where they timed out (after 1 hour) so I set up my censor column for that. I then tried plotting the full data set (~200 trials), and got a nice downward stepped curve. , however this doesnt seem to tell me anything, the model used was just survfit(my.surv ~ 1). If I replace the ~ 1 with one of the variables, for example survfit(mysurv ~ x1), I get a total mess. Lines all over the place (mostly vertical).

So, my questions are these:

1) Is it possible to gain useful information, other than saying that 80% of the trials finished before 2000 seconds, which is what I think this plot shows

2) Is it possible to include factors in the model, for example time ~ x1 + x2, etc. If not how does one ascertain which variables are important by this method.

3) How can I use this to make predictions about new data, not included in this model?

Answer 1

OK, here goes!

Firstly, when you try doing survfit(mysurv ~ x1) and get a total mess it is probably because x1 is continuous or at least takes lots of different values. You are telling R to calculate separate Kaplan-Meier curves for each value of x1 and if you have multiple variables, e.g., mysurv ~ x1 + x2 + x3 it will attempt to make a separate curve for each possible combination of different values.

The most common approach to investigating whether a variable has an important effect on the survival time is through the Cox proportional hazards model. This assumes that each dependent variable scales the hazard rate by a constant amount for all time. It is important to realise that the proportional hazards assumption doesn't always hold, or sometimes it holds but only with time-varying covariates, etc.

The sort of command you want to do is: coxph(my.surv ~ x1 + x2 + as.factor(x3)) where x1 and x2 are metric variables (discrete or continuous) and x3 is a categorical variable.

You will then get an output similar to (I took an example dataset bfeed from library(KMsurv)):

Call:
coxph(formula = my.surv ~ agemth + as.factor(race))


                     coef exp(coef) se(coef)      z      p
agemth           -0.00196     0.998   0.0133 -0.148 0.8800
as.factor(race)2  0.10881     1.115   0.1031  1.056 0.2900
as.factor(race)3  0.25295     1.288   0.0930  2.721 0.0065

Likelihood ratio test=7.65  on 3 df, p=0.0538  n= 927, number of events= 892

You could interpret these results in the following way:

• The coefficient for agemth (mother's age in months) is only slightly negative and its z-score is very close to zero, therefore mother's age in months is not a good linear predictor of hazard (we might want to investigate further for other possible relationships)
• Race appears to be a plausible predictor of hazard, such that the hazard rate for ceasing breast feeding is higher for women with race coded 2 or 3

The Akaike Information Criterion is calculated as $AIC = -2 \log{L} + 2p$ where $L$ is the likelihood and $p$ is the number of parameters in the model. A lower AIC suggests a better model. For the above model the AIC is $-2\times (-5187.290) + 2\times 3 = 10380.58$.

You could then look at removing agemth as a predictor:

Call:
coxph(formula = my.surv ~ as.factor(race))


                  coef exp(coef) se(coef)    z      p
as.factor(race)2 0.111      1.12   0.1024 1.08 0.2800
as.factor(race)3 0.254      1.29   0.0925 2.75 0.0059

Likelihood ratio test=7.63  on 2 df, p=0.0221  n= 927, number of events= 892

Now we obtain an AIC of $-2\times(-5187.301) + 2\times 2 = 10378.602$, which is lower than the AIC with agemth as a predictor. This suggests that the model without agemth is superior, but obviously there would be much more investigation to be done!

Once you have settled on a model you should really check whether proportional hazards is a valid assumption: I point you to http://www.ats.ucla.edu/stat/examples/asa/test_proportionality.htm

Finally you asked about predictions about new data. Assuming a Cox proportional hazards model you can make predictions about median survival, or probability of survival at a certain point in time, or the restricted mean survival. You cannot make predictions about mean survival without making some assumption about how survival extrapolates beyond the data you already have. If you do make some assumption along those lines then you can use predict.survreg but assuming you do not want to make such an assumption...

You can construct Kaplan-Meier curves based on the proportional hazards model by using survfit(my.coxph, newdata=data.frame(race=c(1,2,3))), which will give you predicted median survival plus summary(survfit(...)) for the Kaplan-Meier tables and plot(survfit(...)) for Kaplan-Meier graphs.

If you want any more detail or help I will need to see your data and understand your question more, but I hope this has been helpful!

(Original response)

You will probably want to look at Cox proportional hazards models, which are included in the survival library you have already loaded.

You will need to do something like Surv(data) ~ x1 + x2 but i can't remember syntax off the top of my head.

EDIT+ You test for significance by comparing proportional hazards coefficients to zero or hazard ratios (the former exponentiated) to one, or by comparing models with different numbers of parameters by Akaike Information Criterion.

There is also prediction functionality if you create a new data frame with desired covariates.

If you want to extrapolate to mean survival you will need to make parametric assumptions about survival, e.g., exponential survival.

If you want the answer fleshed out and tidied up you will have to wait till tomorrow :)