Bayesian Correlations: Let’s Talk Options.

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 * s),        (s * s * r), (s * s * r),
(s * s * r), (s * s),        (s * s * r),
(s * s * r), (s * s * r), (s * s)),
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))

## # 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.11     0.12 1.00     3782     2599
## x_s           0.91      0.06     0.79     1.02 1.00     3847     2946
##
## Family Specific Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.42      0.04     0.35     0.52 1.00     3811     2729
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, 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)
##  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.5836596 0.1155676 0.3569717 0.8123137
fixef(f3)[2, ]
##   Estimate  Est.Error       Q2.5      Q97.5
## 0.31047431 0.13742697 0.03672921 0.57820500

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.001     0.060    -0.114     0.119
fixef(f2)["Intercept", ] %>% round(3)
##  Estimate Est.Error      Q2.5     Q97.5
##     0.002     0.117    -0.226     0.233
fixef(f3)["Intercept", ] %>% round(3)
##  Estimate Est.Error      Q2.5     Q97.5
##     0.000     0.134    -0.261     0.266

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     2390     2083
##
## Family Specific Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.42      0.04     0.35     0.51 1.00     2791     2916
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, 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.909      0.06 0.789 1.033
fixef(f5) %>% round(3)
##     Estimate Est.Error  Q2.5 Q97.5
## x_s    0.581     0.117 0.356 0.809
fixef(f6) %>% round(3)
##     Estimate Est.Error  Q2.5 Q97.5
## y_s    0.311     0.137 0.047 0.585

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 cbind().
• Place the cbind() 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).
f7 <-
brm(data = d,
family = gaussian,
cbind(x_s, y_s, z_s) ~ 1,
prior = c(prior(normal(0, 1), class = Intercept),
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.01      0.13    -0.27     0.24 1.00     2359     2462
## ys_Intercept    -0.01      0.13    -0.28     0.25 1.00     2538     2458
## zs_Intercept    -0.00      0.14    -0.28     0.28 1.00     3096     2643
##
## 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     2273     2300
## sigma_ys     1.00      0.10     0.83     1.23 1.00     2294     2296
## sigma_zs     1.02      0.10     0.84     1.25 1.00     3161     2751
##
## 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     2432     2598
## rescor(xs,zs)     0.55      0.09     0.35     0.71 1.00     3190     2990
## rescor(ys,zs)     0.25      0.12     0.01     0.48 1.00     2919     2686
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, 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.

f8 <-
brm(data = d,
family = gaussian,
cbind(x_s, y_s, z_s) ~ 0,
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.19 1.00     1906     1751
## sigma_ys     0.99      0.10     0.82     1.20 1.00     1948     1975
## sigma_zs     1.02      0.10     0.84     1.24 1.00     2842     2253
##
## 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     2399     2234
## rescor(xs,zs)     0.55      0.10     0.35     0.71 1.00     2631     2324
## rescor(ys,zs)     0.26      0.13     0.01     0.50 1.00     2567     1730
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, 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)
##  0

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

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

Other than the sigma ~ 0 syntax, the main thing to notice is we’ve wrapped the entire model formula into the bf() function. 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.91      0.02     0.87     0.93 1.00     2621     2710
## rescor(xs,zs)     0.57      0.07     0.42     0.69 1.00     3307     2727
## rescor(ys,zs)     0.29      0.09     0.11     0.47 1.00     3095     2393
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, 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,
cbind(x, y, z) ~ 1,
prior = c(prior(normal(0, 10), class = Intercept),
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.64      1.20     7.29    12.00 1.00     1864     2582
## y_Intercept    15.57      2.47    10.69    20.34 1.00     2079     2631
## z_Intercept    14.85      3.38     8.00    21.42 1.00     2780     2773
##
## Family Specific Parameters:
##         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_x     8.99      0.87     7.50    10.91 1.00     2108     2278
## sigma_y    18.24      1.81    15.18    22.05 1.00     2114     2240
## sigma_z    26.09      2.62    21.58    31.77 1.00     3099     2736
##
## Residual Correlations:
##             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## rescor(x,y)     0.89      0.03     0.83     0.94 1.00     2252     2314
## rescor(x,z)     0.54      0.09     0.34     0.70 1.00     3067     2982
## rescor(y,z)     0.24      0.12    -0.00     0.47 1.00     2826     2709
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, 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") +
geom_halfeyeh(.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 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
##
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## Random number generation:
##  RNG:     Mersenne-Twister
##  Normal:  Inversion
##  Sample:  Rounding
##
## locale:
##  en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
##  stats     graphics  grDevices utils     datasets  methods   base
##
## other attached packages:
##   tidybayes_1.1.0 GGally_1.4.0    forcats_0.4.0   stringr_1.4.0
##   dplyr_0.8.3     purrr_0.3.2     readr_1.3.1     tidyr_1.0.0
##   tibble_2.1.3    ggplot2_3.2.1   tidyverse_1.2.1 brms_2.10.0
##  Rcpp_1.0.2      mvtnorm_1.0-11
##
## loaded via a namespace (and not attached):
##    colorspace_1.4-1          ellipsis_0.2.0.1
##    ggridges_0.5.1            rsconnect_0.8.13
##    ggstance_0.3.2            markdown_1.0
##    base64enc_0.1-3           rstudioapi_0.10
##    rstan_2.19.2              svUnit_0.7-12
##   DT_0.7                    fansi_0.4.0
##   lubridate_1.7.4           xml2_1.2.0
##   bridgesampling_0.6-0      codetools_0.2-16
##   knitr_1.23                shinythemes_1.1.2
##   zeallot_0.1.0             bayesplot_1.7.0
##   jsonlite_1.6              broom_0.5.2
##   shiny_1.3.2               compiler_3.6.0
##   httr_1.4.0                backports_1.1.4
##   assertthat_0.2.1          Matrix_1.2-17
##   lazyeval_0.2.2            cli_1.1.0
##   later_0.8.0               htmltools_0.3.6
##   prettyunits_1.0.2         tools_3.6.0
##   igraph_1.2.4.1            coda_0.19-2
##   gtable_0.3.0              glue_1.3.1
##   reshape2_1.4.3            cellranger_1.1.0
##   vctrs_0.2.0               nlme_3.1-139
##   blogdown_0.14             crosstalk_1.0.0
##   xfun_0.8                  ps_1.3.0
##   rvest_0.3.4               mime_0.7
##   miniUI_0.1.1.1            lifecycle_0.1.0
##   gtools_3.8.1              zoo_1.8-6
##   scales_1.0.0              colourpicker_1.0
##   hms_0.4.2                 promises_1.0.1
##   Brobdingnag_1.2-6         parallel_3.6.0
##   inline_0.3.15             shinystan_2.5.0
##   RColorBrewer_1.1-2        yaml_2.2.0
##   gridExtra_2.3             loo_2.1.0
##   StanHeaders_2.18.1-10     reshape_0.8.8
##   stringi_1.4.3             dygraphs_1.1.1.6
##   pkgbuild_1.0.3            rlang_0.4.0
##   pkgconfig_2.0.2           matrixStats_0.54.0
##   evaluate_0.14             lattice_0.20-38
##   rstantools_1.5.1          htmlwidgets_1.3
##   labeling_0.3              processx_3.3.1
##   tidyselect_0.2.5          plyr_1.8.4
##   magrittr_1.5              bookdown_0.12
##   R6_2.4.0                  generics_0.0.2
##   pillar_1.4.2              haven_2.1.0
##   withr_2.1.2               xts_0.11-2
##   abind_1.4-5               modelr_0.1.4
##   crayon_1.3.4              arrayhelpers_1.0-20160527
##   utf8_1.1.4                rmarkdown_1.13
##   grid_3.6.0                readxl_1.3.1
##   callr_3.2.0               threejs_0.3.1
##   digest_0.6.20             xtable_1.8-4
##  httpuv_1.5.1              stats4_3.6.0
##  munsell_0.5.0             shinyjs_1.0