r语言 - 使用modelsummary包中的msummary来显示模型的多个比较，调整p值的最简单方法

Question

使用lmerTest::lmer()对重复测量数据执行线性回归后，我想调整多重比较。
我运行了几个模型，并使用Bonferroni-Holm来调整每个模型，见下面我的方法。
最后，我想用modelsummary生成一个回归表，其中应该包括调整后的p值和额外的拟合优度统计(类似于这篇文章)。

-也许在modelsummary()中还有一种更简单的方法来调整我还没有意识到的p值?(既适用于单个模型，也适用于/跨一组模型)

兆瓦

library("modelsummary")
library("lmerTest") 
library("parameters")
library("performance")
mod1  <- lmer(mpg ~ hp + (1 | cyl), data = mtcars) 
mod2  <- lmer(mpg ~ hp + (1 | gear), data = mtcars) 
l_mod <- list("mod1" = mod1, "mod2" = mod2)
# msummary(l_mod) # works well
adjMC <- function( model ) {
model_glht <- glht(model)
model_mc_adj <- summary(model_glht, test = adjusted('holm')) # Bonferroni-Holm is less conservative and uniformly more powerful than Bonferroni
return(model_mc_adj)
}
mod1_adj <- adjMC(mod1)
mod2_adj <- adjMC(mod2)
l_mod_adj <- list("mod1_adj" = mod1_adj, "mod2_adj" = mod2_adj)
 

然而，在调整p值后，来自模型的类从"lmerModLmerTest";summary.glht"

class(mod1) # => "lmerModLmerTest"
class(mod1_adj) # => "summary.glht"

"summary.glht"是modelsummary的支持模型列表之一，并且我成功地获得了估计和p值:

parameters::model_parameters(mod1_adj)
# Parameter        | Coefficient |   SE |         95% CI | Statistic | df |      p
# --------------------------------------------------------------------------------
# (Intercept) == 0 |       24.71 | 3.13 | [17.84, 31.57] |      7.89 |  0 | < .001
# hp == 0          |       -0.03 | 0.01 | [-0.06,  0.00] |     -2.09 |  0 | 0.064
# e.g. with modelsummary:
modelsummary::get_estimates(mod1_adj) # gives more info than broom::tidy(mod1_adj)

然而，获得拟合优度统计量没有成功:

performance::model_performance(mod1_adj)
# 'r2()' does not support models of class 'summary.glht'.
# Can't extract residuals from model.
# no 'nobs' method is availableModels of class 'summary.glht' are not yet supported.NULL
broom::glance(mod1_adj) # and also for broom.mixed::glance(mod1_adj)
# => Error: No glance method for objects of class summary.glht
# e.g. with modelsummary:
modelsummary::get_gof(mod1_adj) 
# => Cannot extract information from models of class "summary.glht".

为了能够在最终回归表中包含调整后的p值，我尝试为"summary. light "生成一个自定义类使用自定义函数提取估计和拟合优度信息。我扫描了summary(mod1_adj)所需的信息，例如，summary(mod1_adj)$coef，但没有找到创建fcts所需的所有信息。

names(mod1_adj$test)
# [1] "pfunction"    "qfunction"    "coefficients" "sigma"        "tstat"        "pvalues"      "type" 
tidy.summary.glht <- function(x, ...) {
    s <- summary(x, ...)
    ret <- tibble::tibble(term = ...,
                          estimate = s$test$coefficients,
                          ...      = ... ,
                          p-values = s$test$pvalues)
    ret
}
glance.summary.glht <- function(x, ...) {
  data.frame(
    "Model" = "summary.glht",
    ... = ...,
    "nobs" = stats::nobs(x)
  )
}

Answer 1

问题是broom包确实为glht模型提供了tidy方法，但是没有为这些模型提供了glance方法。由于只有部分支持，modelsummary只能提取估计值，而不能提取拟合优度统计数据，因此它会中断。要查看这一点，您可以在ghlt对象上尝试modelsummary中的get_gof和get_estimates。

相比之下，modelsummary可以很容易地从lmerModLmerTest模型中提取估计和拟合优度。因此，一种方法是将lmerModLmerTest对象传递给modelsummary，但是通过定义modelsummary文档中描述的tidy_custom.CLASSNAME方法来动态地修改p值。

估计模型:

library(modelsummary)
library(lmerTest) 
library(multcomp)
mod  <- lmer(mpg ~ hp + (1 | cyl), data = mtcars) 

然后，我们定义一个适合您的模型类的tidy_custom方法(同样，参见上面链接的文档了解详细信息)。

注意，由于某种原因，当从glht对象而不是从lmerModLmerTest模型提取结果时，术语名称略有修改。这是上游包中的一个问题，所以你可能想在那里报告它(不确定它是broom还是performance，但它很容易检查)。在任何情况下，这都很容易修复，只需在新方法中添加一个gsub调用:

tidy_custom.lmerModLmerTest <- function(x, ...) {
  new <- multcomp::glht(x)
  new <- summary(new, test = adjusted('holm'))
  out <- get_estimates(new)
  out$term <- gsub(" == 0", "", out$term)
  out
}
modelsummary(mod, statistic = "p.value")

24.708 (0.000) (0.064) R2玛格。 0.143 R2电导率。 0.674 BIC 187.8 ICC 0.6