Advice on computing contrasts of categorical predictor across range of interacting continuous predictor fitted with a restricted cubic spline #rms

Question

Background

After fitting a logistic regression model, I am trying to produce a series contrasts between levels of a categorical variable, x1, across the values of an interacting continuous variable, x2, fitted with a restricted cubic spline. Using the contrast() function in the rms package in R I did:

con <- contrast(model, 
          list(x1='a', x2=vals),
          list(x1='b', x2=vals),
          type='individual',
          conf.type="simultaneous")

where vals is a vector of values for the continuous variable x2.

Aim

I would like to present the contrasts as a graph with x2 on the x horizontal axis and the contrast (or odds ratio) on the vertical y axis so that we get a smooth curve and confidence bands over the entire range of x (or e.g. 95% of the values). So something that looks like this:

Example of options

Suppose x2 ranges from 0 to 100, we could compute contrasts at: x2 = 0,1,2,...,100. Or we could use less values of x2, but still cover its range e.g. compute contrasts at x2 = 0,10,20,...,100. In both cases we can get the smooth contrast curve I decribed above. In the former, the curve will be "smoother" since we have more x2 values to plot, though confidence bands will be wider since we have adjusted for having done multiple contrasts (conf.type="simultaneous" in the above R code). In the latter, the curve will be less smooth, but the confidence bands will be narrower since we have done fewer contrasts and hence the penalty applied is less.

Question

How do we select which values of x2 to use when computing the contrasts i.e. in the 'vals' vector above? It seems from my description above that there is a tradeoff between the smoothness of the "contrast curve" and the width of the confidence bands (assuming we are adjusting for simultaneous inference). How does one approach this problem?

Update:

Here is an example where we want contrasts over the range of x2 (here called visit_no.) from 0 to 10. In the first dataframe k, we use all integers, in the second, k2 only a few. The point estimates are the same but the confidence intervals are not.

Answer 1

The problem

Here's a reproducible case, following an example from the page linked from a comment. Although this might seem like a coding issue at first, it does raise questions about what multiple comparisons to control for.

library(survival)
library(rms)
data(cancer)
dat <- colon
dat$sexC <-  ifelse(dat$sex==1, 'F', 'M')
fit <- cph(Surv(time, etype) ~ sexC * rcs(age, 3), data = dat)

Start by comparing a single contrast against a contrast at 5 different times (the default is 95% confidence intervals, CI).

MS_1 <- contrast(fit, list(sexC = 'F', age = 60), list(sexC = 'M', age = 60), conf.type = 'simultaneous' , fun = exp)
MS_1
#   age    Contrast       S.E.      Lower      Upper     Z Pr(>|z|)
# 1  60 -0.09192097 0.09335485 -0.2748931 0.09105118 -0.98   0.3248

MS_5 <- contrast(fit, list(sexC = 'F', age = c(40,50,60,70,80)), list(sexC = 'M', age = c(40,50,60,70,80)), conf.type = 'simultaneous')

MS_5
#     age    Contrast       S.E.       Lower     Upper     Z Pr(>|z|)
# 1    40  0.20882617 0.14419424 -0.15240285 0.5700552  1.45   0.1476
# 2    50 -0.01587881 0.09027919 -0.24204223 0.2102846 -0.18   0.8604
# 3    60 -0.09192097 0.09335485 -0.32578942 0.1419475 -0.98   0.3248
# 4*   70  0.12963034 0.08829577 -0.09156433 0.3508250  1.47   0.1421
# 5*   80  0.47848826 0.19727288 -0.01571104 0.9726876  2.43   0.0153

The first is identical to the default "individual" CI. The report for the second set of contrasts has corrected the CI (not the Z-scores or p-values) for multiple comparisons. The asterisks indicate redundant contrasts that are linearly related to the others. The CI correction is done by the glht() function in the R multcomp package. That function uses a "max-t type test statistic" for a one-step correction, determined by integrating a multivariate distribution defined by the coefficient covariance matrix and the set of requested contrasts.

Although there are many redundant contrasts, the multiple comparison correction changes as the number of requested comparisons changes. For illustration, here are 2 more contrast sets and a function to calculate the width of the 95% CI for the corresponding contrast at age = 60.

MS_3 <- contrast(fit, list(sexC = 'F', age = c(40,60,80)), list(sexC = 'M', age = c(40,60,80)), conf.type = 'simultaneous')
MS_100 <- contrast(fit, list(sexC = 'F', age = 1:100), list(sexC = 'M', age = 1:100), conf.type = 'simultaneous')
CIwidth60 <- function(x) x$Upper[which(x$age==60)]-x$Lower[which(x$age==60)]


CIwidth60(MS_1)
# 0.3659443
CIwidth60(MS_3)
# 0.4449882 
CIwidth60(MS_5)
# 0.4677369 
CIwidth60(MS_100)
# 0.4884676

The CI widths do increase somewhat with the number of requested contrasts. There seems to be an asymptote. The following took over a minute on my decade-old laptop, but shows that a 10-fold further increase in the number of points makes almost no difference.

MS_1000 <- contrast(fit, list(sexC = 'F', age = seq(1,100,by=0.1)), list(sexC = 'M', age = seq(1,100,by=0.1)), conf.type = 'simultaneous')

CIwidth60(MS_1000)
# 0.4884937 

What to do

If there is a specific set of contrasts of interest, then a multiple-comparison correction considering all of them is important.

For a continuous display like what's in the question, however, a correction that depends simply on the number of time points at which you happen to interpolate the continuous predictions would seem to be potentially misleading.

I would suggest a simple transparent approach: use the pointwise individual CI and specify that choice in your figure legend. Pointwise CI are the default both for contrast.rms() and related functions like Predict(). That way a reader will know just what you have done and will be able most readily to reproduce your work.

An alternative (assuming that there is always an asymptotic limit to the CI widths) would be to choose a very large number of time points and use simultaneous CI. Again, specify exactly what you did.