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)$?