paul-buerkner / brms

brms R package for Bayesian generalized multivariate non-linear multilevel models using Stan
https://paul-buerkner.github.io/brms/
GNU General Public License v2.0
1.25k stars 177 forks source link

Unexpected results of `emtrends()` with meausurement error model #1620

Closed kthayashi closed 3 months ago

kthayashi commented 3 months ago

I've encountered (what I think might be) unexpected behavior when using emmeans (specifically emtrends() here) with brms measurement error models. I'm not totally sure if this is a brms or emmeans issue, but I'm opening it here since it seems to be associated with multivariate brms models with mi() terms.

The issue arises when trying to calculate the effect of some predictor variable that has measurement error, e.g. emtrends(fit_mi, ~ xobs, var = "xobs", resp = "y"). Some observations are that (1) the result of this function call changes between calls, (2) the estimated trend is often quite off from the corresponding slope coefficient, and (3) the reported lower/upper HPD have very large magnitudes relative to the estimated trend. Perhaps the first point has to do with the measurement error associated with xobs, but I'm not sure if/how this can be controlled when calculating emmeans/emtrends, and I'm really not sure what's going on with the other points.

Here's a reproducible example:

library(brms)
library(emmeans)

set.seed(253)

# Simulate data
n <- 100 # number of observations
a <- 0 # intercept
b <- 0.5 # slope
sigma <- 1 # response scale
xtrue <- runif(n, -3, 3) # true x
xsd <- abs(rnorm(n, 0, 0.5)) # measurement error
xobs <- rnorm(n, xtrue, xsd) # observed x
y <- rnorm(n, a + b*xtrue, sigma) # response
df <- data.frame(y, xobs, xsd)

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

# Specify formula for model without measurement error (for comparison)
form <- bf(y ~ xobs) + gaussian()

# Fit models
fit_mi <- brm(form_mi, df)
fit <- brm(form, df)

# emtrends for model with measurement error
emtrends(fit_mi, ~ xobs, var = "xobs", resp = "y") # results change per call
   xobs xobs.trend lower.HPD upper.HPD
 -0.013       1.41      -154       147

Point estimate displayed: median 
HPD interval probability: 0.95 
# emtrends for model without measurement error
emtrends(fit, ~ xobs, var = "xobs")
   xobs xobs.trend lower.HPD upper.HPD
 -0.013      0.441     0.326     0.555

Point estimate displayed: median 
HPD interval probability: 0.95 

Session info:

``` R version 4.3.1 (2023-06-16) Platform: aarch64-apple-darwin20 (64-bit) Running under: macOS Ventura 13.4 Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0 locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 time zone: America/Los_Angeles tzcode source: internal attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] emmeans_1.10.0 brms_2.20.4 Rcpp_1.0.11 loaded via a namespace (and not attached): [1] tidyselect_1.2.0 dplyr_1.1.3 farver_2.1.1 loo_2.7.0 fastmap_1.1.1 [6] tensorA_0.36.2 shinystan_2.6.0 promises_1.2.1 shinyjs_2.1.0 digest_0.6.33 [11] estimability_1.5 mime_0.12 lifecycle_1.0.4 StanHeaders_2.32.6 ellipsis_0.3.2 [16] processx_3.8.2 magrittr_2.0.3 posterior_1.5.0 compiler_4.3.1 rlang_1.1.3 [21] tools_4.3.1 igraph_1.5.1 utf8_1.2.4 labeling_0.4.3 prettyunits_1.2.0 [26] bridgesampling_1.1-2 htmlwidgets_1.6.2 curl_5.1.0 pkgbuild_1.4.2 plyr_1.8.9 [31] dygraphs_1.1.1.6 abind_1.4-5 miniUI_0.1.1.1 withr_2.5.1 grid_4.3.1 [36] stats4_4.3.1 fansi_1.0.6 xts_0.13.2 xtable_1.8-4 colorspace_2.1-0 [41] inline_0.3.19 ggplot2_3.5.0 scales_1.3.0 gtools_3.9.5 cli_3.6.2 [46] mvtnorm_1.2-4 crayon_1.5.2 generics_0.1.3 RcppParallel_5.1.7 rstudioapi_0.15.0 [51] reshape2_1.4.4 rstan_2.32.6 stringr_1.5.0 shinythemes_1.2.0 bayesplot_1.11.1 [56] parallel_4.3.1 matrixStats_1.0.0 base64enc_0.1-3 vctrs_0.6.5 V8_4.3.3 [61] Matrix_1.5-4.1 jsonlite_1.8.7 callr_3.7.3 crosstalk_1.2.1 glue_1.7.0 [66] codetools_0.2-19 ps_1.7.5 DT_0.32 distributional_0.3.2 stringi_1.7.12 [71] gtable_0.3.4 later_1.3.1 QuickJSR_1.1.3 munsell_0.5.0 tibble_3.2.1 [76] colourpicker_1.3.0 pillar_1.9.0 htmltools_0.5.6 Brobdingnag_1.2-9 R6_2.5.1 [81] shiny_1.7.5 lattice_0.21-8 markdown_1.9 backports_1.4.1 threejs_0.3.3 [86] httpuv_1.6.11 rstantools_2.4.0 coda_0.19-4.1 gridExtra_2.3 nlme_3.1-162 [91] checkmate_2.2.0 zoo_1.8-12 pkgconfig_2.0.3 ```
paul-buerkner commented 3 months ago

I think this question is better asked on stan discourse.

kthayashi commented 3 months ago

Thank you for the suggestion - I've now posted this question on the Stan Discourse: https://discourse.mc-stan.org/t/unexpected-behavior-of-emtrends-with-brms-models-with-me-vs-mi/34525