Residual correlation of multivariate model for binary response variable

General
,
1

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:

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!