Problems specifying bounds on Wiener model parameter in brms model

Interfaces brms
1

I am trying to model accuracies and reaction times in psychophysical data using the Wiener diffusion model in brms. see the excellent tutorial by Henrik Singmann. The Wiener family has four distributional parameters: the drift rate, boundary separation, bias and non-decision time.

I am having trouble sampling from the model (both for sampling convergence and even chain initialization). The problem seems to arise from the non-decision time parameter: the NDT parameter creates a hard lower bound for the lpdf output such that the random variable output by the other model parameters must be larger than the NDT. If this is not the case the family function throws an error. NDT must also be positive.

This last part is handled by using a log link function for that parameter in the family by default. The problem is that very large NDTs can be sampled in chain initialization and exploration, leading to many rejected samples and problems initializing.

Since I don’t have any theoretical interest in the NDT, I would like to fix or strongly bound the parameter to facilitate sampling in the rest of the model, but I’m unsure how to achieve this.

A simplified version of the model with problems is shown below:

model_formula <- brmsformula(rt | dec(response_yes) ~ factor_a +
                               (factor_a | p | subj) +
                               (1 | im_filename),
                             bs ~ factor_b +
                               (factor_b | p | subj),
                             ndt ~ 1,
                             bias ~ factor_b +
                               (factor_b | p | subj),
                             family = wiener(link_bs = "log",
                                             link_ndt = "log",
                                             link_bias = "logit")
                             )

I would like to be able to do something like:

ndt = 0.5

(i.e. fixing the ndt to a single value)

or

ndt ~ inv_logit(factor_b)*2.0

(i.e. bounding the NDT parameter in the model to be within 0–2 by using an inverse logit link).

Neither of these attempts works. Can anyone more familiar with brms or this distributional model suggest how I could achieve these bounds or fixed parameter?

Thanks!

sessionInfo()
[I’m running this in a docker container]:
R version 3.5.2 (2018-12-20)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS: /usr/lib/openblas-base/libblas.so.3
LAPACK: /usr/lib/libopenblasp-r0.2.19.so

locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8
[4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

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

other attached packages:
[1] bindrcpp_0.2.2 foreach_1.4.4 here_0.1 forcats_0.3.0 stringr_1.3.1 dplyr_0.7.8
[7] purrr_0.2.5 readr_1.3.0 tidyr_0.8.2 tibble_1.4.2 tidyverse_1.2.1 brms_2.7.0
[13] ggplot2_3.1.0 Rcpp_1.0.0

loaded via a namespace (and not attached):
[1] nlme_3.1-137 matrixStats_0.54.0 xts_0.11-2 lubridate_1.7.4
[5] threejs_0.3.1 httr_1.4.0 rprojroot_1.3-2 rstan_2.18.2
[9] tools_3.5.2 backports_1.1.3 R6_2.3.0 DT_0.5
[13] lazyeval_0.2.1 colorspace_1.3-2 withr_2.1.2 tidyselect_0.2.5
[17] gridExtra_2.3 prettyunits_1.0.2 processx_3.2.1 Brobdingnag_1.2-6
[21] compiler_3.5.2 rvest_0.3.2 cli_1.0.1 xml2_1.2.0
[25] shinyjs_1.0 colourpicker_1.0 scales_1.0.0 dygraphs_1.1.1.6
[29] mvtnorm_1.0-8 ggridges_0.5.1 callr_3.1.0 digest_0.6.18
[33] StanHeaders_2.18.0-1 base64enc_0.1-3 pkgconfig_2.0.2 htmltools_0.3.6
[37] readxl_1.2.0 htmlwidgets_1.3 rlang_0.3.0.1 rstudioapi_0.8
[41] shiny_1.2.0 bindr_0.1.1 generics_0.0.2 zoo_1.8-4
[45] jsonlite_1.6 crosstalk_1.0.0 gtools_3.8.1 inline_0.3.15
[49] magrittr_1.5 loo_2.0.0 bayesplot_1.6.0 Matrix_1.2-15
[53] munsell_0.5.0 abind_1.4-5 stringi_1.2.4 yaml_2.2.0
[57] pkgbuild_1.0.2 plyr_1.8.4 grid_3.5.2 parallel_3.5.2
[61] promises_1.0.1 crayon_1.3.4 miniUI_0.1.1.1 lattice_0.20-38
[65] haven_2.0.0 hms_0.4.2 ps_1.2.1 pillar_1.3.1
[69] igraph_1.2.2 markdown_0.9 shinystan_2.5.0 codetools_0.2-15
[73] reshape2_1.4.3 stats4_3.5.2 rstantools_1.5.1 glue_1.3.0
[77] modelr_0.1.2 httpuv_1.4.5.1 cellranger_1.1.0 gtable_0.2.0
[81] assertthat_0.2.0 mime_0.6 xtable_1.8-3 broom_0.5.1
[85] coda_0.19-2 later_0.7.5 rsconnect_0.8.12 iterators_1.0.10
[89] shinythemes_1.1.2 bridgesampling_0.6-0

2

Replacing ndt ~ 1 with ndt = 0.5 should be working. What is the error message you receive?

You can also ndt non-linearily by replacing ndt ~ 1 with bf(ndt ~ inv_logit(factor_b)*2.0).

3

Hi Paul, thanks for the quick response. If I replace ndt ~ 1 with ndt = 0.5 I get the following error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

variable "min_Y" does not exist.
  error in 'model12975707253_file129440de5dd' at line 32, column 27
  -------------------------------------------------
    30:   int<lower=1> K_bs;  // number of population-level effects 
    31:   matrix[N, K_bs] X_bs;  // population-level design matrix 
    32:   real<lower=0,upper=min_Y> ndt;  // non-decision time parameter 
                                  ^
    33:   int<lower=1> K_bias;  // number of population-level effects 
  -------------------------------------------------

Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname,  : 
  failed to parse Stan model 'file129440de5dd' due to the above error.
SYNTAX ERROR, MESSAGE(S) FROM PARSER:

variable "min_Y" does not exist.
  error in 'model12945d21090_file129440de5dd' at line 32, column 27
  -------------------------------------------------
    30:   int<lower=1> K_bs;  // number of population-level effects 
    31:   matrix[N, K_bs] X_bs;  // population-level design matrix 
    32:   real<lower=0,upper=min_Y> ndt;  // non-decision time parameter 
                                  ^
    33:   int<lower=1> K_bias;  // number of population-level effects 
  -------------------------------------------------

Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname,  : 
  failed to parse Stan model 'file129440de5dd' due to the above error.

If I try replacing ndt ~ 1 with bf(ndt ~ inv_logit(instruction_condition) * 2.0, nl = TRUE) I get the following error:

Error: Expecting a single value when fixing parameters.

4

The first one was a minor bug I just fixed on github, that is ndt = 0.5 should now work.

The second was my own stupidity. Instead of

bf(ndt ~ inv_logit(instruction_condition) * 2.0)

it should be

nlf(ndt ~ inv_logit(instruction_condition) * 2.0)
1 Like
5

For

nlf(ndt ~ inv_logit(instruction_condition) * 2.0)

I get the error

Error: No non-linear parameters specified.

when calling get_prior.

6

This is because your non-linear formula contains no parameters (just data), which is not allowed and probably
not what you want. I guess you want one parameter per instruction_condition (or their contrasts), right? In this case, add

nlf(ndt ~ inv_logit(bndt) * 2.0)
bndt ~ instruction_condition

to your set of formulas (the name bndt is arbitrary of course).

7

Ah yes of course, sorry Paul. I needed to mentally jump up the hierarchy!