Bayesian Correlations: Let's Talk Options.

By A. Solomon Kurz

February 16, 2019

[edited Dec 11, 2022]

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(upper = list(continuous = wrap("cor", stars = FALSE)))

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 × 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) 
##   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     0.00      0.06    -0.13     0.12 1.00     3718     2422
## x_s           0.91      0.06     0.79     1.03 1.00     3756     2874
## 
## 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     3209     2812
## 
## 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 ‘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.5826984 0.1194922 0.3459339 0.8149767
fixef(f3)[2, ]
##   Estimate  Est.Error       Q2.5      Q97.5 
## 0.30874774 0.13877552 0.03055009 0.57645259

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) %>% 
  data.frame() %>% 
  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))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

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.062    -0.126     0.118
fixef(f2)["Intercept", ] %>% round(3)
##  Estimate Est.Error      Q2.5     Q97.5 
##    -0.001     0.119    -0.234     0.234
fixef(f3)["Intercept", ] %>% round(3)
##  Estimate Est.Error      Q2.5     Q97.5 
##    -0.001     0.137    -0.271     0.269

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) 
##   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
## x_s     0.91      0.06     0.79     1.02 1.00     3286     2475
## 
## 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     3499     2542
## 
## 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).

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) %>% 
  data.frame() %>% 
  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.059 0.789 1.021
fixef(f5) %>% round(3)
##     Estimate Est.Error  Q2.5 Q97.5
## x_s    0.579     0.118 0.341 0.808
fixef(f6) %>% round(3)
##     Estimate Est.Error  Q2.5 Q97.5
## y_s    0.312     0.141 0.044 0.597

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) 
##   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
## xs_Intercept    -0.00      0.13    -0.26     0.26 1.00     2026     2269
## ys_Intercept    -0.00      0.14    -0.27     0.27 1.00     2292     2784
## zs_Intercept    -0.00      0.14    -0.28     0.27 1.00     2966     2695
## 
## 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.20 1.00     1730     2325
## sigma_ys     1.00      0.10     0.82     1.22 1.00     1722     2815
## sigma_zs     1.03      0.11     0.85     1.26 1.00     3260     2682
## 
## 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     2099     2480
## rescor(xs,zs)     0.55      0.10     0.34     0.72 1.00     3089     2798
## rescor(ys,zs)     0.25      0.13    -0.02     0.49 1.00     2781     2438
## 
## 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).

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) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 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     1734     2118
## sigma_ys     0.99      0.10     0.83     1.21 1.00     1900     2102
## sigma_zs     1.02      0.11     0.83     1.25 1.00     2467     2189
## 
## 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     2026     2350
## rescor(xs,zs)     0.55      0.09     0.35     0.72 1.00     2935     2670
## rescor(ys,zs)     0.26      0.13     0.00     0.49 1.00     2509     2408
## 
## 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).

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) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 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.86     0.93 1.00     3003     2877
## rescor(xs,zs)     0.57      0.07     0.43     0.69 1.00     3543     3011
## rescor(ys,zs)     0.29      0.09     0.12     0.46 1.00     2950     2869
## 
## 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 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) 
##   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
## x_Intercept     9.65      1.22     7.18    11.99 1.00     2020     2106
## y_Intercept    15.62      2.49    10.69    20.41 1.00     2194     2113
## z_Intercept    14.77      3.47     7.79    21.35 1.00     3070     2977
## 
## Family Specific Parameters: 
##         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_x     8.98      0.87     7.47    10.92 1.00     2054     2469
## sigma_y    18.17      1.80    15.08    22.19 1.00     2401     2259
## sigma_z    26.17      2.56    21.77    31.70 1.00     3001     2697
## 
## 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     2212     2539
## rescor(x,z)     0.55      0.09     0.36     0.71 1.00     2935     2613
## rescor(y,z)     0.25      0.12     0.01     0.48 1.00     2625     2736
## 
## 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).

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) {
  as_draws_df(fit) %>% 
    select(starts_with("b_"), -contains("Intercept")) %>% 
    set_names("rho") 
}

