Calculating correlation between intercepts and slope magnitude

Interfaces brms
1

Hi everyone,

I’m still fairly new to brms and Bayesian statistics. My question is about intercept-slope correlations as well as slope-slope correlations in brms. I repeatedly measured activity levels and risk-taking behaviour in animals over 8 experimental days. A simplified version of my dataset is as follows:

ID      MASS       DAY      ACTIVITY      RISK
1       0.7         1          25          41
2       0.9         1         106          22
3       0.6         1          81          37
1       0.7         2          41          68
2       0.9         2          93          36
3       0.6         2          57          11

I am interested in:

  1. Seeing if more active individuals are more “behaviourally plastic” (i.e. have steeper slopes) than less active individuals.
  2. Whether individuals that are more “behaviourally plastic” (i.e. have steeper slopes) in their activity are also more behaviourally plastic in their risk-taking behaviour.

I have set my model out as follows:

ACTIVITY <- bf(ACTIVITY ~ MASS + DAY + (DAY|a|ID) , family = gaussian)
RISK <- bf(RISK ~ MASS + DAY + (DAY|a|ID) , family = gaussian)

Model<- brm(ACTIVITY + RISK + set_rescor(FALSE),
                          data = dat,
                          cores = 4,
                          chains = 4,
                          warmup = 1000,
                          iter = 10000)

In the output for this model I get an intercept slope correlation (i.e. cor(ACTIVITY_Intercept,ACTIVITY_DAY)). However, as I understand it, the correlation provided is the rank correlation between individual differences from the average population intercept and the average population slope (i.e. this takes into account both the magnitude AND direction of the slopes). However, I am interested in the correlation between the intercept and the slope magnitude (i.e. do more active individuals have steeper slopes, regardless of direction, than less active individuals). So essentially what I’m looking for is a relationship between the intercept and the absolute value of the slope. Is there any way to estimate this relationship from the model?

Similarly, in regards to my second question, I get a slope-slope correlation from the model output (i.e. cor(ACTIVITY_DAY, RISK_DAY)). Again, if I understand this correctly, this describes whether the rank order of slopes is maintained between the two behaviours. I’m also wondering if there is anyway to estimate the correlation between the absolute magnitude of individuals’ slopes across the two behaviours, so that I can determine whether individuals that are more “behaviourally plastic” in one behaviour are also more “behaviourally plastic” in another behaviour.

Any help would be greatly appreciated!

Thanks in advance.

1 Like
2

Hi,

I have doubts whether this is a sensible quantity to get from the model - the model assumes that the intercept and slope are gaussian distributed and their raw values correlated. So a correlation between absolute value of slope and the intercept (while possible to compute) is IMHO heavily restricted by this assumption. If this is really the quantity of interest, I would consider writing a model that has it as a “first-class citizen” (although this might be challenging)

If you really want to compute this correlation, I think this is relatively straightforward to achieve - you take the samples for the correlation matrix and the standard deviations, use it to get samples of the full covariance matrix. For each posterior sample of the covariance matrix you then draw a large number of draws from the multivariate normal, take the absolute value of the slope and compute their correlation, i.e. for each posterior sample you compute separate value of this transformed correlation. In pseudocode:

samples_cov_orig <- array [n_samples x 2 x 2] of covariance matrices from the mode
for(s in 1:n_samples) {
   draws <- 1000 samples from multivaraite normal with covariance samples_cov_orig[s]
   draws[,2] <- abs(draws[,2)
   cor_abs[s] <- cor(draws)[1,2]
}

If you get stuck in computing or need other help, feel free to ask for clarifications.

Best of luck!

3

OK, I didn’t think my answer through and there is another way to think about the correlation you discussed: one could take the fitted coefficients for your actual subjects, take the absolute value of the slope and compute for each sample a correlation between the intercept and the abs(slope) across the observed subjects. This way you get the correlation for your observed data, without it necessarily saying anything about the generalizability of this correlation. The plus is that this correlation is less affected by the modelling assumptions - but the model still plays a role as the coefficents are squished together assuming the multivariate normal distribution which could distort the correlation of interest.

4

Hi Martin,

Great. Thanks for the input and help.

I also found this preprint (https://doi.org/10.32942/osf.io/bnugw) which outlines ways to compute these correlations (see section 2.1; Equation 9 & section 2.4; Equation 16). There is also supplementary R scripts that compute the correlations in rstan, however, I had some trouble translating this into brms syntax.

But I will give your suggestion a go and let you know if I get stuck computing the correlation.

Thanks again!

1 Like