tests/testthat/test-ConfidenceIntervals.R

In BranchGLM: Efficient Best Subset Selection for GLMs via Branch and Bound Algorithms

# Confidence intervals tests
## Comparing to glm
### Linear
test_that("Confidence intervals tests linear", {
  library(BranchGLM)
  set.seed(8621)
  x <- sapply(rep(0, 10), rnorm, n = 1000, simplify = TRUE)
  x <- cbind(1, x)
  beta <- rnorm(11)
  y <- rnorm(n = 1000, mean = x %*% beta, sd = 1)
  Data <- cbind(y, x[,-1]) |>
    as.data.frame()
  
  ### Getting confidence intervals
  CIs <- confint(BranchGLM(y ~ ., data = Data, family = "gaussian", 
                                       link = "identity"))
  CIsWald <- confint.default(BranchGLM(y ~ ., data = Data, family = "gaussian", 
                                       link = "identity"))
  
  ### Getting glm confints
  LMCIs <- confint(lm(y ~ ., data = Data))
  LMWald <- confint.default(lm(y ~ ., data = Data))
  
  ### Checking results
  expect_equal(unname(CIs$CIs), unname(LMCIs), tolerance = 1e-2)
  expect_equal(CIsWald, LMWald, tolerance = 1e-2)
})

### Gamma
test_that("Confidence intervals tests gamma", {
  library(BranchGLM)
  set.seed(8621)
  x <- sapply(rep(0, 10), rnorm, n = 1000, simplify = TRUE)
  x <- cbind(1, x)
  beta <- rnorm(11)
  y <- rgamma(n = 1000, shape = 1 / 2, scale = exp(x %*% beta) * 2)
  Data <- cbind(y, x[,-1]) |>
    as.data.frame()
  
  ### Getting confidence intervals
  Fish <- BranchGLM(y ~ ., data = Data, family = "gamma", link = "log", 
                            method = "Fisher")
  FishCIs <- confint(Fish)
  FishCIsWald <- confint.default(Fish)
  BFGS <- BranchGLM(y ~ ., data = Data, family = "gamma", link = "log", 
                            method = "BFGS")
  BFGSCIs <- confint(BFGS)
  BFGSCIsWald <- confint.default(BFGS)
  LBFGS <- BranchGLM(y ~ ., data = Data, family = "gamma", link = "log", 
                             method = "LBFGS")
  LBFGSCIs <- confint(LBFGS)
  LBFGSCIsWald <- confint.default(LBFGS)
  
  ### glm
  GLM <- glm(y ~ ., data = Data, family = Gamma(link = "log"))
  GLMCIsWald <- confint.default(GLM)
  
  ### Checking results
  expect_equal(FishCIs, BFGSCIs, tolerance = 1e-2)
  expect_equal(FishCIs, LBFGSCIs, tolerance = 1e-2)
  expect_equal(unname(FishCIs$CIs), unname(GLMCIsWald), tolerance = 1e-2)
  expect_equal(FishCIsWald, BFGSCIsWald, tolerance = 1e-2)
  expect_equal(FishCIsWald, LBFGSCIsWald, tolerance = 1e-2)
  expect_equal(FishCIsWald, GLMCIsWald, tolerance = 1e-2)
})

### Logistic
test_that("Confidence intervals tests logistic", {
  library(BranchGLM)
  set.seed(8621)
  x <- sapply(rep(0, 10), rnorm, n = 1000, simplify = TRUE)
  x <- cbind(1, x)
  beta <- rnorm(11)
  y <- rbinom(n = 1000, size = 1, p = 1 / (1 + exp(-x %*% beta)))
  Data <- cbind(y, x[,-1]) |>
    as.data.frame()
  
  ### Getting confidence intervals
  Fish <- BranchGLM(y ~ ., data = Data, family = "binomial", link = "logit", 
                    method = "Fisher")
  FishCIs <- confint(Fish)
  FishCIsWald <- confint.default(Fish)
  BFGS <- BranchGLM(y ~ ., data = Data, family = "binomial", link = "logit", 
                    method = "BFGS")
  BFGSCIs <- confint(BFGS)
  BFGSCIsWald <- confint.default(BFGS)
  LBFGS <- BranchGLM(y ~ ., data = Data, family = "binomial", link = "logit", 
                     method = "LBFGS")
  LBFGSCIs <- confint(LBFGS)
  LBFGSCIsWald <- confint.default(LBFGS)
  
  ### glm
  GLM <- glm(y ~ ., data = Data, family = binomial)
  GLMCIsWald <- confint.default(GLM)
  
  ### Checking results
  expect_equal(FishCIs, BFGSCIs, tolerance = 1e-2)
  expect_equal(FishCIs, LBFGSCIs, tolerance = 1e-2)
  expect_equal(unname(FishCIs$CIs), unname(GLMCIsWald), tolerance = 1e-2)
  expect_equal(FishCIsWald, BFGSCIsWald, tolerance = 1e-2)
  expect_equal(FishCIsWald, LBFGSCIsWald, tolerance = 1e-2)
  expect_equal(FishCIsWald, GLMCIsWald, tolerance = 1e-2)
})

