tests/testthat/test-estimators.R

In cvCovEst: Cross-Validated Covariance Matrix Estimation

# Linear Shrinkage Estimator ##################################################
test_that("linear shrinkage estimator with no shrinkage is S_n", {
  expect_identical(
    linearShrinkEst(mtcars, alpha = 1),
    coop::covar(mtcars)
  )
})

test_that("linear shrinkage estimator with 100% shrinkage is identity", {
  expect_identical(
    linearShrinkEst(mtcars, alpha = 0) %>% unname(),
    diag(ncol(mtcars))
  )
})

test_that("Linear shrinkage estimator produces shrunken estimates", {

  # aAl absolute elements of the estimate should be smaller than the
  # absolute value of the entries in the sample covariance matrix estimate
  dat <- scale(mtcars, center = TRUE, scale = TRUE)
  abs_est <- abs(linearShrinkEst(dat, alpha = 0.1))
  abs_sample_cov <- abs(cov(dat))
  expect_true(
    sum(round((abs_sample_cov - abs_est), digits = 6) >= 0) == 121
  )
})

test_that("Linear shrinkage estimator centers data internally", {
  expect_equal(
    linearShrinkEst(mtcars, 1),
    linearShrinkEst(mtcars - 1, 1)
  )
})

# Ledoit Wolf Linear Shrinkage Estimator #######################################
test_that("LW LS estimator runs without issue", {
  expect_silent(
    linearShrinkLWEst(mtcars)
  )
})

test_that("LW LS estimator centers data internally", {
  expect_equal(
    linearShrinkLWEst(mtcars),
    linearShrinkLWEst(mtcars - 1)
  )
})

# Simple Thresholding Estimator ###############################################
test_that("simple thresholing estimator without thresholding is S_n", {
  expect_identical(
    thresholdingEst(mtcars, gamma = 0),
    coop::covar(mtcars)
  )
})

test_that("simple thresholing estimator with large threshold is 0 matrix", {
  expect_identical(
    thresholdingEst(mtcars, gamma = 1000000) %>% unname(),
    matrix(data = 0, nrow = ncol(mtcars), ncol = ncol(mtcars))
  )
})

test_that("simple thresholding estimator centers data internally", {
  expect_equal(
    thresholdingEst(mtcars, 0.2),
    thresholdingEst(mtcars - 1, 0.2)
  )
})

# Banding Estimator ###########################################################
test_that("banding estimator with k = 0 is diagonal of S_n", {
  expect_identical(
    bandingEst(mtcars, k = 0L) %>% unname(),
    diag(diag(coop::covar(mtcars)))
  )
})

test_that("banding estimator with k >> 0 is S_n", {
  expect_identical(
    bandingEst(mtcars, k = 1000000L),
    coop::covar(mtcars) %>% unname()
  )
})

test_that("banding estimator centers data internally", {
  expect_equal(
    bandingEst(mtcars, 4),
    bandingEst(mtcars - 1, 4)
  )
})

# Tapering Estimator ##########################################################
test_that("tapering estimator with k = 0 is diagonal of S_n", {
  expect_identical(
    taperingEst(mtcars, k = 0L) %>% unname(),
    diag(diag(coop::covar(mtcars)))
  )
})

test_that("tapering estimator with k = 2*p-2 is S_n", {
  expect_identical(
    taperingEst(mtcars, k = 20L) %>% unname(),
    coop::covar(mtcars) %>% unname()
  )
})

test_that("tapering estimator with k = 2*p-4 is not S_n", {
  expect_false(identical(
    taperingEst(mtcars, k = 18L) %>% unname(),
    coop::covar(mtcars) %>% unname()
  ))
})

test_that("tapering estimator with k >> 0 is S_n", {
  expect_identical(
    taperingEst(mtcars, k = 1000000L) %>% unname(),
    coop::covar(mtcars) %>% unname()
  )
})

test_that("tapering estimator centers data internally", {
  expect_equal(
    taperingEst(mtcars, 4),
    taperingEst(mtcars - 1, 4)
  )
})

