Report the Monte Carlo SE of the mean from the draws themselves

NEWS.md

@@ -202,6 +202,11 @@
 * Improvement: The binary, continuous and count IPD likelihoods use Stan's
   fused GLM densities on their canonical links; results agree with 0.1.0 to
   Monte Carlo error.
+* Fix: With the cmdstanr engine, `fit$summary$se_mean` is the Monte Carlo
+  standard error of the posterior mean from `posterior::mcse_mean()`. It was
+  the posterior SD over the square root of the bulk ESS, which is computed on
+  rank-normalized draws and misstates the error for skewed quantities such as
+  a ratio.
 * Fix: `mlumr()` refuses a normal fit whose outcome is constant or reproduced
   exactly by its covariates, where the posterior for the residual SD is
   improper, and warns when the design is saturated.

R/backend_cmdstanr.R

@@ -128,7 +128,9 @@ fit_cmdstanr <- function(model_name, stan_data, chains, iter, warmup,
     cmdstan_summ <- fit$summary(
       variables = NULL,
       mean = mean,
-      se_mean = function(.x) stats::sd(.x) / sqrt(posterior::ess_bulk(.x)),
+      # The Monte Carlo SE of the mean needs the ESS of the draws themselves;
+      # bulk ESS is computed on rank-normalized draws, a different quantity.
+      se_mean = posterior::mcse_mean,
       sd = stats::sd,
       `2.5%` = function(.x) stats::quantile(.x, 0.025),
       `25%` = function(.x) stats::quantile(.x, 0.25),

tests/testthat/test-engine.R

@@ -75,4 +75,14 @@ test_that("cmdstanr backend fits a model end-to-end", {
   expect_true(is.numeric(fit$diagnostics$n_divergent))
   expect_true(is.numeric(fit$diagnostics$n_max_treedepth))
   expect_false(file.exists(file.path("inst", "stan", "mlumr_binary_spfa")))
+
+  # se_mean is the Monte Carlo SE of the mean of the draws, chain by chain;
+  # a risk ratio is skewed, where a rank-normalized ESS would misstate it.
+  vars <- c("mu_index", "rr_index")
+  arr <- fit$stanfit$draws(variables = vars)
+  expected <- vapply(vars, function(v) {
+    posterior::mcse_mean(posterior::extract_variable_matrix(arr, v))
+  }, numeric(1))
+  got <- fit$summary$se_mean[match(vars, fit$summary$variable)]
+  expect_equal(got, unname(expected))
 })