Why hierarchical models?

Hierarchical models are preferred when dealing with data that have inherent grouping structures, such as repeated measures across subjects or items. Unlike simpler models that ignore or simplify within-group variability, hierarchical models incorporate both fixed effects (e.g., the variable cond in our data) and random effects (e.g., part and item in our data). In the context of data with repeated observations per subject and item, in which the independence assumption is violated, hierarchical models can account for these sources of variability, allowing for more accurate estimates that reflect the true complexity of the data. By including random intercepts or slopes, hierarchical models can capture subject- or item-specific influences, which results in a more nuanced understanding of the effects being studied.

One important aspect of hierarchical models involves a process called shrinkage (also known as partial pooling). Shrinkage adjusts the estimates for random effects by pulling extreme estimates closer to the group mean. This regularization reduces the impact of outliers or overly variable estimates, leading to more stable and conservative results. As will be seen below, shrinkage can also affect fixed effects estimates, a key point in the analysis that follows.

Before hierarchical models became mainstream in second language studies, a popular method of analysis involved running separate ANOVAs on data aggregated by subject and by item. There are, of course, negative consequences of such a method. First, two analyses are needed a priori, one for subject-aggregated data, one for item-aggregated data. As a result, participant- and item-level predictors cannot be included in the same analysis. Second, missing items in one or more condition levels must be removed. Third, ANOVAs do not typically focus on effect sizes, which in turn requires a separate method to quantify the magnitude of effects (e.g., $ {\eta}^2 $ ). We will not discuss ANOVAs in this paper, but one of the linear models we will discuss next (LM) is equivalent to an ANOVA run on participant-aggregated data.

Examine Figure 2, which compares three different models for our data. Along the $ x $ -axis, LM represents a simple linear regression run over the aggregated mean response by participant for each condition. LM+B represents the same model, but with a Bonferroni correction applied to the post hoc comparisons, which affects the confidence intervals of those estimates (but notice that the standard errors are not affected). Finally, HLM shows the estimates from the hierarchical model from Code 2. Notice how the estimate for c4 (highlighted in the figure) is substantially closer to the grand mean, i.e., it has shrunk. In addition, the confidence intervals shown under HLM are wider than the intervals under LM.Footnote ¹ Given these two characteristics of the HLM in question, multiple comparisons stemming from such a model are more reliable and generally more conservative than those performed from LM or LM+B—this will also be clear when we examine our model figure later. It is also worth reiterating that HLM includes both by-participant and by-item random intercepts in the same analysis, and is therefore objectively more complex and more comprehensive than the other two models. The reader can find a similar scenario and a more detailed discussion in Gelman, Hill, and Yajima (Reference Gelman, Hill and Yajima2012).

Figure 2. Estimate comparison across three models ( $ x $ -axis). Standard errors (horizontal lines) and 95% confidence intervals (vertical lines) are also shown for each estimate. The dashed line represents the average across all conditions (panels on top).

Shrinkage is not always visible in the estimates of our fixed effects, since its effect is more directly observable in our random effects. The reason why we see shrinkage in our estimate for c4 here is because our data is imbalanced: condition c4 has only 50 observations (cf. $ n= $ 189–215 for the other conditions). This represents a situation that we often encounter, where a particular group of interest is imbalanced and may have too much influence on the estimates of our models. Our hierarchical model here mitigates the problem by pulling the estimate of c4 towards the grand mean, since it does not have enough empirical evidence (i.e., observations) to grant the estimate provided by LM and LM+B.

In general, then, hierarchical models allow for more reliable confidence intervals even without corrections such as Bonferroni, a desirable characteristic, given that such corrections typically reduce statistical power and tend to modify only confidence intervals (or $ p $ -values), keeping estimates stationary. Hierarchical models also reduce the issue of multiple comparisons, since point estimates can be shifted (shrinkage), as demonstrated above—see Gelman & Tuerlinckx (Reference Gelman and Tuerlinckx2000), Gelman et al. (Reference Gelman, Hill and Yajima2012), and Kruschke (Reference Kruschke2015, pp. 567–568). And because multiple comparisons are derived from the estimates of each level being compared, more reliable estimates in our models lead to more reliable comparisons from said models.

In the sections that follow, we will reproduce the analysis in Code 1 using a Bayesian approach. In doing so, we will also discuss fundamental differences between Frequentist and Bayesian statistics, and how favoring the latter can offer intuitive insight on the estimation of parameters, especially when it comes to multiple comparisons.

