# 25 Tools in the Trunk

“This chapter includes some important topics that apply to many different models throughout the book… The sections can be read independently of each other and at any time” (Kruschke, 2015, p. 721).

## 25.1 Reporting a Bayesian analysis

Bayesian data analyses are not yet standard procedure in many fields of research, and no conventional format for reporting them has been established. Therefore, the researcher who reports a Bayesian analysis must be sensitive to the background knowledge of his or her specific audience, and must frame the description accordingly. (p. 721)

At the time of this writing (mid 2022), this is still generally the case. See Aczel et al. (2020), *Discussion points for Bayesian inference*, for a recent discussion from several Bayesian scholars. For a take from a familiar voice, see Kruschke’s (2021) *Bayesian analysis reporting guidelines*.

### 25.1.1 Essential points.

Recall the basic steps of a Bayesian analysis from Section 2.3 (p. 25): Identify the data, define a descriptive model, specify a prior, compute the posterior distribution, interpret the posterior distribution, and, check that the model is a reasonable description of the data. Those steps are in logical order, with each step building on the previous step. That logical order should be preserved in the report of the analysis. (p. 722)

Kruschke then gave recommendations for motivating Bayesian inference. His (2018) paper with Liddell, *The Bayesian New Statistics: Hypothesis testing, estimation, meta-analysis, and power analysis from a Bayesian perspective*, might be helpful in this regard. Many of the other points Kruschke made in this section (e.g., adequately reporting the data structure, the priors, evidence for convergence) can be handled by adopting open science practices.

If your data and research questions are simple and straightforward, you might find it easy to detail these and other concerns in the primary manuscript. The harsh reality is many journals place tight constraints on word and/or page limits. If your projects are not of the simple and straightforward type, supplemental materials are your friend. Regardless of a journal’s policy on hosting supplemental materials on the official journal website, you can detail your data, priors, MCMC diagnostics, and all the other fine-grained details of your analysis in supplemental documents hosted in publicly-accessible repositories like the Open Science Framework (OSF) or GutHub. If possible, do consider making your data openly available. Regardless of the status of your data, please consider making all your **R** scripts available as supplementary material. To reiterate from Chapter 3, I strongly recommend checking out R Notebooks for that purpose. They are a type of R Markdown document with augmentations that make them more useful for working scientists. You can learn more about them here and here. And for a more comprehensive overview, check out Xie, Allaire, and Grolemund’s (2022) *R markdown: The definitive guide*.

### 25.1.2 Optional points.

For more thoughts on robustness checks, check out a couple Gelman’s blog posts, *What’s the point of a robustness check?* and *Robustness checks are a joke*, along with the action in the comments section.

In addition to posterior predictive checks, which are great (see Kruschke, 2013), consider prior predictive checks, too. For a great introduction to the topic, check out Gabry, Simpson, Vehtari, Betancourt, and Gelman’s (2019) *Visualization in Bayesian workflow*.

### 25.1.3 Helpful points.

For more ideas on open data, check out Rouder’s (2016) *The what, why, and how of born-open data*. You might also check out Klein and colleagues’ (2018) *A practical guide for transparency in psychological science* and Martone, Garcia-Castro, and VandenBos’s (2018) *Data sharing in psychology*.

As to posting your model fits, this could be done in any number of ways, including as official supplemental materials hosted by the journal, on GitHub, or on the OSF. At a base level, this means saving your fits as external files. We’ve already been modeling this with our `brm()`

code throughout this book. With the `save`

argument, we saved the model fits within the `fits`

folder on GitHub. You might adopt a similar approach. But do be warned: **brms** fit objects contain a copy of the data used to create them. For example, here’s how we might reload `fit24.1`

from last chapter.

`.1 <- readRDS("fits/fit24.01.rds") fit24`

By indexing the fit object with `$data`

, you can see the data.

```
library(tidyverse)
library(brms)
.1$data %>%
fit24glimpse()
```

```
## Rows: 16
## Columns: 4
## $ Count <dbl> 20, 68, 5, 15, 94, 7, 16, 10, 84, 119, 29, 54, 17, 26, 14, …
## $ Hair <fct> Black, Black, Black, Black, Blond, Blond, Blond, Blond, Bru…
## $ Eye <chr> "Blue", "Brown", "Green", "Hazel", "Blue", "Brown", "Green"…
## $ `Hair:Eye` <chr> "Black_Blue", "Black_Brown", "Black_Green", "Black_Hazel", …
```

Here’s a quick way to remove the data from the fit object.

`.1$data <- NULL fit24`

Confirm it worked.

`.1$data fit24`

`## NULL`

Happily, the rest of the information is still there for you. E.g., here’s the summary.

`summary(fit24.1)`

