调整模型多重比较的 p 值的最简单方法，使用 modelsummary 包中的 msummary 显示

Question

使用 lmerTest::lmer() 对重复测量数据执行线性回归后，我想针对多重比较进行调整。
我运行几个模型并使用 Bonferroni-Holm 调整每个模型，请参阅下面的方法。
最终，我想用 modelsummary which should include the adjusted p-values and additional goodness-of-fit statistics (alike in ).

生成回归 table

– 也许 modelsummary() 中还有更简单的方法来调整我还不知道的 p 值？ （对于单独的模型以及 for/across 一组模型）

MWE

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 值后，模型中的 class 从“lmerModLmerTest”变为“summary.glht”

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

"summary.glht" 属于 list of supported models of 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 值 table 我尝试为“summary.glht”生成自定义 class，并使用自定义函数来提取估计值和拟合优度信息。我扫描了 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，但通过定义 tidy_custom.CLASSNAME method as described in the modelsummary documentation.

来动态修改 p 值

预估型号：

library(modelsummary)
library(lmerTest) 
library(multcomp)

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

然后，我们定义一个适合您的模型的 tidy_custom 方法 class（同样，请参阅上面链接的文档了解详细信息）。

请注意，由于某些原因，从 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")

	Model 1
(Intercept)	24.708
	(0.000)
hp	-0.030
	(0.064)
Num.Obs.	32
R2 Marg.	0.143
R2 Cond.	0.674
AIC	181.9
BIC	187.8
ICC	0.6

调整模型多重比较的 p 值的最简单方法，使用 modelsummary 包中的 msummary 显示

Easiest way to adjust the p-values for multiple comparisons for model(s) to be shown with msummary from modelsummary package

r

regression

lme4

lmertest

modelsummary

MWE