loo moment_match crashes R if save_all_pars not specified

[n-kall]

As mentioned in #1126 loo with moment matching doesn't work without save_all_pars=save_pars(all = TRUE). But it seems worse than just not having an appropriate warning message, as it can crash R if tried without the parameter.

Example: With save_all_pars

 library(brms)
 m <- brm(yield ~ N*P*K, npk, save_all_pars = save_pars(all = TRUE))
 loo(m, moment_match = TRUE)

Yields:

Computed from 4000 by 24 log-likelihood matrix

         Estimate  SE
elpd_loo    -81.2 3.1
p_loo         7.5 1.3
looic       162.4 6.3
------
Monte Carlo SE of elpd_loo is 0.1.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     20    83.3%   377       
 (0.5, 0.7]   (ok)        4    16.7%   538       
   (0.7, 1]   (bad)       0     0.0%   <NA>      
   (1, Inf)   (very bad)  0     0.0%   <NA>      

All Pareto k estimates are ok (k < 0.7).
See help('pareto-k-diagnostic') for details.

Whereas,

 library(brms)
 m <- brm(yield ~ N*P*K, npk)
 loo(m, moment_match = TRUE)

Yields:

 *** caught segfault ***
address (nil), cause 'memory not mapped'

Traceback:
 1: .External(list(name = "CppMethod__invoke_notvoid", address = <pointer: 0x564515d9c940>,     dll = list(name = "Rcpp", path = "/lib/R/library/Rcpp/libs/Rcpp.so",         dynamicLookup = TRUE, handle = <pointer: 0x564515fcbf70>,         info = <pointer: 0x564513e69b30>), numParameters = -1L),     <\
pointer: 0x564524ed9590>, <pointer: 0x56451dc4f180>, .pointer,     ...)
 2: object@.MISC$stan_fit_instance$unconstrain_pars(pars)
 3: .local(object, ...)
 4: rstan::unconstrain_pars(x, pars = .rstan_relist(theta, skeleton))
 5: rstan::unconstrain_pars(x, pars = .rstan_relist(theta, skeleton))
 6: FUN(newX[, i], ...)
 7: apply(pars, 1, FUN = function(theta) {    rstan::unconstrain_pars(x, pars = .rstan_relist(theta, skeleton))})
 8: unconstrain_pars_stanfit(x$fit, pars = pars, ...)
 9: unconstrain_pars(x, pars = pars, ...)
10: loo::loo_moment_match.default(x, loo = loo, post_draws = as.matrix,     log_lik_i = .log_lik_i, unconstrain_pars = .unconstrain_pars,     log_prob_upars = .log_prob_upars, log_lik_i_upars = .log_lik_i_upars,     k_threshold = k_threshold, newdata = newdata, resp = resp,     ...)
11: doTryCatch(return(expr), name, parentenv, handler)
12: tryCatchOne(expr, names, parentenv, handlers[[1L]])
13: tryCatchList(expr, classes, parentenv, handlers)
14: tryCatch(expr, error = function(e) {    call <- conditionCall(e)    if (!is.null(call)) {        if (identical(call[[1L]], quote(doTryCatch)))             call <- sys.call(-4L)        dcall <- deparse(call)[1L]        prefix <- paste("Error in", dcall, ": ")        LONG <- 75L        sm <- strsplit(conditionMessage(e), "\n")[[1L]] \
       w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w")        if (is.na(w))             w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L],                 type = "b")        if (w > LONG)             prefix <- paste0(prefix, "\n  ")    }    else prefix <- "Error : "    msg <- paste0(prefix, conditionMessage(e), "\n")  \
  .Internal(seterrmessage(msg[1L]))    if (!silent && isTRUE(getOption("show.error.messages"))) {        cat(msg, file = outFile)        .Internal(printDeferredWarnings())    }    invisible(structure(msg, class = "try-error", condition = e))})
15: try(loo::loo_moment_match.default(x, loo = loo, post_draws = as.matrix,     log_lik_i = .log_lik_i, unconstrain_pars = .unconstrain_pars,     log_prob_upars = .log_prob_upars, log_lik_i_upars = .log_lik_i_upars,     k_threshold = k_threshold, newdata = newdata, resp = resp,     ...))
16: loo_moment_match.brmsfit(x = .x1, loo = .x2, newdata = .x3, resp = .x4,     k_threshold = .x5, check = .x6)
17: loo_moment_match(x = .x1, loo = .x2, newdata = .x3, resp = .x4,     k_threshold = .x5, check = .x6)
18: eval(expr, envir, ...)
19: eval(expr, envir, ...)
20: eval2(call, envir = args, enclos = envir)
21: do_call("loo_moment_match", moment_match_args)
22: .loo(x = .x1, newdata = .x2, resp = .x3, model_name = .x4, pointwise = .x5,     k_threshold = .x6, save_psis = .x7, moment_match = .x8, reloo = .x9,     moment_match_args = .x10, reloo_args = .x11)
23: eval(expr, envir, ...)
24: eval(expr, envir, ...)
25: eval2(call, envir = args, enclos = envir)
26: do_call(paste0(".", criterion), args)
27: .fun(criterion = .x1, pointwise = .x2, resp = .x3, k_threshold = .x4,     save_psis = .x5, moment_match = .x6, reloo = .x7, moment_match_args = .x8,     reloo_args = .x9, x = .x10, model_name = .x11, use_stored = .x12)
28: eval(expr, envir, ...)
29: eval(expr, envir, ...)
30: eval2(call, envir = args, enclos = envir)
31: do_call(compute_loo, args)
32: .fun(models = .x1, criterion = .x2, pointwise = .x3, compare = .x4,     resp = .x5, k_threshold = .x6, save_psis = .x7, moment_match = .x8,     reloo = .x9, moment_match_args = .x10, reloo_args = .x11)
33: eval(expr, envir, ...)
34: eval(expr, envir, ...)
35: eval2(call, envir = args, enclos = envir)
36: do_call(compute_loolist, args)
37: loo.brmsfit(m, moment_match = TRUE)
38: loo(m, moment_match = TRUE)
An irrecoverable exception occurred. R is aborting now ...

Session Info:

 R version 4.0.5 (2021-03-31)
 Platform: x86_64-pc-linux-gnu (64-bit)
 Running under: Ubuntu 18.04.5 LTS

 Matrix products: default
 BLAS:   /lib/R/lib/libRblas.so
 LAPACK: /lib/R/lib/libRlapack.so

 locale:
  [1] LC_CTYPE=en_GB.utf8       LC_NUMERIC=C
  [3] LC_TIME=en_DK.utf8        LC_COLLATE=en_GB.utf8
  [5] LC_MONETARY=en_GB.utf8    LC_MESSAGES=en_GB.utf8
  [7] LC_PAPER=fi_FI.utf8       LC_NAME=C
  [9] LC_ADDRESS=C              LC_TELEPHONE=C
 [11] LC_MEASUREMENT=en_GB.utf8 LC_IDENTIFICATION=C

 attached base packages:
 [1] stats     graphics  grDevices utils     datasets  methods   base

 other attached packages:
 [1] brms_2.15.5 Rcpp_1.0.6

 loaded via a namespace (and not attached):
  [1] nlme_3.1-152         matrixStats_0.58.0   xts_0.12.1
  [4] threejs_0.3.3        rstan_2.26.1         backports_1.2.1
  [7] tools_4.0.5          utf8_1.2.1           R6_2.5.0
 [10] DT_0.18              DBI_1.1.1            mgcv_1.8-35
 [13] projpred_2.0.2       colorspace_2.0-1     tidyselect_1.1.1
 [16] gridExtra_2.3        prettyunits_1.1.1    processx_3.5.2
 [19] Brobdingnag_1.2-6    emmeans_1.6.0        curl_4.3.1
 [22] compiler_4.0.5       cli_2.5.0            shinyjs_2.0.0
 [25] colourpicker_1.1.0   scales_1.1.1         dygraphs_1.1.1.6
 [28] mvtnorm_1.1-1        ggridges_0.5.3       callr_3.7.0
 [31] stringr_1.4.0        digest_0.6.27        StanHeaders_2.26.1
 [34] minqa_1.2.4          base64enc_0.1-3      pkgconfig_2.0.3
 [37] htmltools_0.5.1.1    lme4_1.1-27          fastmap_1.1.0
 [40] htmlwidgets_1.5.3    rlang_0.4.11         shiny_1.6.0
 [43] generics_0.1.0       zoo_1.8-9            jsonlite_1.7.2
 [46] crosstalk_1.1.1      gtools_3.8.2         dplyr_1.0.6
 [49] inline_0.3.18        magrittr_2.0.1       loo_2.4.1
 [52] bayesplot_1.8.0      Matrix_1.3-2         munsell_0.5.0
 [55] fansi_0.4.2          abind_1.4-5          lifecycle_1.0.0
 [58] stringi_1.6.2        MASS_7.3-53.1        pkgbuild_1.2.0
 [61] plyr_1.8.6           grid_4.0.5           parallel_4.0.5
 [64] promises_1.2.0.1     crayon_1.4.1         miniUI_0.1.1.1
 [67] lattice_0.20-41      splines_4.0.5        ps_1.6.0
 [70] pillar_1.6.1         igraph_1.2.6         boot_1.3-27
 [73] estimability_1.3     markdown_1.1         shinystan_2.5.0
 [76] codetools_0.2-18     reshape2_1.4.4       stats4_4.0.5
 [79] rstantools_2.1.1     glue_1.4.2           V8_3.4.2
 [82] RcppParallel_5.1.4   vctrs_0.3.8          nloptr_1.2.2.2
 [85] httpuv_1.6.1         gtable_0.3.0         purrr_0.3.4
 [88] assertthat_0.2.1     ggplot2_3.3.3        mime_0.10
 [91] xtable_1.8-4         coda_0.19-4          later_1.2.0
 [94] rsconnect_0.8.17     tibble_3.1.2         shinythemes_1.2.0
 [97] gamm4_0.2-6          ellipsis_0.3.2       bridgesampling_1.1-2

[paul-buerkner]

Thanks for opening this issue. A segfault doesn't sounds that great and should probably be fixed rather from the rstan side, as there seems to be a problem in rstan::unconstrain_pars. Would you mind opening an issue for rstan as well?

[n-kall]

Filed an issue for rstan, stan-dev/rstan#939

[paul-buerkner]

Thanks! So I can add a check for save_pars(all = TRUE) in brms but this may uncessarily stop some models from running that would otherwise have run, because those models do not actually have any parameters that were not saved. So I am not sure there is anything reasonable to do from the brms side, as rstan should capture the problem and throw an informative error that than brms can catch and enrich with the please use save_pars(all = TRUE) info.

[n-kall]

I agree, it seems there's not much to be fixed on the brms side. Also, I think it's the case that the bug is not present in rstan on cran (2.21.1), but is on the latest development version (2.26.1), so it's unlikely this will cause issues for users at the moment. So I think you can probably mark this as closed.

[paul-buerkner]

Ok, thank you!