# Ledoit Wolf Nonlinear Shrinkage Estimator ####################################
test_that("LW NLS estimator runs without issue", {
  expect_silent(
    nlShrinkLWEst(mtcars)
  )

  # make sure that nlShrinkLWEst can handle case where n = p

  library(MASS)
  set.seed(123)
  Sigma <- matrix(0.5, nrow = 50, ncol = 50) + diag(0.5, nrow = 50)
  dat <- mvrnorm(n = 50, mu = rep(0, 50), Sigma = Sigma)
  expect_false(
    any(is.na(nlShrinkLWEst(dat)))
  )
})

# Dense linear Shrinkage Estimator ####################################
test_that("Dense linear shrinkage estimator runs without issue", {
  expect_silent(
    denseLinearShrinkEst(mtcars)
  )
})

test_that("Dense linear shrinkage estimator produces shrunken estimates", {

  # Mean covariance is -0.05 in sample covarian matrix.
  # All absolute covariance values are larger than 0.057.
  # Estimator should therefore produce smaller estimates in all off diagonal
  # entry (given scaling).
  dat <- scale(mtcars, center = TRUE, scale = TRUE)
  abs_est <- abs(denseLinearShrinkEst(dat))
  abs_sample_cov <- abs(cov(dat))
  expect_equal(
    sum(round((abs_sample_cov - abs_est), digits = 6) >= 0), 121
  )
  expect_equal(diag(abs_sample_cov), diag(abs_est))
  expect_true(abs_est[2, 5] != abs_sample_cov[2, 5])
})

test_that("dense linear shrinkage estimator centers data internally", {
  expect_equal(
    denseLinearShrinkEst(mtcars),
    denseLinearShrinkEst(mtcars - 1)
  )
})

# SCAD Thresholding Estimator ##################################################

test_that("SCAD estimator doesn't generate any errors for no reason", {
  expect_silent(
    scadEst(mtcars, lambda = 0.1)
  )
})

test_that("SCAD estimator generates zero matrix for large lambda", {
  dat <- scale(mtcars, center = TRUE, scale = TRUE)
  expect_equal(
    sum(scadEst(dat, lambda = 10) == 0),
    121
  )
})

test_that("SCAD estimator centers data internally", {
  expect_equal(
    scadEst(mtcars, 0.2),
    scadEst(mtcars - 1, 0.2)
  )
})

# POET Estimator ###############################################################

test_that("Verify POET estimator's results", {
  # check that the off diagonal elemets of the POET estimator equal to the
  # off diagonal elements of the rank-1 approximation of the sample covariance
  # matrix

  # compute the rank-1 approx
  dat <- scale(mtcars, center = TRUE, scale = TRUE)
  sample_cov_mat <- coop::covar(dat)
  eig_decomp <- RSpectra::eigs_sym(sample_cov_mat, 1)
  rank_one_sample_cov <- eig_decomp$values *
    eig_decomp$vectors %*% t(eig_decomp$vectors)

  # compute the POET estimate with a large lambda
  poet_estimate <- poetEst(dat, k = 1, lambda = 10)

  # remove the diagonal
  diag(rank_one_sample_cov) <- 0
  diag(poet_estimate) <- 0
  poet_estimate <- poet_estimate %>% unname()

  # compare
  expect_identical(round(rank_one_sample_cov, 10), round(poet_estimate, 10))
})

test_that("POET estimator centers data internally", {
  expect_equal(
    poetEst(mtcars, 2, 3),
    poetEst(mtcars - 1, 2, 3)
  )
})

# Robust POET Estimator ########################################################

test_that("Verify Robust POET estimator's ranks", {
  # check that the rank of robust POET equal to k when lambda is large
  # compute the robust POET estimate with a large lambda
  dat <- scale(mtcars, center = TRUE, scale = TRUE)
  k <- ceiling(ncol(dat) / 5)
  robust_poet_estimate <- robustPoetEst(dat, k,
    lambda = 10,
    var_est = "sample"
  )

  # compare
  library(Matrix)
  expect_equal(Matrix::rankMatrix(robust_poet_estimate)[1], k)
})