### Poisson
test_that("Confidence intervals tests methods", {
  library(BranchGLM)
  set.seed(8621)
  x <- sapply(rep(0, 10), rnorm, n = 1000, simplify = TRUE, sd = 0.1)
  x <- cbind(1, x)
  beta <- rnorm(11)
  y <- rpois(n = 1000, lambda = exp(x %*% beta))
  Data <- cbind(y, x[,-1]) |>
    as.data.frame()
  
  ### Getting confidence intervals
  Fish <- BranchGLM(y ~ ., data = Data, family = "poisson", link = "log", 
                    method = "Fisher")
  FishCIs <- confint(Fish)
  FishCIsWald <- confint.default(Fish)
  BFGS <- BranchGLM(y ~ ., data = Data, family = "poisson", link = "log", 
                    method = "BFGS")
  BFGSCIs <- confint(BFGS)
  BFGSCIsWald <- confint.default(BFGS)
  LBFGS <- BranchGLM(y ~ ., data = Data, family = "poisson", link = "log", 
                     method = "LBFGS")
  LBFGSCIs <- confint(LBFGS)
  LBFGSCIsWald <- confint.default(LBFGS)
  
  ### glm
  GLM <- glm(y ~ ., data = Data, family = poisson)
  GLMCIsWald <- confint.default(GLM)
  
  ### Checking results
  expect_equal(FishCIs, BFGSCIs, tolerance = 1e-2)
  expect_equal(FishCIs, LBFGSCIs, tolerance = 1e-2)
  expect_equal(unname(FishCIs$CIs), unname(GLMCIsWald), tolerance = 1e-2)
  expect_equal(FishCIsWald, BFGSCIsWald, tolerance = 1e-2)
  expect_equal(FishCIsWald, LBFGSCIsWald, tolerance = 1e-2)
  expect_equal(FishCIsWald, GLMCIsWald, tolerance = 1e-2)
})

## Testing which and parm
test_that("Testing which and parm", {
  library(BranchGLM)
  set.seed(8621)
  x <- sapply(rep(0, 10), rnorm, n = 1000, simplify = TRUE)
  x <- cbind(1, x)
  beta <- rnorm(11)
  y <- rnorm(n = 1000, mean = x %*% beta, sd = 1)
  Data <- cbind(y, x[,-1]) |>
    as.data.frame()
  
  ### Getting confidence intervals
  expect_error(confint(BranchGLM(y ~ ., data = Data, family = "gaussian", link = "identity"), 
                       parm = 11:1), NA)
  expect_error(CI1 <- confint(BranchGLM(y ~ ., data = Data, family = "gaussian", link = "identity"), 
                              parm = 1:11), NA)
  expect_error(CI2 <- confint(BranchGLM(y ~ ., data = Data, family = "gaussian", link = "identity"), 
                 parm = c("(Intercept)", colnames(Data)[-1])), NA)
  expect_equal(CI1, CI2)
  
  ### Testing which in plot
  expect_error(plot(CI1, which = c("(Intercept)", colnames(Data)[-1])), NA)
  expect_error(plot(CI1, which = 1:11), NA)
  expect_error(plot(CI1, which = "AlL"), NA)
})

## Testing errors in CI functions
test_that("Testing errors in CI functions", {
  library(BranchGLM)
  set.seed(8621)
  x <- sapply(rep(0, 10), rnorm, n = 1000, simplify = TRUE)
  x <- cbind(1, x)
  beta <- rnorm(11)
  y <- rnorm(n = 1000, mean = x %*% beta, sd = 1)
  Data <- cbind(y, x[,-1]) |>
    as.data.frame()
  
  ### Testing confidence intervals
  expect_error(confint(BranchGLM(y ~ ., data = Data, family = "gaussian", link = "identity"), 
                       parm = 24))
  expect_error(confint(BranchGLM(y ~ ., data = Data, family = "gaussian", link = "identity"), 
                       parm = "apple"))
  CI <- confint(BranchGLM(y ~ ., data = Data, family = "gaussian", link = "identity"))
  expect_error(plot(CI, which = 24))
  expect_error(plot(CI, which = "apple"))
  expect_error(plotCI(t(CI$CIs)))
})