```
## Family: poisson
## Links: mu = log
## Formula: Count ~ 1 + (1 | Hair) + (1 | Eye) + (1 | Hair:Eye)
## Data: NULL (Number of observations: )
## Draws: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
## total post-warmup draws = 8000
##
## Group-Level Effects:
## ~Eye (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.80 1.51 0.30 5.80 1.00 2907 5001
##
## ~Hair (Number of levels: 4)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.89 1.52 0.35 6.10 1.00 3302 5159
##
## ~Hair:Eye (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.95 0.30 0.53 1.67 1.00 2534 4849
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 3.13 1.47 -0.04 6.12 1.00 4494 4726
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

## 25.2 Functions for computing highest density intervals

You can find a copy of Kruschke’s scripts, including `DBDA2E-utilities.R`

, at https://github.com/ASKurz/Doing-Bayesian-Data-Analysis-in-brms-and-the-tidyverse/tree/master/data.R.

### 25.2.1 R code for computing HDI of a grid approximation.

We can imagine the grid approximation of a distribution as a landscape of poles sticking up from each point on the parameter grid, with the height of each pole indicating the probability mass at that discrete point. We can imagine the highest density region by visualizing a rising tide: We gradually flood the landscape, monitoring the total mass of the poles that protrude above water, stopping the flood when 95% (say) of the mass remains protruding. The waterline at that moment defines the highest density region (e.g., Hyndman, 1996). (p. 725)

We can use Kruschke’s `HDIofGrid()`

function for such a task.

```
<- function(probMassVec, credMass = 0.95) {
HDIofGrid
# Arguments:
# probMassVec is a vector of probability masses at each grid point.
# credMass is the desired mass of the HDI region.
# Return value:
# A list with components:
# indices is a vector of indices that are in the HDI
# mass is the total mass of the included indices
# height is the smallest component probability mass in the HDI
# Example of use: For determining HDI of a beta(30,12) distribution
# approximated on a grid:
# > probDensityVec = dbeta( seq(0,1,length=201) , 30 , 12 )
# > probMassVec = probDensityVec / sum( probDensityVec )
# > HDIinfo = HDIofGrid( probMassVec )
# > show( HDIinfo )
<- sort(probMassVec, decreasing = TRUE)
sortedProbMass <- min(which(cumsum(sortedProbMass) >= credMass))
HDIheightIdx <- sortedProbMass[HDIheightIdx]
HDIheight <- sum(probMassVec[probMassVec >= HDIheight])
HDImass
return(list(indices = which(probMassVec >= HDIheight),
mass = HDImass, height = HDIheight))
}
```

I found Kruschke’s description of his `HDIofGrid()`

a bit opaque. Happily, we can understand this function with a little help from an example posted at https://rdrr.io/github/kyusque/DBDA2E-utilities/man/HDIofGrid.html.

```
<- dbeta(seq(0, 1, length = 201), 30, 12)
prob_density_vec <- prob_density_vec / sum(prob_density_vec)
prob_mass_vec <- HDIofGrid(prob_mass_vec)
HDI_info
show(HDI_info)
```

```
## $indices
## [1] 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
## [20] 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
## [39] 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
##
## $mass
## [1] 0.9528232
##
## $height
## [1] 0.004448336
```

To walk that through a bit, `prob_density_vec`

is a vector of density values for \(\operatorname{Beta} (30, 12)\) based on 201 evenly-spaced values spanning the parameter space for \(\theta\) (i.e., from 0 to 1). In the second line, we converted those density values to the probability metric by dividing each by their sum, which we then saved as `prob_mass_vec`

. In the third line we shoved those probability values into Kruschke’s `HDIofGrid()`

and saved the results as `HDI_info`

. The output of the fourth line, `show(HDI_info)`

, showed us the results (i.e., the contents of `HDI_info`

).

As to those results, the values in saved as `$indices`

are the row numbers for all cases in `prob_mass_vec`

that were within the HDI. The value in `$mass`

showed the actual width of the HDI. Because we’re only working with finite samples (i.e., `length = 201`

), we won’t likely get a perfect 95% HDI. The value in `$height`

is the density value for *the waterline that defines the highest density region*. A plot might make that less abstract.

```
library(ggthemes)
theme_set(
theme_base() +
theme(plot.background = element_rect(color = "transparent"))
)
# wrangle
tibble(row = 1:length(prob_density_vec),
theta = seq(0, 1, length = length(prob_density_vec)),
density = prob_mass_vec,
cred = if_else(row %in% HDI_info$indices, 1, 0)) %>%
# plot
ggplot(aes(x = theta, y = density)) +
# HDI
geom_area(aes(fill = cred == 1)) +
# density line
geom_line(color = "skyblue", size = 1) +
# waterline
geom_hline(yintercept = HDI_info$height, linetype = 3, size = 1/4) +
# fluff
annotate(geom = "text", x = .2, y = 0.0046,
label = '"waterline" that defines all points\ninside the highest density region') +
annotate(geom = "text", x = .715, y = 0.01,
label = "95.28% HDI", color = "black", size = 5) +
scale_fill_manual(values = c("transparent", "grey67")) +
xlab(expression(theta))
```

You might have noticed our `theme_set()`

lines at the top of that code block. For the plots in this chapter, we’ll give a nod to the source material and make them with a similar aesthetic to Kruschke’s figures.

### 25.2.2 HDI of unimodal distribution is shortest interval.

The algorithms [in the next sections] find the HDI by searching among candidate intervals of mass \(M\). The shortest one found is declared to be the HDI. It is an approximation, of course. See Chen & Shao (1999) for more details, and Chen, He, Shao, and Xu (2003) for dealing with the unusual situation of multimodal distributions. (p. 727)

### 25.2.3 R code for computing HDI of a MCMC sample.

In this section, Kruschke provided the code for his `HDIofMCMC()`

function. We recreate it, below, with a few mild formatting changes. The `ggthemes::theme_base()`

theme gets us most of the way there.

```
<- function(sampleVec, credMass = .95) {
HDIofMCMC
# Computes highest density interval from a sample of representative values,
# estimated as shortest credible interval.
# Arguments:
# sampleVec
# is a vector of representative values from a probability distribution.
# credMass
# is a scalar between 0 and 1, indicating the mass within the credible
# interval that is to be estimated.
# Value:
# HDIlim is a vector containing the limits of the HDI
<- sort(sampleVec)
sortedPts <- ceiling(credMass * length(sortedPts))
ciIdxInc <- length(sortedPts) - ciIdxInc
nCIs <- rep(0, nCIs)
ciWidth for (i in 1:nCIs) {
<- sortedPts[i + ciIdxInc] - sortedPts[i]
ciWidth[i]
}<- sortedPts[which.min(ciWidth)]
HDImin <- sortedPts[which.min(ciWidth) + ciIdxInc]
HDImax <- c(HDImin, HDImax)
HDIlim
return(HDIlim)
}
```

Let’s continue working with `fit24.1`

to see how Kruschke’s `HDIofMCMC()`

works. First we need to extract the posterior draws.

`<- as_draws_df(fit24.1) draws `

Here’s how you might use the function to get the HDIs for the intercept parameter.

`HDIofMCMC(draws$b_Intercept)`

`## [1] -0.06511808 6.08724849`

Kruschke’s `HDIofMCMC()`

works very much the same as the summary functions from **tidybayes**. For example, here’s good old `tidybayes::mode_hdi()`

.

```
library(tidybayes)
mode_hdi(draws$b_Intercept)
```

```
## y ymin ymax .width .point .interval
## 1 3.207135 -0.06511808 6.087248 0.95 mode hdi
```

If you’d like to use **tidybayes** to just pull the HDIs without the extra information, just use the `hdi()`

function.

`hdi(draws$b_Intercept)`

```
## [,1] [,2]
## [1,] -0.06511808 6.087248
```

