Variance Shrinkage

Overview

In this section, we demonstrate how to use variance shrinkage when fitting a matrix linear model with MatrixLM.jl.

Within the matrix linear model framework,

\[Y = XBZ^T + E,\]

where

• $Y_{n \times m}$ is the response matrix,
• $X_{n \times p}$ is the row-predictor matrix,
• $Z_{m \times q}$ is the response-attribute matrix,
• $B_{p \times q}$ is the coefficient matrix, and
• $E_{n \times m}$ is the error matrix.

The rows of $E$ are assumed to have covariance matrix $\Sigma$ across the $m$ responses. Estimating $\Sigma$ directly from the residuals can be unstable when the number of responses is large relative to the sample size. Variance shrinkage stabilizes this estimate by shrinking noisy sample quantities toward structured targets using the analytic shrinkage procedures proposed by Schäfer and Strimmer (2005)^[1] and Opgen-Rhein and Strimmer (2007)^[2].

MatrixLM.jl provides a simple Boolean interface for this choice:

mlm(data)         # no variance shrinkage
mlm(data, false)  # no variance shrinkage
mlm(data, true)   # use the a shrinkage estimator

When varShrinkage=true, MatrixLM.jl applies a shrinkage estimator to the covariance matrix. The covariance matrix is factored into a product of the variances and the correlation matrix as described by Barnard et al. (2000)^[3]. The variances and the correlations are each shrunk separately. The variances are shrunk towards a common variance target after using the variance stabilizing log transformation. The off-diagonal entries of the correlation matrix are shrunk towards a common correlation after using the variance stabilizing arctanh transformation.

Data Generation

We begin by simulating data from a matrix linear model with correlated, heteroskedastic errors. This setting is useful for illustrating variance shrinkage because the response variables have both unequal variances and nonzero correlations.

using MatrixLM, LinearAlgebra, Random, Statistics
using StableRNGs

rng = StableRNG(2026)

# Matrix dimensions
n = 100
m = 250
p = 5
q = 4

# Row and column predictors
X = randn(rng, n, p)
Z = randn(rng, m, q)

# True coefficient matrix
B = [
     1.5   0.0   0.5  -1.0;
     0.0   1.0   0.0   0.5;
    -0.5   0.0   1.5   0.0;
     0.0  -1.0   0.5   1.0;
     1.0   0.5   0.0   0.0
]

To illustrate variance shrinkage, we generate an error matrix with correlated responses and unequal response variances. Correlation is introduced through a common random effect shared by all responses, while each response is assigned its own variance.

# Unequal response variances
variances = range(0.5, 2.0, length=m)
standard_deviations = sqrt.(variances)

# Common correlation between responses
rho = 0.6

# Generate correlated, heteroskedastic errors
common_noise = randn(rng, n, 1)
independent_noise = randn(rng, n, m)

E = (
    sqrt(rho) .* common_noise .+
    sqrt(1 - rho) .* independent_noise
) .* reshape(standard_deviations, 1, m)

# True covariance matrix
R = fill(rho, m, m)
R[diagind(R)] .= 1.0
Dhalf = Diagonal(standard_deviations)
Sigma_true = Dhalf * R * Dhalf

# Generate the response matrix
Y = X * B * Z' + E

Now construct a RawData object containing Y, X, and Z.

dat = RawData(
    Response(Y),
    Predictors(X, Z)
)

We can use the function show() to respectively display a readable summary of the matrices and dimensions stored in a RawData object.

show(dat)

RawData
  Response matrix Y: 100 × 250
  Design matrix X: 100 × 5
  Design matrix Z: 250 × 4
  X includes intercept: false
  Z includes intercept: false
  Preview of Y first row (first 3 columns): [0.1345, -1.0693, -0.1144, …]
  Preview of X first row (first 3 columns): [0.5803, 0.1856, 0.3414, …]
  Preview of Z first row (first 3 columns): [-0.9384, -0.4889, -0.3803, …]

Model Estimation Without Variance Shrinkage

Set the positional Boolean argument to false to estimate the residual covariance matrix without shrinkage.

est_no_shrinkage = mlm(
    dat,
    false;
    addXIntercept=false,
    addZIntercept=false
)

We can use the function show() to respectively display a readable summary of the matrices and dimensions stored in a Mlm object.

show(est_no_shrinkage)

