Update error term covariance matrix of multiple linear regression

This function updates the error term covariance matrix of a multiple linear regression.

Usage

update_Sigma(kappa, E, N, S)

Arguments

kappa: The degrees of freedom (a natural number greater than J-1) of the Inverse Wishart prior for Sigma. Per default, kappa = J + 1.
E: The scale matrix of dimension J-1 x J-1 of the Inverse Wishart prior for Sigma. Per default, E = diag(J - 1).
N: The draw size.
S: A matrix, the sum over the outer products of the residuals $(\epsilon_n)_{n=1,\dots,N}$.

Value

A matrix, a draw from the Inverse Wishart posterior distribution of the error term covariance matrix in a multiple linear regression.

Details

This function draws from the posterior distribution of the covariance matrix $\Sigma$ in the linear utility equation $$U_n = X_n\beta + \epsilon_n,$$ where $U_n$ is the (latent, but here assumed to be known) utility vector of decider $n = 1,\dots,N$, $X_n$ is the design matrix build from the choice characteristics faced by $n$, $\beta$ is the coefficient vector, and $\epsilon_n$ is the error term assumed to be normally distributed with mean $0$ and unknown covariance matrix $\Sigma$. A priori we assume the (conjugate) Inverse Wishart distribution $$\Sigma \sim W(\kappa,E)$$ with $\kappa$ degrees of freedom and scale matrix $E$. The posterior for $\Sigma$ is the Inverted Wishart distribution with $\kappa + N$ degrees of freedom and scale matrix $E^{-1}+S$, where $S = \sum_{n=1}^{N} \epsilon_n \epsilon_n'$ is the sum over the outer products of the residuals $(\epsilon_n = U_n - X_n\beta)_n$.

Examples

### true error term covariance matrix
(Sigma_true <- matrix(c(1,0.5,0.2,0.5,1,0.2,0.2,0.2,2), ncol=3))
#>      [,1] [,2] [,3]
#> [1,]  1.0  0.5  0.2
#> [2,]  0.5  1.0  0.2
#> [3,]  0.2  0.2  2.0
### coefficient vector
beta <- matrix(c(-1,1), ncol=1)
### draw data
N <- 100
X <- replicate(N, matrix(rnorm(6), ncol=2), simplify = FALSE)
eps <- replicate(N, rmvnorm(mu = c(0,0,0), Sigma = Sigma_true), simplify = FALSE)
U <- mapply(function(X, eps) X %*% beta + eps, X, eps, SIMPLIFY = FALSE)
### prior parameters for covariance matrix
kappa <- 4
E <- diag(3)
### draw from posterior of coefficient vector
outer_prod <- function(X, U) (U - X %*% beta) %*% t(U - X %*% beta)
S <- Reduce(`+`, mapply(outer_prod, X, U, SIMPLIFY = FALSE))
Sigma_draws <- replicate(100, update_Sigma(kappa, E, N, S))
apply(Sigma_draws, 1:2, mean)
#>           [,1]        [,2]        [,3]
#> [1,] 1.2410191  0.61306570  0.16269546
#> [2,] 0.6130657  1.15402615 -0.04130621
#> [3,] 0.1626955 -0.04130621  1.50028635
apply(Sigma_draws, 1:2, stats::sd)
#>           [,1]      [,2]      [,3]
#> [1,] 0.2093024 0.1533339 0.1488460
#> [2,] 0.1533339 0.1634999 0.1488295
#> [3,] 0.1488460 0.1488295 0.2210208