e_plot_model_contrasts: Compute and plot all contrasts and test results from a linear...

In erikerhardt/erikmisc: Erik Erhardt's miscellaneous functions for solving complex data analysis workflows

Compute and plot all contrasts and test results from a linear model by automating the use of emmeans tables and plots.

Description

Variable labels can be provided by labelling your data with the labelled::var_label() command.

Usage

e_plot_model_contrasts(
  fit = NULL,
  dat_cont = NULL,
  choose_contrasts = NULL,
  sw_table_in_plot = TRUE,
  sw_points_in_plot = TRUE,
  sw_ribbon_in_plot = TRUE,
  adjust_method = c("none", "tukey", "scheffe", "sidak", "bonferroni", "dunnettx",
    "mvt")[4],
  CI_level = 0.95,
  sw_glm_scale = c("link", "response")[2],
  sw_print = TRUE,
  sw_marginal_even_if_interaction = FALSE,
  sw_TWI_plots_keep = c("singles", "both", "all")[1],
  sw_TWI_both_orientation = c("wide", "tall")[1],
  sw_plot_quantiles_values = c("quantiles", "values")[1],
  plot_quantiles = c(0.05, 0.25, 0.5, 0.75, 0.95),
  sw_quantile_type = 7,
  plot_values = NULL,
  emmip_rg.limit = 10000
)

Arguments