Just in case you’re curious, Kruschke’s `HDIofMCMC()`

function returns the same information as `tidybayes::hdi()`

. Let’s confirm.

`HDIofMCMC(draws$b_Intercept) == hdi(draws$b_Intercept)`

```
## [,1] [,2]
## [1,] TRUE TRUE
```

Identical.

### 25.2.4 R code for computing HDI of a function.

The function described in this section finds the HDI of a unimodal probability density function that is specified mathematically in R. For example, the function can find HDI’s of normal densities or of beta densities or of gamma densities, because those densities are specified as functions in R. (p. 728)

If you recall, we’ve been using this function off and on since Chapter 4. Here is it, again, with mildly reformatted code and parameter names.

```
<- function(name, width = .95, tol = 1e-8, ... ) {
hdi_of_icdf
# Arguments:
# `name` is R's name for the inverse cumulative density function
# of the distribution.
# `width` is the desired mass of the HDI region.
# `tol` is passed to R's optimize function.
# Return value:
# Highest density iterval (HDI) limits in a vector.
# Example of use: For determining HDI of a beta(30, 12) distribution, type
# `hdi_of_icdf(qbeta, shape1 = 30, shape2 = 12)`
# Notice that the parameters of the `name` must be explicitly stated;
# e.g., `hdi_of_icdf(qbeta, 30, 12)` does not work.
# Adapted and corrected from Greg Snow's TeachingDemos package.
<- 1.0 - width
incredible_mass <- function(low_tail_prob, name, width, ...) {
interval_width name(width + low_tail_prob, ...) - name(low_tail_prob, ...)
}<- optimize(interval_width, c(0, incredible_mass),
opt_info name = name, width = width,
tol = tol, ...)
<- opt_info$minimum
hdi_lower_tail_prob
return(c(name(hdi_lower_tail_prob, ...),
name(width + hdi_lower_tail_prob, ...)))
}
```

Here’s how it works for the standard normal distribution.

`hdi_of_icdf(qnorm, mean = 0, sd = 1)`

`## [1] -1.959964 1.959964`

By default, it returns 95% HDIs. Here’s how it’d work if you wanted the 80% intervals for \(\operatorname{Beta}(2, 2)\).

`hdi_of_icdf(qbeta, shape1 = 2, shape2 = 2, width = .8)`

`## [1] 0.1958001 0.8041999`

## 25.3 Reparameterization

There are situations in which one parameterization is intuitive to express a distribution, but a different parameterization is required for mathematical convenience. For example, we may think intuitively of the standard deviation of a normal distribution, but have to parameterize the distribution in terms of the precision (i.e., reciprocal of the variance). (p. 729)

The details in the rest of this section are beyond the scope of this project.

## 25.4 Censored data in ~~JAGS~~ brms

“In many situations some data are censored, which means that their values are known only within a certain range” (p. 732) Happily, **brms** is capable of handling censored variables. The setup is a little different from how Kruschke described for JAGS. From the `brmsformula`

section of the **brms** reference manual (Bürkner, 2022f), we read:

With the exception of categorical, ordinal, and mixture families, left, right, and interval censoring can be modeled through

`y | cens(censored) ~ predictors`

. The censoring variable (named`censored`

in this example) should contain the values`'left'`

,`'none'`

,`'right'`

, and`'interval'`

(or equivalently`-1`

,`0`

,`1`

, and`2`

) to indicate that the corresponding observation is left censored, not censored, right censored, or interval censored. For interval censored data, a second variable (let’s call it`y2`

) has to be passed to`cens`

. In this case, the formula has the structure`y | cens(censored,y2) ~ predictors`

. While the lower bounds are given in`y`

, the upper bounds are given in`y2`

for interval censored data. Intervals are assumed to be open on the left and closed on the right:`(y,y2]`

.

We’ll make sense of all this in just a moment. First, let’s see how Kruschke described the example in the text.

To illustrate why it is important to include censored data in the analysis, consider a case in which \(N = 500\) values are generated randomly from a normal distribution with \(\mu = 100\) and \(\sigma = 15\). Suppose that values above \(106\) are censored, as are values in the interval between \(94\) and \(100\). For the censored values, all we know is the interval in which they occurred, but not their exact value. (p. 732)

I’m not aware that we have access to Kruschke’s censored data, so we’ll just make our own based on his description. We’ll start off by simulating the idealized uncensored data, `y`

, based on \(\operatorname{Normal}(100, 15)\).

```
<- 500
n
set.seed(25)
<- tibble(y = rnorm(n, mean = 100, sd = 15)) d
```

To repeat, Kruschke described two kinds of censoring:

- “values above 106 are censored,”
- “as are values in the interval between 94 and 100.”

⚠️ Rather than examine both kinds of censoring at once, as in the text, I’m going to slow down and break this section up. First, we’ll explore right censoring. Second, we’ll explore interval censoring. For our grand finale, we’ll combine the two, as in the text.

### 25.4.1 Right censoring: When you’re uncertain beyond a single threshold.

In the case where the “values above 106 are censored,” we need to save a single threshold value. We’ll call it `t1`

.

`<- 106 t1 `

Now we can use our `t1`

threshold value to make a right censored version of the `y`

column called `y1`

. We’ll also make a character value listing out a row’s censoring status in nominal terms.

```
<-
d %>%
d mutate(y1 = if_else(y > t1, t1, y),
cen1 = if_else(y > t1, "right", "none"))
d
```

```
## # A tibble: 500 × 3
## y y1 cen1
## <dbl> <dbl> <chr>
## 1 96.8 96.8 none
## 2 84.4 84.4 none
## 3 82.7 82.7 none
## 4 105. 105. none
## 5 77.5 77.5 none
## 6 93.3 93.3 none
## 7 126. 106 right
## 8 108. 106 right
## 9 101. 101. none
## 10 99.1 99.1 none
## # … with 490 more rows
```

When the values in `y1`

are not censored, we see `cen1 == "none"`

. For all cases where the original `y`

value exceeded the `t1`

threshold (i.e., 106), `y1 == 106`

and `cen1 == "right"`

. We are using the terms `"none"`

and `"right"`

because they are based on the block quote from the **brms** reference manual. For plotting purposes, we’ll make a new variable, `y_na`

, that only has values for which `cen1 == "none"`

.

