sjPlot
sjPlot copied to clipboard
What is the predicted "risk score" for cox models?
I ran the following time-varying covariate cox model:
model_int <- survival::coxph(Surv(tstart, tstop, status) ~ nat + dist +
sex + sf + popDens + medInc +
popDens*nat + popDens*dist + popDens*medInc,
data = model_df)
and I get the following summary
Call:
survival::coxph(formula = Surv(tstart, tstop, status) ~ nat +
dist + sex + sf + popDens + medInc + popDens * nat + popDens *
dist + popDens * medInc, data = model_df)
n= 43846, number of events= 154
(251 observations deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
nat -0.121507 0.885585 0.124637 -0.975 0.329614
dist -0.195417 0.822492 0.123732 -1.579 0.114255
sexm -0.021949 0.978290 0.168305 -0.130 0.896240
sftransient 0.659347 1.933530 0.182534 3.612 0.000304 ***
popDens 0.552056 1.736821 0.178899 3.086 0.002030 **
medInc 0.006984 1.007008 0.074640 0.094 0.925454
nat:popDens 0.645145 1.906263 0.199470 3.234 0.001219 **
dist:popDens 0.190398 1.209730 0.151516 1.257 0.208894
popDens:medInc 0.019754 1.019950 0.073993 0.267 0.789491
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
I plotted the interaction effect by using plot_model to generate predicted values for 5 quantiles of nat across a range of popDens values (please ignore the title):
sjPlot::plot_model(model_int, type = "emm",
terms = c("popDens [-0.369923, 10]", "nat [-1.44193, -1.0567426, 0.0988196,0.9012933, 1.551297]"),
title = "Predicted Hazard Ratios",
legend.title = "Proportion Natural \nHabitat (quantiles)") +
scale_color_discrete(labels = c("0%", "25%", "50%", "75%", "100%")) +
theme_bw()
And I get the following plot:
So, how does plot_model generate risk scores? Is it the hazard ratio? If so is it the ratio of hazard rates between the new data and sample data? I haven't found anything definitive in the sjPlot or ggemeans documentation.