Thank you so much for the detailed response!! I have two follow-ups and hope you can help me clarify this method: 1) In the complete output you shared above (fit), how come there are 3 coefficients for percentage vs only 1 in my output? 2) In your final output by anova( ), does it mean the spline-transformed one (non-linear) did not differ from the non-transformed one (percentage), correct? Also, why the df is different between percentage & non-linear?

@RBeginner the 3 coefficients represent details of the spline fit. The summary(fit) shows the model's prediction of the net effect of percentage between its 1st and 3rd quartiles, taking the entire spline fit into account. It's a prediction from the model, rather than a comprehensive summary of it.

@RBeginner For anova(fit), there is 1 linear term (1 df) and 2 nonlinear terms (2 df) for percentage when you have 4 knots. The report for percentage includes all 3 terms. The report for Nonlinear is for the the 2 nonlinear terms themselves. In this case the high p-value for Nonlinear means the nonlinear terms didn't contribute significantly to the model beyond the linear term. In this model, however, including percentage didn't improve significantly over a null model with no predictors, as indicated by the high p value for Total in that display of Wald test results.

Thank you for the great response. So, can I simplify the anova ( ) output as follows: the p-value from percentage tells us the difference between percentage-included model vs percentage-excluded model whereas the p-value from non-linear informs us of the difference between non-transformed (percentage-included) model vs spline-transformed (percentage-included) model?

@RBeginner yes, with some caveats. Here there's only the one continuous predictor percentage, so the percentage report in anova() is for including percentage (with its spline fit) versus the no-predictor model. In a model with additional predictors, that would be for models with versus without percentage but including the other predictors. Also, the best test for the overall model versus null model in general is the "Model Likelihood Ratio Test" included in the output from print(fit).

Got it!! Can you please explain why the odds ratio did not match from the output between summary(fit) vs fit? For example, the 1st quartile's odds ratio based on summary(fit) is 15.3 but it is only 1.002 (e^0.0028) according to the output from fit ?

@RBeginner look at the column headings of summary(fit). An odds ratio is between two defined situations, here "High" and "Low." 15.3 is the "Low" value of percentage that was used. The odds ratio for disease between "High" and "Low" is shown under the "Effect" column, 0.74273. The coef value of 0.0028 is for the change in log-odds per 1 unit change in the linear spline term alone. The odds ratio in the summary(fit) report is for a 22.7-unit difference in percentage and it includes the nonlinear spline terms.

So, the summary(fit) can be summarized as follows: for a 22.7-unit difference in percentage (spline-term included), the odds of having a disease will be 0.74273 (i.e. less likely)?

That's only true for that particular 22.7-unit difference, between percentage values of 15.3 and and 38. With a spline, the effect of a unit change in a continuous predictor on log-odds can vary over the range of predictor values.

So what can we infer from the terms percentage' and percentage'' from the fit output? I know both of them are non-linear spline terms but I can't figure out how to report them in my paper.

Also, I tried to get the odds ratio and 95% CI from '''fit''' output by using the following code '''exp(cbind(OR = coef(fit), confint(fit)))''' but it won't work as I got the following message: $ operator is invalid for atomic vectors. This is strange as this code always works for the regular logistic regression. Can you pls help me spot where it went wrong? Thanks!!

I added to the answer to handle the dependence of odds ratio on specific predictor values with a spline, and to add a link to how Harrell recommends to report spline fits in logistic regression. Use plot(Predict(fit)) to get the partial-effect plot with respect to percentage. Use rms functions like summary and Predict (capital "P") and contrast for getting predictions rather than using your own code, as the structure of rms objects isn't always the same as for other R packages.

Thank you so much! I have now understood it much better. For the partial-effect plot, do you know how I can change the Y-axis from log-odds to the odds ratio?

@RBeginner an odds ratio comes from a log-odds difference between two situations. You'd have to specify a reference value of percentage and then use the fun argument of Predict to specify the calculation you want. The display would then depend heavily on your choice of reference value and would emphasize extreme odds ratios. I'd recommend sticking with the log-odds versus percentage plot.

I see - but the Fig 2 of the article (you shared you with me in another post) used hazard ratio, that's why I also want to use odds ratio in my graphs. I thought we can simply convert log-odds to odds ratio by taking exponential (I don't know how to do that with the graph tho)? Are you able to provide me with the complete code for it?

@RBeginner the help page for Predict describes a fun argument that lets you transform predictions however you wish before plotting. With the argument ref.zero=TRUE to Predict(), the median of the percentage values becomes the reference (log-odds = 0) for your odds ratios. Adding the further argument fun=function(x) exp(x) to Predict() will then do what you wish. Be careful what you wish for, however.

Thank you! If I want to use percentage =10 as the reference, do you know why the following code is wrong? final$limits["Adjust to","percentage"] <- 10 .... I got the error message of Error in final$limits["Adjust to", "percentage"] <- 10 : incorrect number of subscripts on matrix

@RBeginner rms uses a datadist object for such things, not the data frame itself (which is what you called final). For the benefit of future visitors, please edit your question to start your code with library(rms) and to specify the steps for running datadist() and setting the options accordingly after you present the final data frame, so that those less familiar with the package can reproduce what you did. Note that you will have to run fit<-update(fit) after you modify the datadist to get the plot you want.

Will do!! I want to ask one final question about both graphs: I notice the x-axis (predictor) stops at around a value of 50. Do you know what this is based on? Can we extrapolate the x-axis beyond our real data? The max value of percentage is 80.6 but the 3rd Qu is 38 (only three points beyond a value of 50).

The first paragraph of the Predict help page explains that the display is between the "High:prediction" and "Low:prediction" values in the datadist. You can tweak them as you did the "Adjust to" value, if you wish. But don't go too far. The "restricted" part of "restricted cubic spline" in your case is a linear association of log-odds with percentage beyond the outer knots, for exponential growth in your odds-ratio scale. So the extrapolation won't add much useful.

You shouldn't put much faith in extrapolations far beyond the outer knots in any event. Also, think about what will be the most useful display for your audience. The odds ratio with respect to the value at a percentage of 10 puts a lot of emphasis on that choice of reference. A simple display of predicted odds or probability of disease versus percentage might be better.