test_that("Verify Robust POET estimator's results", {
  # In this specific example, robust POET estimator is supposed to return a
  # matrix of all 13 or all 19.78292 dependent on the variance estimation method
  Y <- matrix(c(1:12, 2:13, 3:14, 4:15, 5:16), 12, 5)
  est <- matrix(rep(13L, 25), 5, 5)
  est_mad <- matrix(rep((3 * 1.4826)**2, 25), 5, 5)
  robust_poet_estimate <- robustPoetEst(Y, 1,
    lambda = 10,
    var_est = "sample"
  )
  robust_poet_estimate_mad <- robustPoetEst(Y, 1,
    lambda = 10,
    var_est = "mad"
  )
  # compare
  expect_equal(est, robust_poet_estimate)
  expect_equal(est_mad, robust_poet_estimate_mad)
})

test_that("robust POET estimator centers data internally", {
  expect_equal(
    robustPoetEst(mtcars, 2, 3, "sample"),
    robustPoetEst(mtcars - 1, 2, 3, "sample")
  )
})

# Adaptive Lasso Estimator #####################################################

test_that("adaptive Lasso estimator with no penalty is Sn", {
  expect_identical(
    adaptiveLassoEst(mtcars, lambda = 0, n = 0) %>% unname(),
    coop::covar(mtcars) %>% unname()
  )
})

test_that("adaptive Lasso estimator with no penalty and non-zero n is Sn", {
  expect_identical(
    adaptiveLassoEst(mtcars, lambda = 0, n = 1) %>% unname(),
    coop::covar(mtcars) %>% unname()
  )
})

test_that("adaptive Lasso estimator with large threshold is 0 matrix", {
  expect_identical(
    adaptiveLassoEst(mtcars, lambda = 1000000, n = 0) %>% unname(),
    matrix(data = 0, nrow = ncol(mtcars), ncol = ncol(mtcars))
  )
})

test_that("adaptive LASSO estimator centers data internally", {
  expect_equal(
    adaptiveLassoEst(mtcars, 0.2, 0.8),
    adaptiveLassoEst(mtcars - 1, 0.2, 0.8)
  )
})

# Operator Norm Shrinkage Estimator, Spiked Covariance Model ###################

test_that(paste("returns the sample covariance matrix when the number of",
                "spikes is zero"), {
  library(MASS)
  set.seed(123)
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(1, nrow = 10))
  expect_equal(spikedOperatorShrinkEst(dat, p_n_ratio = 0.1), sampleCovEst(dat))
})

test_that(paste("estimated eigenvalues are approximately equal to population",
                "eigenvalues in spike model with one spike"), {
  library(MASS)
  set.seed(62342)
  eig_vals <- c(10, rep(0.5, 9))
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(eig_vals, nrow = 10))
  expect_equal(eigen(spikedOperatorShrinkEst(dat, p_n_ratio = 0.1))$values,
               eig_vals, tolerance = 0.5)
})

test_that(paste("noise terms are exactly equal to corresponding eigenvalues of",
                "estimate when specified"), {
  library(MASS)
  set.seed(1231)
  eig_vals <- c(10, rep(0.5, 9))
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(eig_vals, nrow = 10))
  estimated_eig_vals <- eigen(
    spikedOperatorShrinkEst(dat, p_n_ratio = 0.1, noise = 0.5)
  )$values
  expect_equal(estimated_eig_vals[2:10], eig_vals[2:10])
})

test_that(paste("number of spikes in estimates equals to num_spikes when, even",
                "if specified number is wrong"), {
  library(MASS)
  set.seed(1231)
  eig_vals <- c(50, 45, 40, rep(1, 7))
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(eig_vals, nrow = 10))
  estimated_eig_vals <- eigen(
    spikedOperatorShrinkEst(dat, p_n_ratio = 0.1, num_spikes = 2)
  )$values
  expect_equal(sum(estimated_eig_vals > 2), 2)
})

test_that("operator norm shrinkage estimator centers data internally", {
  expect_equal(
    spikedOperatorShrinkEst(mtcars, p_n_ratio = 0.1, 2),
    spikedOperatorShrinkEst(mtcars - 1, p_n_ratio = 0.1, 2)
  )
})

# Frobenius Norm Shrinkage Estimator, Spiked Covariance Model ##################

