Calculating robust SEs following full matching (MatchIt) with mixed-effect models (lme4)

I am attempting to estimate a treatment effect for a study that features both (a) hierarchical data (students within classrooms; implemented with lme4) and (b) matching between treatment and control students (implemented with the "full" method of MatchIt).

I have been following the excellent MatchIt vignette "Estimating Effects After Matching," which details how to estimate treatment effects after matching using robust SEs, relying on marginaleffects to estimate effects and sandwich to calculate robust SEs. However, while the vignette shows how to do this for ordinary linear regression models (lm), I'm running into problems when trying to apply this method to our mixed-effect model (lmer).

marginaleffects::avg_comparisons appears to work with lmer model objects (with the caveat that "marginaleffects only takes into account the uncertainty in fixed-effect parameters"):

library(lme4)
library(marginaleffects)

data("sleepstudy")

model <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

avg_comparisons(model)

However, when I request cluster-robust SEs (by using the vcov argument), I get a message saying that sandwich isn't compatible with lme4 model objects.

library(sandwich)

# Creating an arbitrary cluster ID for the sake of this example
sleepstudy$Cluster = as.numeric(sleepstudy$Subject) %% 10 

avg_comparisons(model, vcov = ~ Cluster)

Message:

Unable to extract a variance-covariance matrix using this vcov argument. Standard errors are computed using the default variance instead. Perhaps the model or argument is not supported by the sandwich ('HC0', 'HC3', ~clusterid, etc.) or clubSandwich ('CR0', etc.) packages. If you believe that the model is supported by one of these two packages, you can open a feature request on Github.

I've been Googling, and I can't seem to find a way to estimate the treatment effect that uses both (a) the mixed-effects model appropriate to the data structure and (b) robust standard errors appropriate for the matching method. Any suggestions?

You don't need to fit a mixed effects model. Fit a linear model with treatment as a covariates and use a two-way cluster-robust standard error with subclass (i.e., subclass membership from the full matching) and Subject as the clustering variables. That is, you should follow the usual steps for estimating effects after matching but have vcov = ~subclass + Subject in the call to avg_comparisons(). There is no reason to use a mixed effects model if your goal is to estimate an average treatment effect.