How can I calculate the odds of exposure for >1 outcome group, combined, using the rms::lrm function. Predict? Contrast?

Question

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
library("rms")
library("dplyr")

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

#run logistic regression with 5 restricted cubic splines
dd <- datadist(ToothGrowth)
options(datadist='dd')
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) 

and 

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.

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