mavericks
Reputation: 1165
Adjust p-values obtained with lmerTest::lmer() for multiple comparisons
I want to fit linear mixed models using lmerTest::lmer() and gradually add random and fixed effects (see code below).
Later I aim to compile a regression table including all models with jtools::export_summs() or huxtable::huxreg().
Before this step, I would like to adjust the p-values obtained in the regressions for multiple comparisons using the Bonferroni-Holm (BH) approach.
I stored each of the adjusted models in a list and wrote a function to apply BH to my models as follows:
summary(glht(model), test = adjusted('holm'))
However, when I compile the regression table of the list with adjusted models via huxreg(list_lm_models_adj) or export_summs(list_lm_models_adj) I get the following error message:
"Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column"
Checking the summary of the adjusted and the unadjusted model reveals that the structure seems to change when applying summary(glht(model), test = adjusted('holm')).
Comparing the output of summary(model_lm2) and summary(model2_adjusted) it seems like the random effects are lost in transition.
# Define models
# ------------------------------------------------------------------
# base model: fixed effect: cat1
model_lm0 <- lm(likertscore ~ cat1, data = df_long)
# + random effect: subject => (1 | subject)
model_lm <- lmer(likertscore ~ cat1 + (1 | subject), data = df_long)
# + fixed effect: index => + index
model_lm1 <- lmer(likertscore ~ cat1 + index + (1 | subject), data = df_long)
# full model
# + random effect: group => (1 | group)
model_lm2 <- lmer(likertscore ~ cat1 + index + (1 | subject) + (1 | group), data = df_long)
# 1) unadjusted models => regression table
# ------------------------------------------------------------------
# Store models in list and output regression table
list_lm_models <- list()
list_lm_models[["model_lm"]] <- model_lm
list_lm_models[["model_lm1"]] <- model_lm1
list_lm_models[["model_lm2"]] <- model_lm2
huxreg(list_lm_models)
# ==> provides regression table with unadjusted p-values
# 2) adjusted models => regression table
# ------------------------------------------------------------------
# Function to adjust p-values
adjMC <- function( model_name ) {
model_glht <- glht(model_name)
model_MCadj <- summary(model_glht, test = adjusted('holm')) # Bonferroni-Holm
return(model_MCadj)
}
# Apply function to models
model_lm_adj <- adjMC( model_name = model_lm )
model_lm1_adj <- adjMC( model_name = model_lm1 )
model_lm2_adj <- adjMC( model_name = model_lm2 )
# Store adjusted models in list and output regression table
list_lm_models_adj <- list()
list_lm_models_adj[["model_lm"]] <- model_adjusted
list_lm_models_adj[["model_lm1"]] <- model_lm1_adj
list_lm_models_adj[["model_lm2"]] <- model_lm2_adj
huxreg(list_lm_models_adj)
Any help is greatly appreciated!
AddOn1:
The error occurs when calling huxreg(list_lm_models_adj) or export_summs(list_lm_models_adj) respectively.
df_long looks as follows:
'data.frame': 1715 obs. of 5 variables:
$ subject : Factor w/ 245 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ task : Factor w/ 7 levels "Q1_Level1",..: 1 1 1 1 1 1 1 1 1 1 ...
..- attr(*, "contrasts")= num [1:7, 1:6] 0.25 0.25 0.25 0.25 -0.333 ...
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr "Q1_Level1" "Q1_Level2" "Q1_Level3" "Q1_Level4" ...
.. .. ..$ : chr "c1_CogMem_vs_MechFun" "c2_Cog_vs_Mem" "c3_CogOnly_Math_vs_Words" "c4_MemOnly_Codes_vs_Encrypt" ...
$ likertscore : int 4 3 4 7 4 7 4 7 7 2 ...
$ index : int 7 7 7 7 7 7 7 7 7 7 ...
$ session_code: Factor w/ 24 levels "1t75nw8b","2wwkn7pm",..: 15 15 15 15 15 15 15 15 15 15 ...
> headTail(df_long,8,8)
subject task likertscore index session_code
1 1 Q1_Level1 4 7 lo0h31ts
2 2 Q1_Level1 3 7 lo0h31ts
3 3 Q1_Level1 4 7 lo0h31ts
4 4 Q1_Level1 7 7 lo0h31ts
5 5 Q1_Level1 4 7 lo0h31ts
6 6 Q1_Level1 7 7 lo0h31ts
7 7 Q1_Level1 4 7 lo0h31ts
8 8 Q1_Level1 7 7 lo0h31ts
... <NA> <NA> ... ... <NA>
1708 238 Q1_Level7 1 2 5tc0tw92
1709 239 Q1_Level7 3 2 5tc0tw92
1710 240 Q1_Level7 3 5 v9z7sllr
1711 241 Q1_Level7 4 5 v9z7sllr
1712 242 Q1_Level7 2 5 v9z7sllr
1713 243 Q1_Level7 1 5 v9z7sllr
1714 244 Q1_Level7 4 5 v9z7sllr
1715 245 Q1_Level7 3 5 v9z7sllr
AddOn2: Minimal Working Example
# MWE
# ------------------------------------------------------------------
library("tidyverse")
library("lmerTest")
library("multcomp")
library("huxtable") # or alternatively
# library("jtools")
states <- as.data.frame(state.x77)
df_wide <- states[, c("Frost", "Area")]
colnames(df_wide) <- c("cat1_level1", "cat1_level2")
# add column with "SubjectIDs":
df_wide$subject <- c(paste0("S", 1:(nrow(df_wide))))
df_long <- df_wide %>%
gather(cat1, likertscore, -subject)
# Define models
# ------------------------------------------------------------------
# base model: fixed effect: cat1
model_lm0 <- lm(likertscore ~ cat1, data = df_long)
# + random effect: subject => (1 | subject)
model_lm <- lmer(likertscore ~ cat1 + (1 | subject), data = df_long)
# 1) unadjusted models => regression table
# ------------------------------------------------------------------
# Store models in list and output regression table
list_lm_models <- list()
list_lm_models[["model_lm0"]] <- model_lm0
list_lm_models[["model_lm"]] <- model_lm
huxreg(list_lm_models)
# ==> provides regression table with unadjusted p-values
# 2) adjusted models => regression table
# ------------------------------------------------------------------
# Function to adjust p-values
adjMC <- function( model_name ) {
model_glht <- glht(model_name)
model_MCadj <- summary(model_glht, test = adjusted('holm')) # Bonferroni-Holm
return(model_MCadj)
}
# Apply function to models
model_lm0_adj <- adjMC( model_name = model_lm0 )
model_lm_adj <- adjMC( model_name = model_lm )
# Store adjusted models in list and output regression table
list_lm_models_adj <- list()
list_lm_models_adj[["model_lm0"]] <- model_lm0_adj
list_lm_models_adj[["model_lm"]] <- model_lm_adj
huxreg(list_lm_models_adj) # huxtable
# export_summs(list_lm_models_adj) # jtools wrapper for huxtable::huxreg
# ==> Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column
AddOn3: REPREX
library("tidyverse")
library("lmerTest")
#> Loading required package: lme4
#> Loading required package: Matrix
#>
#> Attaching package: 'Matrix'
#> The following object is masked from 'package:tidyr':
#>
#> expand
#>
#> Attaching package: 'lmerTest'
#> The following object is masked from 'package:lme4':
#>
#> lmer
#> The following object is masked from 'package:stats':
#>
#> step
library("multcomp")
#> Loading required package: mvtnorm
#> Loading required package: survival
#> Loading required package: TH.data
#> Loading required package: MASS
#>
#> Attaching package: 'MASS'
#> The following object is masked from 'package:dplyr':
#>
#> select
#>
#> Attaching package: 'TH.data'
#> The following object is masked from 'package:MASS':
#>
#> geyser
library("huxtable") # or alternatively
#>
#> Attaching package: 'huxtable'
#> The following object is masked from 'package:dplyr':
#>
#> add_rownames
#> The following object is masked from 'package:purrr':
#>
#> every
#> The following object is masked from 'package:ggplot2':
#>
#> theme_grey
# library("jtools")
states <- as.data.frame(state.x77)
df_wide <- states[, c("Frost", "Area")]
colnames(df_wide) <- c("cat1_level1", "cat1_level2")
# add column with "SubjectIDs":
df_wide$subject <- c(paste0("S", 1:(nrow(df_wide))))
df_long <- df_wide %>%
gather(cat1, likertscore,-subject)
# Define models
# base model: fixed effect: cat1
model_lm0 <- lm(likertscore ~ cat1, data = df_long)
# + random effect: subject => (1 | subject)
model_lm <- lmer(likertscore ~ cat1 + (1 | subject), data = df_long)
# 1) unadjusted models => regression table
# Store models in list and output regression table
list_lm_models <- list()
list_lm_models[["model_lm0"]] <- model_lm0
list_lm_models[["model_lm"]] <- model_lm
huxreg(list_lm_models)
#> Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
#> TMB was built with Matrix version 1.2.15
#> Current Matrix version is 1.2.17
#> Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
#> Warning in knit_print.huxtable(x, ...): Unrecognized output format "markdown". Using `to_screen` to print huxtables.
#> Set options("huxtable.knitr_output_format") manually to "latex", "html", "rtf", "docx", "pptx", "md" or "screen".
───────────────────────────────────────────────────── model_lm0 model_lm
─────────────────────────────────── (Intercept) 104.460 104.460
(8532.732) (8532.732)
cat1cat1_level2 70631.420 *** 70631.420 ***
(12067.105) (12064.550)
sd__(Intercept) 1241.576
(NA)
sd__Observation 60322.752
(NA)
─────────────────────────────────── N 100 100
R2 0.259
logLik -1241.651 -1221.720
AIC 2489.303 2451.441
───────────────────────────────────────────────────── *** p < 0.001; ** p < 0.01; * p < 0.05.
Column names: names, model_lm0, model_lm
# ==> provides regression table with unadjusted p-values
# 2) adjusted models => regression table
# Function to adjust p-values
adjMC <- function( model_name ) {
model_glht <- glht(model_name)
model_MCadj <- summary(model_glht, test = adjusted('holm')) # Bonferroni-Holm
return(model_MCadj)
}
# Apply function to models
model_lm0_adj <- adjMC( model_name = model_lm0 )
model_lm_adj <- adjMC( model_name = model_lm )
# Store adjusted models in list and output regression table
list_lm_models_adj <- list()
list_lm_models_adj[["model_lm0"]] <- model_lm0_adj
list_lm_models_adj[["model_lm"]] <- model_lm_adj
huxreg(list_lm_models_adj) # huxtable
#> Warning: Unknown or uninitialised column: 'term'.
#> Warning: Unknown or uninitialised column: 'term'.
#> Error in fix.by(by.x, x): 'by' must specify a uniquely valid column
# export_summs(list_lm_models_adj) # jtools wrapper for huxtable::huxreg
# ==> Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column
Created on 2019-08-24 by the reprex package (v0.3.0)
Upvotes: 0
Views: 4161
Answers (1)
dash2
Reputation: 2262
Your issue can be seen by
tidy(model_lm_adj)
# A tibble: 2 x 6
lhs rhs estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0 104. 8533. 0.0122 0.990
2 cat1cat1_level2 0 70631. 12067. 5.85 0.00000000963
From ?huxreg:
Models must have a generics::tidy() method defined, which should return "term", "estimate", "std.error", "statistic" and "p.value".
The summary.glht class has a tidy method, but it doesn't return a "term" column. So huxreg gets confused. The broom package treads a fine line between "anything goes" and "whips and chains" when it comes to enforcing standards on stats packages that use it.
I will try and improve the error reporting code. Meanwhile, you probably want to use tidy_override:
adj_override <- function( model ) {
model_glht <- glht(model)
model_mc_adj <- summary(model_glht, test = adjusted('holm')) # Bonferroni-Holm
pvals <- tidy(model_mc_adj)$p.value
return(tidy_override(model, p.value = pvals))
}
overridden_models <- lapply(list_lm_models, adj_override)
huxreg(overridden_models)
By the way, you probably don't need all the machinery with glht here. You could just use stats::p.adjust.
Upvotes: 3
Related Questions
- Adjust p-values for multiple comparisons in Matlab
- increase iterations for new version of lmer?
- Easiest way to adjust the p-values for multiple comparisons for model(s) to be shown with msummary from modelsummary package
- FDR correction - extracting p-values from lmer() and creating vectors for use in p.adjust in R
- Pairwise comparisons on lmer using lsmeans or difflsmeans
- Computing p-values for a null random effect model in lm4/lmerTest
- lmerTest::anova not showing p-values
- lmer output very low P values
- lmer Get p values from anova
- anova() does not display p-value when used with lmerTest