Statistical model

In this subsection, we will run our hierarchical Bayesian model, which will have a single fixed predictor (cond). I will assume that you have never run a Bayesian model or used the brms package. While specific details of our models are outside the scope of this paper, I will attempt to provide all the necessary information to reproduce and apply these methods to your data. In the next subsection, we will work on our multiple comparisons based on the model bFit1 that we run and discuss below.

In Code 5, we define our model as resp ~ cond using the same familiar syntax from the lm() function. Notice that this is the most complex model thus far, as it includes random slopes for cond by participant—this model specification yields a singular fit using a Frequentist approach with lmer(). Because we are running a linear regression, we assume that the response variable is normally distributed, hence family = gaussian() in the code. We could also use another family value if we wanted to better accommodate outliers. For example, we could specify family = student(), which would assume a student $ t $ distribution for our response variable resp. This can be useful in situations where we have smaller sample sizes (higher uncertainty) and/or potential outliers.

By default, our model will use four chains to sample from the posterior. To speed up the process, we can use four processing cores. The argument save_pars is necessary if the reader wishes to calculate the Bayes Factor later on using the bayes_factor() function in brms. If that is the case, a null model is also needed: that is why bFitNull is also present in the code. The reader can run both models and subsequently run bayes_factor (bFit1, bFitNull) to confirm that the Bayes Factor will favor the null model (this is in line with Code 1), i.e., that cond has no effect on respFootnote ³ —the line is commented out in Code 5 because we will not explore those results, as mentioned above. Finally, we set our priors (see below) and save the model (optionally) to have access to the compiled Stan code later. While this last step is not necessary, if you want to better understand the specifics of our model, it is useful to examine the raw Stan code later on.

The model in question, bFit1, will not run until we create the variable priors1 referred to in the code. First, it should be clarified that we do not need to specify priors ourselves: brm() will use its own set of default priors if we do not provide ours. Here, however, we will assume some mildly informative priors. The default priors used by brm() will often be flat, which means the model starts by assuming that values are equally plausible. This, however, is neither realistic nor ideal: given the data at hand, we know that no responses are lower than 4 or higher than 8, given Figure 1.Footnote ⁴ We also know that extreme reaction times are unlikely in the literature of lexical decision tasks. Importantly, we know that not all values for our condition levels are equally likely. Imagine that we were estimating the height of a given person. A height of seven meters is clearly impossible, and our model should not start by assuming that such an extraordinary height is just as likely as, say, 1.75 meters. Likewise, a reaction time of 6 (approximately 400 ms) of much more likely than a reaction time of 10 (approximately 22,000 ms) in such tasks, and reaction times under 2 (approximately 100ms) are considered to be unreliable (Luce, Reference Luce1991; Whelan, Reference Whelan2008).

In Code 6, we first use the function get_prior() to inspect the priors assumed by the model (default priors). The reader can inspect the output of get_prior() in Code 7. Note that we are feeding get_prior() with the nonhierarchical version of the model to simplify our analysis here, but you could certainly extract and modify priors on the mixed effects as well. We then set our priors with the function set_prior(). Here, we are assuming that all of the slopes (c2 to c5) follow a Gaussian (normal) distribution centered around zero with a standard deviation of 1: $ \mathcal{N}\left(\mu =0,\sigma =1\right) $ . Thus, the priors in question contain no directional bias, and they also introduce a relatively high degree of uncertainty, given the standard deviation we assume for each distribution. For comparison, in our actual data, the standard deviation of our response variable, resp, ranges from $ .16 $ to $ .21 $ for all five conditions, cond, as can be verified in Table 2.Footnote ⁵ These priors tell the model that absurd values are simply not plausible. Finally, our intercept ( $ {\hat{\beta}}_0 $ ) is assumed to have a mean of 7 (which approximates the observed mean of c1 already visualized).

In summary, we are telling our model that there should be no differences between the conditions, and we are constraining the range of potential differences by saying that values near zero are more plausible than values far away from zero. This starting point is very distinct from assuming flat priors. We are merely telling the model that certain values are much more plausible than others, and we are being conservative by assuming priors centered at zero with relatively wide distributions.

