# Bayesian Correlations: Let's Talk Options.

[edited Apr 21, 2021]

## tl;dr

There’s more than one way to fit a Bayesian correlation in brms.

## Here’s the deal.

In the last post, we considered how we might estimate correlations when our data contain influential outlier values. Our big insight was that if we use variants of Student’s \(t\)-distribution as the likelihood rather than the conventional normal distribution, our correlation estimates were less influenced by those outliers. And we mainly did that as Bayesians using the brms package. Click here for a refresher.

Since the brms package is designed to fit regression models, it can be surprising when you discover it’s handy for correlations, too. In short, you can fit them using a few tricks based on the multivariate syntax.

Shortly after uploading the post, it occurred to me we had more options and it might be useful to walk through them a bit.

## I assume things.

For this post, I’m presuming you are vaguely familiar with linear regression–both univariate and multivariate–, have a little background with Bayesian statistics, and have used Paul Bürkner’s brms packge. As you might imagine, all code in is R, with a heavy use of the tidyverse.

## We need data.

First, we’ll load our main packages.

```
library(mvtnorm)
library(brms)
library(tidyverse)
```

We’ll use the mvtnorm package to simulate three positively correlated variables.

```
m <- c(10, 15, 20) # the means
s <- c(10, 20, 30) # the sigmas
r <- c(.9, .6, .3) # the correlations
# here's the variance/covariance matrix
v <-
matrix(c((s[1] * s[1]), (s[2] * s[1] * r[1]), (s[3] * s[1] * r[2]),
(s[2] * s[1] * r[1]), (s[2] * s[2]), (s[3] * s[2] * r[3]),
(s[3] * s[1] * r[2]), (s[3] * s[2] * r[3]), (s[3] * s[3])),
nrow = 3, ncol = 3)
# after setting our seed, we're ready to simulate with `rmvnorm()`
set.seed(1)
d <-
rmvnorm(n = 50, mean = m, sigma = v) %>%
as_tibble() %>%
set_names("x", "y", "z")
```

Our data look like so.

```
library(GGally)
theme_set(theme_gray() +
theme(panel.grid = element_blank()))
d %>%
ggpairs()
```

Do note the Pearson’s correlation coefficients in the upper triangle.

In order to exploit all the methods we’ll cover in this post, we need to standardize our data. Here we do so by hand using the typical formula

\[z_{x_i} = \frac{x_i - \overline x}{s_x}\]

where \(\overline x\) is the observed mean and \(s_x\) is the observed standard deviation.

```
d <-
d %>%
mutate(x_s = (x - mean(x)) / sd(x),
y_s = (y - mean(y)) / sd(y),
z_s = (z - mean(z)) / sd(z))
head(d)
```

```
## # A tibble: 6 x 6
## x y z x_s y_s z_s
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3.90 11.5 -6.90 -0.723 -0.308 -0.928
## 2 17.7 29.5 4.01 0.758 0.653 -0.512
## 3 20.4 33.8 41.5 1.05 0.886 0.917
## 4 20.3 42.1 34.8 1.04 1.33 0.663
## 5 -3.64 -26.8 43.5 -1.53 -2.36 0.994
## 6 13.9 17.3 47.6 0.347 0.00255 1.15
```

There are at least two broad ways to get correlations out of standardized data in brms. One way uses the typical univariate syntax. The other way is an extension of the multivariate `cbind()`

approach. Let’s start univariate.

And for a point of clarification, we’re presuming the Gaussian likelihood for all the examples in this post.

## Univariate

If you fit a simple univariate model with standardized data and a single predictor, the coefficient for the slope will be in a correlation-like metric. Happily, since the data are all standardized, it’s easy to use regularizing priors.

```
f1 <-
brm(data = d,
family = gaussian,
y_s ~ 1 + x_s,
prior = c(prior(normal(0, 1), class = Intercept),
prior(normal(0, 1), class = b),
prior(normal(0, 1), class = sigma)),
chains = 4, cores = 4,
seed = 1)
```

Take a look at the model summary.