```
<-
d %>%
d mutate(y_na = ifelse(cen1 == "none", y, NA))
d
```

```
## # A tibble: 500 × 4
## y y1 cen1 y_na
## <dbl> <dbl> <chr> <dbl>
## 1 96.8 96.8 none 96.8
## 2 84.4 84.4 none 84.4
## 3 82.7 82.7 none 82.7
## 4 105. 105. none 105.
## 5 77.5 77.5 none 77.5
## 6 93.3 93.3 none 93.3
## 7 126. 106 right NA
## 8 108. 106 right NA
## 9 101. 101. none 101.
## 10 99.1 99.1 none 99.1
## # … with 490 more rows
```

Now let’s get a better sense of the data with a few histograms.

```
%>%
d pivot_longer(-cen1) %>%
mutate(name = factor(name, levels = c("y", "y1", "y_na"))) %>%
ggplot(aes(x = value)) +
geom_histogram(size = .25, binwidth = 2, boundary = 106,
fill = "skyblue4", color = "white") +
xlab(NULL) +
scale_y_continuous(NULL, breaks = NULL) +
facet_wrap(~ name, ncol = 3)
```

The original un-censored data are in `y`

. The censored data, for which all values above 106 were simply recoded 106, is `y1`

. The `y_na`

column shows what the distribution would look like if the censored values were omitted.

Here’s how we might fit a model which only uses the uncensored values, those in `y_na`

.

```
# define the stanvars
<- mean(d$y_na, na.rm = T)
mean_y <- sd(d$y_na, na.rm = T)
sd_y
<-
stanvars stanvar(mean_y, name = "mean_y") +
stanvar(sd_y, name = "sd_y")
# fit the model
.1a <-
fit25brm(data = d,
family = gaussian,
~ 1,
y_na prior = c(prior(normal(mean_y, sd_y * 100), class = Intercept),
prior(normal(0, sd_y), class = sigma)),
chains = 4, cores = 4,
stanvars = stanvars,
file = "fits/fit25.01a")
```

Check the summary for the naïve model.

`print(fit25.1a)`

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y_na ~ 1
## Data: d (Number of observations: 333)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 91.43 0.51 90.40 92.41 1.00 3410 2678
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 9.39 0.36 8.70 10.13 1.00 3708 2881
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

Relative to the true data-generating process for the original variable `y`

, \(\operatorname{Normal}(100, 15)\), those parameters look pretty biased. Now let’s practice fitting our first censored model. Here we use the `y1 | cens(cen1)`

syntax in the left side of the model `formula`

to indicate that the criterion variable, `y1`

, has been censored according to the process as defined in the `cens()`

function. Within `cens()`

we inserted our nominal variable, `cens1`

, which indicates which cases are censored (`"right"`

) or not (`"none"`

).

This model is one of the rare occasions where we’ll set our initial values for the model intercept. In my first few attempts, `brm()`

had great difficulty initializing the chains using the default initial values. We’ll help it out by setting them at `mean_y`

. Recall that when you set custom initial values in **brms**, you save them in a list with the number of lists equaling the number of HMC chains. Because we’re using the default `chains = 4`

, well need four lists of intercept start values, `mean_y`

. You can set them to different values, if you’d like.

```
<- list(Intercept = mean_y)
inits
<- list(inits, inits, inits, inits)
inits_list
.2a <-
fit25brm(data = d,
family = gaussian,
| cens(cen1) ~ 1,
y1 prior = c(prior(normal(mean_y, sd_y * 100), class = Intercept),
prior(normal(0, sd_y), class = sigma)),
chains = 4, cores = 4,
stanvars = stanvars,
# here we insert our start values for the intercept
inits = inits_list,
file = "fits/fit25.02a")
```

Now check the summary for the model accounting for the censoring.

`print(fit25.2a)`

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y1 | cens(cen1) ~ 1
## Data: d (Number of observations: 500)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 99.43 0.71 98.09 100.90 1.00 2939 2197
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 14.29 0.61 13.17 15.56 1.00 2918 2440
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

To get a better sense of what these models are doing, we might do a posterior predictive check with `predict()`

.

```
# original
<-
p1 %>%
d pivot_longer(cols = c(y, y_na)) %>%
ggplot(aes(x = value)) +
geom_histogram(size = .25, binwidth = 2, boundary = 106,
fill = "skyblue4", color = "white") +
labs(subtitle = "Data",
x = NULL) +
scale_x_continuous(breaks = 3:5 * 25) +
scale_y_continuous(NULL, breaks = NULL) +
facet_wrap(~ name, ncol = 1)
# posterior simulations
set.seed(25)
<-
p2 predict(fit25.1a,
summary = F,
ndraws = 4) %>%
data.frame() %>%
mutate(sim = 1:n()) %>%
pivot_longer(-sim) %>%
ggplot(aes(x = value)) +
geom_histogram(size = .25, binwidth = 2, boundary = 106,
fill = "skyblue", color = "white") +
geom_vline(xintercept = t1, linetype = 2) +
labs(subtitle = "Simulations",
x = NULL) +
scale_x_continuous(breaks = 3:5 * 25) +
scale_y_continuous(NULL, breaks = NULL) +
coord_cartesian(xlim = range(d$y)) +
facet_wrap(~ sim, ncol = 2, labeller = label_both)
# combine
library(patchwork)
+ p2 +
p1 plot_layout(widths = c(1, 2)) +
plot_annotation(title = "PP check for the naïve model, fit25.1a")
```

On the left column, the have the original uncensored data `y`

and the `y_na`

data where the censored values were removed. The four lighter histograms on the right are individual simulations based on our naïve model, `fit25.1a`

, which only used the `y_na`

data. In all cases, the simulated data are biased to be too small compared with the original data, `y`

. They also have the undesirable quality that they don’t even match up with the `y_na`

data, which drops sharply off at the 106 threshold, whereas the simulated data all contain values exceeding that point (depicted by the dashed vertical lines). Now let’s see what happens when we execute a few posterior predictive checks with our first censored model, `fit25.2a`

.

