more flexible residual correlation structures in multivariate models

[wpetry]

In multivariate brms models, the correlation structure between the residuals of the response variables is currently (i) estimated using rescor = TRUE or (ii) fully zeroed out using rescor = FALSE. Would it be possible to allow users to specify residual correlations between a subset of response variables rather than all-or-nothing? To pose the question in another way, could the user specify which elements of the residual correlation matrix should be zero (i.e., making that matrix sparse)?

To illustrate my motivation, we can imagine adding a new, uncorrelated response variable (beak) to the blue tit data used in vignette("brms_multivariate").

library(brms)
data("BTdata", package = "MCMCglmm")
set.seed(3240)

# make a new uncorrelated response variable
BTdata$beak <- rnorm(nrow(BTdata), mean = 0, sd = 1)

Suppose the user has prior information that beak will not be correlated with the other two responses and wishes to constrain the model accordingly. The desired correlation matrix would therefore be:

I attempted to fit this model by setting up one formula for the responses with correlated residuals, one formula for the response to be left out of the residual correlation, then combining them together while specifying no residual correlation between the two formulas:

form_tars_back <- bf(mvbind(tarsus, back) ~ sex + hatchdate + (1|p|fosternest) + (1|q|dam)) +
  set_rescor(TRUE)

form_beak <- bf(beak ~ sex + hatchdate + (1|p|fosternest) + (1|q|dam))

form_all <- mvbf(
  form_tars_back,
  form_beak,
  rescor = FALSE
)

fit <- brm(form_all, BTdata, chains = 2, cores = 2)
fit

No residual correlations are estimated in this model, presumably because the mvbf(..., rescor = FALSE) overrides the correlation in form_tars_back.

Another example of when this could be useful arises when responses come from different distributions, and the user wants to specify residual correlation between the the Gaussian responses. Building on the bird data example, imagine we also counted parasites on each bird:

BTdata$parasites <- rpois(nrow(BTdata), lambda = 8)

I understand that brms has not implemented residual correlations except for multivariate Gaussian and t-distributions—I am not asking for that feature.

The bird example is bit contrived to illustrate what I think is the broader value of more flexible specification of residual correlations. Implementing finer control over the residual correlation structure would open up a broader range of models that use latent variables. Ultimately, it would be nice to be able to fit an instrumental variable model to non-Gaussian response data. brms already can fit IV models where the observed responses are Gaussian (see this blog post). It'd be great to be able to treat the Gaussian part as latent, just adding a non-Gaussian observation process (e.g., Poisson) as a higher level in the hierarchical model.

Please forgive me if this issue is already being addressed via #304. I see the connections between what I'm asking & latent variables, but it's not clear to me that the fix would allow for residual correlations between a subset of responses when there is a mixture of response distributions.

Thank you for your work on this fantastic package.

[paul-buerkner]

This feature is indeed related to #304 and will be implemented in brms 3.0

[mvuorre]

The real question on everyone's minds is "When do I get to try brms 3.0?" 😀

[paul-buerkner]

Once I got a specific project funded and other ones in review that I am currently writing on. Whenever this will be :-D

[djmit91]

Will the update also have the option to separate the residual covariances by a grouping factor, equivalent to the option to do it on the random effects with gr(ID, by = GROUP)?

I thought maybe I could do it like this bf(mvbind(y1,y2) ~ 1 + (1|2|obs), sigma ~0) + set_rescor(FALSE) Where obs is just a unique identifier for each row of the dataframe. But don't think that's doing what I thought it would.

Fixing sigma to a very small number (e.g. sigma = 0.01) gives better results, but mixing is really poor.

Sorry if this is the wrong thread, couldn't see this elsewhere, but it's related. Cheers, David

set.seed(42)

vcov <- matrix(c(1  , 0.5, 
                 0.5,   1), , nrow = 2, ncol = 2)
y <- mvrnorm(n = 100, mu = c(0,1), Sigma = vcov)
ds <- list(y1=y[,1], y2=y[,2], obs=1:100 )

library(brms)

model <- bf(mvbind(y1,y2) ~ 1 + (1|2|obs), sigma ~0)
fit <- brm(model, gaussian(), data = ds, chains = 4, warmup = 500,
           iter = 1500, cores = 4)

#gives
#~obs (Number of levels: 100) 
#                              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#sd(y1_Intercept)                   0.18      0.13     0.01     0.49 1.00      669      652
#sd(y2_Intercept)                   0.22      0.15     0.01     0.54 1.00      748     1933
#cor(y1_Intercept,y2_Intercept)     0.03      0.58    -0.96     0.96 1.00      955     1141

model2 <- bf(mvbind(y1,y2) ~ 1)+ set_rescor(T)
fit2 <- brm(model2, data = ds, chains = 4, warmup = 500,
           iter = 1500, cores = 4)

# Family Specific Parameters: 
#         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#sigma_y1     0.95      0.07     0.83     1.09 1.00     3458     3160
#sigma_y2     0.99      0.07     0.86     1.14 1.00     3773     3281

# Residual Correlations: 
#               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
# rescor(y1,y2)     0.56      0.07     0.42     0.68 1.00     3494     2900

[paul-buerkner]

That is a reasonable request. In fact, residual correlations varying by group would be covered by issue #805 I think.

[djmit91]

Thanks Paul, Ah, yes! Bit more complex, but that'd also do what I am looking for. Hope you have time/funds/inclination to implement this sometime.

Awesome package that's quickly becoming a one-stop-shop for most of my analyses.