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