我有一个对比使用lrm在rms封装.如何在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)

推荐答案

是的,你可以使用emmeans来计算47岁和33岁女性的比值比.(我把年龄四舍五入了,因为47.356和32.634似乎对一个人的年龄来说很特殊.

为了清晰起见,我重新调整了模型,没有rms个铃铛和口哨.我也会计算emmeansmarginaleffects的对比度.

# ... 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)

对于emmeans,首先我们指定参考网格(在本例中,性别=两个不同年龄的女性),然后我们计算成对对比度.

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

我们用marginaleffects来计算对比度(即.比较)一步.

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

R相关问答推荐

在R中使用自定义函数时如何删除该函数的一部分?

使用sensemakr和fixest feols模型(R)

提取R中值和列名的所有可能组合

管道末端运行功能

使用R中的Shapetime裁剪格栅文件

如何 bootstrap glm回归、估计95%置信区间并绘制它?

在垂直轴中包含多个ggplot2图中的平均值

多重RHS固定估计

使用R中相同值创建分组观测指标

用预测NLS处理R中生物学假设之上的误差传播

在某些栏和某些条件下,替换dfs列表中的NA

提取第一个下划线和最后一个下划线之间的任何内容,例外情况除外

当我们有多个反斜杠和/特殊字符时使用Gsubing

有没有可能用shiny 的书签恢复手风琴面板?

如何筛选截止年份之前最后一个测量年度的所有观测值以及截止年份之后所有年份的所有观测值

如何显示准确的p值而不是<;0.001*?

我需要使用ggplot2制作堆叠条形图

R-使用stri_trans_General()将其音译为德语字母

根据排名的顶点属性调整曲线图布局(&Q)

使用显式二元谓词子集化sfc对象时出错