```
<-
p1 %>%
d pivot_longer(cols = c(y, y1)) %>%
ggplot(aes(x = value)) +
geom_histogram(size = .25, binwidth = 2, boundary = 106,
fill = "skyblue4", color = "white") +
labs(subtitle = "Data",
x = NULL) +
scale_x_continuous(breaks = 3:5 * 25) +
scale_y_continuous(NULL, breaks = NULL) +
facet_wrap(~ name, ncol = 1, scales = "free_y")
set.seed(25)
<-
p2 predict(fit25.2a,
summary = F,
ndraws = 4) %>%
data.frame() %>%
mutate(sim = 1:n()) %>%
pivot_longer(-sim) %>%
ggplot(aes(x = value)) +
geom_histogram(size = .25, binwidth = 2, boundary = 106,
fill = "skyblue", color = "white") +
labs(subtitle = "Simulations",
x = NULL) +
scale_x_continuous(breaks = 3:5 * 25) +
scale_y_continuous(NULL, breaks = NULL) +
coord_cartesian(xlim = range(d$y)) +
facet_wrap(~ sim, ncol = 2, labeller = label_both)
# combine
+ p2 +
p1 plot_layout(widths = c(1, 2)) +
plot_annotation(title = "PP check for the censored model, fit25.2a")
```

Now the posterior predictive draws match up nicely from the original uncensored data `y`

in terms of central tendency, dispersion, and their overall shapes. This is the magic of a model that accounts for right (or left) censoring. It can take all those stacked-up 106 values and approximate the underlying distribution, had they been more accurately recorded.

### 25.4.2 Interval censoring: When your values are somewhere in the middle.

Time to move onto interval censoring. Kruschke’s example in the text included values censored in the interval between 94 and 100. This will require two more thresholds, which we’ll call `t2`

and `t3`

.

```
<- 94
t2 <- 100 t3
```

We should revisit part of the block quote from the **brms** reference manual, above:

For interval censored data, a second variable (let’s call it

`y2`

) has to be passed to`cens`

. In this case, the formula has the structure`y | cens(censored,y2) ~ predictors`

. While the lower bounds are given in`y`

, the upper bounds are given in`y2`

for interval censored data. Intervals are assumed to be open on the left and closed on the right:`(y,y2]`

It’s a little unclear, to me, if this is how Kruschke defined his intervals, but since we’re working with **brms** we’ll just use this convention. Thus, we will define “values in the interval between 94 and 100” as `y >= t2 & y < t3`

. Here we’ll follow the **brms** convention and define a new variable `y2`

for our lower bound and `y3`

for our upper bound. Additionally, the new `cen2`

variable will tell us whether a case is interval censored (`"interval"`

) or not (`"none"`

). We have also redefined `y_na`

in terms of our `cen2`

variable.

```
<-
d %>%
d mutate(y2 = if_else(y >= t2 & y < t3, t2, y),
y3 = if_else(y >= t2 & y < t3, t3, y),
cen2 = if_else(y >= t2 & y < t3, "interval", "none")) %>%
mutate(y_na = ifelse(cen2 == "none", y, NA))
d
```

```
## # A tibble: 500 × 7
## y y1 cen1 y_na y2 y3 cen2
## <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <chr>
## 1 96.8 96.8 none NA 94 100 interval
## 2 84.4 84.4 none 84.4 84.4 84.4 none
## 3 82.7 82.7 none 82.7 82.7 82.7 none
## 4 105. 105. none 105. 105. 105. none
## 5 77.5 77.5 none 77.5 77.5 77.5 none
## 6 93.3 93.3 none 93.3 93.3 93.3 none
## 7 126. 106 right 126. 126. 126. none
## 8 108. 106 right 108. 108. 108. none
## 9 101. 101. none 101. 101. 101. none
## 10 99.1 99.1 none NA 94 100 interval
## # … with 490 more rows
```

Now let’s get a better sense of the new data with a few histograms.

```
%>%
d pivot_longer(cols = c(y, y2, y3, y_na)) %>%
mutate(name = factor(name, levels = c("y", "y2", "y3", "y_na"))) %>%
ggplot(aes(x = value)) +
geom_histogram(size = .25, binwidth = 2, boundary = 106,
fill = "skyblue4", color = "white") +
labs(subtitle = "Our data have been updated",
x = NULL) +
scale_y_continuous(NULL, breaks = NULL) +
facet_wrap(~ name, ncol = 2)
```

Now we fit two models. The first, `fit25.1b`

, will only used the values not interval censored, `y_na`

. The second model, `fit25.2b`

, will use the `cens()`

function to fit account for the interval censored data with the `y2`

, `y3`

, and `cen2`

columns. As before, we’ll want to use `inits`

to help the censored model run smoothly.

```
# naïve
.1b <-
fit25brm(data = d,
family = gaussian,
~ 1,
y_na prior = c(prior(normal(mean_y, sd_y * 100), class = Intercept),
prior(normal(0, sd_y), class = sigma)),
chains = 4, cores = 4,
stanvars = stanvars,
file = "fits/fit25.01b")
# censored
.2b <-
fit25brm(data = d,
family = gaussian,
| cens(cen2, y3) ~ 1,
y2 prior = c(prior(normal(mean_y, sd_y * 100), class = Intercept),
prior(normal(0, sd_y), class = sigma)),
chains = 4, cores = 4,
stanvars = stanvars,
# here we insert our start values for the intercept
inits = inits_list,
file = "fits/fit25.02b")
```

Review the model summaries.

`print(fit25.1b)`

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y_na ~ 1
## Data: d (Number of observations: 424)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 100.24 0.77 98.74 101.73 1.00 3302 2636
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 15.79 0.53 14.81 16.87 1.00 3801 2459
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

`print(fit25.2b)`

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y2 | cens(cen2, y3) ~ 1
## Data: d (Number of observations: 500)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 99.75 0.65 98.49 101.07 1.00 3475 2566
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 14.61 0.46 13.74 15.56 1.00 2822 2343
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

In this case, both models did a reasonable job approximating the original data-generating parameters, \(\operatorname{Normal}(100, 15)\). However, the naïve model `fit25.1b`

did so at a considerable loss of precision, as indicated by the posterior standard deviations. To give a better sense, here are the model parameters in an interval plot.

