2
$\begingroup$

I fitted a Cox PH model in R with the survival package and the coxph function. I get the beta estimates from this model. How can I use these coefficients to manually predict on new data, like the predict function does.

In a linear regression this is just the matrix multiplication X %*% beta if $X$ is the data and $beta$ is the vector of coefficients.

How is this in the Cox model? I also see that predict has several options for types of predictions.

here is a minimal example:

library(survival)
data("ovarian")
m <- coxph(formula = Surv(futime, fustat) ~., data=ovarian)

these two give different results:

head(as.matrix(ovarian[, -c(1:2)]) %*% m$coefficients)

      [,1]
1 10.102002
2 10.371810
3  9.706097
4  6.820160
5  7.357138
6  7.627324

head(predict(m, ovarian))
          1           2           3           4           5           6 
 2.66935119  2.93915962  2.27344680 -0.61249088 -0.07551308  0.19467374 
$\endgroup$

2 Answers 2

3
$\begingroup$

I think the discrepancy is due to the fact that predict function in package survival does not yield x*beta but (x-mean)*beta.

$\endgroup$
2
$\begingroup$

It would be good to understand the model from an algebraic standpoint. To help in understanding and to handle complexities such as categorical predictors, nonlinearity and interaction, the R rms package provides two ways to represent a fitted model algebraically. The big step is to get the linear predictor $X\hat{\beta}$ then you need to turn that into things like median survival time and Prob$(T > t | X)$, i.e., survival probabilities.

The cph function is a front-end to coxph, and if you have $\LaTeX$ on your system and you type latex(f) where f is the fit object result from cph, you'll get the full form of the fitted model in terms of survival probabilities. If you don't have $\LaTeX$ you can use Function(f) to get the full algebraic form in R notation, but just for the linear predictor.

$\endgroup$
10
  • $\begingroup$ thanks. When I fit the model in cph and use the latex command I get an output that only prints the linear predictor $X\beta$, the same I cant get with coef(fit) and does not help me. $\endgroup$
    – spore234
    Apr 24, 2016 at 12:53
  • 1
    $\begingroup$ Tell me how you installed $\LaTeX$ and how did you make the system aware of where the $\LaTeX$ executables are. $\endgroup$ Apr 24, 2016 at 12:56
  • $\begingroup$ I use TeXlive on a linux system. Here is what I get for the ovarian dataset from survival package: i.imgur.com/R5asEzF.png $\endgroup$
    – spore234
    Apr 24, 2016 at 13:10
  • $\begingroup$ Excellent. You may have too specity surv=TRUE to cph to get $S_{0}(t)$ defined in the output. Also see the cph time.inc argument. $\endgroup$ Apr 24, 2016 at 17:06
  • $\begingroup$ thanks, that works. I now have $S(t)$ for the intervals and can calculate the survival probabilities at time $t$ for each person. However, that is not exactly what I want, and it's not what predict does (I do not specify any $t$s). I just want to manually reproduce these values: fit <- coxph(...); preds <- predict(fit, newdata) $\endgroup$
    – spore234
    Apr 24, 2016 at 17:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.