我经常使用emmeans
各种统计模型来计算自定义对比。它的优势之一是它的多功能性:它与大量的软件包兼容。我最近发现它emmeans
与该brms
软件包兼容,但无法使其正常工作。我将使用此处提供的数据集进行示例多项逻辑回归分析。我还将在另一个包 ( nnet
) 中进行相同的分析,以证明我需要什么。
library(brms)
library(nnet)
library(emmeans)
# read in data
ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
该数据集包含 200 名学生的变量。结果变量是 prog、程序类型、三级分类变量(一般、学术、职业)。预测变量是社会经济地位,ses,一个三级分类变量。现在通过 nnet 包进行分析nnet
# first relevel so 'academic' is the reference level
ml$prog2 <- relevel(ml$prog, ref = "academic")
# run test in nnet
test_nnet <- multinom(prog2 ~ ses,
data = ml)
现在在 brms 中运行相同的测试
# run test in brms (note: will take 30 - 60 seconds)
test_brm <- brm(prog2 ~ ses,
data = ml,
family = "categorical")
我不会打印两个模型的输出,但两者的系数大致相等
现在创建一个 emmeans 对象,它允许我们进行 pariwise 测试
# pass into emmeans
rg_nnet <- ref_grid(test_nnet)
em_nnet <- emmeans(rg_nnet,
specs = ~prog2|ses)
# regrid to get coefficients as logit
em_nnet_logit <- regrid(em_nnet,
transform = "logit")
em_nnet_logit
# output
# ses = low:
# prog2 prob SE df lower.CL upper.CL
# academic -0.388 0.297 6 -1.115 0.3395
# general -0.661 0.308 6 -1.415 0.0918
# vocation -1.070 0.335 6 -1.889 -0.2519
#
# ses = middle:
# prog2 prob SE df lower.CL upper.CL
# academic -0.148 0.206 6 -0.651 0.3558
# general -1.322 0.252 6 -1.938 -0.7060
# vocation -0.725 0.219 6 -1.260 -0.1895
#
# ses = high:
# prog2 prob SE df lower.CL upper.CL
# academic 0.965 0.294 6 0.246 1.6839
# general -1.695 0.363 6 -2.582 -0.8072
# vocation -1.986 0.403 6 -2.972 -0.9997
#
# Results are given on the logit (not the response) scale.
# Confidence level used: 0.95
所以现在我们有了我们可爱的emmeans()
对象,我们可以用它来执行大量不同的比较。
但是,当我尝试对brms
对象做同样的事情时,我什至没有通过将 brms 对象转换为参考网格的第一步,然后我收到一条错误消息
# do the same for brm
rg_brm <- ref_grid(test_brm)
Error : The select parameter is not predicted by a linear formula. Use the 'dpar' and 'nlpar' arguments to select the parameter for which marginal means should be computed.
Predicted distributional parameters are: 'mugeneral', 'muvocation'
Predicted non-linear parameters are: ''
Error in ref_grid(test_brm) :
Perhaps a 'data' or 'params' argument is needed
显然,不出所料,我不知道有一些步骤可以让贝叶斯软件与emmeans
. 显然,我需要在流程的某个阶段指定一些额外的参数,但我不确定这些参数是否需要brms
在emmeans
. 我在网上搜索过,但找不到简单但详尽的指南。
任何知道如何的人都可以帮我把brms
模型变成一个emmeans
对象吗?