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.28k stars 184 forks source link

Error in prepare predictions, non-numeric argument to pnorm #1467

Closed jflournoy closed 1 year ago

jflournoy commented 1 year ago

I have a model with a predictor modeled shifted lognormal with measurement error that samples fine using rstan.

When I ask for pp_check or bayes_R2 I get the following error that seems to come from prepare_predictions:

Error in pnorm(.x1, mean = .x2, sd = .x3) : 
  Non-numeric argument to mathematical function

Full traceback, below.

Reproducible example:

library(brms)

b_data <- data.frame(y = rnorm(500),
                     x_se = rgamma(500, 1, 1))
b_data$x <- brms::rshifted_lnorm(500, meanlog = 0, sdlog = b_data$x_se, shift = .5)
b_formula <- bf(y ~ 1 + mi(x), family = gaussian()) + 
  bf(x | mi(sdy = x_se) ~ 1, family = shifted_lognormal()) + 
  set_rescor(rescor = FALSE)

b_standata <- brms::make_standata(b_formula, data = b_data)
b_standata$lbmi_x
b_fit <- brm(formula = b_formula, data = b_data,
             cores = 4, chains = 4, iter = 2000, warmup = 1000)

bayes_R2(b_fit)
traceback()
pp_check(b_fit, resp = 'y')
traceback()

packageVersion('brms')
packageVersion('rstan')
> packageVersion('brms')
[1] ‘2.17.7’
> packageVersion('rstan')
[1] ‘2.26.13’

Traceback:

> bayes_R2(b_fit)
Error in pnorm(.x1, mean = .x2, sd = .x3) : 
  Non-numeric argument to mathematical function
> traceback()
23: pnorm(.x1, mean = .x2, sd = .x3)
22: eval(expr, envir, ...)
21: eval(expr, envir, ...)
20: eval2(call, envir = args, enclos = envir)
19: do_call(pdist, c(list(lb), args))
18: rcontinuous(n = prod(dim), dist = "norm", mean = Y, sd = sdy, 
        lb = sdata[[paste0("lbmi_", vmi)]], ub = sdata[[paste0("ubmi_", 
            vmi)]])
17: prepare_predictions_sp(x, draws, sdata, ...)
16: prepare_predictions.btl(x$dpars[[dp]], draws = draws, sdata = sdata, 
        data = data, ...)
15: prepare_predictions(x$dpars[[dp]], draws = draws, sdata = sdata, 
        data = data, ...)
14: prepare_predictions.brmsterms(x$terms[[resp]], draws = draws, 
        sdata = sdata, ...)
13: prepare_predictions(x$terms[[resp]], draws = draws, sdata = sdata, 
        ...)
12: prepare_predictions.mvbrmsterms(bterms, draws = draws, sdata = sdata, 
        data = x$data, prep_ranef = prep_ranef, meef = meef, resp = resp, 
        sample_new_levels = sample_new_levels, nug = nug, new = new, 
        oos = oos, stanvars = x$stanvars)
11: prepare_predictions(bterms, draws = draws, sdata = sdata, data = x$data, 
        prep_ranef = prep_ranef, meef = meef, resp = resp, sample_new_levels = sample_new_levels, 
        nug = nug, new = new, oos = oos, stanvars = x$stanvars)
10: prepare_predictions.brmsfit(object, newdata = newdata, re_formula = re_formula, 
        resp = resp, ndraws = ndraws, draw_ids = draw_ids, check_response = FALSE, 
        ...)
9: prepare_predictions(object, newdata = newdata, re_formula = re_formula, 
       resp = resp, ndraws = ndraws, draw_ids = draw_ids, check_response = FALSE, 
       ...)
8: posterior_epred.brmsfit(.x1, sort = .x2, resp = .x3)
7: .fun(.x1, sort = .x2, resp = .x3)
6: eval(expr, envir, ...)
5: eval(expr, envir, ...)
4: eval2(call, envir = args, enclos = envir)
3: do_call(posterior_epred, args_ypred)
2: bayes_R2.brmsfit(b_fit)
1: bayes_R2(b_fit)
> pp_check(b_fit, resp = 'y')
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Error in pnorm(.x1, mean = .x2, sd = .x3) : 
  Non-numeric argument to mathematical function
> traceback()
23: pnorm(.x1, mean = .x2, sd = .x3)
22: eval(expr, envir, ...)
21: eval(expr, envir, ...)
20: eval2(call, envir = args, enclos = envir)
19: do_call(pdist, c(list(lb), args))
18: rcontinuous(n = prod(dim), dist = "norm", mean = Y, sd = sdy, 
        lb = sdata[[paste0("lbmi_", vmi)]], ub = sdata[[paste0("ubmi_", 
            vmi)]])
17: prepare_predictions_sp(x, draws, sdata, ...)
16: prepare_predictions.btl(x$dpars[[dp]], draws = draws, sdata = sdata, 
        data = data, ...)
15: prepare_predictions(x$dpars[[dp]], draws = draws, sdata = sdata, 
        data = data, ...)
14: prepare_predictions.brmsterms(x$terms[[resp]], draws = draws, 
        sdata = sdata, ...)
13: prepare_predictions(x$terms[[resp]], draws = draws, sdata = sdata, 
        ...)
12: prepare_predictions.mvbrmsterms(bterms, draws = draws, sdata = sdata, 
        data = x$data, prep_ranef = prep_ranef, meef = meef, resp = resp, 
        sample_new_levels = sample_new_levels, nug = nug, new = new, 
        oos = oos, stanvars = x$stanvars)
11: prepare_predictions(bterms, draws = draws, sdata = sdata, data = x$data, 
        prep_ranef = prep_ranef, meef = meef, resp = resp, sample_new_levels = sample_new_levels, 
        nug = nug, new = new, oos = oos, stanvars = x$stanvars)
10: prepare_predictions.brmsfit(object, newdata = newdata, re_formula = re_formula, 
        resp = resp, ndraws = ndraws, draw_ids = draw_ids, check_response = FALSE, 
        ...)
9: prepare_predictions(object, newdata = newdata, re_formula = re_formula, 
       resp = resp, ndraws = ndraws, draw_ids = draw_ids, check_response = FALSE, 
       ...)
8: posterior_predict.brmsfit(.x1, newdata = .x2, resp = .x3, draw_ids = .x4)
7: posterior_predict(.x1, newdata = .x2, resp = .x3, draw_ids = .x4)
6: eval(expr, envir, ...)
5: eval(expr, envir, ...)
4: eval2(call, envir = args, enclos = envir)
3: do_call(method, pred_args)
2: pp_check.brmsfit(b_fit, resp = "y")
1: pp_check(b_fit, resp = "y")
paul-buerkner commented 1 year ago

Thanks for opening this issue! Is this problem still present in the latest github version of brms?

jflournoy commented 1 year ago

Yep, same behavior after installing from github

> bayes_R2(b_fit)
Error in pnorm(.x1, mean = .x2, sd = .x3) : 
  Non-numeric argument to mathematical function
> packageVersion('brms')
[1] ‘2.18.8
paul-buerkner commented 1 year ago

Thanks. Will try to fix this in the next brms release.

paul-buerkner commented 1 year ago

I think this issue was the same as for the other one that shifted_lognormal is not compatible with missing values. I am sorry that I don't have a fix for this issue in the sense that these two features work well together.