Since we have now defined our priors, we can finally run our model. Running a Bayesian model takes much longer than running an equivalent Frequentist model. For that reason, it is a good idea to save the model for later use: save(bFit1, file = “bFit1.RData”) will save the model in RData format. You can later read it by using load("bFit1.RData"). Importantly, we can save however many objects in a single RData file. The population-level effects of our model, bFit1, are shown in Code 8—the reader should run bFit1 to explore the entire output of the model.

In the Estimate column, we can see the mean estimate for each condition (recall that each estimate represents an entire posterior distribution). We also see the uncertainty for each estimate, quantified as the standard deviation of the posterior distribution,Footnote ⁶ and the two-sided 95% credible intervals (l–95% CrI and u–95% CrI) based on posterior quantiles Bürkner (Reference Bürkner2017, p. 11). For symmetrical posterior distributions (e.g., Gaussian distributions), this interval coincides with the highest density interval (HDI). To calculate the exact 95% HDI of a posterior distribution, we can use the hdi() function from the bayestestR package. Following Kruschke (Reference Kruschke2015, p. 342), I favor HDIs over CrIs, as HDIs represent the narrowest interval containing the most probable values, offering a more intuitive interpretation, especially when dealing with skewed distributions—although such skewness is not a concern in the present paper.

Next, we see Rhat ( $ \hat{R} $ should be $ <1.01 $ ),Footnote ⁷ which tells us that the model has converged. Finally, we see two columns showing us the effective sample size (ESS) of our posterior samples: The bulk ESS corresponds to the main body of the posterior distribution, while the tail ESS corresponds to the tails (i.e., more extreme values) in the posterior distribution. By default, each of the four chains draws 2,000 samples from the posterior. However, we use a warmup of 1,000 samples. As a result, only 1,000 samples from each chain are actually used here, totaling 4,000 samples. Of these 4,000, however, some “steps” in our random walk will be correlated with one another, which reduces the amount of information they provide about the parameter space. Therefore, we should prioritize uncorrelated steps, which results in a proper subset of the 4,000 steps in question. That is why all the numbers in both ESS columns are less than 4,000. There is no universal rule of thumb as to how large the ESS should be (higher is always better, of course), but 100 samples per chain is a common recommendation. If the ESS is lower than 400, we could increase the number of steps by adding iter = 4000, for example, which would give us 3,000 post-warmup steps for each chain.Footnote ⁸ We can also change the number of warmup steps in the model, allowing it to “settle” in more informative regions of parameter values before extracting the samples we will actually use.

The reader will notice that the estimates in our Bayesian model are not very different from the estimates in our Frequentist model, run in Code 1. This makes sense, since our priors are only mildly informative. Ultimately, the posterior distribution in our model is (mostly) driven by the patterns in the data. In a more realistic application, we would attempt to set more precise priors given what we know from the literature on the topic of interest.

While tables can be useful to inspect the results from our model, figures are often a better option as they show the actual posterior distribution for each parameter in our model. In this paper, I will employ a specific type of figure, shown in Figure 4. The figure contains different pieces of information, some of which are standard in Bayesian models. First, note that a dashed line marks zero. Around that line, we find the ROPE. Now, we can assess the posterior distribution and the HDI of each parameter (our primary goal) as well as their location relative to the ROPE (our secondary goal). As usual, all conditions in the figure must be interpreted relative to c1 (our intercept; not shown). The central portion of each distribution represents the most plausible values given the data (and our priors). It is important to reiterate that the HDIs in the figure are a summary statistic from actual probability distributions, unlike Frequentist confidence intervals. As such, values at the tails of the HDIs are less plausible than values at the center of the HDIs. The figure therefore gives us a comprehensive view of effect sizes and parameter estimation.

Figure 4. Posterior distributions from hierarchical model with associated 95% HDIs. Region of practical equivalence is represented by shaded area around zero.

Next, we can move to our secondary goal. The reader will notice that all of our posterior distributions include zero in their 95% HDI. Logically, this means that every posterior in Figure 4 is at least in part within the ROPE. Simply put, zero is a credible parameter value representing the difference between each of c2–c5 and c1, hence the result of the Bayes Factor alluded to earlier, which favored the null model. That being said, we should also observe how much of each distribution is within the ROPE. For example, as we compare the distributions of c5 and c4, we can see that a null effect is more plausible for c4 than for c5, even though we cannot accept any null effect for the slopes in question, given that no HDI is entirely within the ROPE.

