Hierarchical brms model: interpretation of cor(Intercept, predictor)

Modeling
1

I’m trying to figure out the interpretation of a group level correlation named as follows in the model summary: cor(Intercept,year)

Could you please explain this simple words? For example, how should I simulate my data to see a positive or negative correlation.

Simulated data

6 groups, each next has higher mean “y” with higher SD, association between “y” and year is positive within each group

set.seed(1)
n = 7
df2 = bind_rows(
  tibble(y = rnorm(n, 10, 2), group = rep(times = n, "A")) %>% arrange(y) %>% mutate(year = seq(1, n, 1)), 
  tibble(y = rnorm(n, 20, 6), group = rep(times = n, "B")) %>% arrange(y) %>% mutate(year = seq(1, n, 1)),
  tibble(y = rnorm(n, 30, 10), group = rep(times = n, "C")) %>% arrange(y) %>% mutate(year = seq(1, n, 1)),
  tibble(y = rnorm(n, 40, 14), group = rep(times = n, "D")) %>% arrange(y) %>% mutate(year = seq(1, n, 1)),
  tibble(y = rnorm(n, 50, 16), group = rep(times = n, "E")) %>% arrange(y) %>% mutate(year = seq(1, n, 1)),
  tibble(y = rnorm(n, 60, 20), group = rep(times = n, "F")) %>% arrange(y) %>% mutate(year = seq(1, n, 1)))

Model

`m= brm(y ~ year + (year |d| group), df2)

m`

>  Family: gaussian 
>   Links: mu = identity; sigma = identity 
> Formula: y ~ year + (year | d | group) 
>    Data: df2 (Number of observations: 42) 
> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
>          total post-warmup samples = 4000
> 
> Group-Level Effects: 
> ~group (Number of levels: 6) 
>                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
> sd(Intercept)          17.53      6.63     8.74    34.65 1.00     1371     1891
> sd(year)                3.05      1.41     1.35     6.73 1.00     1115     1843
> cor(Intercept,year)     0.25      0.39    -0.58     0.84 1.00     1586     2178
> 
> Population-Level Effects: 
>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
> Intercept    20.12      6.95     6.39    34.02 1.00     1441     1835
> year          3.93      1.28     1.27     6.50 1.00     1446     1298
> 
> Family Specific Parameters: 
>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
> sigma     3.99      0.56     3.07     5.25 1.00     2108     2274
> 
> 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).

Group level effects

conditions = tibble(group = unique(df2$group),
                    row.names = unique(df2$group))
conditional_effects(m_default2, effects = "year", re_formula = NULL, conditions = conditions)

image
image700×432 5.48 KB

1 Like
2

Hi,
it means that the varying intercept and the varying effect of year are assumed to be jointly distributed as multivariate normal.

So to simulate the data, you would do:

library(MASS)

# Currently ignoring the LKJ prior on the correlation matrix and "faking it" - for larger correlation matrices
# this will not work
cor_intercept_year <- runif(1, -1, 1)
corM <- matrix(c(1, cor_intercept_year, cor_intercept_year, 1), nrow = 2, ncol = 2)
sd_intercept <- abs(rnorm(1))
sd_year <- abs(rnorm(1))

all_sds <- c(sd_intercept, sd_year)
covM <- diag(all_sds) %*% corM %*% diag(all_sds) 

n_groups <- 13

group_terms <- MASS::mvrnorm(n = n_groups, mu = c(0,0), Sigma = covM)
r_Intercept_group <- group_terms[,1]
r_year_group <- group_terms[,2]

as in “standard” vayring effects/intercepts, the parameters of the multivariate normal (the marginal standard deviations and the correlations) are esimated from data.

Does that make sense?

1 Like