CODEWITHSUNDEEP
rregressionlme4lmertestmodelsummary
mavericks
mavericks

Reputation: 1165

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

Having used lmerTest::lmer() to perform linear regression with repeated measures data I would like to adjust for multiple comparisons.
I ran several models and used Bonferroni-Holm to adjust each of them, see my approach below.
Eventually, I would like to generate a regression table with modelsummary which should include the adjusted p-values and additional goodness-of-fit statistics (alike in this SO post).

– Maybe there is also a much easier way in modelsummary() to adjust the p-values that I am not aware of yet? (both for models individually as well as for/across a group of models altogether)

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)

However, upon adjusting the p-values the class from the model changes from "lmerModLmerTest" to "summary.glht"

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

"summary.glht" is among the list of supported models of modelsummary, and I succeed in obtaining the estimates and p-values:

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)

However, obtaining the goodness-of-fit statistics did not succeed:

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".

To be able to include the adjusted p-values in the final regression table I tried to generate a custom class for "summary.glht" with custom functions to extract estimates and goodness-of-fit information. I scanned summary(mod1_adj) for the required information, e.g., summary(mod1_adj)$coef but did not find all the info needed to create the 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)
  )
}

Upvotes: 1

Views: 1674

Answers (1)

Vincent
Vincent

Reputation: 17725

The problem is that the broom package does have a tidy method for glht models, but does not have a glance method for such models. Since it is only partially supported, modelsummary can only extract estimates, but not goodness-of-fit statistics, so it breaks. To see this, you could try the get_gof and get_estimates from modelsummary on a ghlt object.

In contrast, modelsummary can easily extract both estimates and goodness-of-fit from lmerModLmerTest models. One approach would thus be to pass a lmerModLmerTest object to modelsummary, but to modify the p values on the fly by defining a tidy_custom.CLASSNAME method as described in the modelsummary documentation.

Estimate model:

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

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

Then, we define a tidy_custom method appropriate for your model class (again, see the documentation linked above for details).

Note that, for some reason, the term names are modified slightly when extracting results from a glht object instead of a lmerModLmerTest model. This is a problem in the upstream package, so you may want to report it there (not sure if it's broom or performance, but it would be easy to check). In any case, this is easy to fix for our purposes by simply adding a gsub call to our new method:

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

Upvotes: 1

Related Questions