If we have a multivariate distribution we can easily estimate the correlations between all our variables by specifying a model like the one here below:

fit <- brm(
bf(mvbind(x, y, z, w) ~ 1) + set_rescor(TRUE),
data = df
)

However, what are the options if we would like to estimate partial correlations? For example, the relationship between x and y by controlling for their relationships with z and w?

Is there a way to specify the brms model such that we can directly estimate the partial correlations? Or, is the only feasible option to use the correlation estimates to calculate the partial correlations by hand using matrix algebra?

In fact, Stan works with canonical partial correlations under the hood: 10.9 Correlation matrices | Stan Reference Manual , so you can extract the unconstrained parameters and those will represent the partial correlatoins. You should be able to access the unconstrained values by using the diagnostic_file option.

Constructing the partial correlations you need yourself is also an option (and might be easier). Note that you can compute posterior samples of any quantity (including partial correlations) by computing the quantity separately for each sample. So the only difference from directly working with partial correlations is that brms puts prior on correlations, not partial correlations.

Perfect, thanks for the clarification! My plan was indeed to compute the posterior samples for the partial correlations by computing them separately for each sample.

I know this has been solved. I just want to add that you can calculate the partial correlations from the inverse of the correlation matrix. See wikipedia. This would have to be done for each sample but might be worthwhile for some applications.