`print(f1)`

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y_s ~ 1 + x_s
## Data: d (Number of observations: 50)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.00 0.06 -0.12 0.12 1.00 3689 2507
## x_s 0.91 0.06 0.79 1.03 1.00 3385 2456
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.42 0.05 0.35 0.52 1.00 3199 2787
##
## Samples were drawn 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 ‘Population-Level Effects’ has the summary information for our intercept and slope. Notice how our `x_s`

slope is the same as the Pearson’s correlation.

`cor(d$x, d$y)`

`## [1] 0.9119708`

Since this approach only yields one correlation at a time, we have to fit two more models to get the other two correlations. To do so with haste, we can use the `update()`

syntax.

```
f2 <-
update(f1,
newdata = d,
formula = z_s ~ 1 + x_s)
f3 <-
update(f2,
newdata = d,
formula = z_s ~ 1 + y_s)
```

With the `fixef()`

function, we can easily isolate the \(\beta\) estimates.

`fixef(f2)[2, ]`

```
## Estimate Est.Error Q2.5 Q97.5
## 0.5829389 0.1205587 0.3448574 0.8235006
```

`fixef(f3)[2, ]`

```
## Estimate Est.Error Q2.5 Q97.5
## 0.31107268 0.14237961 0.02815067 0.58747426
```

There’s another thing I’d like to point out. Plotting the model results will help make the point.

```
# define the predictor values you'd like the fitted values for
nd <- tibble(x_s = seq(from = -3, to = 3, length.out = d %>% nrow()))
# wrangle
fitted(f1,
newdata = nd) %>%
as_tibble() %>%
bind_cols(nd) %>%
# plot
ggplot(aes(x_s)) +
geom_vline(xintercept = 0, color = "white") +
geom_hline(yintercept = 0, color = "white") +
geom_point(data = d,
aes(y = y_s)) +
geom_smooth(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = "identity",
alpha = 1/4, size = 1/2) +
coord_cartesian(xlim = range(d$x_s),
ylim = range(d$y_s))
```

The blue line is the posterior mean and the surrounding gray ribbon depicts the 95% posterior interval. Notice how the data and their respective fitted lines pass through [0, 0]? This is a consequence of modeling standardized data. We should always expect the intercept of a model like this to be 0. Here are the intercept summaries for all three models.

`fixef(f1)["Intercept", ] %>% round(3)`

```
## Estimate Est.Error Q2.5 Q97.5
## 0.000 0.062 -0.123 0.121
```

`fixef(f2)["Intercept", ] %>% round(3)`

```
## Estimate Est.Error Q2.5 Q97.5
## 0.000 0.117 -0.234 0.230
```

`fixef(f3)["Intercept", ] %>% round(3)`

```
## Estimate Est.Error Q2.5 Q97.5
## 0.002 0.138 -0.267 0.270
```

Within simulation error, they’re all centered on zero. So instead of estimating the intercept, why not just bake that into the models? Here we refit the models by fixing the intercept for each to zero.

```
f4 <-
update(f1,
formula = y_s ~ 0 + x_s)
f5 <-
update(f4,
newdata = d,
formula = z_s ~ 0 + x_s)
f6 <-
update(f4,
newdata = d,
formula = z_s ~ 0 + y_s)
```

Let’s take a look at the summary for the first.

