`marginaleffects` error: Matrix columns are not supported

Summary

I am trying to using marginaleffects::plot_predictions() to plot predictions from a GLMM. A minimal working example is provided below.

The aim

I'm trying to plot the predictions only for the four-way interaction scaled_exp * Agreement * Condition * Group in the following GLMM:

glmer(Acceptance ~  scaled_exp * Agreement * Condition * Group +
                    scaled_age * Agreement * Group +
                    (1 | Item) +
                    (1 + Condition | Subject)

However, it seems marginaleffects::plot_predictions() runs into problems (only with GLMMs, not with LMMs) if not every single predictor is included in the conditions list (i.e., we're not interested in plotting the effect of scaled_age, it's only included to control for aging effects on the studys' results), but if I try to run it without scaled_age, I get the following error.

Warning: Matrix columns are not supported as predictors and are therefore omitted. This may prevent computation of the quantities of interest. You can construct your own prediction dataset and supply it explicitly to the newdata argument.Warning: Some columns are a multi-column type (such as a matrix column): [2, 6]. setDT will retain these columns as-is but subsequent operations like grouping and joining may fail. Please consider as.data.table() instead which will create a new column for each embedded column.Error: Unable to compute predicted values with this model. You can try to supply a different dataset to the newdata argument. This error was also raised: object 'scaled_age' not found. Bug Tracker: https://github.com/vincentarelbundock/marginaleffects/issues

I don't run into these problems with corresponding LMMs with similar fixed and random effects structures, so I suspect it has something to do with that.

MWE

# Load necessary libraries
library(lme4)
library(tidyverse)
library(marginaleffects)

# Create a sample dataset
set.seed(123)
n <- 5000

sample_data <- tibble(
  Subject = as.factor(rep(1:500, each = 10)),  # Increase the number of subjects
  Group = factor(rep(c("Group1", "Group2"), each = n / 2)),
  Item = sample(1:120, n, replace = TRUE),
  Condition = factor(rep(c("Cond1", "Cond2"), times = n / 2)),
  Agreement = factor(sample(c("Agree", "Disagree"), n, replace = TRUE)),
  Acceptance = rbinom(n, 1, 0.5),
  scaled_age = scale(rnorm(n, 40, 15)),
  scaled_exp = scale(rnorm(n, 0.3, 1.4))
)

# Fit the generalized linear mixed model
glmm_mod <- glmer(Acceptance ~ 
                    scaled_exp * Agreement * Condition * Group +
                    scaled_age * Agreement * Group +
                    (1 | Item) +
                    (1 + Condition | Subject), 
                  data = sample_data, 
                  family = binomial, 
                  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

# Generate predictions for the four-way interaction using plot_predictions()
plot_predictions(
  glmm_mod,
  condition = list("scaled_exp", "Agreement", "Condition", "Group"),
  re.form = NA
) +
  aes(linetype = Agreement) +  # Map linetype to Agreement
  facet_wrap(~ Group + Condition, ncol = 4) +
  scale_color_manual(values = c("Agree" = "blue", "Disagree" = "red")) +
  scale_fill_manual(values = c("Agree" = "blue", "Disagree" = "red")) +
  scale_linetype_manual(values = c("Agree" = "solid", "Disagree" = "longdash")) +  
  theme_classic() +
  labs(x = "Scaled Experience", y = "Predicted Probability of Acceptance")

If you read the error message again carefully, you'll see that a lot of it has to do with "Matrix Columns" not being supported by marginaleffects. And if you inspect your dataset, you'll notice that some of your columns are indeed of class "matrix":

> class(sample_data$scaled_age)
[1] "matrix" "array" 

Turns out that the scale() function returns a matrix by default, and not a numeric vector. To fix this, you can use the drop() function:

library(lme4)
library(tidyverse)
library(marginaleffects)

set.seed(123)
n <- 5000

sample_data <- tibble(
  Subject = as.factor(rep(1:500, each = 10)),  # Increase the number of subjects
  Group = factor(rep(c("Group1", "Group2"), each = n / 2)),
  Item = sample(1:120, n, replace = TRUE),
  Condition = factor(rep(c("Cond1", "Cond2"), times = n / 2)),
  Agreement = factor(sample(c("Agree", "Disagree"), n, replace = TRUE)),
  Acceptance = rbinom(n, 1, 0.5),
  scaled_age = drop(scale(rnorm(n, 40, 15))),
  scaled_exp = drop(scale(rnorm(n, 0.3, 1.4)))
)

glmm_mod <- glmer(Acceptance ~ 
                    scaled_exp * Agreement * Condition * Group +
                    scaled_age * Agreement * Group +
                    (1 | Item) +
                    (1 + Condition | Subject), 
                  data = sample_data, 
                  family = binomial, 
                  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

plot_predictions(
  glmm_mod,
  condition = list("scaled_exp", "Agreement", "Condition", "Group"),
  re.form = NA
)