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

Exgaussian() fitting, bug in pp_check() #1046

Closed rosskempner closed 3 years ago

rosskempner commented 3 years ago

Hello, I am fitting an exgaussian() model and experiencing an issue when it comes to using pp_check() on the brmsfit.

formula_exgaussian1 = bf(formula = reaction_time | trunc(ub = 550) ~ 
                        1 + 
                        all_random_mode +
                        bonus_for_speed + 
                        all_random_mode * bonus_for_speed +
                        block_num + 
                        block_num * all_random_mode + 
                        block_num * bonus_for_speed + 
                        block_num * bonus_for_speed * all_random_mode +
                        bigram_ideal_first_surprisal  +  
                        bigram_ideal_first_surprisal * all_random_mode + 
                        bigram_ideal_first_surprisal * bonus_for_speed + 
                        bigram_ideal_first_surprisal * bonus_for_speed * all_random_mode +
                        finger_transition +
                        (1 + block_num + bigram_ideal_first_surprisal | participant_id)
                      )

get_prior(formula = formula_exgaussian1, data = experiment_df, family = exgaussian())
priors_exgaussian1 = c(
    set_prior("normal(400,0.3)", class = "Intercept"),
    set_prior("normal(0.5,0.2)", class = "sigma"),
    set_prior("normal(-0.01,0.01)", class = "b", coef = "block_num"),
    set_prior("normal(1,0.4)", class = "b", coef = "bigram_ideal_first_surprisal"),
    set_prior("normal(0.1,0.1)", class = "sd", coef = "bigram_ideal_first_surprisal", group = "participant_id"),
    set_prior("normal(0,0.01)", class = "sd", coef = "block_num", group = "participant_id"),
    set_prior("normal(0.4,1)", class = "sd", coef = "Intercept", group = "participant_id"),
    set_prior("normal(0,0.2)", class = "b", coef = "all_random_mode"),
    set_prior("normal(0,0.2)", class = "b", coef = "bonus_for_speed"),
    set_prior("normal(0,0.1)", class = "b", coef = "all_random_mode:bonus_for_speed"),
    set_prior("normal(0,0.0001)", class = "b", coef = "all_random_mode:block_num"),
    set_prior("normal(0,0.0001)", class = "b", coef = "bonus_for_speed:block_num"),
    set_prior("normal(0,0.0001)", class = "b", coef = "all_random_mode:bonus_for_speed:block_num"),
    set_prior("normal(0,0.0001)", class = "b", coef = "all_random_mode:bigram_ideal_first_surprisal"),
    set_prior("normal(0,0.001)", class = "b", coef = "bonus_for_speed:bigram_ideal_first_surprisal"),
    set_prior("normal(0,0.001)", class = "b", coef = "all_random_mode:bonus_for_speed:bigram_ideal_first_surprisal"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitiondj"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitiondk"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitionfd"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitionfj"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitionfk"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitionjd"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitionjf"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitionjk"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitionkd"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitionkf"),
    set_prior("normal(0,0.1)", class = "b", coef = "finger_transitionkj"),
    set_prior("lkj(2)", class = "cor")
)        
fit_exgaussian1 <- brm(formula = formula_exgaussian1,
                        data = experiment_df,
                        family = exgaussian(),
                        prior = priors_exgaussian1,
  sample_prior = "only")
pp_check(fit_exgaussian1)

The fitting process works, plotting works and summary on the brms fit object works but an error comes at pp_check, Below is the error message :

Error in qexgaussian(.x1, mu = .x2, sigma = .x3, beta = .x4) : could not find function "qexgaussian"
18.
eval(expr, envir, ...)
17.
eval(expr, envir, ...)
16.
eval2(call, envir = args, enclos = parent.frame())
15.
do_call(qdist, c(list(out), args))
14.
rcontinuous(n = prep$nsamples, dist = "exgaussian", mu = get_dpar(prep, "mu", i = i), sigma = get_dpar(prep, "sigma", i = i), beta = get_dpar(prep, "beta", i = i), lb = prep$data$lb[i], ub = prep$data$ub[i])
13.
FUN(X[[i]], ...)
12.
lapply(X, FUN, ...)
11.
plapply(seq_len(N), pp_fun, cores = cores, prep = object, ...)
10.
posterior_predict.brmsprep(prep, transform = transform, sort = sort, ntrys = ntrys, negative_rt = negative_rt, cores = cores, summary = FALSE)
9.
posterior_predict(prep, transform = transform, sort = sort, ntrys = ntrys, negative_rt = negative_rt, cores = cores, summary = FALSE)
8.
posterior_predict.brmsfit(.x1, newdata = .x2, resp = .x3, subset = .x4)
7.
posterior_predict(.x1, newdata = .x2, resp = .x3, subset = .x4)
6.
eval(expr, envir, ...)
5.
eval(expr, envir, ...)
4.
eval2(call, envir = args, enclos = parent.frame())
3.
do_call(method, pred_args)
2.
pp_check.brmsfit(fit_exgaussian1)
1.
pp_check(fit_exgaussian1)

Operating System: MAC OS Catalina brms Version: 2.14.0++

rosskempner commented 3 years ago

hi @paul-buerkner, do you need more information for this post to be helpful? If so please let me know.

Also, please note that I sent this over to the bayesplot issue forum and they sent me to you, since they believe this was an issue in brms before it gets to bayesplot.

paul-buerkner commented 3 years ago

Thank you! I get what the problem is. The error is not a bug but the error message is not helpful. Basically, it should say, that posterior predictions for this model cannot be evaluated.

rosskempner commented 3 years ago

Got it thanks! might that be due to NA values due to too wide parameters or values less than zero.

What are some things we need to watch out for with the exguassian()?

paul-buerkner commented 3 years ago

The reason is simply that the quantile function of the exgaussian is non-analytical (to my current knowledge) and hence we have no qexgaussian function, which is required for truncation evaluation.

rosskempner commented 3 years ago

This makes sense, since without truncation pp_check() works. Does this imply that we can not use truncation feature with exgaussian()?

paul-buerkner commented 3 years ago

Just not posterior_predict in the truncated exgaussian case.

Am Mo., 23. Nov. 2020 um 14:27 Uhr schrieb rosskempner < notifications@github.com>:

This makes sense, since without truncation pp_check() works. Does this imply that we can not use truncation feature with exgaussian()?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/paul-buerkner/brms/issues/1046#issuecomment-732161191, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADCW2AB4OT356RBPCFIYM6TSRJPL3ANCNFSM4T6ED4TQ .

rosskempner commented 3 years ago

thanks!

paul-buerkner commented 3 years ago

posterior_predict with exgaussian should now work as brms will now use (slow but working) rejection sampling when the quantile function is unavailable.