Plot the odds ratios (ORs) for a specific covariate from different logistic regression models

Perhaps you are interested in the effect of a specific covariate from several different fitted models. For example, we might be interested in the effect of gender (sex) from several logistic regression models (all with the same outcome). An illustrative graph would be to plot the OR of gender from each of the different fitted logistic regression models in the same plot.

To do this, assign a similar object for each model (e.g., summfit1, summfit2, etc.) using the summary.Design function of a lrm model fit. The needed information is the Effect, Lower 0.95, and Upper 0.95 columns of the output for the Odds Ratio row (2nd row). Then create a data frame from all three model summaries (plotinfo). Then generate the plot.

For this example, let's use our samplefile.txt data file --- not the best example, but a good illustration of the process.

• samplefile.txt

library(Hmisc) library(Design)

x<-read.table("samplefile.txt", header=TRUE) x<-upData(x, labels=c(age="Age", race="Race", sex="Sex", weight="Weight", visits="No. of Visits", tx="Treatment"), levels=list(sex=c("Female", "Male"), race=c("Black", "Caucasian", "Other"), tx=c("Drug X", "Placebo")), units=c(age="years", weight="lbs.")) contents(x)

# Generate an outcome variable x<-upData(x, # Specify population model for log odds that Y=1 L = .4*(sex=='Male') + .045*(age-50) + (log(weight - 10)-5.2)*(-2*(sex=='Female') +

1. 1. *(sex=='Male')), # Simulate binary y to have Prob(y=1) = 1/[1+exp(-L)] y = ifelse(runif(n) < plogis(L), 1, 0))

ddist <- datadist(x) ; options(datadist='ddist')

fit1 <- lrm(y ~ visits + sex * (age + rcs(weight,4)), data = x, x=TRUE, y=TRUE) fit2 <- lrm(y ~ visits + tx + sex * (age + rcs(weight,4)), data = x, x=TRUE, y=TRUE) fit3 <- lrm(y ~ sex * (age + rcs(weight,4)), data = x, x=TRUE, y=TRUE)

summfit1 <- summary(fit1, sex = NA, est.all = FALSE)[2, c("Effect", "Lower 0.95", "Upper 0.95")] summfit2 <- summary(fit2, sex = NA, est.all = FALSE)[2, c("Effect", "Lower 0.95", "Upper 0.95")] summfit3 <- summary(fit3, sex = NA, est.all = FALSE)[2, c("Effect", "Lower 0.95", "Upper 0.95")] plotinfo <- data.frame(rbind(summfit1, summfit2, summfit3), model = 1:3)

with(plotinfo, { plot(Effect ~ model, axes = FALSE, xlab = "Model", ylab = "OR", ylim = c(min(Lower.0.95), max(Upper.0.95))) axis(side = 2) axis(side = 1, at = 1:3, labels = paste("Model", 1:3)) box() arrows(x0 = 1:3, x1 = 1:3, y0 = Lower.0.95, y1 = Upper.0.95, code = 3, angle = 90) })