I have a contrast using lrm in the rms package. How can I recreate this in emmeans?
# example from rms package
library("rms")
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)
myData <- data.frame(y=y,sex=sex,age=age)
dd <- datadist(myData)
options(datadist="dd")
f <- lrm(y ~ age*sex,data=myData)
# what is the difference in log odds for two ages for females
k <- contrast(f, list(sex='female', age=47.356), list(sex='female', age=32.634))
# Can this be done in emmeans
print(k,fun=exp)
Yes, you can use emmeans to compute the odds ratio for females at the ages of 47 and 33. (I've rounded the ages as 47.356 and 32.634 seem oddly specific for a person's age.)
For clarity I refit the model without the rms bells and whistles. I'll also calculate the same contrast with both emmeans and marginaleffects.
# ... code to generate myData ...
fit.lrm <- lrm(y ~ age * sex, data = myData)
# Log odds ratio
k <- rms::contrast(
fit.lrm,
list(sex = "female", age = 47),
list(sex = "female", age = 33)
)
# Odds ratio
print(k, fun = exp)
#> sex Contrast S.E. Lower Upper Z Pr(>|z|)
#> 1 female 12.70547 NA 4.68864 34.42982 5 0
#>
#> Confidence intervals are 0.95 individual intervals
fit.glm <- glm(y ~ age * sex, family = binomial(), data = myData)
With emmeans first we specify the reference grid (in this case, sex = female at two different ages) and then we calculate the pairwise contrast.
emm <- emmeans(fit.glm, ~ age + sex, at = list(sex = "female", age = c(47, 33)))
# Use `type = "response"` to get the odds ratio (rather than the log odds ratio)
emmeans::contrast(emm, method = "pairwise", type = "response")
#> contrast odds.ratio SE df null z.ratio p.value
#> age47 female / age33 female 12.7 6.46 Inf 1 4.998 <.0001
#>
#> Tests are performed on the log odds ratio scale
# `pairs` is short-hand for `contrast(emm, method = "pairwise")`
confint(pairs(emm, type = "response"))
#> contrast odds.ratio SE df asymp.LCL asymp.UCL
#> age47 female / age33 female 12.7 6.46 Inf 4.69 34.4
#>
#> Confidence level used: 0.95
#> Intervals are back-transformed from the log odds ratio scale
With marginaleffects we calculate the contrast (ie. the comparison) in one step.
marginaleffects::comparisons(
fit.glm,
variables = list(age = c(47, 33)),
newdata = datagrid(sex = "female"),
type = "link",
transform = exp
)
#>
#> Term Contrast sex Estimate Pr(>|z|) S 2.5 % 97.5 % age
#> age 47 - 33 female 12.7 <0.001 20.7 4.69 34.4 40.4
#>
#> Columns: rowid, term, contrast, estimate, p.value, s.value, conf.low, conf.high, sex, predicted_lo, predicted_hi, predicted, y, age
#> Type: link