title: Easiest way to adjust the p-values for multiple comparisons for model(s) to be shown with msummary from modelsummary package
body: <p>Having used <code>lmerTest::lmer()</code> to perform linear regression with repeated measures data I would like to adjust for multiple comparisons.<br />
  I ran several models and used Bonferroni-Holm to adjust each of them, see my approach below.<br />
  Eventually, I would like to generate a regression table with <a href="https://vincentarelbundock.github.io/modelsummary/articles/modelsummary.html" rel="nofollow noreferrer">modelsummary</a> which should include the adjusted p-values and additional goodness-of-fit statistics (alike in <a href="https://stackoverflow.com/questions/65723430/how-to-extract-the-goodness-of-fit-statistics-from-lmer-model-for-msummary-fro/65725091#65725091">this SO post</a>).</p>
  <p><em>– Maybe there is also a much easier way in <code>modelsummary()</code> to adjust the p-values that I am not aware of yet?</em>
  (both for models individually as well as for/across a group of models altogether)</p>
  <h3>MWE</h3>
  <pre><code>library(&quot;modelsummary&quot;)
  library(&quot;lmerTest&quot;) 
  library(&quot;parameters&quot;)
  library(&quot;performance&quot;)
  
  mod1  &lt;- lmer(mpg ~ hp + (1 | cyl), data = mtcars) 
  mod2  &lt;- lmer(mpg ~ hp + (1 | gear), data = mtcars) 
  
  l_mod &lt;- list(&quot;mod1&quot; = mod1, &quot;mod2&quot; = mod2)
  # msummary(l_mod) # works well
  
  adjMC &lt;- function( model ) {
  model_glht &lt;- glht(model)
  model_mc_adj &lt;- summary(model_glht, test = adjusted('holm')) # Bonferroni-Holm is less conservative and uniformly more powerful than Bonferroni
  return(model_mc_adj)
  }
  mod1_adj &lt;- adjMC(mod1)
  mod2_adj &lt;- adjMC(mod2)
  
  l_mod_adj &lt;- list(&quot;mod1_adj&quot; = mod1_adj, &quot;mod2_adj&quot; = mod2_adj)
   
  </code></pre>
  <p>However, upon adjusting the p-values the class from the model changes from &quot;lmerModLmerTest&quot; to &quot;summary.glht&quot;</p>
  <pre><code>class(mod1) # =&gt; &quot;lmerModLmerTest&quot;
  class(mod1_adj) # =&gt; &quot;summary.glht&quot;
  </code></pre>
  <p>&quot;summary.glht&quot; is among the <a href="https://vincentarelbundock.github.io/modelsummary/articles/modelsummary.html#supported-models-1" rel="nofollow noreferrer">list of supported models of modelsummary</a>,
  and I succeed in obtaining the estimates and p-values:</p>
  <pre><code>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 | &lt; .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)
  </code></pre>
  <p>However, obtaining the goodness-of-fit statistics did not succeed:</p>
  <pre><code>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)
  # =&gt; Error: No glance method for objects of class summary.glht
  
  # e.g. with modelsummary:
  modelsummary::get_gof(mod1_adj) 
  # =&gt; Cannot extract information from models of class &quot;summary.glht&quot;.
  </code></pre>
  <p>To be able to include the adjusted p-values in the final regression table I tried to generate a custom class for &quot;summary.glht&quot; with custom functions to extract estimates and goodness-of-fit information.
  I scanned <code>summary(mod1_adj)</code> for the required information, e.g., <code>summary(mod1_adj)$coef</code> but did not find all the info needed to create the fcts.</p>
  <pre><code>names(mod1_adj$test)
  # [1] &quot;pfunction&quot;    &quot;qfunction&quot;    &quot;coefficients&quot; &quot;sigma&quot;        &quot;tstat&quot;        &quot;pvalues&quot;      &quot;type&quot; 
  
  tidy.summary.glht &lt;- function(x, ...) {
      s &lt;- summary(x, ...)
      ret &lt;- tibble::tibble(term = ...,
                            estimate = s$test$coefficients,
                            ...      = ... ,
                            p-values = s$test$pvalues)
      ret
  }
  
  glance.summary.glht &lt;- function(x, ...) {
    data.frame(
      &quot;Model&quot; = &quot;summary.glht&quot;,
      ... = ...,
      &quot;nobs&quot; = stats::nobs(x)
    )
  }
  </code></pre>