```
# wrangle
bind_rows(as_draws_df(fit25.1b),
as_draws_df(fit25.2b)) %>%
mutate(fit = rep(c("fit25.1b", "fit25.2b"), each = n() / 2)) %>%
pivot_longer(b_Intercept:sigma) %>%
mutate(parameter = if_else(name == "b_Intercept", "mu", "sigma")) %>%
# plot
ggplot(aes(x = value, y = fit)) +
stat_interval(.width = c(.5, .95)) +
scale_color_manual("ETI",
values = c("skyblue2", "skyblue4"),
labels = c("95%", "50%")) +
labs(x = "marginal posterior",
y = NULL) +
theme(axis.ticks.y = element_blank()) +
facet_wrap(~ parameter, scales = "free_x", labeller = label_parsed)
```

One could also argue the parameters for the censored model `fit25.2b`

were a little less biased. However, I’d be leery of making that as a general claim based on this example, alone.

### 25.4.3 Complex censoring: When the going gets tough…

Okay, let’s put what we’ve learned together and practice with data that are both right AND interval censored. This is going to require some tricky coding on our part. Since interval censoring is in play, we’ll need two `y`

variables, again. When the data are indeed interval censored, their lower and upper bounds will be depicted by `y4`

and `y5`

, respectively. When the data are right censored, `y4`

will contain the threshold 106 and `y5`

will contain the original value from `y`

. All this tricky information will be indexed in our nominal variable `cen3`

. Once again, we update our `y_na`

variable, too.

```
<-
d %>%
d mutate(y4 = if_else(y >= t2 & y < t3, t2,
if_else(y > t1, t1, y)),
y5 = if_else(y >= t2 & y < t3, t3, y),
cen3 = if_else(y >= t2 & y < t3, "interval",
if_else(y > t1, "right", "none"))) %>%
mutate(y_na = ifelse(cen3 == "none", y, NA))
d
```

```
## # A tibble: 500 × 10
## y y1 cen1 y_na y2 y3 cen2 y4 y5 cen3
## <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <chr>
## 1 96.8 96.8 none NA 94 100 interval 94 100 interval
## 2 84.4 84.4 none 84.4 84.4 84.4 none 84.4 84.4 none
## 3 82.7 82.7 none 82.7 82.7 82.7 none 82.7 82.7 none
## 4 105. 105. none 105. 105. 105. none 105. 105. none
## 5 77.5 77.5 none 77.5 77.5 77.5 none 77.5 77.5 none
## 6 93.3 93.3 none 93.3 93.3 93.3 none 93.3 93.3 none
## 7 126. 106 right NA 126. 126. none 106 126. right
## 8 108. 106 right NA 108. 108. none 106 108. right
## 9 101. 101. none 101. 101. 101. none 101. 101. none
## 10 99.1 99.1 none NA 94 100 interval 94 100 interval
## # … with 490 more rows
```

Explore the new data with histograms.

```
%>%
d pivot_longer(cols = c(y, y4, y5, y_na)) %>%
mutate(name = factor(name, levels = c("y", "y4", "y5", "y_na"))) %>%
ggplot(aes(x = value)) +
geom_histogram(size = .25, binwidth = 2, boundary = 106,
fill = "skyblue4", color = "white") +
labs(subtitle = "Our data have been updated, again",
x = NULL) +
scale_y_continuous(NULL, breaks = NULL) +
facet_wrap(~ name, ncol = 2)
```

In the text, Kruschke reported he had 255 uncensored values (p. 732). Here’s the breakdown of our data.

```
%>%
d count(cen3)
```

```
## # A tibble: 3 × 2
## cen3 n
## <chr> <int>
## 1 interval 76
## 2 none 257
## 3 right 167
```

We got really close! Now we fit our two models: `fit25.1c`

, which ignores the censoring problem and just drops those cases, and `fit25.2c`

, which makes slick use of the `cens()`

function.

```
# naïve
.1c <-
fit25brm(data = d,
family = gaussian,
~ 1,
y_na prior = c(prior(normal(mean_y, sd_y * 100), class = Intercept),
prior(normal(0, sd_y), class = sigma)),
chains = 4, cores = 4,
stanvars = stanvars,
file = "fits/fit25.01c")
# censored
.2c <-
fit25brm(data = d,
family = gaussian,
| cens(cen3, y5) ~ 1,
y4 prior = c(prior(normal(mean_y, sd_y * 100), class = Intercept),
prior(normal(0, sd_y), class = sigma)),
chains = 4, cores = 4,
stanvars = stanvars,
# here we insert our start values for the intercept
inits = inits_list,
file = "fits/fit25.02c")
```

Compare the model summaries.

`print(fit25.1c)`

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y_na ~ 1
## Data: d (Number of observations: 257)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 89.87 0.64 88.62 91.12 1.00 3532 2900
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 10.12 0.44 9.30 11.01 1.00 4041 2475
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

`print(fit25.2c)`

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y4 | cens(cen3, y5) ~ 1
## Data: d (Number of observations: 500)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 99.48 0.70 98.13 100.89 1.00 3074 2240
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 14.27 0.58 13.21 15.45 1.00 2792 2545
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

The results for the naïve model are a mess and those from the censored model look great! Also note how the latter used all avaliable cases, whereas the former dropped nearly half of them.

Before we can make our version of Figure 25.4, we’ll need to extract the posterior draws. We’ll start with `fit25.1c`

.

```
<-
draws as_draws_df(fit25.1c) %>%
mutate(mu = b_Intercept,
`(mu-100)/sigma` = (b_Intercept - 100) / sigma)
head(draws)
```

```
## # A draws_df: 6 iterations, 1 chains, and 6 variables
## b_Intercept sigma lprior lp__ mu (mu-100)/sigma
## 1 89 9.9 -11 -968 89 -1.10
## 2 90 10.6 -11 -968 90 -0.98
## 3 91 9.7 -11 -968 91 -0.98
## 4 90 10.3 -11 -967 90 -0.99
## 5 90 9.8 -11 -968 90 -0.97
## 6 90 9.6 -11 -968 90 -1.08
## # ... hidden reserved variables {'.chain', '.iteration', '.draw'}
```

These subplots look a lot like those from back in Section 16.2. Since this is the last plot from the book, it seems like we should go all out. To reduce some of the code redundancy with the six subplots of the marginal posteriors, we’ll make a custom geom, `geom_hist()`

.

