Main function to fit the Bayesian Multivariate Linear Model

This function fits a Bayesian multivariate linear regression model using conjugate priors and returns a detailed output, including summaries and simulations of fixed effects and hyperparameters, as well as analytic posterior distributions.

Usage

mlvr(formula, data, priors = list(), n_sims = 1000)

Arguments

formula: Formula to specify the model.
data: A data frame containing the covariates and response variables.
priors: A list specifying the prior parameters. Default is `list()`.
n_sims: Number of posterior simulations. Default is 1000.

Value

A list containing summaries, simulations of fixed effects and hyperparameters, and the analytic posterior distributions.

Details

## Model and Priors: The model assumes a multivariate normal distribution for the response variables: $$Y = XB + E$$ where $Y$ is the matrix of response variables, $X$ is the matrix of covariates, $B$ is the matrix of coefficients, and $E$ is the error matrix.

The following prior distributions are used: - $B$ | $\Sigma$ follows a matrix-variate normal distribution: $$B | \Sigma \sim \mathcal{N}(B_0, \Sigma \otimes A^{-1})$$ - $\Sigma$ follows an inverse-Wishart distribution: $$\Sigma \sim \text{IW}(\nu_0, V_0)$$

## Posterior Distributions: - The posterior of $B$ given $\Sigma$ is: $$B | \Sigma, Y, X \sim \mathcal{N}(B_n, \Sigma \otimes (X'X + A)^{-1})$$ - The marginal posterior of $\Sigma$ given $Y, X$ is: $$\Sigma | Y, X \sim \text{IW}(\nu_0 + n, V_0 + S)$$

The analytic forms of the posterior parameters are computed as follows: - $B_n = (X'X + A)^{-1}(X'Y + A B_0)$ - $V_n = V_0 + S$, where $S = (Y - X B_n)'(Y - X B_n) + (B_n - B_0)' A (B_n - B_0)$

## Output: The function returns a list containing: - `summary.fixed`: A summary of the marginal posterior distributions of the fixed effects. - `marginals.fixed`: A list with the simulations of the fixed effects. - `summary.hyperpar`: A summary of the posterior distributions of the hyperparameters (elements of $\Sigma$). - `marginals.hyperpar`: A list with the simulations of the hyperparameters. - `analytic`: A list containing the following analytic posterior parameters: - `B_n`: The posterior mean of the coefficient matrix $B$, calculated as $(X'X + A)^{-1}(X'Y + A B_0)$. - `V_n`: The posterior scale matrix for the inverse-Wishart distribution of $\Sigma$, given by $V_0 + S$. - `df`: The degrees of freedom for the marginal posterior distribution of $\Sigma$, equal to $\nu_0 + n - m + 1$. - `cov_matrix`: The covariance matrix for the marginal posterior distribution of $B | X, Y$, computed as $V_n / \text{df} \otimes (X'X + A)^{-1}$. - `A`: The precision matrix used for $B$. - `formula`: The formula specifying the model. - `X`: The design matrix of covariates. - `Y`: The response matrix. - `call`: The original function call.

## Details: The analytic posterior distributions are computed using the following theoretical results: - The posterior distribution of the coefficient matrix $B$ is derived using the matrix-variate normal distribution: $$B | \Sigma, Y, X \sim \mathcal{N}(B_n, \Sigma \otimes (X'X + A)^{-1})$$ where $B_n$ is the posterior mean and $\Sigma \otimes (X'X + A)^{-1}$ is the Kronecker product of the covariance matrix $\Sigma$ and the precision matrix of the covariates.

- The posterior distribution of the covariance matrix $\Sigma$ follows an inverse-Wishart distribution: $$\Sigma | Y, X \sim \text{IW}(\nu_0 + n, V_0 + S)$$ where $S$ is the sum of the squared residuals and the deviation of $B_n$ from its prior mean $B_0$.

Author

Joaquín Martínez-Minaya jmarmin@eio.upv.es

Examples

if (FALSE) { # \dontrun{
data <- data.frame(X1 = rnorm(100), X2 = rnorm(100), Y1 = rnorm(100), Y2 = rnorm(100))
formula <- as.matrix(data[, c("Y1", "Y2")]) ~ X1 + X2
fit <- mlvr(formula, data)
summary(fit)
plot(fit)
} # }