Matrix linear model fit
  Observations: 100
  Responses: 250
  Row predictors: 5
  Column predictors: 4
  Coefficient matrix B: 5 × 4
  Residual covariance matrix sigma(Σ): 250 × 250
  Weighted fit: false
  Preview of B first row (first 3 columns): [1.5029, 0.0057, 0.4933, …]
  Preview of sigma(Σ) first row (first 3 columns): [0.4968, 0.3239, 0.3731, …]
  Covariance shrinkage: none

The fitted object contains the coefficient estimates in B, the estimated error covariance matrix in sigma, and the coefficient variance estimates in varB.

size(est_no_shrinkage.B), size(est_no_shrinkage.sigma),
size(est_no_shrinkage.varB)

((5, 4), (250, 250), (5, 4))

Model Estimation With Variance Shrinkage

Set the positional Boolean argument to true to use the recommended variance shrinkage estimator.

est_shrinkage = mlm(
    dat,
    true;
    addXIntercept=false,
    addZIntercept=false
)

show(est_shrinkage)

Matrix linear model fit
  Observations: 100
  Responses: 250
  Row predictors: 5
  Column predictors: 4
  Coefficient matrix B: 5 × 4
  Residual covariance matrix sigma(Σ): 250 × 250
  Weighted fit: false
  Preview of B first row (first 3 columns): [1.5029, 0.0057, 0.4933, …]
  Preview of sigma(Σ) first row (first 3 columns): [0.5658, 0.332, 0.3739, …]
  Covariance shrinkage target: R
  Shrinkage coefficient: (correlation = 1.0, variance = 0.1292)

The estimator separately shrinks

1. arctanh-transformed sample correlations toward their common mean, and
2. log-transformed sample variances toward their common mean.

The fitted object records the estimated shrinkage intensities.

est_shrinkage.lambda

(correlation = 1.0, variance = 0.12917425638429167)

The returned lambda object contains separate shrinkage intensities for the correlation and variance components.

Comparing the Two Fits

Variance shrinkage changes the estimated error covariance matrix and, consequently, the estimated variances and standard errors of the coefficients. It does not change the least-squares coefficient estimates.

maximum(abs.(est_no_shrinkage.B - est_shrinkage.B))

0.0

The covariance estimates generally differ:

covariance_difference =
    norm(est_no_shrinkage.sigma - est_shrinkage.sigma)

covariance_difference

18.828114017531252

The coefficient variance estimates can also differ because they depend on the estimated error covariance matrix:

coefficient_variance_difference =
    norm(est_no_shrinkage.varB - est_shrinkage.varB)

coefficient_variance_difference

1.0564812118034594e-5

For this simulated example, we can compare both covariance estimators with the known data-generating covariance matrix.

error_no_shrinkage =
    norm(est_no_shrinkage.sigma - Sigma_true)

error_shrinkage =
    norm(est_shrinkage.sigma - Sigma_true)

(
    no_shrinkage = error_no_shrinkage,
    R_shrinkage = error_shrinkage
)

(no_shrinkage = 28.345336570959645, R_shrinkage = 17.11807784892948)

Standard Errors and T-statistics

The coefficient standard errors are obtained from varB.

se_no_shrinkage = sqrt.(est_no_shrinkage.varB)
se_shrinkage = sqrt.(est_shrinkage.varB)

5×4 Matrix{Float64}:
 0.00644423  0.00564131  0.00585292  0.00584954
 0.0069296   0.00606621  0.00629376  0.00629012
 0.00645633  0.00565191  0.00586391  0.00586052
 0.00590891  0.00517269  0.00536672  0.00536362
 0.0062813   0.00549868  0.00570494  0.00570164

The corresponding t-statistics can be obtained with t_stat.

tstats_no_shrinkage = t_stat(est_no_shrinkage)
tstats_shrinkage = t_stat(est_shrinkage)

Because variance shrinkage affects varB, it may also affect t-statistics, confidence intervals, and hypothesis-testing results.

Summary

Variance shrinkage can be enabled in MatrixLM.jl by passing a Boolean as the second positional argument to mlm:

est = mlm(data, true)

The two principal calls are:

mlm(data, false)  # no variance shrinkage
mlm(data, true)   # variance shrinkage

The estimator stabilizes the residual covariance estimate by separately shrinking transformed correlations and transformed variances. The resulting covariance estimate is stored in est.sigma, the coefficient variance estimates are stored in est.varB, and the estimated shrinkage intensities are stored in est.lambda.

References