tests/testthat/test-splicing.R

In tramvs: Optimal Subset Selection for Transformation Models

### Test splicing / regression tests
test_that("aliases pass arguments by name", {
  set.seed(24101968)
  library("tramvs")

  N <- 1e2
  P <- 5
  nz <- 3
  beta <- rep(c(3, 0), c(nz, P - nz))
  X <- matrix(rnorm(N * P), nrow = N, ncol = P)
  Y <- 1 + X %*% beta + rnorm(N)

  dat <- data.frame(y = Y, x = X)

  res <- LmVS(y ~ ., data = dat, supp_max = 2)
  expect_equal(nrow(SIC(res)), 2)

  res <- LmVS(y ~ ., data = dat, mandatory = y ~ x.5)
  expect_true("x.5" %in% support(res))

  # Arguments to modFUN
  expect_no_error(BoxCoxVS(y ~ ., data = dat, supp_max = 2, order = 4))

  # Splicing size larger than smallest support
  res <- LmVS(y ~ ., data = dat, k_max = 2)
  expect_equal(support(res), c("x.1", "x.2", "x.3"))
})

test_that("scaling the SIC penalty", {
  set.seed(24101968)
  library("tramvs")

  N <- 1e2
  P <- 5
  nz <- 3
  beta <- rep(c(3, 0), c(nz, P - nz))
  X <- matrix(rnorm(N * P), nrow = N, ncol = P)
  Y <- 1 + X %*% beta + rnorm(N)

  dat <- data.frame(y = Y, x = X)

  res <- LmVS(y ~ ., data = dat)
  res_abess <- LmVS(y ~ ., data = dat, ic_scale = 0.5)
  nll <- sapply(res$all_fits, \(fit) -c(logLik(fit$mod)))
  expect_equal(SIC(res_abess)$SIC, (SIC(res)$SIC + nll) / 2)
})

test_that("splicing from random initialization", {
  set.seed(24101968)
  library("tramvs")

  N <- 2e2
  P <- 12
  nz <- 3
  beta <- rep(c(2, 0), c(nz, P - nz))
  X <- matrix(rnorm(N * P), nrow = N, ncol = P)
  Y <- 1 + X %*% beta + rnorm(N)

  dat <- data.frame(y = Y, x = X)

  for (seed in 1:3) {
    set.seed(seed)
    res <- abess_tram(y ~ ., dat, modFUN = Lm, supp = 3, init = FALSE)
    expect_equal(res$A, c("x.1", "x.2", "x.3"))
    expect_equal(res$I, setdiff(names(coef(res)), res$A))

    # Splicing size smaller than support
    set.seed(seed)
    res <- abess_tram(y ~ ., dat, modFUN = Lm, supp = 3, init = FALSE,
                      k_max = 1)
    expect_equal(res$A, c("x.1", "x.2", "x.3"))
  }
})

test_that("mandatory covariates stay active", {
  set.seed(24101968)
  library("tramvs")

  N <- 1e2
  P <- 6
  nz <- 3
  beta <- rep(c(2, 0), c(nz, P - nz))
  X <- matrix(rnorm(N * P), nrow = N, ncol = P)
  Y <- 1 + X %*% beta + rnorm(N)

  dat <- data.frame(y = Y, x = X)

  for (supp in 1:3) {
    res <- abess_tram(y ~ ., dat, modFUN = Lm, supp = supp, init = FALSE,
                      mandatory = y ~ x.6)
    expect_true("x.6" %in% res$A)
    expect_length(res$A, supp + 1)
  }
})

test_that("ties in initialization and splicing", {
  set.seed(24101968)
  library("tramvs")

  expect_equal(tramvs:::.a0_init(c(0.5, 0.3, 0.3, 0.1), 2), c(1, 2))
  expect_equal(tramvs:::.ak_compute(c(0, 0, -1), 1), 1)
  expect_equal(tramvs:::.ik_compute(c(-1, -1, 0), 1), 1)

  # Duplicated covariate
  N <- 1e2
  X <- matrix(rnorm(N * 3), nrow = N, ncol = 3)
  Y <- 1 + X %*% c(2, 1, 1) + rnorm(N)

  dat <- data.frame(y = Y, x = X[, c(1, 2, 2, 3)])
  res <- tramvs(y ~ ., data = dat, modFUN = Lm)
  expect_equal(sapply(res$all_fits, \(fit) length(fit$A)), 1:4)
})

test_that("mmlt", {
  set.seed(24101968)
  library("tramvs")

  N <- 3e2
  S <- matrix(c(1, 0.6, 0, 0.6, 1, 0, 0, 0, 1), nrow = 3)
  Z <- matrix(rnorm(N * 3), nrow = N) %*% chol(S)
  dat <- data.frame(y1 = Z[, 1], y2 = Z[, 2], y3 = Z[, 3])

  m1 <- BoxCox(y1 ~ 1, data = dat)
  m2 <- BoxCox(y2 ~ 1, data = dat)
  m3 <- BoxCox(y3 ~ 1, data = dat)
  args <- list(m1, m2, m3, formula = ~ 1, data = dat)

  res <- mmltVS(args)
  expect_equal(res$all_fits[[1]]$A, "y2.y1.(Intercept)")
  expect_true("y2.y1.(Intercept)" %in% support(res))
  expect_equal(dim(coef(res)), c(3, 3))

  # Supplying m0
  m0 <- do.call("mmlt", args)
  expect_no_error(mmltVS(args, m0 = m0, supp_max = 1))
})

test_that("plot restores par", {
  set.seed(24101968)
  library("tramvs")

  N <- 1e2
  P <- 5
  nz <- 3
  beta <- rep(c(3, 0), c(nz, P - nz))
  X <- matrix(rnorm(N * P), nrow = N, ncol = P)
  Y <- 1 + X %*% beta + rnorm(N)

  dat <- data.frame(y = Y, x = X)
  res <- tramvs(y ~ ., data = dat, modFUN = Lm)

  opar <- par("mar")
  plot(res)
  plot(res, which = "path")
  expect_equal(par("mar"), opar)
})