Multiple comparisons

In this section, we will use our hierarchical model bFit1 to extract all the multiple comparisons we wish. There are at least two ways of accomplishing this task. The easy (automatic) way involves a function called hypothesis() from brms. While this function does provide us with estimates on any contrast of interest, its main application is in hypothesis testing, i.e., our secondary goal here—see example in Veríssimo (Reference Veríssimo2024) (p. 14). Still, because at times this is what is needed, and because it is an easy approach to implement, we will start by exploring the function in question. Then, we will work our way through a more manual approach that will focus on our primary goal here, namely, parameter estimation.

Comparisons via hypothesis testing

Given the levels of our factor cond, we will begin by defining a null hypothesis, represented in Code 9 by the contrasts variable. Here, we simulate three null hypotheses by stating that the posterior distributions of c3 and c2 are equal, as are the posteriors of c5 and c4, and of c4 and c3. We can then use hypothesis() to test our model against said hypotheses using a given posterior probability threshold, the alpha Footnote ⁹ argument in the function hypothesis().Footnote ¹⁰ We could also include directional hypotheses such as condc4 > condc2.

To understand what the function hypothesis() is doing, let us pick the hypothesis condc3 = condc2 as an example. The function takes the posterior distribution of each parameter (c3 and c2) and generates the posterior distribution for their difference. Then it calculates the credible interval for this new posterior distribution (here, the interval is 95% given our threshold defined by alpha). If zero is not found within that interval, we will have strong evidence that condition c3 is not equal to condition c2 (the column Star in our output will have an asterisk). If zero is found within that interval, we would not have enough evidence to conclude that a difference exists between the two conditions.

The output generated by hypothesis() is shown in Table 3 with resulting estimates and credible intervals. Because we have assigned the result to a variable (h), we can also plot the resulting posterior distributions simply by running plot(h).

Table 3. Simplified output of the hypothesis() function from brms

In our comparisons table, we can see that two of the three hypotheses can be rejected (assuming a posterior threshold of $ .05 $ ), namely, c3–c2 ( $ \hat{\beta}=.10 $ ) and c5–c4 ( $ \hat{\beta}=-.09 $ ): in both cases, the credible interval does not include zero. To be clear, the estimates in our table represent the means of the posterior distributions that result from each comparison of distributions. These are therefore single-point estimates that summarize entire posterior distributions (as is always the case in Bayesian models). For example, we saw in Code 8 that the mean estimate for c2 is $ \hat{\beta}=-.06 $ , and the mean estimate for c3 is $ \hat{\beta}=.05 $ —both of which represent the difference between each condition and c1, our intercept. Thus, it is not surprising that the posterior distribution of c2–c3 has a mean estimate of $ \hat{\beta}=-.10 $ (estimates have been rounded to two decimal places).

You may be wondering how c1 could be included in the multiple comparisons. It is true that all of the comparisons involving c1 are by definition in the result of bFit1, since c1 is our intercept. However, if we want to list all comparisons, we will have to include those comparisons as well. There are two important details in this scenario: first, c1 does not exist in the model as it is called Intercept.Footnote ¹¹ Thus, we cannot simply add condc1 to our contrasts variable as the function hypothesis() will not be able to find such a parameter in the model. Rather, you’d simply type Intercept. Second, in a standard model, slopes encompass only their differences relative to the intercept. As a result, if we want to test the hypothesis that c1 and c3 are identical, for example, we cannot simply write Intercept = c3. Instead, you’d type c3 = 0, which is essentially the same as Intercept = c3 + Intercept. Indeed, c3 + Intercept gives you the actual posterior draws for the absolute values of c3. Note that this has nothing to do with the Bayesian framework that we are employing. This is just how a typical regression model works.

How could we report our results thus far, assuming only the comparisons in Code 9? We can be more or less comprehensive on the amount of detail we provide (see Kruschke, Reference Kruschke2021). Here is a comprehensive example that emphasizes the hypothesis testing approach just described and adds some information on the parameter estimation from Figure 4:

Comparisons via parameter estimation

