Unexpected behavior of emtrends() with brms models with me() vs. mi()

Having read that modeling measurement error in `brms` via `me()` terms are soft-deprecated in favor of `mi()`, I’ve been trying to understand how each of these terms work in practice. In doing so, I’ve stumbled across some unexpected behavior (at least to me) of `emmeans` (specifically `emtrends()`) when using `mi()` terms.

The issue arises when trying to calculate the effect of some predictor variable specified via `mi()` using `emtrends()`. Some observations are that (1) the result of `emtrends()` changes every time, (2) the estimated trend deviates from the trend implied by the coefficient(s), and (3) the lower/upper HPD have an unexpectedly large range and magnitude. This does not seem to be an issue when the predictor is specified via `me()`. Does someone have an explanation for this behavior?

Here’s a reproducible example:

``````library(brms)
library(emmeans)

set.seed(359)

# Simulate data
n <- 100
a <- 0
b <- 1
sigma <- 1
x <- rnorm(n, 0, 1)
xsd <- abs(rnorm(n, 0, 1))
xobs <- rnorm(n, x, xsd)
mu <- a + b*x
y <- rnorm(n, mu, 1)
df <- data.frame(y, xobs, xsd)

##### Model with me() #####

# Specify model formula
form_me <- bf(y ~ me(xobs, xsd))

# Fit model
fit_me <- brm(form_me, df, chains = 4, cores = 4)

# Model summary
summary(fit_me)

# emtrends (gives expected results)
# Specify arbitrarily small xsd to get effects without measurement error
emtrends(fit_me, ~ xobs, var = "xobs", at = list(xsd = 1e-10))

##### Model with mi() #####

# Specify model formula
form_mi <- bf(y ~ mi(xobs)) + bf(xobs | mi(xsd) ~ 1) + set_rescor(FALSE)

# Fit model
fit_mi <- brm(form_mi, df, chains = 4, cores = 4)

# Model summary
summary(fit_mi)

# emtrends (gives unexpected results)
emtrends(fit_mi, ~ xobs, var = "xobs", resp = "y", at = list(xsd = 1e-10))
``````

`summary(fit_me)` returns:

`````` Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ me(xobs, xsd)
Data: df (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -0.03      0.12    -0.27     0.20 1.00     2620     2714
mexobsxsd     1.02      0.12     0.77     1.26 1.00     2407     2529

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.00      0.09     0.83     1.19 1.00     1965     2570

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
``````

`summary(fit_mi)` returns:

`````` Family: MV(gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: y ~ mi(xobs)
xobs | mi(xsd) ~ 1
Data: df (Number of observations: 100)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
y_Intercept       -0.03      0.12    -0.26     0.20 1.00     5834     3015
xobs_Intercept     0.07      0.12    -0.17     0.31 1.00     5731     3294
y_mixobs           1.03      0.12     0.79     1.27 1.00     2903     3079

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_y        1.00      0.09     0.84     1.20 1.00     3787     2917
sigma_xobs     1.05      0.09     0.89     1.26 1.00     3657     3437

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
``````

The estimates from the two models are largely equivalent (except for the differences inherent to the different model specifications). However, using `emtrends()` with the two models seems to produce different results:

`emtrends(fit_me, ~ xobs, var = "xobs", at = list(xsd = 1e-10))` returns expected estimates:

``````   xobs xobs.trend lower.HPD upper.HPD
0.0223       1.02     0.774      1.26

Point estimate displayed: median
HPD interval probability: 0.95
``````

`emtrends(fit_mi, ~ xobs, var = "xobs", resp = "y", at = list(xsd = 1e-10))` returns unexpected estimates:

``````   xobs xobs.trend lower.HPD upper.HPD
0.0223         15     -1107      1202

Point estimate displayed: median
HPD interval probability: 0.95
``````

System info:

• Operating System: macOS Ventura 13.4
• brms Version: 2.20.4
• emmeans Version: 1.10.0