```
<- function(xintercept = xintercept, binwidth = binwidth, ...) {
geom_hist
list(
geom_histogram(fill = "skyblue", color = "white", size = .2, binwidth = binwidth, boundary = 106),
geom_vline(xintercept = xintercept, color = "skyblue3", size = 1/2, linetype = 3),
stat_pointinterval(aes(y = 0), point_interval = mode_hdi, .width = .95),
scale_y_continuous(NULL, breaks = NULL)
)
}
```

Now we have our `geom_hist()`

, here are the first three histograms for the marginal posteriors from `fit25.1c`

.

```
<-
p1 %>%
draws ggplot(aes(x = mu)) +
geom_hist(xintercept = 100, binwidth = 0.1) +
xlab(expression(mu))
<-
p3 %>%
draws ggplot(aes(x = sigma)) +
geom_hist(xintercept = 15, binwidth = 0.08) +
xlab(expression(sigma))
<-
p4 %>%
draws ggplot(aes(x = `(mu-100)/sigma`)) +
geom_hist(xintercept = 0, binwidth = 0.015) +
xlab(expression((mu-100)/sigma))
```

The histogram of the censored data with the posterior predictive density curves superimposed will take a little more work.

```
<- 50
n_lines
<-
p2 %>%
draws slice(1:n_lines) %>%
expand(nesting(.draw, mu, sigma),
y_na = seq(from = 40, to = 120, by = 1)) %>%
mutate(density = dnorm(x = y_na, mean = mu, sd = sigma)) %>%
ggplot(aes(x = y_na)) +
geom_histogram(data = d,
aes(y = stat(density)),
size = .25, binwidth = 2, boundary = 106,
fill = "skyblue4", color = "white") +
geom_line(aes(y = density, group = .draw),
size = 1/4, alpha = 1/3, color = "skyblue") +
scale_x_continuous("y", limits = c(40, 110)) +
scale_y_continuous(NULL, breaks = NULL)
```

Now extract the posterior draws from our censored model, `fit25.2`

, and repeat the process.

```
<-
draws as_draws_df(fit25.2c) %>%
mutate(mu = b_Intercept,
`(mu-100)/sigma` = (b_Intercept - 100) / sigma)
<-
p5 %>%
draws ggplot(aes(x = mu)) +
geom_hist(xintercept = 100, binwidth = 0.1) +
xlab(expression(mu))
<-
p7 %>%
draws ggplot(aes(x = sigma)) +
geom_hist(xintercept = 15, binwidth = 0.1) +
xlab(expression(sigma))
<-
p8 %>%
draws ggplot(aes(x = `(mu-100)/sigma`)) +
geom_hist(xintercept = 0, binwidth = 0.01) +
xlab(expression((mu-100)/sigma))
<-
p6 %>%
draws slice(1:n_lines) %>%
expand(nesting(.draw, mu, sigma),
y_na = seq(from = 40, to = 120, by = 1)) %>%
mutate(density = dnorm(x = y_na, mean = mu, sd = sigma)) %>%
ggplot(aes(x = y_na)) +
geom_histogram(data = d,
aes(y = stat(density)),
size = .25, binwidth = 2, boundary = 106,
fill = "skyblue4", color = "white") +
geom_line(aes(y = density, group = .draw),
size = 1/4, alpha = 1/3, color = "skyblue") +
scale_x_continuous("y", limits = c(40, 110)) +
scale_y_continuous(NULL, breaks = NULL)
```

Finally, combine the subplots, and annotate a bit.

```
| p2) / (p3 | p4) / (p5 | p6) / (p7 | p8)) +
((p1 plot_annotation(title = "This is our final plot, together.",
caption = expression(atop(italic("Upper quartet")*": Censored data omitted from analysis; parameter estimates are too small. ", italic("Lower quartet")*": Censored data imputed in known bins; parameter estimates are accurate."))) &
theme(plot.caption = element_text(hjust = 0))
```

To learn more about censored data, check out the nice lecture by Gordon Fox, *Introduction to analysis of censored and truncated data*.

### 25.4.4 Bonus: Truncation.

Though Kruschke didn’t cover it, truncation is often confused with censoring. When you have actually collected all your values, but some are just less certain than you’d like, you have a censoring issue. That is, you can think of censoring as an issue with *measurement precision*. Truncation has to do with *data collection*. When there’s some level in your criterion variable that you haven’t collected data on, that’s a truncation issue. For example, imagine you wanted to predict scores on an IQ test with socioeconomic status (SES). There might be a relation, there. But now imagine you only collected data on people with an IQ above 110. Those IQ data are severely left truncated. If you only care about the relation of SES and IQ for those with high IQ scores, your data are fine. But if you want to make general statements across the full range of IQ values, standard regression methods will likely produce biased results. Happily, **brms** allows users to accommodate truncated criterion variables with the `trunc()`

function, which works in a similar way to the `cens()`

function. For details on fitting truncated regression models with `trunc()`

, check out the *Additional response information* subsection of the `brmsformula`

section of the **brms** reference manual. The end of the Fox lecture, above, briefly covers truncated regression. You can also find a brief vignette on truncation by the good folks at the UCLA Institute for Digital Research and Education.

## 25.5 What Next?

“If you have made it this far and you are looking for more, you might peruse posts at [Kruschke’s] blog, https://doingbayesiandataanalysis.blogspot.com/, and search there for topics that interest you.” In addition to the other references Kruschke mentioned, you might also check out McElreath’s *Statistical rethinking*, both first (2015) and second (2020) editions. Much like this ebook, I have recoded both editions of *Statistical rethinking* in a **bookdown** form (Kurz, 2021, 2020). Click here for the first and here for the second. You can find other pedagogical material at my academic blog, https://solomonkurz.netlify.com/post/. For assistance on **brms**-related issues, check out the **brms** section on the Stan forums at https://discourse.mc-stan.org/c/interfaces/brms/36. Nicenboim, Schad, and Vasishth also have a nice new (2021) ebook called *An introduction to Bayesian data analysis for cognitive science*, which highlights **brms** in a way that could compliment the material we have practiced together.

Happy modeling, friends!

## Session info

`sessionInfo()`

