Cox Proportional Hazards - Survival Function Estimation (on new data)?

Question

I understand that the relationship between the hazard function, the baseline hazard function and the covariates is the following:

$$ ln(\frac{h(T)}{h_0(T)}) = \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p $$

Where my $\hat{\beta_j}$'s are given by coef in the code below:

library(survival)  

cox <- coxph(Surv(time, status) ~ trt + celltype + karno + diagtime + age + prior, data = veteran)  
summary(cox)

What I am having trouble understanding is the next step - how exactly is R estimating the survival function for new observations?:

cox_fit_newdata <- survfit(cox, newdata = veteran[c(1, 18, 77), ])  
plot(cox_fit_newdata)

I found this document, and page 7 seems to describe what i'm trying to do:

Once the maximum likelihood estimates have been obtained, it may be of interest to estimate the survival probability of a new or existing individual with specific covariate settings at a particular point in time

I'm just having real difficulty connecting the math in here with the survival function produced by R. I'm able to extract the predicted survival rates for these three observations, as seen on the graph, with the code: summary(cox_fit_newdata)[6]

I would like to reproduce the steps taken - how is R going from $\hat{\beta_1}, \dots, \hat{\beta_p}$, and a set of new observations $x_{i1}, \dots, x_{ip}$, to the estimated survival function for that observation displayed on the graph: $\hat{S}_i(T)$?

Answer 1

The model estimates the coefficients for each variable and a baseline hazard with covariate values set to zero.

The values contained in new data are then multiplied by the obtained coefficients and summed. The summed value is exponentiated giving the hazard ratio. The baseline hazard at each time point is then multiplied by the hazard ratio, transformed to the cumulative hazard and then transformed to survival.