get_rescor <- function(fit) {
  as_draws_df(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.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/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_3.0.2 GGally_2.1.2    forcats_0.5.1   stringr_1.4.1  
##  [5] dplyr_1.0.10    purrr_0.3.4     readr_2.1.2     tidyr_1.2.1    
##  [9] tibble_3.1.8    ggplot2_3.4.0   tidyverse_1.3.2 brms_2.18.0    
## [13] Rcpp_1.0.9      mvtnorm_1.1-3  
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.1         backports_1.4.1      plyr_1.8.7          
##   [4] igraph_1.3.4         svUnit_1.0.6         splines_4.2.0       
##   [7] crosstalk_1.2.0      TH.data_1.1-1        rstantools_2.2.0    
##  [10] inline_0.3.19        digest_0.6.30        htmltools_0.5.3     
##  [13] fansi_1.0.3          magrittr_2.0.3       checkmate_2.1.0     
##  [16] googlesheets4_1.0.1  tzdb_0.3.0           modelr_0.1.8        
##  [19] RcppParallel_5.1.5   matrixStats_0.62.0   xts_0.12.1          
##  [22] sandwich_3.0-2       prettyunits_1.1.1    colorspace_2.0-3    
##  [25] rvest_1.0.2          ggdist_3.2.0         haven_2.5.1         
##  [28] xfun_0.35            callr_3.7.3          crayon_1.5.2        
##  [31] jsonlite_1.8.3       lme4_1.1-31          survival_3.4-0      
##  [34] zoo_1.8-10           glue_1.6.2           gtable_0.3.1        
##  [37] gargle_1.2.0         emmeans_1.8.0        distributional_0.3.1
##  [40] pkgbuild_1.3.1       rstan_2.21.7         abind_1.4-5         
##  [43] scales_1.2.1         DBI_1.1.3            miniUI_0.1.1.1      
##  [46] xtable_1.8-4         stats4_4.2.0         StanHeaders_2.21.0-7
##  [49] DT_0.24              htmlwidgets_1.5.4    httr_1.4.4          
##  [52] threejs_0.3.3        arrayhelpers_1.1-0   RColorBrewer_1.1-3  
##  [55] posterior_1.3.1      ellipsis_0.3.2       reshape_0.8.9       
##  [58] pkgconfig_2.0.3      loo_2.5.1            farver_2.1.1        
##  [61] sass_0.4.2           dbplyr_2.2.1         utf8_1.2.2          
##  [64] labeling_0.4.2       tidyselect_1.1.2     rlang_1.0.6         
##  [67] reshape2_1.4.4       later_1.3.0          munsell_0.5.0       
##  [70] cellranger_1.1.0     tools_4.2.0          cachem_1.0.6        
##  [73] cli_3.4.1            generics_0.1.3       broom_1.0.1         
##  [76] ggridges_0.5.3       evaluate_0.18        fastmap_1.1.0       
##  [79] yaml_2.3.5           processx_3.8.0       knitr_1.40          
##  [82] fs_1.5.2             nlme_3.1-159         mime_0.12           
##  [85] projpred_2.2.1       xml2_1.3.3           compiler_4.2.0      
##  [88] bayesplot_1.9.0      shinythemes_1.2.0    rstudioapi_0.13     
##  [91] gamm4_0.2-6          reprex_2.0.2         bslib_0.4.0         
##  [94] stringi_1.7.8        highr_0.9            ps_1.7.2            
##  [97] blogdown_1.15        Brobdingnag_1.2-8    lattice_0.20-45     
## [100] Matrix_1.4-1         nloptr_2.0.3         markdown_1.1        
## [103] shinyjs_2.1.0        tensorA_0.36.2       vctrs_0.5.0         
## [106] pillar_1.8.1         lifecycle_1.0.3      jquerylib_0.1.4     
## [109] bridgesampling_1.1-2 estimability_1.4.1   httpuv_1.6.5        
## [112] R6_2.5.1             bookdown_0.28        promises_1.2.0.1    
## [115] gridExtra_2.3        codetools_0.2-18     boot_1.3-28         
## [118] colourpicker_1.1.1   MASS_7.3-58.1        gtools_3.9.3        
## [121] assertthat_0.2.1     withr_2.5.0          shinystan_2.6.0     
## [124] multcomp_1.4-20      mgcv_1.8-40          parallel_4.2.0      
## [127] hms_1.1.1            grid_4.2.0           coda_0.19-4         
## [130] minqa_1.2.5          rmarkdown_2.16       googledrive_2.0.0   
## [133] shiny_1.7.2          lubridate_1.8.0      base64enc_0.1-3     
## [136] dygraphs_1.1.1.6
Posted on:
February 16, 2019
Length:
19 minute read, 3848 words
Tags:
Bayesian brms R tutorial
See Also:
Yes, you can compute standardized regression coefficients with multiple imputation
Sum-score effect sizes for multilevel Bayesian cumulative probit models
Just use multilevel models for your pre/post RCT data