```
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.1.1 tidybayes_3.0.2 ggthemes_4.2.4 brms_2.17.0
## [5] Rcpp_1.0.8 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.8
## [9] purrr_0.3.4 readr_2.0.1 tidyr_1.2.0 tibble_3.1.6
## [13] ggplot2_3.3.5 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.4.1 plyr_1.8.6
## [4] igraph_1.2.6 svUnit_1.0.6 sp_1.4-5
## [7] splines_4.1.2 crosstalk_1.1.1 TH.data_1.0-10
## [10] rstantools_2.1.1 inline_0.3.19 digest_0.6.29
## [13] htmltools_0.5.2 rsconnect_0.8.24 fansi_1.0.3
## [16] magrittr_2.0.3 checkmate_2.0.0 tzdb_0.1.2
## [19] modelr_0.1.8 RcppParallel_5.1.4 matrixStats_0.61.0
## [22] xts_0.12.1 sandwich_3.0-1 prettyunits_1.1.1
## [25] colorspace_2.0-3 rvest_1.0.1 ggdist_3.1.1
## [28] haven_2.4.3 xfun_0.30 callr_3.7.0
## [31] crayon_1.5.1 jsonlite_1.8.0 lme4_1.1-27.1
## [34] survival_3.2-13 zoo_1.8-9 glue_1.6.2
## [37] gtable_0.3.0 emmeans_1.7.1-1 distributional_0.3.0
## [40] pkgbuild_1.2.0 rstan_2.21.5 abind_1.4-5
## [43] scales_1.1.1 mvtnorm_1.1-2 emo_0.0.0.9000
## [46] DBI_1.1.1 miniUI_0.1.1.1 xtable_1.8-4
## [49] HDInterval_0.2.2 stats4_4.1.2 StanHeaders_2.21.0-7
## [52] DT_0.19 htmlwidgets_1.5.3 httr_1.4.2
## [55] threejs_0.3.3 arrayhelpers_1.1-0 posterior_1.1.0.9000
## [58] ellipsis_0.3.2 pkgconfig_2.0.3 loo_2.5.1
## [61] farver_2.1.0 sass_0.4.1 dbplyr_2.1.1
## [64] utf8_1.2.2 labeling_0.4.2 tidyselect_1.1.2
## [67] rlang_1.0.2 reshape2_1.4.4 later_1.3.0
## [70] munsell_0.5.0 cellranger_1.1.0 tools_4.1.2
## [73] cli_3.2.0 generics_0.1.2 broom_0.7.12
## [76] ggridges_0.5.3 evaluate_0.15 fastmap_1.1.0
## [79] processx_3.5.2 knitr_1.38 fs_1.5.2
## [82] nlme_3.1-153 mime_0.11 projpred_2.0.2
## [85] xml2_1.3.2 compiler_4.1.2 bayesplot_1.9.0
## [88] shinythemes_1.2.0 rstudioapi_0.13 gamm4_0.2-6
## [91] reprex_2.0.1 bslib_0.3.1 stringi_1.7.6
## [94] highr_0.9 ps_1.6.0 Brobdingnag_1.2-6
## [97] lattice_0.20-45 Matrix_1.3-4 nloptr_1.2.2.2
## [100] markdown_1.1 shinyjs_2.0.0 tensorA_0.36.2
## [103] vctrs_0.4.0 pillar_1.7.0 lifecycle_1.0.1
## [106] jquerylib_0.1.4 bridgesampling_1.1-2 estimability_1.3
## [109] raster_3.4-13 httpuv_1.6.2 R6_2.5.1
## [112] bookdown_0.26 promises_1.2.0.1 gridExtra_2.3
## [115] codetools_0.2-18 boot_1.3-28 colourpicker_1.1.0
## [118] MASS_7.3-54 gtools_3.9.2 assertthat_0.2.1
## [121] withr_2.5.0 shinystan_2.5.0 multcomp_1.4-17
## [124] mgcv_1.8-38 parallel_4.1.2 hms_1.1.0
## [127] grid_4.1.2 coda_0.19-4 minqa_1.2.4
## [130] rmarkdown_2.13 shiny_1.6.0 lubridate_1.7.10
## [133] base64enc_0.1-3 dygraphs_1.1.1.6
```

### References

*Nature Human Behaviour*, 1–3. https://doi.org/10.1038/s41562-019-0807-z

*brms reference manual, Version 2.17.0*. https://CRAN.R-project.org/package=brms/brms.pdf

*Development of Modern Statistics and Related Topics: Vols. Volume 1*(pp. 38–52). World Scientific. https://doi.org/10.1142/9789812796707_0004

*Journal of Computational and Graphical Statistics*,

*8*(1), 69–92. https://doi.org/10.1080/10618600.1999.10474802

*Journal of the Royal Statistical Society: Series A (Statistics in Society)*,

*182*(2), 389–402. https://doi.org/10.1111/rssa.12378

*The American Statistician*,

*50*(2), 120–126. https://doi.org/10.1080/00031305.1996.10474359

*Collabra: Psychology*,

*4*(1), 1–15. https://doi.org/10.1525/collabra.158

*British Journal of Mathematical and Statistical Psychology*,

*66*(1), 45–56. https://doi.org/10.1111/j.2044-8317.2012.02063.x

*Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan*. Academic Press. https://sites.google.com/site/doingbayesiandataanalysis/

*Nature Human Behaviour*,

*5*(10), 1282–1291. https://doi.org/10.1038/s41562-021-01177-7

*Psychonomic Bulletin & Review*,

*25*(1), 178–206. https://doi.org/10.3758/s13423-016-1221-4

*Statistical rethinking with brms, ggplot2, and the tidyverse: Second Edition*(version 0.2.0). https://bookdown.org/content/4857/

*Statistical rethinking with brms, ggplot2, and the tidyverse*(version 1.2.0). https://doi.org/10.5281/zenodo.3693202

*The American Psychologist*,

*73*(2), 111–125. https://doi.org/10.1037/amp0000242

*Statistical rethinking: A Bayesian course with examples in R and Stan*(Second Edition). CRC Press. https://xcelab.net/rm/statistical-rethinking/

*Statistical rethinking: A Bayesian course with examples in R and Stan*. CRC press. https://xcelab.net/rm/statistical-rethinking/

*An introduction to Bayesian data analysis for cognitive science*. https://vasishth.github.io/bayescogsci/book/

*Behavior Research Methods*,

*48*(3), 1062–1069. https://doi.org/10.3758/s13428-015-0630-z

*R Markdown: The definitive guide*. Chapman and Hall/CRC. https://bookdown.org/yihui/rmarkdown/