Compare KM curve to Cox Proportional Hazard model with multiple variables

I have the cox proportional hazard model that has 4 explanatory variables.

newm3<-coxph(Surv(HoldTime, Indicator) ~ 1 + HIC+SH+DeckDur+Bushels_Scallops, df_model)

I have used the predict function to produce a plot of the model:

#identify variables for new data
HIC<-sort(unique(df_model$HIC)) SH<-mean(df_model$SH)
DeckDur<-mean(df_model$DeckDur) Bushels_Scallops=mean(df_model$Bushels_Scallops)

#use expand.grid to make new dataset for predict
newdata<-expand.grid(HIC,SH,DeckDur,Bushels_Scallops)
names(newdata)<-c("HIC","SH","DeckDur","Bushels_Scallops")

#predict and plot
predict_fit <- survfit(newm3,  newdata= newdata)
plot(predict_fit,ylab="Probability of Survival",xlab="Holding Time (hrs)")

I have been searching around trying to figure out how to compare the coxph fit to a KM curve to assess the model goodness-of-fit. I have only come across examples with one explanatory variable. Does anyone have suggestion for an example with multiple variables or advice on how to proceed?

I did not provide an example because the dataset has ~2000 rows and 19 columns.

• There is no multivariate version of Kaplan-Meier estimates. You can compute and then compare separate estimates across mutually exclusive strata, but that only allows you to compare categorically on one dimension at a time. So you could compare your multivariate model to a K-M curve, but the K-M curve can't account for your covariates. May 6, 2017 at 15:44
• Thanks for your response. That is what I figured, but thought it would be helpful to see if anyone had a different answer. Do you think it is appropriate to compare a multivariate Cox model to a KM curve to assess goodness of fit? I have tested for Cox model assumptions. May 6, 2017 at 16:40
• If you are confident the covariates are relevant, it's probably more informative to compare the Cox model to other multivariate models involving different parametric versions of duration dependency. If you want to check the assumption that the parameters are informative, then comparing multivariate Cox to K-M would be interesting. May 6, 2017 at 17:05
• @ulfelder: Isn't a no covariate model fit the same as a KM curve?
– DWin
May 7, 2017 at 16:53

(I don't understand why you didn't pose the question using one of the available datasets in the 'survival' package.) Generally you can get a no-covariate models with code like:

newm3<-coxph(Surv(HoldTime, Indicator) ~ 1 , df_model)

The coxph function will return a deviance estimate and I'm reasonably sure you can use anova (which is really analysis of deviance).

bladder1 <- bladder[bladder$enum < 5, ] mod_base <- coxph(Surv(stop, event) ~ 1 , bladder1) bladder1 <- bladder[bladder$enum < 5, ]
mod1 <- coxph(Surv(stop, event) ~ (rx + size + number)  , bladder1)
anova(mod1, mod_base)
#----
Analysis of Deviance Table
Cox model: response is  Surv(stop, event)
Model 1: ~ (rx + size + number)
Model 2: ~ 1
loglik  Chisq Df P(>|Chi|)
1 -588.10
2 -599.26 22.321  3 5.594e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This is a fairly generalizable strategy. In the survival package there is also a survfit function that can generate KM-like curves from model-objects that with no newdata argument will be the estimated survival function determined at the mean of the covariates.