Restricted Cubic Spline Function Summary Intepretation

Question

I am trying to incorporate spline transformation into my logistic regression and finally piece together the following (working) R code (pls see it below). However, I have no idea how to interpret this summary. My outcome is a binary variable (disease; yes/no) and my predictor is a spline-transformed continuous variable (percentage).

Can someone please walk me through the output below in detail (what are low, high, diff, effect)? Which one is my odds ratio?

Also, what does the number "4" do? I tried to change it to other numbers like 5 or 6 but the output stays the same.

fit <- lrm(formula = disease ~ rcs(percentage, 4), data = final)
summary (fit)  
             Effects              Response : disease 

 Factor      Low  High Diff. Effect   S.E.    Lower 0.95 Upper 0.95
 percentage  15.3 38   22.7  -0.29743 0.77288 -1.81220   1.2174    
  Odds Ratio 15.3 38   22.7   0.74273      NA  0.16329   3.3784    

Here is my data:

structure(list(percentage = c(5.5, 72.1, 7.9, 80.6, 56.3, 11.5, 
15.3, 12.3, 30.9, 27.5, 0.3, 5.3, 19.6, 19.8, 0.3, 40.5, 16.8, 
38, 13.8, 29.9, 15.8, 15.3, 22.8, 17.2, 41.2, 17.2, 31.6, 41.2, 
19.6, 38, 41.2, 29.9, 15.3, 29.9, 38, 30.9, 31.6, 15.3, 15.3, 
38, 31.6, 41.3, 21.4, 0.4, 41.2, 7.6, 29.9), 
    disease = structure(c(1L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 2L, 
    1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L, 
    2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L), .Label = c("none", "disease"), class = "factor")), row.names = c(NA, 
-47L), class = "data.frame")

Answer 1

The summary() function works somewhat differently for objects from the rms package than they do for other R objects. From the help page for summary.rms:

By default, inter-quartile range effects (odds ratios, hazards ratios, etc.) are printed for continuous factors, ...

so what you have in the summary() is for model-prediction differences between the 1st (Low) and the 3rd (High) quartile of your percentage predictor variable, which differ by 22.7 units:

quantile(final$percentage,c(0.25,0.75))
#  25%  75% 
# 15.3 38.0 

Those values of "Low" and "High" are based solely on the data and have nothing to do with the model itself. They are chosen as one example of how the predictor is associated with outcome.

In your case, the reported "Effect" for "percentage" is the log-odds difference between those cases, with the corresponding odds ratio and odds-ratio confidence interval shown on the second line. See the help page of predict.lrm for details.

That's a value of odds ratio, but with a spline the log-odds difference and odds ratio for a given change in the predictor will depend on the predictor values. With the rms package you can get a summary for different values of the predictor. Here are two examples for a 22.7-unit difference in percentage, one starting from percentage = 10 and the other from percentage = 25.

summary(fit,percentage=c(10,10+22.7))
#              Effects              Response : disease 
# 
#  Factor      Low High Diff. Effect   S.E.    Lower 0.95 Upper 0.95
#  percentage  10  32.7 22.7  -0.73109 0.77538 -2.25080   0.78862   
#   Odds Ratio 10  32.7 22.7   0.48138      NA  0.10531   2.20040   

summary(fit,percentage=c(30,30+22.7))
#              Effects              Response : disease 
# 
#  Factor      Low High Diff. Effect S.E.   Lower 0.95 Upper 0.95
#  percentage  30  52.7 22.7  1.7632 1.5807 -1.33490     4.8612  
#   Odds Ratio 30  52.7 22.7  5.8308     NA  0.26319   129.1800  

So with a spline the odds ratio for a specified difference in a predictor can switch from < 1 to > 1, depending on the specific values.

Just typing fit in your case gives a more complete report of coefficients, measures of goodness-of-fit, etc.:

fit
# Logistic Regression Model
#  
#  lrm(formula = disease ~ rcs(percentage, 4), data = final)
#  
#                        Model Likelihood    Discrimination    Rank Discrim.    
#                              Ratio Test           Indexes          Indexes    
#  Obs            47    LR chi2      2.75    R2       0.078    C       0.608    
#   none          17    d.f.            3    g        0.631    Dxy     0.216    
#   disease       30    Pr(> chi2) 0.4315    gr       1.880    gamma   0.222    
#  max |deriv| 9e-06                         gp       0.107    tau-a   0.102    
#                                            Brier    0.221                     
#  
#               Coef    S.E.   Wald Z Pr(>|Z|)
#  Intercept     0.7586 1.0646  0.71  0.4761  
#  percentage    0.0028 0.0927  0.03  0.9758  
#  percentage'  -0.1540 0.2829 -0.54  0.5860  
#  percentage''  0.5206 0.7076  0.74  0.4619  

Harrell's book explains what's shown in this display. The individual terms for percentage, percentage' and percentage'' represent details of the spline fit but should not be interpreted individually. See this page with Frank Harrell's explanation in an answer.

The "4" in rcs(percentage, 4) specifies the number of knots that are used in fitting the continuous predictor. Although more knots give you more flexibility for fitting details of the relationship, changing the number of knots might not necessarily make a big difference in the inter-quartile range effect that is reported. You can see this by calling plot(Predict(fit)) for this model and for models with correspondingly fewer or more knots. The "High," "Low," and "Diff" values are fixed by the data, but I did find numeric differences in the "Effect" etc. with different numbers of knots.

fit3 <- lrm(formula = disease ~ rcs(percentage, 3), data = final)
summary(fit3)
#              Effects              Response : # disease 
# 
#  Factor      Low  High Diff. Effect  S.E.    Lower 0.95 Upper 0.95
#  percentage  15.3 38   22.7  0.20266 0.50987 -0.79666   1.2020    
#   Odds Ratio 15.3 38   22.7  1.22470      NA  0.45083   3.3267    

 
fit5 <- lrm(formula = disease ~ rcs(percentage, 5), data = final)
summary(fit5)
#              Effects              Response : # disease 
# 
#  Factor      Low  High Diff. Effect   S.E.    Lower 0.95 Upper 0.95
#  percentage  15.3 38   22.7  -0.47287 0.79507 -2.03120   1.0854    
#   Odds Ratio 15.3 38   22.7   0.62321      NA  0.13118   2.9607

You also can test for the significance of the predictor and of the nonlinear spline terms with the anova() function (which reports Wald tests on rms objects, again unlike what you might find for other anova() functions in R).

anova(fit)
#                 Wald Statistics          Response: disease 
# 
#  Factor     Chi-Square d.f. P     
#  percentage 1.57       3    0.6670
#   Nonlinear 1.54       2    0.4621
#  TOTAL      1.57       3    0.6670

This particular data set shows no significant evidence of an association between percentage and disease. Note that the number of minority-class cases in the data set is only 17, which in this type of study can probably only support a single degree of freedom in the fit without overfitting. Your 4-knot restricted spline uses up 3 degrees of freedom.