Im trying to replicate the rms::contrast function using the marginaleffects package in R, but am having difficulties. Particularly, I would like to request simultaneous confidence intervals for conditional contrasts in marginal effects, but cannot see how this is done. My initial attempt doesn't even seem to get the same results as rms::contrast when just simply computing point wise confidence intervals so I must be doing quite a bit wrong!
Attempted R code:
library(rms)
library(marginaleffects)
set.seed(1)
age <- rnorm(200,40,12)
sex <- factor(sample(c('female','male'),200,TRUE))
logit <- (sex=='male') + (age-40)/5
y <- ifelse(runif(200) <= plogis(logit), 1, 0)
f <- lrm(y ~ rcs(age,3)*sex)
age_vals <- c(30,40,50)
k <- rms::contrast(f, list(sex='male', age=age_vals), list(sex='female', age=age_vals), type="individual")
# conf.type="simultaneous")
print(k, fun=exp)
comparisons(
f,
comparison = "difference",
variables = "sex",
newdata = data.frame(age=age_vals),
transform = exp
)
First, your examples do not match because the default scale of predicted values is “fitted” for marginaleffects::comparisons() but “lp” for rms::contrast(). You can change the type argument to get matching results:
library(rms)
library(marginaleffects)
set.seed(1)
age <- rnorm(200,40,12)
sex <- factor(sample(c('female','male'),200,TRUE))
logit <- (sex=='male') + (age-40)/5
y <- ifelse(runif(200) <= plogis(logit), 1, 0)
f <- lrm(y ~ rcs(age,3)*sex)
age_vals <- c(30,40,50)
k <- rms::contrast(f, list(sex='male', age=age_vals), list(sex='female', age=age_vals), conf.type="individual")
print(k, fun = exp)
# age Contrast S.E. Lower Upper Z Pr(>|z|)
# 1 30 1.29132 NA 0.2526292 6.600606 0.31 0.7587
# 2 40 2.73344 NA 1.0020803 7.456181 1.96 0.0495
# 3 50 9.67841 NA 0.6565387 142.674938 1.65 0.0982
#
# Confidence intervals are 0.95 individual intervals
comparisons(
f,
variables = "sex",
newdata = data.frame(age=age_vals),
type = "lp",
transform = exp
)
#
#
# Term Contrast Estimate Pr(>|z|) S 2.5 % 97.5 % age
# sex male - female 1.29 0.7587 0.4 0.253 6.60 30
# sex male - female 2.73 0.0495 4.3 1.002 7.46 40
# sex male - female 9.68 0.0982 3.3 0.657 142.67 50
#
# Columns: rowid, term, contrast, estimate, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, age, y
# Type: lp
Second, unfortunately there is currently no way to adjust confidence intervals for multiple comparisons automatically in marginaleffects. However, there is a p_adjust argument which allows you to adjust p values instead. But note that the CIs will remain untouched, even if you use p_adjust.