# Residual correlation of multivariate model for binary response variable

**URL:** <https://discourse.mc-stan.org/t/residual-correlation-of-multivariate-model-for-binary-response-variable/36745>\
**Category:** General\
**Tags:** multivariate-normal, brms\
**Created:** [September 27, 2024, 11:00am UTC](https://discourse.mc-stan.org/t/residual-correlation-of-multivariate-model-for-binary-response-variable/36745 "2024-09-27T11:00:24Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Emma0305](https://avatars.discourse-cdn.com/v4/letter/e/b487fb/32.png) [@Emma0305](https://discourse.mc-stan.org/u/Emma0305)\
**Post date:** [September 27, 2024, 11:00am UTC](https://discourse.mc-stan.org/t/residual-correlation-of-multivariate-model-for-binary-response-variable/36745/1 "2024-09-27T11:00:24Z")

</div>

Hi all,

I’m working on a multivariate model that involves a continuous response variable (`nicheScore`) and a binary response variable (`PLBin`). My goal is to estimate the **residual correlation** between these two responses. However, I understand that in `brms`, residual correlations (`set_rescor(TRUE)`) currently can only be directly estimated in models where all responses follow **Gaussian** or **Student-t** distributions. Therefore, I wonder if I can calculate the correlation and its 95% credible interval (CI) manually from the posterior samples.

Here’s the approach suggested by ChatGPT:

```no-highlight
A <- ape::vcv.phylo(targetTree) # Phylogenetic covariance matrix
bf_PL <- bf(PLBin ~ (1|d|gr(phylo, cov = A)), family = bernoulli()) # Binary response
bf_NicheScore <- bf(nicheScore ~ (1|d|gr(phylo, cov = A)), family = gaussian()) # Continuous response

bform2 <- bf_PL + bf_NicheScore + set_rescor(FALSE) # Multivariate model without residual correlation

# Fit the model
fit2 <- brm(bform2, data = allTogether, data2 = list(A = A), chains = 2, warmup = warmup, iter = iter, thin = thin)
fit2 <- add_criterion(fit2, "loo")

# Extract posterior predictions for both responses
fitted_PLBin <- posterior_epred(fit2, resp = "PLBin") # Posterior predictions for PLBin (binary)
fitted_NicheScore <- posterior_epred(fit2, resp = "nicheScore") # Posterior predictions for nicheScore (continuous)

# Calculate the correlation between predictions for each posterior sample
n_samples <- dim(fitted_PLBin)[1] # Number of posterior samples
correlations <- numeric(n_samples)

for (i in 1:n_samples) {
  correlations[i] <- cor(fitted_PLBin[i,], fitted_NicheScore[i,])
}

# Compute the 95% credible interval (CI) for the correlation
ci_95 <- quantile(correlations, probs = c(0.025, 0.975))
print(ci_95)

# Plot the posterior distribution of the correlations
hist(correlations, main = "Posterior Correlation", xlab = "Correlation", freq = FALSE)
abline(v = ci_95, col = "red", lwd = 2) # Add lines for the credible interval

```

Would this be an appropriate way to estimate the residual correlation between a continuous and a binary response? Thanks in advance for your help!