Function for safe use of hypothesis argument

[AKrishna-psy]

Hi all,

as documented in the package and in the warning provided when calling avg_predictions using the hypothesis argument below, the b* wildcard may not be a stable parameter when code is rerun.

It is essential to check the order of estimates when specifying hypothesis tests using positional indices like b1, b2, etc. The indices of estimates can change depending on the order of rows in the original dataset, user-supplied arguments, model-fitting package, and version of marginaleffects.

I do a lot of analyses of multifactorial designs where third- or fourth-order interactions are of interest, often with specific contrasts, so I rely on this function a lot. To keep my analyses reproducible, I recently wrote a function as_marginaleffects_hypothesis that takes a user's input based on term values and translates it to row numbers, given a base avg_predictions output to work from. This allows a user to save the output of avg_predictions without the hypothesis argument, then call as_marginaleffects_hypothesis later on to get a b*-style hypothesis argument robustly.

Would there be any interest in my sharing this function? I'm not familiar with github (neither etiquette nor typical procedures), so apologies if this is impertinent or asked in the wrong place. If there is interest, I can share my code here or elsewhere!

Best,
A. Krishna

[vincentarelbundock]

Thanks for posting. This sounds very interesting!

I do use a function internally to do something like that, but I've found the problem to be very difficult. The issue is that unique row identifiers depend on columns that vary depending on the function and on user input. There is something like this here:

marginaleffects/R/get_labels.R

Line 1 in 21ab438

