Dear all,
I am using R 3.4.3 on Windows 10. I am preparing some teaching materials
and I'm having trouble matching the by-hand version with the R code.
I have fitted a Cox model - let's use the ovarian data as an example:
library(survival)
data(ovarian)
ova_mod <- coxph(Surv(futime,fustat)~age+rx,data=ovarian)
If I want to make predict survival for a new set of individuals at 100
days then that is trivial using predict.coxph e.g.:
newdata <- data.frame(futime=rep(100,5),fustat=rep(1,5),age=c(45,50,
55,60,65),rx=c(1,2,1,2,1))
preds <- predict(ova_mod,newdata,type="expected")
survpreds <- exp(-1*preds)
survpreds
[1] 0.9967635 0.9969739 0.9859543 0.9868629 0.9401437
However, due to centering I believe, I am finding this a bit difficult to
replicate by hand.
To replicate the analysis I need the baseline survival at 100 days, the
regression coefficients, and the mean/proportions.
Baseline survival at 100 days: summary(survfit(ova_mod),time=100)$surv =
0.9888
Regression coefficient: ova_mod$coef = 0.1473 (age) & -0.8040
Mean age: mean(ovarian$age) = 56.1654
Proportions with rx: prop.table(table(ovarian$rx)) = 0.5 0.5
So, to recreate the predicted survival probability I would work out the
linear predictor for my first individual:
LP1 = ((45-56.1654)*0.1437)+((1*0.5)*-0.8040) = -2.0470
However, predict(ova_mod,newdata,type="lp") suggests that it should be
-1.2430. Therefore have I misinterpreted the centering? I.e that we take
off the mean value from continuous variables, and multiply by the
proportion with that response for binary/categorical variables?
Then, assuming the I have the correct LP of -1.2430, I need to raise the
baseline survival estimate at 100 days to the exp(LP1). This gives 0.9968
which is almost the predicted value of 0.9968 of above. However, I need to
get the linear predictors to agree with the output first!
Many thanks for your help with this.
Kind regards,
Laura
[[alternative HTML version deleted]]