test_that(paste("returns the sample covariance matrix when the number of",
                "spikes is zero"), {
  library(MASS)
  set.seed(123)
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(1, nrow = 10))
  expect_equal(spikedFrobeniusShrinkEst(dat, p_n_ratio = 0.1),
               sampleCovEst(dat))
})

test_that(paste("estimated eigenvalues are approximately equal to population",
                "eigenvalues in spike model with one spike"), {
  library(MASS)
  set.seed(62342)
  eig_vals <- c(10, rep(0.5, 9))
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(eig_vals, nrow = 10))
  expect_equal(eigen(spikedFrobeniusShrinkEst(dat, p_n_ratio = 0.1))$values,
               eig_vals, tolerance = 0.5)
})

test_that(paste("noise terms are exactly equal to corresponding eigenvalues of",
                "estimate when specified"), {
  library(MASS)
  set.seed(1231)
  eig_vals <- c(10, rep(0.5, 9))
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(eig_vals, nrow = 10))
  estimated_eig_vals <- eigen(
    spikedFrobeniusShrinkEst(dat, p_n_ratio = 0.1, noise = 0.5)
  )$values
  expect_equal(estimated_eig_vals[2:10], eig_vals[2:10])
})

test_that(paste("number of spikes in estimates equals to num_spikes when, even",
                "if specified number is wrong"), {
  library(MASS)
  set.seed(1231)
  eig_vals <- c(50, 45, 40, rep(1, 7))
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(eig_vals, nrow = 10))
  estimated_eig_vals <- eigen(
    spikedFrobeniusShrinkEst(dat, p_n_ratio = 0.1, num_spikes = 2)
  )$values
  expect_equal(sum(estimated_eig_vals > 2), 2)
})

test_that("Frobenius norm shrinkage estimator centers data internally", {
  expect_equal(
    spikedFrobeniusShrinkEst(mtcars, p_n_ratio = 0.1, 2),
    spikedFrobeniusShrinkEst(mtcars - 1, p_n_ratio = 0.1, 2)
  )
})

# Stein Loss Shrinkage Estimator, Spiked Covariance Model ##################

test_that(paste("returns the sample covariance matrix when the number of",
                "spikes is zero"), {
  library(MASS)
  set.seed(123)
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(1, nrow = 10))
  expect_equal(spikedSteinShrinkEst(dat, p_n_ratio = 0.1),
               sampleCovEst(dat))
})

test_that(paste("estimated eigenvalues are approximately equal to population",
                "eigenvalues in spike model with one spike"), {
  library(MASS)
  set.seed(62342)
  eig_vals <- c(10, rep(0.5, 9))
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(eig_vals, nrow = 10))
  expect_equal(eigen(spikedSteinShrinkEst(dat, p_n_ratio = 0.1))$values,
               eig_vals, tolerance = 0.5)
})

test_that(paste("noise terms are exactly equal to corresponding eigenvalues of",
                "estimate when specified"), {
  library(MASS)
  set.seed(1231)
  eig_vals <- c(10, rep(0.5, 9))
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(eig_vals, nrow = 10))
  estimated_eig_vals <- eigen(
    spikedSteinShrinkEst(dat, p_n_ratio = 0.1, noise = 0.5)
  )$values
  expect_equal(estimated_eig_vals[2:10], eig_vals[2:10])
})

test_that(paste("number of spikes in estimates equals to num_spikes when, even",
                "if specified number is wrong"), {
  library(MASS)
  set.seed(1231)
  eig_vals <- c(50, 45, 40, rep(1, 7))
  dat <- mvrnorm(n = 100, mu = rep(0, 10), Sigma = diag(eig_vals, nrow = 10))
  estimated_eig_vals <- eigen(
    spikedSteinShrinkEst(dat, p_n_ratio = 0.1, num_spikes = 2)
  )$values
  expect_equal(sum(estimated_eig_vals > 2), 2)
})

test_that("Stein loss norm shrinkage estimator centers data internally", {
  expect_equal(
    spikedSteinShrinkEst(mtcars, p_n_ratio = 0.1, 2),
    spikedSteinShrinkEst(mtcars - 1, p_n_ratio = 0.1, 2)
  )
})