get_labels <- function(x, idx = NULL, by = NULL, wrap_parens = FALSE, hypothesis_by = NULL) {

But it's not perfect and does not guarantee uniques. Frankly, I wasn't convinced this could be done in a general predictable way, so I gave up on this idea. But I'd be curious to know what you are using.

[AKrishna-psy]

Yeah, it's a knotty problem - I freely admit that my solution hasn't been tested for edge cases and is built on my usage of the package (lmerTest objects with mostly factorial predictors and multiple random intercepts).

as_marginaleffects_hypothesis <- function(
  base_marginaleffects,
  hypothesis_lhs,
  hypothesis_rhs = NULL
) {
  # expected input:
  # base_marginaleffects: the result of a call to avg_predictions
  # or similar functions without the hypothesis argument
  # hypothesis_lhs: a list containing a set of two element types, namely:
  # text chunks - unnamed (!) character strings, such as "+", "-", "(", and ")"
  # row identifiers - character vectors with one element per term column 
  # in base_marginaleffects; each value should be named by column name and 
  # hold the respective value of the desired row
  # hypothesis_rhs (optional): a similarly structured list defining the 
  # right-hand side of the hypothesis formula
  
  # output:
  # a string consisting of the concatenated content of hypothesis_lhs, 
  # ending on "=0" if hypothesis_rhs is not specified or "=" + 
  # concatenated content of hypothesis_rhs if it is.
  
  # in concatenation, row identifiers in the hypothesis_* arguments
  # are referenced against the rows of base_marginaleffects in order 
  # to identify the unique (!) row they describe. The b* wildcard for 
  # this row is then placed in the output at the respective position.
  
  # add a rowname column to read from later
  base_marginaleffects <-
    base_marginaleffects |>
    mutate(
      rowname = str_c("b", row_number())
    )
  hypothesis_lhs <-
    # take the input list and iterate through each element
    # to translate row identifiers into b* wildcards
    modify_if(
      hypothesis_lhs,
      # row identifiers have names, other elements do not
      ~ !is.null(names(.x)),
      \(x) {
        filtered <- 
          # for each element in the row identifier vector...
          imap(
          x, 
          \(y, y_name) {
            # ...output a logical vector that checks whether
            # the respective column in base_marginaleffects
            # has the desired value
            pluck(base_marginaleffects, y_name) == y
          }
          ) |>
          # reduce to a single logical vector that is 
          # TRUE when all columns have the desired values
          reduce(\(a,b) {a & b})
        # output the rowname of that respective row
        # could throw an error here if length != 1
        base_marginaleffects$rowname[filtered]
      }
      ) |>
    # having translated rownames, compress list to a 
    # single character vector
    list_c() |>
    # make the vector into a single text string
    str_flatten()
  # if no rhs argument was supplied, the hypothesis is 
  # lhs = 0
  if (is.null(hypothesis_rhs)) {hypothesis_rhs <- "0"} else {
    # same procedure as for lhs above
    hypothesis_rhs <-
      modify_if(
        hypothesis_rhs,
        ~ !is.null(names(.x)),
        \(x) {
          filtered <- 
            imap(
              x, 
              \(y, y_name) {
                pluck(base_marginaleffects, y_name) == y
              }
            ) |>
            reduce(\(a,b) {a & b})
          base_marginaleffects$rowname[filtered]
        }
      ) |>
      list_c() |>
      str_flatten()
  }
  # put both sides together
  output <-
    str_c(hypothesis_lhs, "=", hypothesis_rhs)
  return(output)
}

An example:

mod <- 
  lmer(
    mpg ~ cyl*am + (1|hp), 
    data = mtcars |> 
      mutate(cyl = factor(cyl), am = factor(am))
    )

# let the hypothesis be:
# the effect of am (transmission) on mpg is larger for 
# six-cylinder cars than for four-cylinder cars

# version 1:

ver1 <- c("cyl","am")
avg_predictions(
  mod,
  variables = ver1
)

# test is (b4 - b3) - (b2 - b1) = 0

avg_predictions(
  mod,
  variables = ver1,
  hypothesis = "(b4 - b3) - (b2 - b1) = 0"
)

# version 2:

ver2 <- c("am","cyl")
avg_predictions(
  mod,
  variables = ver2
)

# test is (b5 - b2) - (b4 - b1) = 0

avg_predictions(
  mod,
  variables = ver2,
  hypothesis = "(b5 - b2) - (b4 - b1) = 0"
)

# using as_marginaleffects_hypothesis

avg_predictions(
  mod,
  variables = ver1,
  hypothesis = as_marginaleffects_hypothesis(
    base_marginaleffects = avg_predictions(mod, variables = ver1),
    hypothesis_lhs = 
      list(
        "(", 
        c(cyl = "6", am = "1"), 
        "-", 
        c(cyl = "6", am = "0"),
        ")-(",
        c(cyl = "4", am = "1"), 
        "-", 
        c(cyl = "4", am = "0"), 
        ")"
      )
  )
)


avg_predictions(
  mod,
  variables = ver2,
  hypothesis = as_marginaleffects_hypothesis(
    base_marginaleffects = avg_predictions(mod, variables = ver2),
    hypothesis_lhs = 
      list(
        "(", 
        c(cyl = "6", am = "1"), 
        "-", 
        c(cyl = "6", am = "0"),
        ")-(",
        c(cyl = "4", am = "1"), 
        "-", 
        c(cyl = "4", am = "0"), 
        ")"
      )
  )
)

I will say I'm not very happy with the input format of the hypothesis_lhs/_rhs arguments. It seems nonintuitive. Perhaps it would be better to have a string argument that takes the form "c1 + c2 = c3 - c4", i.e., parallel to the normal hypotheses argument, and then a list that requires a definition for each parameter in the string argument, like list(c1 = list(factor1 = "factor1valueA", factor2 = "factor2valueA"), c2 = list(factor1 = "factor1valueB", factor2 = "factor2valueA"), c3 = list(factor1 = "factor1valueA", factor2 = "factor2valueB"), c4 = list(factor1 = "factor1valueB", factor2 = "factor2valueB")).

[vincentarelbundock]

On your last point, that's exactly where I ended up everytime. A fully general and robust approach is no longer convenient.

In the end, I kind of think they the easiest approach is to use the b# syntax and to write a couple assertions with stopifnot(), as explained in the warning message.

It feels clunky, but it's easy and quite safe.

[AKrishna-psy]

Very fair. I think I might clean up my function and change its input format as I described anyway, if only for personal use. I'll post it here when I'm done in case anyone else wants to use it. If you do end up revisiting this, feel free to use any portion or the entirety of what I post here with no need for attribution!

[vincentarelbundock]

Thanks, I appreciate all that. And yes, please do post back when you're happy with the results!

[AKrishna-psy]

Done. This is my final version for now; I think it actually ends up handling fairly well. The clunky aspect of needing the output of the marginaleffects function without the hypothesis argument is unavoidable, as is the necessity of using exactly the same variables (or newdata) argument for the subsequent call with the hypothesis argument. For the latter, it might make sense for users to explicitly save their variables/newdata argument as its own object before calling the function for the first time; this would ensure consistency.

library(tidyverse)
library(lmerTest)
library(marginaleffects)

# Made with: 
# marginaleffects 0.25.0
# lmerTest 3.1-3
# lme4 1.1-35.5
# Matrix 1.7-2
# lubridate 1.9.4
# forcats 1.0.0
# stringr 1.5.1
# dplyr 1.1.4
# purrr 1.0.4
# readr 2.1.5
# tidyr 1.3.1
# tibble 3.2.1
# ggplot2 3.5.1
# tidyverse 2.0.0

# Function ----------------------------------------------------------------

as_marginaleffects_hypothesis <- function(
    base_marginaleffects,
    hypothesis,
    hypothesis_terms
) {
  # expected input:
  
  # base_marginaleffects: the result of a call to avg_predictions
  # or similar functions without the hypothesis argument
  
  # hypothesis: a hypothesis of the form "term1 + term2 = 0" or similar;
  # any format that would be valid for a b* wildcard-style hypothesis argument
  # in marginaleffects is acceptable; the terms need NOT be named b*
  
  # hypothesis_terms: a list with elements named after the terms specified in 
  # the hypothesis argument - each element should be a named character vector 
  # containing the column names and respective values in base_marginaleffects; 
  # the combination must uniquely identify one row
  
  # output:
  # a string of the form of a b* wildcard hypothesis term for marginaleffects.
  
  # add a rowname column to read from later
  base_marginaleffects <-
    base_marginaleffects |>
    mutate(
      rowname = str_c("b", row_number())
    )
  # for each element of hypothesis_terms, get the respective rowname
  hypothesis_terms <-
    map(
      hypothesis_terms,
      \(current_term) {
        filtered <- 
          imap(
            current_term,
          \(row_value, column_name) {
            # get a logical vector that checks whether
            # the respective column in base_marginaleffects
            # has the desired value
            pluck(base_marginaleffects, column_name) == row_value
            }
          ) |>
          # reduce to a single logical vector that is 
          # TRUE when all columns have the desired values
          reduce(\(a,b) {a & b})
        # output the rowname of that respective row
        if (length(base_marginaleffects$rowname[filtered]) > 1) {
          msg <-
            str_c(
              "Error in as_marginaleffects_hypothesis:\nThe row values ",
              str_remove(
                reduce2(
                current_term, 
                names(current_term), 
                \(x,y,z) {
                str_c(
                  x,
                  str_c(
                    z, " = ", y
                  ),
                  sep = "; "
                )
              },
              .init = ""
              ),
              "; "
              ),
              " do not identify unique rows in base_marginaleffects."
            )
          stop(msg)
        }
        base_marginaleffects$rowname[filtered]
      }
      )
  output <-
    reduce2(
      hypothesis_terms,
      names(hypothesis_terms),
      \(hypothesis_string, b_wildcard, term_to_replace) {
        str_replace_all(hypothesis_string, term_to_replace, b_wildcard)
      },
      .init = hypothesis
    ) |>
    set_names(names(hypothesis))
  return(output)
}

# Example ---------------------------------------------

mod <- 
  lmer(
    mpg ~ cyl*am + (1|hp), 
    data = mtcars |> 
      mutate(cyl = factor(cyl), am = factor(am))
    )

# let the hypothesis be:
# the effect of am (transmission) on mpg is larger for 
# six-cylinder cars than for four-cylinder cars

# version 1:

ver1 <- c("cyl","am")
avg_predictions(
  mod,
  variables = ver1
)

# test is (b4 - b3) - (b2 - b1) = 0

# version 2:

ver2 <- c("am","cyl")
avg_predictions(
  mod,
  variables = ver2
)

# test is (b5 - b2) - (b4 - b1) = 0

# using as_marginaleffects_hypothesis

avg_predictions(
  mod,
  variables = ver1,
  hypothesis = as_marginaleffects_hypothesis(
    base_marginaleffects = avg_predictions(mod, variables = ver1),
    hypothesis = "(t1-t2)-(t3-t4)=0",
    hypothesis_terms = 
      list(
        t1 = c(cyl = "6", am = "1"), 
        t2 = c(cyl = "6", am = "0"),
        t3 = c(cyl = "4", am = "1"), 
        t4 = c(cyl = "4", am = "0")
      )
  )
)


avg_predictions(
  mod,
  variables = ver2,
  hypothesis = as_marginaleffects_hypothesis(
    base_marginaleffects = avg_predictions(mod, variables = ver2),
    hypothesis = "(t1-t2)-(t3-t4)=0",
    hypothesis_terms = 
      list(
        t1 = c(cyl = "6", am = "1"), 
        t2 = c(cyl = "6", am = "0"),
        t3 = c(cyl = "4", am = "1"), 
        t4 = c(cyl = "4", am = "0")
      )
  )
)

The function retains names in the hypothesis argument, so the output can still be made more easily legible by naming elements of the hypothesis character vector.

By the way, thanks for developing this package, it's been very useful to me!

[vincentarelbundock]

This is very cool, thanks for posting! I re-opened the discussion. Hopefully this gives this a bit more visibility.