<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/diffusion-wiener-model-analysis-with-brms-part-iii-hypothesis-tests-of-parameter-estimates/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://singmann.org/wiener-model-analysis-with-brms-part-iii/"> Computational Psychology - Henrik Singmann</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This is the third part of my blog series on fitting the 4-parameter Wiener model with <a rel="nofollow" href="https://cran.r-project.org/package=brms" target="_blank"><code>brms</code></a>. The <a rel="nofollow" href="http://singmann.org/wiener-model-analysis-with-brms-part-i/" target="_blank">first part</a> discussed how to set up the data and model. The <a rel="nofollow" href="http://singmann.org/wiener-model-analysis-with-brms-part-ii/" target="_blank">second part</a> was concerned with (mostly graphical) model diagnostics and the assessment of the adequacy (i.e., the fit) of the model. This third part will inspect the parameter estimates of the model with the goal of determining whether there is any evidence for differences between the conditions. As before, this part is completely self sufficient and can be run without running the code of Parts I or II.</p>
<p>As I promised in the second part of this series of blog posts, the third part did not take another two months to appear. No, this time it took almost eight month. I apologize for this, but we all know the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Planning_fallacy" target="_blank">planning fallacy</a> and a lot of more important things got into the way (e.g., teaching).</p>
<p>As this part is relatively long, I will provide a brief overview. The next section contains a short explanation for the way in which we will perform hypothesis testing. This is followed by a short section loading some packages and the fitted model object and giving a small recap of the model. After this comes one relatively long section looking at the drift rate parameters in various ways. Then we will take look at each of the other three parameters in turn. Of especial importance will be the subsection on the non-decision time. As described in more detail below, I believe that this parameter cannot be interpreted. In the end, I give a brief overview of some of the limits of the present model and how it could be improved upon.</p>
<h1>Bayesian Hypothesis Testing</h1>
<p>The goal of this post is to provide evidence for differences in parameter estimates between conditions. This posts will present different ways to do so. Importantly, different ways of how to produce such evidence is only meant in the technical sense. In statistical terms we will always do basically the same thing: inspect difference distributions resulting from linear combinations of cell-wise posterior distributions of the group-level model parameter estimates. The somewhat technical phrase “linear combinations of cell-wise posterior distributions” often simply means the difference between two distributions. For example, the difference distribution resulting from subtracting the posterior of the speed condition from the posterior of the accuracy condition.</p>
<p>As a reminder, a posterior distribution is the probability distribution of the parameter conditional on data and model (where the latter includes the parameter priors). It answers the question which parameters are likely given our prior knowledge and the data. Therefore, the posterior distribution of the difference answers, for example, which difference values between two conditions are likely or not. With such a difference distribution we can then do two things.</p>
<p>First, we can check whether the x%-highest posterior density (HPD) or credibility interval of this difference distribution includes 0. If 0 is within the 95% HPD interval it could be seen a plausible value. If 0 is outside the 95% interval we could regard it as not plausible enough and would conclude that there is evidence for a difference.</p>
<p>Second, we can evaluate how much of the difference distribution is on one side of 0. If this value is considerably away from 50%, this constitutes evidence for a difference. For example, if all of the posterior samples for a specific difference are larger than zero, this provides considerable evidence that the difference is above 0.</p>
<p>The approach of investigating posterior distributions to gauge differences between conditions is only one approach for hypothesis testing in a Bayesian setting. And, at least in the psychological literature, it is not the most popular one. More specifically, many of the more vocal proponents of Bayesian statistics in the psychological literature advocate hypothesis testing using Bayes factors (e.g., <span rel="{ 'api_user_id': '881472', 'pages': '', 'items': '{8DDMVQM7},{6KZWXEP9}', 'format': '%a%, %d%, %p%', 'brackets': '', 'etal': '', 'separator': '', 'and': '' }"></span>). One prominent exception to this rule in psychology is maybe John <span rel="{ 'api_user_id': '881472', 'pages': '', 'items': '{298PA4EJ}', 'format': '%a% (%d%)', 'brackets': '', 'etal': '', 'separator': '', 'and': '' }"></span>. However, he proposes yet another approach of inference based on posterior distributions as used here. In general, I agree with many of the argument pro Bayes factor, especially in cases as the current one in which all relevant hypothesis or competing models are nested within one large (super) model.</p>
<p>The main difficulty when using Bayes factors is their extreme sensitivity to the parameter priors. In a situation with nested models, this is in principle not such a big problem, because one could use Jeffrey’s default prior approach (e.g., <span rel="{ 'api_user_id': '881472', 'pages': '', 'items': '{WIJMV4QQ}', 'format': '%a%, %d%, %p%', 'brackets': '', 'etal': '', 'separator': '', 'and': '' }"></span>). <span rel="{ 'api_user_id': '881472', 'pages': '', 'items': '{45ZMK47X}', 'format': '%a% (%d%, %p%)', 'brackets': '', 'etal': '', 'separator': '', 'and': '' }"></span> have extended this approach to general ANOVA designs (I am sure they were not the first to have this idea, but they were at least the first to popularize this idea in psychology). <a rel="nofollow" href="https://psyarxiv.com/9g4et/" target="_blank">Quentin Gronau and colleagues</a> have applied it to accumulator models, including the diffusion model. The general idea is to reparameterize the model using effect parameters which are normalized using, for example, the residual variance. For example, for a two sample design parameterize the model using a standardized difference such as Cohen’s d. Then it is comparatively easy and uncontroversial to put a prior on the standardized effect size measure. In the present case, in which the model does not contain a residual variance parameter, one could use the variance estimate of the group-level distribution for each parameter for such a normalization.</p>
<p>Unfortunately, <code>brms</code> does to the best of my knowledge not contain the ability to specify a parameterization and prior distribution in line with Jeffrey’s default Bayes factor. And as far as I remember a discussion I had on this topic with Paul Bürkner some time ago, it is also unlikely <code>brms</code> will ever get this ability. Consequently, I feel that <code>brms</code> is not the right tool for model selection using Bayes factors. Whereas it offers this ability now from a technical side (using our <a rel="nofollow" href="https://cran.r-project.org/package=bridgesampling" target="_blank"><code>bridgesampling</code></a> package), it only allows models with an unnormalized parameterization. I believe that such a parameterization is in most cases not appropriate for Bayes factors based model selection as the priors cannot be specified in a ‘default’ manner. Thus, I cannot recommend <code>brms</code> for Bayes factor based model selection at the moment. In sum, the reason for basing our inferences solely on posterior distributions in the present case is practical constraints and not philosophical considerations.</p>
<p>One final word of caution for the psychological readership. Whereas Bayes factors are clearly extremely popular in psychology, this is not the case in many other scientific disciplines. For example, the patron saint of applied Bayesian statistics, Andrew Gelman, is a <a rel="nofollow" href="https://andrewgelman.com/2017/07/21/bayes-factor-term-came-references-generally-hate/" target="_blank">self-declared opponent of Bayes factors</a>: “I generally hate Bayes factors myself”. As far as I can see, this disagreement comes from the different type of data different people work with. When working with observational (or correlational) data, as Andrew Gelman usually does, tests of the presence of effects (or of nullity) are either a big no-no (e.g., when wanting to do causal inference) or simply not interesting. We know that the real world is full of relationships, especially small ones, between arbitrary things. So getting effects simply by increasing N is just not interesting and estimation is the more interesting approach. In contrast, for experimental data, we often have true null hypothesis and testing of those makes a lot of sense. For example, <a rel="nofollow" href="https://slate.com/health-and-science/2017/06/daryl-bem-proved-esp-is-real-showed-science-is-broken.html" target="_blank">if Bem was right and there truly were PSI</a>, we could surely exploit this somehow. But as far as we can tell, the effect is truly null. In this case we really need ypothesis testing.</p>
<h1>Getting Started</h1>
<p>We start with loading some packages for analyzing the posterior. Since the beginning of this series, I have more and more become a fan of the whole <code>tidyverse</code>, so we import it completely. We of course also need <code>brms</code>. As shown below, we will need a few more packages (especially <code>emmeans</code> and <code>tidybayes</code>), but these are only loaded when needed.</p>
<pre data-enlighter-language="null">library(&quot;brms&quot;)
library(&quot;tidyverse&quot;)
theme_set(theme_classic()) # theme for ggplot2
options(digits = 3)</pre>
<p>Then we will also need the posterior samples, we can load them in the same way as before from my github page. Note that we neither need the data nor the posterior predictive distribution this time.</p>
<pre data-enlighter-language="null">tmp &lt;- tempdir()
download.file(&quot;https://singmann.github.io/files/brms_wiener_example_fit.rda&quot;,
file.path(tmp, &quot;brms_wiener_example_fit.rda&quot;))
load(file.path(tmp, &quot;brms_wiener_example_fit.rda&quot;))</pre>
<p>We begin with looking at the group-level posteriors. An overview of their posterior distributions can be obtained using the <code>summary</code> function.</p>
<pre data-enlighter-language="null">#                                    Estimate Est.Error l-95% CI u-95% CI
# conditionaccuracy:frequencyhigh      -2.944    0.1971   -3.345   -2.562
# conditionspeed:frequencyhigh         -2.716    0.2135   -3.125   -2.299
# conditionaccuracy:frequencynw_high    2.238    0.1429    1.965    2.511
# conditionspeed:frequencynw_high       1.989    0.1785    1.626    2.332
# bs_conditionaccuracy                  1.898    0.1448    1.610    2.186
# bs_conditionspeed                     1.357    0.0813    1.200    1.525
# ndt_conditionaccuracy                 0.323    0.0173    0.289    0.358
# ndt_conditionspeed                    0.262    0.0154    0.232    0.293
# bias_conditionaccuracy                0.471    0.0107    0.449    0.491
# bias_conditionspeed                   0.499    0.0127    0.474    0.524
# Warning message:
# There were 7 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
# See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup</pre>
<p>As a reminder, we have data from a lexical decision task (i.e., participants have to decide whether presented strings are a word or not) and <code>frequency</code> is the factor determining the true status of a string, with <code>high</code> referring to words and <code>nw_high</code> to non-words. Consequently, for the drift rate (the first four rows in the results table) the frequency factor determines the sign of the parameter estimates with the drift rate for words (rows 1 and 2) being clearly negative (i.e., those trials mostly hit the lower boundary for the word decision) and the drift rate for non-words (rows 3 and 4) being clearly positive (i.e., those trials mostly hit the upper boundary for non-word decisions). Furthermore, there could be differences between the drift rates in the accuracy or speed conditions. Specifically, in the speed conditions drift rates seem to be less extreme (i.e., nearer to 0) compared to the accuracy conditions.</p>
<p>The other three parameters, only differ between the <code>condition</code> factor. Given the experimental manipulation of accuracy versus speed condition, we expect differences for the boundary separation, parameters starting with <code>bs_</code>. For the non-decision time, parameters starting with <code>ndt_</code>, there also appears to be a small effect as the 95% only overlap slightly. However, as discussed in detail below, we should be careful in interpreting this difference. Finally, for bias, parameters starting with <code>bias_</code>, there might be a difference or not. Furthermore, at least in the accuracy condition there appears to be a bias for “word” responses.</p>
<p>One way to test differences between conditions is using the <code>hypothesis</code> function in <code>brms</code>. However, I was not able to get it to work with the current model. I suspect the reason for this is the somewhat unconventional parameterizations where each cell gets one parameter (in some sense each cell has its own intercept, but there is no overall intercept). This contrasts with a more “standard” parameterization in which there is one intercept (for either the unweighted means or one of the cells) and the remaining parameters capture the differences between the intercept and the cell means. As a reminder, I chose this unconventional parameterization in the first place to make the specification of the parameters priors easier. Additionally, this is a common parameterization when programming cognitive models by hand.</p>
<h1><code>emmeans</code> and <code>tidybayes</code>: Differences in the Drift Rate</h1>
<p>An alternative is to use the great <a rel="nofollow" href="https://cran.r-project.org/package=emmeans" target="_blank"><code>emmeans</code></a> package by Russel Lenth. I am a huge fan of <code>emmeans</code> and use it all the time when using “normal” statistical models (e.g., ANOVAs, mixed models), independent of whether I use frequentist methods (e.g., via <code>afex</code>) or Bayesian methods (e.g., <a rel="nofollow" href="https://cran.r-project.org/package=rstanarm" target="_blank"><code>rstanarm</code></a> or <code>brms</code>). Unfortunately, it appears as if <code>emmeans</code> at the moment only allows an analysis of the main parameter of the response distribution for models estimated with <code>brms</code>, which in our case is the drift rate. If someone were to extend <code>emmeans</code> to allow using <code>brms</code> models with all parameters, I would be very happy and thankful. In any case, I highly recommend to check out the <a rel="nofollow" href="https://cran.r-project.org/package=emmeans" target="_blank"><code>emmeans</code> vignettes</a> to get an overview of what type of follow-up tests are all possible with this great package.</p>
<p>As I <a rel="nofollow" href="https://discourse.mc-stan.org/t/how-do-i-get-marginal-effects-for-categorical-variables-to-condition-on-an-average-rather-than-a-category/5323/6?u=henrik_singmann" target="_blank">recently learned</a>, <code>emmeans</code> works quite nicely together with <a rel="nofollow" href="https://cran.r-project.org/package=tidybayes" target="_blank"><code>tidybayes</code></a>, a package that enables working with posterior draws within the <code>tidyverse</code>. <code>tidybayes</code> has a surprisingly large package footprint (i.e., it imports quite a lot of other packages) for a package with a comparatively small functionality. I guess this is a consequence of being embedded within the <code>tidyverse</code>. In any case, many of the imported packages are already in the search path thanks to loading the <code>tidyverse</code> above and attaching should not take that long here.</p>
<pre data-enlighter-language="null">library(&quot;emmeans&quot;)
library(&quot;tidybayes&quot;)</pre>
<p>We begin with <code>emmeans</code> only to assure ourselves that it works as expected. For this, we get the estimated marginal means plus 95%-highest posterior density (HPD) intervals which match the output of the fixed effects for the estimate of the central tendency (which is the median of the posterior samples in both cases). As a reminder, the fact that the cell estimates match the parameter estimates is of course a consequence of the unusual parameterization which is picked up correctly by <code>emmeans</code>. The lower and upper bounds of the intervals differ slightly between the <code>summary</code> output from <code>brms</code> and <code>emmeans</code>, a consequence of using different ways of calculating the intervals (i.e., quantiles versus HPD intervals).</p>
<pre data-enlighter-language="null">fit_wiener %>%
  emmeans( ~ condition*frequency) 