An alternative to the method explored above involves the actual estimation of effect sizes from multiple comparisons. For example, if we wish to compare c3 and c4, we could simply subtract their posterior distributions and generate a figure with the resulting distribution. While this is also accomplished when we run hypothesis(), we can directly extract all pairwise comparisons using emmeans() as seen earlier in Code 2. This is shown in Code 10.Footnote ¹⁴ Note that here emmeans() outputs the median and the HPD (highest posterior density), which is the smallest interval in the posterior that contains a given proportion of the total probability (e.g., 95% HPD). Although not identical in their meanings, HPDs and HDIs are often used interchangeably in distributions that are approximately Gaussian, our case here.

While functions such as emmeans() are extremely useful, there are two advantages of doing things manually at least once: first, it forces us to understand exactly what is being done, which is especially useful to those new to Bayesian models; second, it allows us to customize our comparisons.

In what follows, we will visualize all the comparisons for our variable cond. We will, however, extract posterior samples manually to generate a custom figure that is informative and comprehensive. The goal here is not to suggest that the same type of figure should necessarily be used in similar analyses. Instead, our goal will be to better understand multiple comparisons in the process of generating a complete figure of what we can extract from our models. Using functions such as hypothesis(), which automatically accomplish a task for us, can be extremely practical, but such functions can also conceal what is actually happening in the background, which can be a problem if we want to better understand the specifics. We will start with our final product, i.e., a figure with all ten multiple comparisons for cond. As previously mentioned, we will also add a ROPE to our figure: even though the figure focuses on parameter estimation, having the ROPE in it allows the reader to reject or accept null effects as well. In simple terms, this makes our figure more flexible and comprehensive. Here are the elements shown in Figure 5:

1. 1. Posterior distributions for all multiple comparisons based on bFit1. The figure also displays the same comparisons from a nonhierarchical model (analogous to model LM in Figure 2; not run here), bFit0, also included in the bFits.RData file. These posteriors are shown with a dark gray border.

2. 2. Means and HDIs of our hierarchical model bFit1 (right).

3. 3. Percentage of each HDI from our hierarchical model contained within the ROPE (left).

4. 4. Random effects by participant, represented with grey circles around the means of the posterior distributions.

Figure 5. A complete figure containing posterior distributions of multiple comparisons using our hierarchical model. Posterior distributions from a nonhierarchical model (analogous to LM in Figure 2) are shown with dark gray borders.

At this point, the reader should immediately see the shrinkage caused by bFit1 in Figure 5 relative to the nonhierarchical model bFit0. The priors on fixed effects used in both models are identical, which means any shrinkage due to priors can be ignored in this comparison. Thus, the only difference between the two Bayesian models is that one is hierarchical, and the other is not. Every posterior distribution involving c4, the condition with the smallest number of observations, shows considerable shrinkage: our hierarchical model is more conservative and estimates such comparisons closer to zero, which is exactly what we would expect given the shrinkage observed in the estimates of our fixed effects for c4.

By design, Figure 5 is comprehensive in what it shows. As such, it requires multiple steps to be generated. In general, if the reader is looking for a quick method, I would recommend using the much simpler approach in Code 9, which creates a similar result in very few steps, although that approach will be more focused on hypothesis testing. Nevertheless, as mentioned earlier, reproducing Figure 5 can be useful to fully understand what we are doing exactly, which in turn can help us better understand the models we are using. Ultimately, it also focuses on the magnitude of effects.

Below, I list the necessary steps involved in creating Figure 5, some of which are required for both models. The reader can find a complete script to reproduce all the steps and the figure on the OSF repository mentioned at the beginning of the paper.

1. 5. Extract draws using the as_draws_df() function

2. 6. Select and rename variables (e.g., from b_Intercept to simply c1)

3. 7. Create comparisons

4. 8. Long-transform the data

5. 9. Create statistical summaries for means and HDIs using the mean_hdi() function

6. 10. Generate ROPE using the rope() function

7. 11. Extract random effects from bFit1

8. 12. Create a summary for mean (random) effects by participant

9. 13. Generate a figure

Heteroscedasticity

In our analysis above, we assume that the variance ( $ {\sigma}^2 $ ), and therefore the standard deviation ( $ \sigma $ ), is the same across all of our conditions, i.e., they are homoscedastic. This is indeed supported by the data insofar as our conditions have very similar variances, as shown in Table 2, and thus the variance in our residuals will be constant across all levels of cond. Sometimes, however, our conditions will not necessarily meet the assumption of equal variances. Indeed, as pointed out by Birdsong (Reference Birdsong2018, p. 1), “non-uniformity is an inherent characteristic of both early and late bilingualism.” For example, heteroscedastic data are frequently found in age effects, as older learners show higher variance than younger learners when it comes to ultimate attainment (Vanhove, Reference Vanhove2013, p. 11).

