I have been using a logistic model from rms::lrm to estimate the odds ratio of a binary exposure on an outcome, using splines [time] as Frank Harrell recommended. I am using Predict and Contrast to estimate time-varying coefficients which can be further manipulated.

An example using the inbuilt ToothGrowth dataset:

#Install and load packages

#develop a binary predictor variable with values A and B 
ToothGrowth <- ToothGrowth %>% 
  mutate(dose_binary = case_when(       
      dose >1 ~ "A", TRUE   ~ "B")). 

#run logistic regression with 5 restricted cubic splines
dd <- datadist(ToothGrowth)
mod <- lrm(supp ~ dose_binary + rcs(len, 5), data = ToothGrowth)

#predict odds of "supp" in group A and group B, for values of len 1:25
p1 <- Predict(mod, dose_binary = "A", len = 1:25)
p2 <- Predict(mod, dose_binary = "B", len = 1:25)

This works very well, but I would like to estimate the odds (not odds ratio) of exposure in both outcome groups, combined (e.g., the average of p1 and p2); along the time-varying spline function (and I have additional covariates in the final model too).

The underlying problem is to estimate the association between therapy A vs therapy B on the probability of a severe outcome: 

OR (severe) = odds(severe)/odds(control) 

But I lack data on controls. I wish to estimate the odds in controls, borrowing probabilities I DO have from different but related sources (bold indicates data I have access to):  

A) OR(s/m) = odds(severe)/odds(mild) 


B) OR(severe+mild) = odds(severe + mild)/odds(control) 

My logic was to calculate the joint probabilities in both the “severe” and “mild” groups (from A) and estimate the odds in the control group that I lack by: 

Odds(control) = odds(severe + mild, from A)/OR(severe + mild, from B) 

All of these values are modeled along spline functions which is why your package rms has been so useful.  

I’m aware this seems a complex way of doing something simple but the problem is generalizable to other outcomes, exploring duration of protection.

I spent a long time looking at the documentation of Contrast and Predict to see, for example, whether it's possible to incorporate a time-varying adjustment to the predictions directly through weights or another option.


1 Answer 1


Your code brings in an amazing amount of dependencies just in order to handle a trivial recode (see alternate below). And you have not demonstrated equivalences of doses vs. outcomes in order to pool unequal doses. You will find that the outcome depends on how much greater than 1 is the dose.

To your question, you refer to outcome groups but then say you want to combine groups on the basis of baseline variables. p1 and p2 use the same outcome, supp.

Decide on which probability you'd like to estimate, and estimate that, optionally converting it to odds after the fact. Are you wanting to estimate the average probability over the two improper dose groups? (Why not estimate a dose-response curve with continuous dose?). If you want the (almost uninterpretable) average probability you can average the predicted values. You'd have to do a little programming (perhaps using the bootstrap) to get confidence limits.

ToothGrowth <- transform(ToothGrowth, dose_binary = dose > 1)  # base R
  • $\begingroup$ Thank you for much for feedback and apologies for dependencies - you'll see I'm new to R and spend a long time doing things which should be simple. I expanded the query to explain why I wish to perform this calculation. If it still seems illogical, I'll explore other options. Thanks again for your inputs, and also for rms your support for it, which are exceptional resources. $\endgroup$
    – josh
    Aug 8 at 5:23
  • $\begingroup$ Thanks for expanding the OP. But now I'm still not clear on the ultimate goal. If you start with probabilities for now (ignore odds) which events do you seek the probabilities of? $\endgroup$ Aug 8 at 10:53
  • $\begingroup$ I seek the adjusted probability of the entire population (exposed and non-exposed, combined) of experiencing the outcome; on each day, modelled over a spline. Does that make sense or am I trying to calculate something which doesn't (or shouldn't) exist? $\endgroup$
    – josh
    Aug 8 at 13:45
  • $\begingroup$ That exists, but is arbitrary because of the choice of the proportion in the population exposed vs proportion not exposed. The result you seek is a weighted average with those proportions as weights. And you've already made "exposed" arbitrary by dichotomizing dose. $\endgroup$ Aug 8 at 20:24
  • $\begingroup$ Ok, got it. My real problem does use a dichotomous exposure - I only dichotomised using the ToothGrowth dataset as an example - apologies if that wasn't clear. Thank you very much for helpful feedback. $\endgroup$
    – josh
    Aug 10 at 4:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.