#  condition frequency emmean lower.HPD upper.HPD
#  accuracy  high       -2.94     -3.34     -2.56
#  speed     high       -2.72     -3.10     -2.28
#  accuracy  nw_high     2.24      1.96      2.50
#  speed     nw_high     1.99      1.64      2.34
# 
# HPD interval probability: 0.95</pre>
<h2>Using HPD Intervals And Histograms</h2>
<p>As a first test, we are interested in assessing whether there is evidence for a difference between speed and accuracy conditions for both words (i.e., <code>frequency</code> = <code>high</code>) and non-words (i.e., <code>frequency</code> = <code>nw_high</code>). There are many ways to do this with <code>emmeans</code> one of them is via the <code>by</code> argument and the <code>pairs</code> function.</p>
<p> </p>
<pre data-enlighter-language="null">fit_wiener %>%
  emmeans(&quot;condition&quot;, by = &quot;frequency&quot;) %>% 
  pairs
# frequency = high:
#  contrast         estimate lower.HPD upper.HPD
#  accuracy - speed   -0.225   -0.6964     0.256
# 
# frequency = nw_high:
#  contrast         estimate lower.HPD upper.HPD
#  accuracy - speed    0.249   -0.0647     0.550
# 
# HPD interval probability: 0.95</pre>
<p>Here, we do not have a lot of evidence that there is a difference for either stimulus type, as both HPD intervals include 0.</p>
<p>Instead of getting the summary of the distribution via <code>emmeans</code>, we can also use the capabilities of <code>tidybayes</code> and extract the samples in a tidy way. Then we use one of the convenient aggregation functions coming with <code>tidybayes</code> and aggregate the samples based on the same conditioning variable. After trying a few different options, I have the feeling that <code>emmeans</code>‘ <code>hpd.summary()</code> function uses the same approach for calculating HPD intervals as <code>tidybayes</code>, as both results match.</p>
<pre data-enlighter-language="null">samp1 &lt;- fit_wiener %>%
  emmeans(&quot;condition&quot;, by = &quot;frequency&quot;) %>% 
  pairs %>% 
  gather_emmeans_draws()
samp1 %>% 
  median_hdi()
# # A tibble: 2 x 8
# # Groups:   contrast [1]
#   contrast         frequency .value  .lower .upper .width .point .interval
#   </pre>