The situation just described brings us to one additional advantage of a Bayesian approach, namely, that we can directly deal with heteroscedastic data by estimating the variance itself. This provides the researcher with a higher degree of flexibility than typical Frequentist methods (e.g., robust regression models) and results in more reliable posterior distributions of effect sizes, as will be shown below. Here, again, it is clear that more reliable posterior distributions will directly impact the reliability of multiple comparisons derived from said distributions. In this section, we will briefly see (a) how to estimate the standard deviation ( $ \sigma $ ) across the conditions in cond and (b) how this can affect our estimates.

The dataset d_het contains a modified version of d in which the condition c4 has a considerably lower variance than all the other conditions. This can be verified in Code 11 (cf. Table 2).

The reader can also see how the new variance of c4 affects the spread of the data, and thus the height of the box plot for c4, in Figure 6. The figure in question can be directly compared to Figure 1 presented at the beginning of this paper.

Figure 6. Overall patterns in the data: box plots and associated means (orange dots) and standard errors (not visible). Notice the variance of the c4 condition.

We will run two models. First, bFit2, which represents the same model as bFit1 run earlier. This will be a typical linear regression where we do not estimate varying standard deviations, i.e., a model that assumes homoscedasticity. We will then run bFit3, a model that estimates standard deviations (sigma) as well as the fixed and random effects estimated by bFit2. While our new model requires a new function in its specification (bf() within brm()), its interpretation is straightforward. These models will be directly comparable as they will share the same priors as well (with the exception of the prior set for sigma in bFit3). Both models, as well as the priors for bFit3, are shown in Code 12.

The reader can see the simplified output of bFit3 in Code 13. In addition to the estimates for cond, we also obtain estimates for $ \sigma $ for each level of cond. Recall that we assume a Gaussian distribution for the response variable. This distribution has two parameters: $ \mu $ (the mean) and $ \sigma $ (the standard deviation). In our model, the fixed effects ( $ \beta $ ) represent $ \mu $ , while $ \sigma $ is now estimated explicitly for each level of cond. The key distinction here is that both parameters of the Gaussian distribution ( $ \mu $ and $ \sigma $ ) are being estimated by our model. In brms, sigma is modeled on the log scale, as reflected in Code 13.

To appreciate the heteroscedasticity in our data, notice how each sigma_condc… deviates from sigma_Intercept. For c1 (our intercept), the estimated log-standard deviation is $ -1.91 $ , which corresponds to $ {\sigma}_1={e}^{-1.91}\approx .15 $ (cf. $ .20 $ in Code 11). For c4, the estimated log-standard deviation includes a more substantial adjustment of $ -.87 $ , making $ {\sigma}_4={e}^{-1.91-.87}\approx .06 $ (cf. $ .09 $ in Code 11). We can therefore see that the estimated sigma for c4 deviates considerably from that of the intercept (c1). This reflects what we already know, given Code 11 and Figure 6. The reader should print the outputs of both bFit2 and bFit3 to verify that the credible intervals of c4 are indeed estimated to be different in both models. Let us now examine how these differences affect our multiple comparisons. This is the goal of Figure 7, where black point ranges represent comparisons drawn from bFit3, and orange point ranges are drawn from bFit2—mixed effects and posterior distributions are all from bFit3.

Figure 7. Multiple comparisons from a model where sigma is also estimated (black). Notice the different HDIs in comparisons involving c4 relative to a model where sigma is not estimated (orange).

In Figure 7, we notice that comparisons involving c4 differ in their HDIs across both models. By estimating sigma, bFit3 generates narrower posteriors for c4, which the reader should verify by plotting the main effects of both models (without comparisons). This, in turn, affects the width of the posterior distributions of our comparisons in Figure 7. For instance, in the comparison between c4 and c2, the mean effect is slightly stronger in bFit3, and its HDI now barely touches the ROPE. Another comparison worth mentioning is c4 and c5, whose HDI has clearly shifted away from the ROPE. This demonstrates how accounting for heteroscedasticity can affect both parameter estimation and hypothesis testing, emphasizing the importance of accurate variance modeling.