fit	(required) model object: lm, glm, lmerModLmerTest (from lme4::lmer, or lmerMod (from lmerTest::as_lmerModLmerTest)
dat_cont	(required) data used for the lm object (only used for variable labels using labelled::var_label()
choose_contrasts	is a list of effects to plot, such as c("hp", "vs:wt"); NULL does all in model.
sw_table_in_plot	T/F put table of results in caption of plot
sw_points_in_plot	T/F include marginal points in numeric plot
sw_ribbon_in_plot	T/F include error bands in numeric plot
adjust_method	see ?emmeans::summary.emmGrid
CI_level	level from ?emmeans::emmeans
sw_glm_scale	for glm fit, choose results on the "link" or "response" (default) scale (e.g., for logistic regression, link is logit scale and response is probability scale)
sw_print	T/F whether to print results as this function runs
sw_marginal_even_if_interaction	T/F whether to also calculate marginal results when involved in interaction(s)
sw_TWI_plots_keep	two-way interaction plots are plotted for each variable conditional on the other. Plots are created separately ("singles", default) or together in a grid ("both"), and "all" keeps the singles and the grid version.
sw_TWI_both_orientation	"tall" or "wide" orientation for when both two-way interaction plots are combined in a grid
sw_plot_quantiles_values	"quantiles" or "values" to specify whether to plot quantiles of the numeric variable or specified values
plot_quantiles	quantiles plotted for numeric:numeric interaction plots, if sw_plot_quantiles_values is "quantiles"
sw_quantile_type	quantile type as specified in ?stats::quantile.
plot_values	a named list for values plotted for a single specified numeric:numeric interaction plot, if sw_plot_quantiles_values is "values". Specify a specific contrast, for example: choose_contrasts = "disp:hp". Then specify the values for each variable, for example: plot_values = list(hp = c(75, 100, 150, 200, 250), disp = c(80, 120, 200, 350, 450))
emmip_rg.limit	from emmeans package, increase from 10000 if "Error: The rows of your requested reference grid would be 10000, which exceeds the limit of XXXXX (not including any multivariate responses)."

Details

Plot interpretation: This EMM plot (Estimated Marginal Means, aka Least-Squares Means) is only available when conditioning on one variable. The blue bars are confidence intervals for the EMMs; don't ever use confidence intervals for EMMs to perform comparisons — they can be very misleading. The red arrows are for the comparisons between means; the degree to which the "comparison arrows" overlap reflects as much as possible the significance of the comparison of the two estimates. If an arrow from one mean overlaps an arrow from another group, the difference is not significant, based on the adjust setting (which defaults to "tukey").

Value

out a list of three lists: "plots", "tables" , and "text", each have results for each contrast that was computed. "plots" is a list of ggplot objects to plot separately or arrange in a grid. "tables" is a list of emmeans tables to print. "text" is the caption text for the plots.

Examples

# Data for testing
dat_cont <-
  dat_mtcars_e

# Set specific model with some interactions
form_model <-
  mpg ~ cyl + disp + hp + wt + vs + am + cyl:vs + disp:hp + hp:vs

fit <-
  lm(
    formula = form_model
  , data    = dat_cont
  )

car::Anova(fit)  #, type = 3)
summary(fit)

# all contrasts from model
fit_contrasts <-
  e_plot_model_contrasts(
    fit                = fit
  , dat_cont           = dat_cont
  , sw_print           = FALSE
  , sw_table_in_plot   = FALSE
  , sw_TWI_plots_keep  = "singles"
  , sw_quantile_type   = 1
  )
fit_contrasts$tables  # to print tables
fit_contrasts$plots   # to print plots
fit_contrasts$text    # to print caption text
fit_contrasts$interp  # to print interpretation text


## To plot all contrasts in a plot grid:
# Since plot interactions have sublists of plots, we want to pull those out
#   into a one-level plot list.
# The code here works for sw_TWI_plots_keep = "singles"
#   which will make each plot the same size in the plot_grid() below.
# For a publications, you'll want to manually choose which plots to show.

# index for plot list,
#   needed since interactions add 2 plots to the list, so the number of terms
#   is not necessarily the same as the number of plots.
i_list <- 0
# initialize a list of plots
p_list <- list()

for (i_term in 1:length(fit_contrasts$plots)) {
  ## i_term = 1

  # extract the name of the plot
  n_list <- names(fit_contrasts$plots)[i_term]

  # test whether the name has a colon ":"; if so, it's an interaction
  if (stringr::str_detect(string = n_list, pattern = stringr::fixed(":"))) {
    # an two-way interaction has two plots

    # first plot
    i_list <- i_list + 1
    p_list[[ i_list ]] <- fit_contrasts$plots[[ i_term ]][[ 1 ]]

    # second plot
    i_list <- i_list + 1
    p_list[[ i_list ]] <- fit_contrasts$plots[[ i_term ]][[ 2 ]]

  } else {
    # not an interaction, only one plot

    i_list <- i_list + 1
    p_list[[ i_list ]] <- fit_contrasts$plots[[ i_term ]]

  } # if
} # for

p_arranged <-
  cowplot::plot_grid(
    plotlist  = p_list
  , nrow      = NULL
  , ncol      = 3
  , labels    = "AUTO"
  )

p_arranged |> print()





# one specific numeric:numeric contrast with specified values
fit_contrasts <-
  e_plot_model_contrasts(
    fit                     = fit
  , dat_cont                = dat_cont
  , choose_contrasts        = "disp:hp"
  , sw_table_in_plot        = TRUE
  , adjust_method           = c("none", "tukey", "scheffe", "sidak", "bonferroni"
                                , "dunnettx", "mvt")[4]  # see ?emmeans::summary.emmGrid
  , CI_level                = 0.95
  , sw_glm_scale            = c("link", "response")[1]
  , sw_print                = FALSE
  , sw_marginal_even_if_interaction = TRUE
  , sw_TWI_plots_keep       = c("singles", "both", "all")[3]
  , sw_TWI_both_orientation = c("tall", "wide")[1]
  , sw_plot_quantiles_values = c("quantiles", "values")[2]    # for numeric:numeric plots
  , plot_quantiles          = c(0.05, 0.25, 0.50, 0.75, 0.95) # for numeric:numeric plots
  , sw_quantile_type        = 1
  , plot_values             = list(hp = c(75, 100, 150, 200, 250)
                                 , disp = c(80, 120, 200, 350, 450)) # for numeric:numeric plots
  )
fit_contrasts$plots$`disp:hp`$both


## GLM on logit and probability scales

dat_cont <-
  dat_cont |>
  dplyr::mutate(
    am_01 =
      dplyr::case_when(
        am == "manual"    ~ 0
      , am == "automatic" ~ 1
      )
  )
labelled::var_label(dat_cont[["am_01"]]) <- labelled::var_label(dat_cont[["am"]])

# numeric:factor interaction
fit_glm <-
  glm(
    cbind(am_01, 1 - am_01) ~
      vs + hp + vs:hp
  , family  = binomial
  , data    = dat_cont
  )

car::Anova(fit_glm, type = 3)
summary(fit_glm)

# all contrasts from model, logit scale
fit_contrasts <-
  e_plot_model_contrasts(
    fit                = fit_glm
  , dat_cont           = dat_cont
  , sw_glm_scale       = c("link", "response")[1]
  , sw_print           = FALSE
  , sw_marginal_even_if_interaction = TRUE
  , sw_TWI_plots_keep  = "both"
  )
fit_contrasts

# all contrasts from model, probability scale
fit_contrasts <-
  e_plot_model_contrasts(
    fit                = fit_glm
  , dat_cont           = dat_cont
  , sw_glm_scale       = c("link", "response")[2]
  , sw_print           = FALSE
  , sw_marginal_even_if_interaction = TRUE
  , sw_TWI_plots_keep  = "both"
  )
fit_contrasts


# numeric:numeric interaction
fit_glm <-
  glm(
    cbind(am_01, 1 - am_01) ~
      disp + hp + disp:hp
  , family  = binomial
  , data    = dat_cont
  )

car::Anova(fit_glm, type = 3)
summary(fit_glm)

fit_contrasts <-
  e_plot_model_contrasts(
    fit                = fit_glm
  , dat_cont           = dat_cont
  , sw_glm_scale       = c("link", "response")[2]
  , sw_print           = FALSE
  , sw_TWI_plots_keep  = "both"
  )
fit_contrasts


# factor:factor interaction
fit_glm <-
  glm(
    cbind(am_01, 1 - am_01) ~
      vs + cyl + vs:cyl
  , family  = binomial
  , data    = dat_cont
  )

car::Anova(fit_glm) #, type = 3)
summary(fit_glm)

fit_contrasts <-
  e_plot_model_contrasts(
    fit                = fit_glm
  , dat_cont           = dat_cont
  , sw_glm_scale       = c("link", "response")[2]
  , sw_print           = FALSE
  , sw_TWI_plots_keep  = "both"
  )
fit_contrasts


## Not run: 
## lme4::lmer mixed-effects model
fit_lmer <-
  lme4::lmer(
    formula = Reaction ~ Days + (Days | Subject)
  , data    = lme4::sleepstudy
  )

fit_lmer_contrasts <-
  e_plot_model_contrasts(
    fit                = fit_lmer
  , dat_cont           = lme4::sleepstudy
  , sw_print           = FALSE
  , sw_table_in_plot   = FALSE
  )
fit_lmer_contrasts


## lmerTest::as_lmerModLmerTest mixed-effects model

fit_lmer <-
  lmerTest::as_lmerModLmerTest(fit_lmer)

fit_lmer_contrasts <-
  e_plot_model_contrasts(
    fit                = fit_lmer
  , dat_cont           = lme4::sleepstudy
  , sw_print           = FALSE
  , sw_table_in_plot   = FALSE
  )
fit_lmer_contrasts

## End(Not run)

erikerhardt/erikmisc documentation built on April 17, 2025, 10:48 a.m.