Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

What is the predicted "risk score" for cox models? #919

Open
zepedae opened this issue Sep 28, 2023 · 1 comment
Open

What is the predicted "risk score" for cox models? #919

zepedae opened this issue Sep 28, 2023 · 1 comment

Comments

@zepedae
Copy link

zepedae commented Sep 28, 2023

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:
image

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.

@strengejacke
Copy link
Owner

Do you have a reproducible example?

Maybe it's a bug in labelling the y-axis. You have three options for survival-models:
https://strengejacke.github.io/ggeffects/articles/introduction_plotmethod.html#survival-models

The above plot looks like that are indeed risk-scores, but it's strange that these vary between 0 and 1.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants