We are getting an intermittent error in calling add_criterion from loo with a brmsfit object fitted from a particular set of data using the cmdstanr backend and a custom beta-binomial family. We have successfully used similar code for fitting many other datasets using cmdstanr with no problem. In addition, it also appears to work fine when rstan is used as the back end instead of cmdstanr. It is this that makes us wonder if there is a bug?
The error is:
Error in validate_ll(log_ratios) : NAs not allowed in input.
The error is intermittent, so to reproduce it reliably we’ve put together a reprex using a while loop. The full reprex is below. Note that I can confirm the error is occurring on both mac and windows 10 machines.
library(brms)
library(cmdstanr)
library(extraDistr)
dat <- dat <- data.frame(log.x = c(-1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -1.00000000, -0.09691001, -0.09691001, -0.09691001, -0.09691001, -0.09691001, -0.09691001, -0.09691001, -0.09691001, -0.09691001, -0.09691001, 0.69019608, 0.69019608, 0.69019608, 0.69019608, 0.69019608, 0.69019608, 0.69019608, 0.69019608, 0.69019608, 0.69019608, 0.99122608, 0.99122608, 0.99122608, 0.99122608, 0.99122608, 0.99122608, 0.99122608, 0.99122608, 0.99122608, 0.99122608, 1.38738983, 1.38738983, 1.38738983, 1.38738983, 1.38738983, 1.38738983, 1.38738983, 1.38738983, 1.38738983, 1.38738983, 1.68930886, 1.68930886, 1.68930886, 1.68930886, 1.68930886, 1.68930886, 1.68930886, 1.68930886, 1.68930886, 1.91062441, 1.91062441, 1.91062441, 1.91062441, 1.91062441, 1.91062441, 1.91062441, 1.91062441, 1.91062441, 1.91062441), suc = c(20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 19, 20, 20, 20, 20, 15, 20, 20, 19, 20, 12, 20, 18, 17, 11, 19, 19, 15, 18, 11, 17, 14, 12, 18, 11, 17, 12, 12), tot = 20)
## Custom family Beta Binomial
beta_binomial2 <- brms::custom_family(
"beta_binomial2", dpars = c("mu", "phi"),
links = c("identity", "log"), lb = c(NA, 0),
type = "int", vars = "trials[n]"
)
stan_funs <- "
real beta_binomial2_lpmf(int y, real mu, real phi, int T) {
return beta_binomial_lpmf(y | T, mu * phi, (1 - mu) * phi);
}
int beta_binomial2_rng(real mu, real phi, int T) {
return beta_binomial_rng(T, mu * phi, (1 - mu) * phi);
}
"
my_stanvars <- brms::stanvar(scode = stan_funs, block = "functions")
beta_binomial2_lpmf <- function(y, mu, phi, trials) {
a <- mu * phi
b <- (1 - mu) * phi
extraDistr::dbbinom(y, trials, alpha = a, beta = b, log = TRUE)
}
posterior_epred_beta_binomial2 <- function(prep) {
mu <- prep$dpars$mu
trials <- prep$data$trials
trials <- matrix(trials, nrow = nrow(mu), ncol = ncol(mu),
byrow = TRUE)
mu * trials
}
beta_binomial2_rng <- function(mu, phi, trials) {
a <- mu * phi
b <- (1 - mu) * phi
out <- numeric(length = length(a))
for (i in seq_along(out)) {
out[i] <- extraDistr::rbbinom(1, trials, alpha = a[i], beta = b[i])
}
out
}
posterior_predict_beta_binomial2 <- function(i, prep, ...) {
mu <- prep$dpars$mu[, i]
phi <- prep$dpars$phi
trials <- prep$data$trials[i]
beta_binomial2_rng(mu, phi, trials)
}
log_lik_beta_binomial2 <- function(i, prep) {
mu <- prep$dpars$mu[, i]
phi <- prep$dpars$phi
trials <- prep$data$trials[i]
y <- prep$data$Y[i]
beta_binomial2_lpmf(y, mu, phi, trials)
}
my_formula <- brms::bf(suc | trials(tot) ~ top * exp(-exp(beta) *
(log.x - nec) * step(log.x - nec)),
top + beta + nec ~ 1, nl = TRUE)
my_priors <- prior("normal(0, 2)", nlpar = "beta") +
prior("normal(0.7, 2.7)", nlpar = "nec", lb = -1, ub = 1.91) +
prior("beta(5, 1)", nlpar = "top", lb = 0, ub = 1)
my_inits <- list(
list(b_beta = -0.10, b_top = 0.93, b_nec = 1.48),
list(b_beta = -1.51, b_top = 0.98, b_nec = -0.14),
list(b_beta = -1.05, b_top = 0.92, b_nec = 1.24),
list(b_beta = -3.46, b_top = 0.99, b_nec = 0.64)
)
set_cmdstan_path("C:/cmdstan")
#options(brms.backend = "cmdstanr")
n <- 0
check <- TRUE
while (check) {
n <- n + 1
print(n)
my_model <- brms::brm(my_formula, family = binomial,
prior = my_priors, chains = 4, cores = 4,
data = dat, inits = my_inits, backend = "cmdstanr",
stanvars = my_stanvars)
loo_pass <- try(add_criterion(my_model, "loo"))
check <- !inherits(loo_pass, "try-error")
}
sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)
Matrix products: default
locale:
[1] LC_COLLATE=English_Australia.1252 LC_CTYPE=English_Australia.1252 LC_MONETARY=English_Australia.1252
[4] LC_NUMERIC=C LC_TIME=English_Australia.1252
system code page: 65001
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] extraDistr_1.9.1 cmdstanr_0.4.0.9001 brms_2.16.3 Rcpp_1.0.8
loaded via a namespace (and not attached):
[1] nlme_3.1-155 matrixStats_0.61.0 xts_0.12.1 threejs_0.3.3 rstan_2.21.3
[6] tensorA_0.36.2 tools_4.1.2 backports_1.4.1 utf8_1.2.2 R6_2.5.1
[11] DT_0.21 DBI_1.1.2 colorspace_2.0-3 tidyselect_1.1.2 gridExtra_2.3
[16] prettyunits_1.1.1 processx_3.5.2 Brobdingnag_1.2-7 compiler_4.1.2 cli_3.2.0
[21] shinyjs_2.1.0 labeling_0.4.2 colourpicker_1.1.1 posterior_1.2.1 scales_1.1.1
[26] dygraphs_1.1.1.6 checkmate_2.0.0 mvtnorm_1.1-3 readr_2.1.2 ggridges_0.5.3
[31] callr_3.7.0 stringr_1.4.0 digest_0.6.29 StanHeaders_2.21.0-7 base64enc_0.1-3
[36] pkgconfig_2.0.3 htmltools_0.5.2 fastmap_1.1.0 htmlwidgets_1.5.4 rlang_1.0.2
[41] shiny_1.7.1 farver_2.1.0 generics_0.1.2 jsonlite_1.8.0 zoo_1.8-9
[46] crosstalk_1.2.0 gtools_3.9.2 dplyr_1.0.8 distributional_0.3.0 inline_0.3.19
[51] magrittr_2.0.2 bayesnec_2.0.2.3 loo_2.4.1 bayesplot_1.8.1 Matrix_1.4-0
[56] munsell_0.5.0 fansi_1.0.2 abind_1.4-5 lifecycle_1.0.1 stringi_1.7.6
[61] pkgbuild_1.3.1 plyr_1.8.6 grid_4.1.2 formula.tools_1.7.1 parallel_4.1.2
[66] promises_1.2.0.1 crayon_1.5.0 miniUI_0.1.1.1 lattice_0.20-45 hms_1.1.1
[71] knitr_1.37 ps_1.6.0 pillar_1.7.0 igraph_1.2.11 markdown_1.1
[76] shinystan_2.6.0 reshape2_1.4.4 codetools_0.2-18 stats4_4.1.2 rstantools_2.1.1
[81] glue_1.6.2 evaluate_0.15 data.table_1.14.2 RcppParallel_5.1.5 operator.tools_1.6.3
[86] vctrs_0.3.8 tzdb_0.2.0 httpuv_1.6.5 gtable_0.3.0 purrr_0.3.4
[91] tidyr_1.2.0 assertthat_0.2.1 ggplot2_3.3.5 xfun_0.30 mime_0.12
[96] xtable_1.8-4 coda_0.19-4 later_1.3.0 diffobj_0.3.5 tibble_3.1.6
[101] shinythemes_1.2.0 ellipsis_0.3.2 bridgesampling_1.1-2