UPDATED: My first post contained unnecessary code parts in the stan program. These parts have been removed now.

In case it’s of interest, I have taken a (probably completely wrong) shot at the bivariate case of this problem, and I’d be curious to hear whether what I’m doing makes sense. Specifically, I took a bivariate distributional brms model and modified the stan code to predict the correlations between F1 and F2.

The model is meant to analyze vowel formants F1 and F2 (the two acoustic dimensions along which vowels are distinguished), for which we are aiming to predict whether accent (whether the talker is a native or non-native speaker of the language) affects the distribution of the vowel in F1-F2 space. For simplicity’s sake I’m assuming vowels have bivariate Gaussian distributions in F1-F2 space, and I’ll focus on one vowel.

The bivariate model of F1 and F2 with linear predictors for the mean and standard deviations, but not the correlations between F1 and F2 (Talker is between Accent since each talker either is a native OR a non-native speaker of the language):

```
bf_F1 =
bf(F1 ~ 1 + Accent + (1 | p | Talker)) +
lf(sigma ~ 1 + Accent + (1 | q | Talker)) + gaussian()
bf_F2 =
bf(F2 ~ 1 + Accent + (1 | p | Talker)) +
lf(sigma ~ 1 + Accent + (1 | q | Talker)) + gaussian()
formula = bf_F1 + bf_F2
```

I basically want to extent the model to also model:

`rho ~ 1 + Accent + (1 | Talker)`

Here’s the modified stancode and some simulated data:

Distributional-Lobanov F1F2-cor-single vowel.stan (11.9 KB)

SimulatedF1F2.csv (81.1 KB)

Following Bloome and Schrage’s (ordinary; not multilevel) model, the linear predictor for correlations models logits (converted into correlations via inv_logit(x) * 2 - 1), though I am wondering whether has unwanted downsides (feedback welcome). Another note that might be helpful is that I stuck as close as possible to the original brms-generated stancode, including some redundancies that I assume have their origin in the goal of brms to generate the code in ways that is maximally general.

And here is some example code to hijack brms’ make_standata() function to prep the input for the model and run it (for this particular example):

```
brm_cor = function(formula, data, filename, ...) {
if (file.exists(filename)) {
m = readRDS(filename)
} else {
standata.new = make_standata(formula, data)
# remove residual variance related input
standata.new$nrescor = NULL
# copy over everything from F1 (same as F2) to the input for the correlations
standata.new$N_F1F2 = standata.new$N_F1
standata.new$K_cor_F1F2 = standata.new$K_F1
standata.new$X_cor_F1F2 = standata.new$X_F1
standata.new$N_3 = standata.new$N_1
standata.new$M_3 = standata.new$M_1
standata.new$J_3_F1F2 = standata.new$J_1_F1
standata.new$Z_3_cor_F1F2_1 = standata.new$Z_1_F1_1
standata.new$NC_3 = standata.new$NC_1
m = stan(
file = "Distributional-Lobanov F1F2-cor-i.stan",
data = standata.new,
chains = 4,
warmup = 1000,
iter = 2000,
cores = 4,
...
)
saveRDS(m, filename)
}
return(m)
}
bf_F1 =
bf(F1 ~ 1 + Accent + (1 | p | Talker)) +
lf(sigma ~ 1 + Accent + (1 | q | Talker)) + gaussian()
bf_F2 =
bf(F2 ~ 1 + Accent + (1 | p | Talker)) +
lf(sigma ~ 1 + Accent + (1 | q | Talker)) + gaussian()
m.simulated.cov = brm_cor(
formula = bf_F1 + bf_F2,
data = d,
filename = "../models/Distributional-SimulatedData-cor.rds"
)
```

The model infers the same means and variances as the brms model without correlations (so far so good). Qualitatively, the correlations also seem to be recovered correctly (with shrinkage to 0, it seems). I’m still running simulations to test whether what I did makes sense.