`print(f4)`

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: y_s ~ x_s - 1
## Data: d (Number of observations: 50)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## x_s 0.91 0.06 0.79 1.03 1.00 3389 2550
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.42 0.04 0.34 0.52 1.00 3421 2600
##
## Samples were drawn 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).
```

Even though it may have seemed like we substantially changed the models by fixing the intercepts to 0, the summaries are essentially the same as when we estimated the intercepts. Here we’ll confirm the summaries with a plot, like above.

```
# wrangle
fitted(f4,
newdata = nd) %>%
as_tibble() %>%
bind_cols(nd) %>%
# plot
ggplot(aes(x_s)) +
geom_vline(xintercept = 0, color = "white") +
geom_hline(yintercept = 0, color = "white") +
geom_point(data = d,
aes(y = y_s)) +
geom_smooth(aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
stat = "identity",
alpha = 1/4, size = 1/2) +
coord_cartesian(xlim = range(d$x_s),
ylim = range(d$y_s))
```

The difference is subtle. By fixing the intercepts at 0, we estimated the slopes (i.e., the correlations) with increased precision as demonstrated by the slightly smaller posterior standard deviations (i.e., the values in the ‘Est.Error’ columns).

Here are the correlation summaries for those last three models.

`fixef(f4) %>% round(3)`

```
## Estimate Est.Error Q2.5 Q97.5
## x_s 0.908 0.061 0.788 1.031
```

`fixef(f5) %>% round(3)`

```
## Estimate Est.Error Q2.5 Q97.5
## x_s 0.581 0.116 0.355 0.801
```

`fixef(f6) %>% round(3)`

```
## Estimate Est.Error Q2.5 Q97.5
## y_s 0.311 0.135 0.049 0.569
```

But anyway, you get the idea. If you want to estimate a correlation in brms using simple univariate syntax, just (a) standardize the data and (b) fit a univariate model with or without an intercept. The slop will be in a correlation-like metric.

## Let’s go multivariate.

If you don’t recall the steps to fit correlations in brms with the multivariate syntax, here they are:

- List the variables you’d like correlations for within
`mvbind()`

. - Place the
`mvbind()`

function within the left side of the model formula. - On the right side of the model formula, indicate you only want intercepts (i.e.,
`~ 1`

). - Wrap that whole formula within
`bf()`

. - Then use the
`+`

operator to append`set_rescor(TRUE)`

, which will ensure brms fits a model with residual correlations.

In addition, you you want to use non-default priors, you’ll want to use the `resp`

argument to specify which prior is associated with which criterion variable. Here’s what that all looks like:

```
f7 <-
brm(data = d,
family = gaussian,
bf(mvbind(x_s, y_s, z_s) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 1), class = Intercept, resp = xs),
prior(normal(0, 1), class = Intercept, resp = ys),
prior(normal(0, 1), class = Intercept, resp = zs),
prior(normal(1, 1), class = sigma, resp = xs),
prior(normal(1, 1), class = sigma, resp = ys),
prior(normal(1, 1), class = sigma, resp = zs),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4,
seed = 1)
```

Behold the summary.

`print(f7)`

```
## Family: MV(gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: x_s ~ 1
## y_s ~ 1
## z_s ~ 1
## Data: d (Number of observations: 50)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## xs_Intercept 0.00 0.14 -0.27 0.27 1.00 2031 2446
## ys_Intercept 0.00 0.14 -0.27 0.28 1.00 2288 2677
## zs_Intercept 0.00 0.15 -0.28 0.29 1.00 2697 2728
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_xs 0.99 0.10 0.82 1.19 1.00 1999 2074
## sigma_ys 1.00 0.10 0.83 1.21 1.00 2265 2463
## sigma_zs 1.03 0.11 0.85 1.26 1.00 3137 2530
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(xs,ys) 0.89 0.03 0.83 0.94 1.00 2458 2815
## rescor(xs,zs) 0.55 0.09 0.35 0.72 1.00 3039 2828
## rescor(ys,zs) 0.26 0.13 -0.00 0.49 1.00 2733 2530
##
## Samples were drawn 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).
```

Look at the ‘Residual Correlations:’ section at the bottom of the output. Since there are no predictors in the model, the residual correlations are just correlations. Now notice how the intercepts in this model are also hovering around 0, just like in our univariate models. Yep, we can fix those, too. We do this by changing our formula to `mvbind(x_s, y_s, z_s) ~ 0`

.

```
f8 <-
brm(data = d,
family = gaussian,
bf(mvbind(x_s, y_s, z_s) ~ 0) + set_rescor(TRUE),
prior = c(prior(normal(1, 1), class = sigma, resp = xs),
prior(normal(1, 1), class = sigma, resp = ys),
prior(normal(1, 1), class = sigma, resp = zs),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4,
seed = 1)
```

Without the intercepts, the rest of the model is the same within simulation variance.

`print(f8)`

```
## Family: MV(gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: x_s ~ 0
## y_s ~ 0
## z_s ~ 0
## Data: d (Number of observations: 50)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_xs 0.98 0.09 0.81 1.18 1.00 1967 2063
## sigma_ys 0.99 0.10 0.82 1.19 1.00 2204 1917
## sigma_zs 1.01 0.10 0.83 1.24 1.00 2690 2817
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(xs,ys) 0.90 0.03 0.83 0.94 1.00 2441 2328
## rescor(xs,zs) 0.55 0.09 0.35 0.71 1.00 2739 2438
## rescor(ys,zs) 0.26 0.12 0.01 0.48 1.00 2608 2468
##
## Samples were drawn 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).
```

If you wanna get silly, we can prune even further. Did you notice how the estimates for \(\sigma\) are all hovering around 1? Since we have no predictors, \(\sigma\) is just an estimate of the population standard deviation. And since we’re working with standardized data, the population standard deviation has to be 1. Any other estimate would be nonsensical. So why not fix it to 1?

With brms, we can fix those \(\sigma\)’s to 1 with a trick of the nonlinear distributional modeling syntax. Recall when you model \(\sigma\), the brms default is to actually model its log. As is turns out, the log of 1 is zero.

`log(1)`

`## [1] 0`

Here’s how to make use of that within `brm()`

.

```
f9 <-
brm(data = d,
family = gaussian,
bf(mvbind(x_s, y_s, z_s) ~ 0,
sigma ~ 0) +
set_rescor(TRUE),
prior(lkj(2), class = rescor),
chains = 4, cores = 4,
seed = 1)
```

Here are the results.

`print(f9)`

```
## Family: MV(gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = log
## mu = identity; sigma = log
## mu = identity; sigma = log
## Formula: x_s ~ 0
## sigma ~ 0
## y_s ~ 0
## sigma ~ 0
## z_s ~ 0
## sigma ~ 0
## Data: d (Number of observations: 50)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(xs,ys) 0.90 0.02 0.87 0.93 1.00 3719 3031
## rescor(xs,zs) 0.57 0.07 0.42 0.69 1.00 3047 2773
## rescor(ys,zs) 0.29 0.09 0.11 0.46 1.00 2839 2615
##
## Samples were drawn 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 correlations are the only things left in the model.

Just to be clear, the multivariate approach does not require standardized data. To demonstrate, here we refit `f7`

, but with the unstandardized variables. And, since we’re no longer in the standardized metric, we’ll be less certain with our priors.

```
f10 <-
brm(data = d,
family = gaussian,
bf(mvbind(x, y, z) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 10), class = Intercept, resp = x),
prior(normal(0, 10), class = Intercept, resp = y),
prior(normal(0, 10), class = Intercept, resp = z),
prior(student_t(3, 0, 10), class = sigma, resp = x),
prior(student_t(3, 0, 10), class = sigma, resp = y),
prior(student_t(3, 0, 10), class = sigma, resp = z),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4,
seed = 1)
```

See, the ‘rescor()’ results are about the same as with `f7`

.

`print(f10)`

```
## Family: MV(gaussian, gaussian, gaussian)
## Links: mu = identity; sigma = identity
## mu = identity; sigma = identity
## mu = identity; sigma = identity
## Formula: x ~ 1
## y ~ 1
## z ~ 1
## Data: d (Number of observations: 50)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## x_Intercept 9.63 1.22 7.14 12.02 1.00 1941 2016
## y_Intercept 15.60 2.49 10.63 20.30 1.00 2238 2441
## z_Intercept 14.71 3.58 7.65 21.51 1.00 3021 2316
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_x 8.93 0.84 7.44 10.70 1.00 2263 2515
## sigma_y 18.13 1.73 15.13 21.77 1.00 2553 2793
## sigma_z 26.15 2.58 21.78 31.92 1.00 2626 2198
##
## Residual Correlations:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(x,y) 0.89 0.03 0.82 0.94 1.00 2540 2630
## rescor(x,z) 0.54 0.09 0.34 0.70 1.00 3122 3224
## rescor(y,z) 0.24 0.12 -0.01 0.47 1.00 2689 2829
##
## Samples were drawn 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).
```

## It’s time to compare methods.

To recap, we’ve compared several ways to fit correlations in brms. Some of the methods were with univariate syntax, others were with the multivariate syntax. Some of the models had all free parameters, others included fixed intercepts and sigmas. Whereas all the univariate models required standardized data, the multivariate approach can work with unstandardized data, too.

Now it might be of help to compare the results from each of the methods to get a sense of which ones you might prefer. Before we do so, we’ll define a couple custom functions to streamline the data wrangling.

```
get_rho <- function(fit) {
posterior_samples(fit) %>%
select(starts_with("b_"), -contains("Intercept")) %>%
set_names("rho")
}
get_rescor <- function(fit) {
posterior_samples(fit) %>%
select(starts_with("rescor")) %>%
set_names("x with y", "x with z", "y with z") %>%
gather(label, rho) %>%
select(rho, label)
}
```

Now let’s put those functions to work and plot.

```
library(tidybayes)
# collect the posteriors from the univariate models
tibble(name = str_c("f", 1:6)) %>%
mutate(fit = map(name, get)) %>%
mutate(rho = map(fit, get_rho)) %>%
unnest(rho) %>%
mutate(predictor = rep(c("x", "x", "y"), each = 4000) %>% rep(., times = 2),
criterion = rep(c("y", "z", "z"), each = 4000) %>% rep(., times = 2)) %>%
mutate(label = str_c(predictor, " with ", criterion)) %>%
select(-c(predictor:criterion)) %>%
# add in the posteriors from the multivariate models
bind_rows(
tibble(name = str_c("f", 7:10)) %>%
mutate(fit = map(name, get)) %>%
mutate(post = map(fit, get_rescor)) %>%
unnest(post)
) %>%
# wrangle a bit just to make the y axis easier to understand
mutate(name = factor(name,
levels = c(str_c("f", 1:10)),
labels = c("1. standardized, univariate",
"2. standardized, univariate",
"3. standardized, univariate",
"4. standardized, univariate, fixed intercepts",
"5. standardized, univariate, fixed intercepts",
"6. standardized, univariate, fixed intercepts",
"7. standardized, multivariate, fixed intercepts",
"8. standardized, multivariate, fixed intercepts",
"9. standardized, multivariate, fixed intercepts/sigmas",
"10. unstandardized, multivariate"))) %>%
# plot
ggplot(aes(x = rho, y = name)) +
geom_vline(data = tibble(label = c("x with y", "x with z", "y with z"),
rho = r),
aes(xintercept = rho), color = "white") +
stat_halfeye(.width = .95, size = 5/4) +
scale_x_continuous(breaks = c(0, r)) +
labs(x = expression(rho),
y = NULL) +
coord_cartesian(0:1) +
theme(axis.ticks.y = element_blank(),
axis.text.y = element_text(hjust = 0)) +
facet_wrap(~ label, ncol = 3)
```

To my eye, a few patterns emerged. First, the point estimates were about the same across methods. Second, fixing the intercepts didn’t seem to effect things, much. But, third, it appears that fixing the sigmas in the multivariate models did narrow the posteriors a bit.

Fourth, and perhaps most importantly, notice how the posteriors for the multivariate models were more asymmetric when they approached 1. Hopefully this makes intuitive sense. Correlations are bound between -1 and 1. However, standardized regression coefficients are not so bound. Accordingly, notice how the posteriors from the univariate models stayed symmetric when approaching 1 and some of their right tails even crossed over 1. So while the univariate approach did a reasonable job capturing the correlation point estimates, their posteriors weren’t quite in a correlation metric. Alternately, the univariate approach did make it convenient to express the correlations with fitted regression lines in scatter plots.

Both univariate and multivariate approaches appear to have their strengths and weaknesses. Choose which methods seems most appropriate for your correlation needs.

Happy modeling.

`sessionInfo()`

```
## R version 4.0.4 (2021-02-15)
## 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.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/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] tidybayes_2.3.1 GGally_2.1.1 forcats_0.5.1 stringr_1.4.0
## [5] dplyr_1.0.5 purrr_0.3.4 readr_1.4.0 tidyr_1.1.3
## [9] tibble_3.1.0 ggplot2_3.3.3 tidyverse_1.3.0 brms_2.15.0
## [13] Rcpp_1.0.6 mvtnorm_1.1-1
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.2.1 plyr_1.8.6
## [4] igraph_1.2.6 svUnit_1.0.3 splines_4.0.4
## [7] crosstalk_1.1.0.1 TH.data_1.0-10 rstantools_2.1.1
## [10] inline_0.3.17 digest_0.6.27 htmltools_0.5.1.1
## [13] rsconnect_0.8.16 fansi_0.4.2 magrittr_2.0.1
## [16] modelr_0.1.8 RcppParallel_5.0.2 matrixStats_0.57.0
## [19] xts_0.12.1 sandwich_3.0-0 prettyunits_1.1.1
## [22] colorspace_2.0-0 rvest_0.3.6 ggdist_2.4.0.9000
## [25] haven_2.3.1 xfun_0.22 callr_3.5.1
## [28] crayon_1.4.1 jsonlite_1.7.2 lme4_1.1-25
## [31] survival_3.2-10 zoo_1.8-8 glue_1.4.2
## [34] gtable_0.3.0 emmeans_1.5.2-1 V8_3.4.0
## [37] distributional_0.2.2 pkgbuild_1.2.0 rstan_2.21.2
## [40] abind_1.4-5 scales_1.1.1 DBI_1.1.0
## [43] miniUI_0.1.1.1 xtable_1.8-4 stats4_4.0.4
## [46] StanHeaders_2.21.0-7 DT_0.16 htmlwidgets_1.5.2
## [49] httr_1.4.2 threejs_0.3.3 arrayhelpers_1.1-0
## [52] RColorBrewer_1.1-2 ellipsis_0.3.1 farver_2.0.3
## [55] reshape_0.8.8 pkgconfig_2.0.3 loo_2.4.1
## [58] dbplyr_2.0.0 utf8_1.1.4 labeling_0.4.2
## [61] tidyselect_1.1.0 rlang_0.4.10 reshape2_1.4.4
## [64] later_1.1.0.1 munsell_0.5.0 cellranger_1.1.0
## [67] tools_4.0.4 cli_2.3.1 generics_0.1.0
## [70] broom_0.7.5 ggridges_0.5.2 evaluate_0.14
## [73] fastmap_1.0.1 yaml_2.2.1 processx_3.4.5
## [76] knitr_1.31 fs_1.5.0 nlme_3.1-152
## [79] mime_0.10 projpred_2.0.2 xml2_1.3.2
## [82] compiler_4.0.4 bayesplot_1.8.0 shinythemes_1.1.2
## [85] rstudioapi_0.13 curl_4.3 gamm4_0.2-6
## [88] reprex_0.3.0 statmod_1.4.35 stringi_1.5.3
## [91] highr_0.8 ps_1.6.0 blogdown_1.3
## [94] Brobdingnag_1.2-6 lattice_0.20-41 Matrix_1.3-2
## [97] nloptr_1.2.2.2 markdown_1.1 shinyjs_2.0.0
## [100] vctrs_0.3.6 pillar_1.5.1 lifecycle_1.0.0
## [103] bridgesampling_1.0-0 estimability_1.3 httpuv_1.5.4
## [106] R6_2.5.0 bookdown_0.21 promises_1.1.1
## [109] gridExtra_2.3 codetools_0.2-18 boot_1.3-26
## [112] colourpicker_1.1.0 MASS_7.3-53 gtools_3.8.2
## [115] assertthat_0.2.1 withr_2.4.1 shinystan_2.5.0
## [118] multcomp_1.4-16 mgcv_1.8-33 parallel_4.0.4
## [121] hms_0.5.3 grid_4.0.4 coda_0.19-4
## [124] minqa_1.2.4 rmarkdown_2.7 shiny_1.5.0
## [127] lubridate_1.7.9.2 base64enc_0.1-3 dygraphs_1.1.1.6
```