tests/testthat/test-gauge.R

In robust2sls: Outlier Robust Two-Stage Least Squares Inference and Testing

test_that("outliers() works correctly", {

  data <- datasets::mtcars
  formula <- mpg ~ cyl + disp | cyl + wt
  data[1, "mpg"] <- NA
  data[2, "cyl"] <- NA
  data[3, "disp"] <- NA
  data[4, "wt"] <- NA

  model1 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "robustified", iterations = 3)

  expect_equal(outliers(model1, 0), 1)
  expect_equal(outliers(model1, 1), 1)
  expect_equal(outliers(model1, 2), 1)
  expect_equal(outliers(model1, 3), 1)

  model2 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "saturated", iterations = 3, shuffle = TRUE,
              shuffle_seed = 42, split = 0.5)

  expect_equal(outliers(model2, 0), 16)
  expect_equal(outliers(model2, 1), 12)
  expect_equal(outliers(model2, 2), 12)
  expect_equal(outliers(model2, 3), 12)

  model3 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "robustified", iterations = "convergence",
              convergence_criterion = 1)

  expect_equal(outliers(model3, 0), 1)
  expect_equal(outliers(model3, 1), 1)
  expect_error(outliers(model3, 2), "fewer iterations than argument")

  model4 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "robustified", iterations = "convergence",
              convergence_criterion = 0.1)

  expect_equal(outliers(model4, 0), 1)
  expect_equal(outliers(model4, 1), 1)
  expect_equal(outliers(model4, 2), 1)
  expect_error(outliers(model4, 3), "fewer iterations than argument")

  model5 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "saturated", iterations = "convergence",
              convergence_criterion = 0.1, shuffle = TRUE,
              shuffle_seed = 42, split = 0.5)

  expect_equal(outliers(model5, 0), 16)
  expect_equal(outliers(model5, 1), 12)
  expect_equal(outliers(model5, 2), 12)
  expect_equal(outliers(model5, 3), 12)
  expect_error(outliers(model5, 4), "fewer iterations than argument")

  model6 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "saturated", iterations = "convergence",
              convergence_criterion = 0.1, shuffle = TRUE,
              shuffle_seed = 24, split = 0.5)

  expect_equal(outliers(model6, 0), 2)
  expect_equal(outliers(model6, 1), 1)
  expect_equal(outliers(model6, 2), 1)
  expect_equal(outliers(model6, 3), 1)
  expect_error(outliers(model6, 4), "fewer iterations than argument")

  model7 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.01,
              initial_est = "saturated", iterations = "convergence",
              convergence_criterion = 0.1, shuffle = TRUE,
              shuffle_seed = 24, split = 0.5)

  expect_equal(outliers(model7, 0), 0)
  expect_equal(outliers(model7, 1), 0)
  expect_equal(outliers(model7, 2), 0)

  # test input error messages

  expect_error(outliers(model1, 1.3), "has to be an integer")
  expect_error(outliers(model2, 1.3), "has to be an integer")
  expect_error(outliers(model1, -1), "has to be >= 0")
  expect_error(outliers(model2, -1), "has to be >= 0")
  expect_error(outliers(model1, "m0"), "has to be numeric")
  expect_error(outliers(model2, "m1"), "has to be numeric")

  class(model1) <- NULL
  class(model2) <- NULL

  expect_error(outliers(model1, 1), "does not have the correct class")
  expect_error(outliers(model2, 1), "does not have the correct class")

})


test_that("outliers_prop() works correctly", {

  data <- datasets::mtcars
  formula <- mpg ~ cyl + disp | cyl + wt
  data[1, "mpg"] <- NA
  data[2, "cyl"] <- NA
  data[3, "disp"] <- NA
  data[4, "wt"] <- NA

  model1 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "robustified", iterations = 3)

  expect_equal(outliers_prop(model1, 0), 1/28)
  expect_equal(outliers_prop(model1, 1), 1/28)
  expect_equal(outliers_prop(model1, 2), 1/28)
  expect_equal(outliers_prop(model1, 3), 1/28)

  model2 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "saturated", iterations = 3, shuffle = TRUE,
              shuffle_seed = 42, split = 0.5)

  expect_equal(outliers_prop(model2, 0), 16/28)
  expect_equal(outliers_prop(model2, 1), 12/28)
  expect_equal(outliers_prop(model2, 2), 12/28)
  expect_equal(outliers_prop(model2, 3), 12/28)

  model3 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "robustified", iterations = "convergence",
              convergence_criterion = 1)

  expect_equal(outliers_prop(model3, 0), 1/28)
  expect_equal(outliers_prop(model3, 1), 1/28)
  expect_error(outliers_prop(model3, 2), "fewer iterations than argument")

  model4 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "robustified", iterations = "convergence",
              convergence_criterion = 0.1)

  expect_equal(outliers_prop(model4, 0), 1/28)
  expect_equal(outliers_prop(model4, 1), 1/28)
  expect_equal(outliers_prop(model4, 2), 1/28)
  expect_error(outliers_prop(model4, 3), "fewer iterations than argument")

  model5 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "saturated", iterations = "convergence",
              convergence_criterion = 0.1, shuffle = TRUE,
              shuffle_seed = 42, split = 0.5)

  expect_equal(outliers_prop(model5, 0), 16/28)
  expect_equal(outliers_prop(model5, 1), 12/28)
  expect_equal(outliers_prop(model5, 2), 12/28)
  expect_equal(outliers_prop(model5, 3), 12/28)
  expect_error(outliers_prop(model5, 4), "fewer iterations than argument")

  model6 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "saturated", iterations = "convergence",
              convergence_criterion = 0.1, shuffle = TRUE,
              shuffle_seed = 24, split = 0.5)

  expect_equal(outliers_prop(model6, 0), 2/28)
  expect_equal(outliers_prop(model6, 1), 1/28)
  expect_equal(outliers_prop(model6, 2), 1/28)
  expect_equal(outliers_prop(model6, 3), 1/28)
  expect_error(outliers_prop(model6, 4), "fewer iterations than argument")

  model7 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.01,
              initial_est = "saturated", iterations = "convergence",
              convergence_criterion = 0.1, shuffle = TRUE,
              shuffle_seed = 24, split = 0.5)

  expect_equal(outliers_prop(model7, 0), 0)
  expect_equal(outliers_prop(model7, 1), 0)
  expect_equal(outliers_prop(model7, 2), 0)

  # test input error messages

  expect_error(outliers_prop(model1, 1.3), "has to be an integer")
  expect_error(outliers_prop(model2, 1.3), "has to be an integer")
  expect_error(outliers_prop(model1, -1), "has to be >= 0")
  expect_error(outliers_prop(model2, -1), "has to be >= 0")
  expect_error(outliers_prop(model1, "m0"), "has to be numeric")
  expect_error(outliers_prop(model2, "m1"), "has to be numeric")

  class(model1) <- NULL
  class(model2) <- NULL

  expect_error(outliers_prop(model1, 1), "does not have the correct class")
  expect_error(outliers_prop(model2, 1), "does not have the correct class")

})


test_that("outlier() works correctly", {

  data <- datasets::mtcars
  formula <- mpg ~ cyl + disp | cyl + wt
  data[1, "mpg"] <- NA
  data[2, "cyl"] <- NA
  data[3, "disp"] <- NA
  data[4, "wt"] <- NA

  model1 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "robustified", iterations = 3)

  model2 <- outlier_detection(data = data, formula = formula,
              ref_dist = "normal", sign_level = 0.05,
              initial_est = "saturated", iterations = 3, shuffle = TRUE,
              shuffle_seed = 42, split = 0.5)

  m1o1 <- matrix(data = c(-1, -1, -1, -1), nrow = 1,
                dimnames = list("Mazda RX4", c("m0", "m1", "m2", "m3")))
  m1o5 <- matrix(data = c(1, 1, 1, 1), nrow = 1,
                dimnames = list("Hornet Sportabout", c("m0", "m1", "m2", "m3")))
  m1o25 <- matrix(data = c(0, 0, 0, 0), nrow = 1,
                dimnames = list("Pontiac Firebird", c("m0", "m1", "m2", "m3")))

  m2o1 <- matrix(data = c(-1, -1, -1, -1), nrow = 1,
                dimnames = list("Mazda RX4", c("m0", "m1", "m2", "m3")))
  m2o5 <- matrix(data = c(0, 0, 0, 0), nrow = 1,
                dimnames = list("Hornet Sportabout", c("m0", "m1", "m2", "m3")))
  m2o6 <- matrix(data = c(1, 1, 1, 1), nrow = 1,
                dimnames = list("Valiant", c("m0", "m1", "m2", "m3")))
  m2o18 <- matrix(data = c(0, 1, 1, 1), nrow = 1,
                dimnames = list("Fiat 128", c("m0", "m1", "m2", "m3")))
  m2o25 <- matrix(data = c(0, 0, 0, 0), nrow = 1,
                dimnames = list("Pontiac Firebird", c("m0", "m1", "m2", "m3")))

  expect_equal(outlier(model1, 1), m1o1)
  expect_equal(outlier(model1, 5), m1o5)
  expect_equal(outlier(model1, 25), m1o25)
  expect_equal(outlier(model2, 1), m2o1)
  expect_equal(outlier(model2, 5), m2o5)
  expect_equal(outlier(model2, 6), m2o6)
  expect_equal(outlier(model2, 18), m2o18)
  expect_equal(outlier(model2, 25), m2o25)

})


test_that("gauge_avar() gives correct errors", {

  p <- generate_param(3, 2, 3, sigma = 2, intercept = TRUE, seed = 42)

  expect_error(gauge_avar(1, 0.05, "robustified", 0, p, 0.5),
               "'ref_dist' must be a character vector of length 1")
  expect_error(gauge_avar(c("n", "n"), 0.05, "robustified", 0, p, 0.5),
               "'ref_dist' must be a character vector of length 1")
  expect_error(gauge_avar("nonexist", 0.05, "robustified", 1, p, 0.4),
               "'ref_dist' must be one of the available reference")
  expect_error(gauge_avar("normal", "a", "saturated", 0, p, 0.4),
               "'sign_level' must be a numeric vector of length 1")
  expect_error(gauge_avar("normal", c(0.01, 0.05), "saturated", 0, p, 0.4),
               "'sign_level' must be a numeric vector of length 1")
  expect_error(gauge_avar("normal", -0.3, "robustified", 1, p, 0.5),
               "'sign_level' must lie strictly between 0 and 1")
  expect_error(gauge_avar("normal", 1.2, "robustified", 1, p, 0.5),
               "'sign_level' must lie strictly between 0 and 1")
  expect_error(gauge_avar("normal", 0.05, 2, 1, p, 0.5),
               "'initial_est' must be a character vector of length 1")
  expect_error(gauge_avar("normal", 0.05, c("a", "b"), 1, p, 0.5),
               "'initial_est' must be a character vector of length 1")
  expect_error(gauge_avar("normal", 0.05, "nonexist", 1, p, 0.5),
               "'initial_est' must be one of the available initial estimators")
  expect_error(gauge_avar("normal", 0.05, "robustified", "1", p, 0.5),
               "'iteration' must either be numeric or 'convergence'")
  expect_error(gauge_avar("normal", 0.05, "robustified", c(0,1), p, 0.5),
               "'iteration' must be a numeric vector of length 1")
  expect_error(gauge_avar("normal", 0.05, "saturated", -2, p, 0.5),
               "'iteration' must be weakly larger than 0")
  expect_error(gauge_avar("normal", 0.05, "saturated", 1.5, p, 0.5),
               "'iteration' must be an integer")
  expect_error(gauge_avar("normal", 0.01, "saturated", 0, p, "0.5"),
               "'split' must be a numeric vector of length 1")
  expect_error(gauge_avar("normal", 0.01, "saturated", 0, p, c(0.2, 0.5)),
               "'split' must be a numeric vector of length 1")
  expect_error(gauge_avar("normal", 0.01, "saturated", 0, p, -1.5),
               "'split' must lie strictly between 0 and 1")
  expect_error(gauge_avar("normal", 0.01, "saturated", 0, p, 2.4),
               "'split' must lie strictly between 0 and 1")
  expect_error(gauge_avar("normal", 0.01, "saturated", 1, p, NULL),
               "'split' cannot be NULL unless initial estimator is 'robustified'")
  expect_silent(gauge_avar("normal", 0.01, "robustified", 1, p, NULL))


  # expect error if specify saturated, split != 0.5, m >= 1
  expect_error(gauge_avar("normal", 0.01, "saturated", iteration = 1,
                          split = 0.3, parameters = p),
               "No theory for m >= 1 if saturated initial estimator is not split in half")
  expect_error(gauge_avar("normal", 0.05, "saturated", iteration = 10,
                          split = 0.3, parameters = p),
               "No theory for m >= 1 if saturated initial estimator is not split in half")

})


test_that("gauge_covar() gives correct errors", {

  p <- generate_param(3, 2, 3, sigma = 2, intercept = TRUE, seed = 42)

  expect_error(gauge_covar(1, 0.01, 0.05, "robustified", 0, p, 0.5),
               "'ref_dist' must be a character vector of length 1")
  expect_error(gauge_covar(c("n", "n"), 0.01, 0.05, "robustified", 0, p, 0.5),
               "'ref_dist' must be a character vector of length 1")
  expect_error(gauge_covar("nonexist", 0.01, 0.05, "robustified", 1, p, 0.4),
               "'ref_dist' must be one of the available reference")
  expect_error(gauge_covar("normal", "a", 0.01, "saturated", 0, p, 0.4),
               "'sign_level1' must be a numeric vector of length 1")
  expect_error(gauge_covar("normal", c(0.01, 0.05), 0.01, "saturated", 0, p, 0.4),
               "'sign_level1' must be a numeric vector of length 1")
  expect_error(gauge_covar("normal", -0.3, 0.01, "robustified", 1, p, 0.5),
               "'sign_level1' must lie strictly between 0 and 1")
  expect_error(gauge_covar("normal", 1.2, 0.01, "robustified", 1, p, 0.5),
               "'sign_level1' must lie strictly between 0 and 1")
  expect_error(gauge_covar("normal", 0.01, "a", "saturated", 0, p, 0.4),
               "'sign_level2' must be a numeric vector of length 1")
  expect_error(gauge_covar("normal", 0.01, c(0.01, 0.05), "saturated", 0, p, 0.4),
               "'sign_level2' must be a numeric vector of length 1")
  expect_error(gauge_covar("normal", 0.01, -0.3, "robustified", 1, p, 0.5),
               "'sign_level2' must lie strictly between 0 and 1")
  expect_error(gauge_covar("normal", 0.01, 1.2, "robustified", 1, p, 0.5),
               "'sign_level2' must lie strictly between 0 and 1")
  expect_error(gauge_covar("normal", 0.05, 0.01, 2, 1, p, 0.5),
               "'initial_est' must be a character vector of length 1")
  expect_error(gauge_covar("normal", 0.05, 0.01, c("a", "b"), 1, p, 0.5),
               "'initial_est' must be a character vector of length 1")
  expect_error(gauge_covar("normal", 0.05, 0.01, "nonexist", 1, p, 0.5),
               "'initial_est' must be one of the available initial estimators")
  expect_error(gauge_covar("normal", 0.05, 0.01, "robustified", "1", p, 0.5),
               "'iteration' must either be numeric or 'convergence'")
  expect_error(gauge_covar("normal", 0.05, 0.01, "robustified", c(0,1), p, 0.5),
               "'iteration' must be a numeric vector of length 1")
  expect_error(gauge_covar("normal", 0.05, 0.01, "saturated", -2, p, 0.5),
               "'iteration' must be weakly larger than 0")
  expect_error(gauge_covar("normal", 0.05, 0.01, "saturated", 1.5, p, 0.5),
               "'iteration' must be an integer")
  expect_error(gauge_covar("normal", 0.01, 0.01, "saturated", 0, p, "0.5"),
               "'split' must be a numeric vector of length 1")
  expect_error(gauge_covar("normal", 0.01, 0.01, "saturated", 0, p, c(0.2, 0.5)),
               "'split' must be a numeric vector of length 1")
  expect_error(gauge_covar("normal", 0.01, 0.01, "saturated", 0, p, -1.5),
               "'split' must lie strictly between 0 and 1")
  expect_error(gauge_covar("normal", 0.01, 0.01, "saturated", 0, p, 2.4),
               "'split' must lie strictly between 0 and 1")
  expect_error(gauge_covar("normal", 0.01, 0.01, "saturated", 1, p, NULL),
               "'split' cannot be NULL unless initial estimator is 'robustified'")
  expect_silent(gauge_covar("normal", 0.01, 0.05, "robustified", 1, p, NULL))

  # expect error if specify saturated, split != 0.5, m >= 1
  expect_error(gauge_covar("normal", 0.01, 0.05, "saturated", iteration = 1,
                           split = 0.3, parameters = p),
                           "No theory for m >= 1 if saturated initial estimator is not split in half")
  expect_error(gauge_covar("normal", 0.01, 0.05, "saturated", iteration = 10,
                           split = 0.3, parameters = p),
               "No theory for m >= 1 if saturated initial estimator is not split in half")



})

test_that("gauge_covar() and gauge_avar() coincide for s=t=c", {

  p <- generate_param(3, 2, 3, sigma = 2, intercept = TRUE, seed = 42)

  avar1 <- gauge_avar(ref_dist = "normal", sign_level = 0.01,
                      initial_est = "robustified", iteration = 0,
                      parameters = p, split = NULL)
  covar1 <- gauge_covar(ref_dist = "normal", sign_level = 0.01,
                        sign_level2 = 0.01, initial_est = "robustified",
                        iteration = 0, parameters = p, split = NULL)
  expect_identical(avar1, covar1)
  avar2 <- gauge_avar(ref_dist = "normal", sign_level = 0.05,
                      initial_est = "robustified", iteration = 0,
                      parameters = p, split = NULL)
  covar2 <- gauge_covar(ref_dist = "normal", sign_level = 0.05,
                        sign_level2 = 0.05, initial_est = "robustified",
                        iteration = 0, parameters = p, split = NULL)
  expect_identical(avar2, covar2)
  avar3 <- gauge_avar(ref_dist = "normal", sign_level = 0.01,
                      initial_est = "robustified", iteration = "convergence",
                      parameters = p, split = NULL)
  covar3 <- gauge_covar(ref_dist = "normal", sign_level = 0.01,
                        sign_level2 = 0.01, initial_est = "robustified",
                        iteration = "convergence", parameters = p, split = NULL)
  expect_identical(avar3, covar3)
  avar4 <- gauge_avar(ref_dist = "normal", sign_level = 0.01,
                      initial_est = "robustified", iteration = 5,
                      parameters = p, split = NULL)
  covar4 <- gauge_covar(ref_dist = "normal", sign_level = 0.01,
                        sign_level2 = 0.01, initial_est = "robustified",
                        iteration = 5, parameters = p, split = NULL)
  expect_identical(avar4, covar4)
  avar5 <- gauge_avar(ref_dist = "normal", sign_level = 0.01,
                      initial_est = "saturated", iteration = 0,
                      parameters = p, split = 0.5)
  covar5 <- gauge_covar(ref_dist = "normal", sign_level = 0.01,
                        sign_level2 = 0.01, initial_est = "saturated",
                        iteration = 0, parameters = p, split = 0.5)
  expect_identical(avar5, covar5)
  avar6 <- gauge_avar(ref_dist = "normal", sign_level = 0.01,
                      initial_est = "saturated", iteration = 0,
                      parameters = p, split = 0.25)
  covar6 <- gauge_covar(ref_dist = "normal", sign_level = 0.01,
                        sign_level2 = 0.01, initial_est = "saturated",
                        iteration = 0, parameters = p, split = 0.25)
  expect_identical(avar6, covar6)
  avar7 <- gauge_avar(ref_dist = "normal", sign_level = 0.01,
                      initial_est = "saturated", iteration = 0,
                      parameters = p, split = 0.33)
  covar7 <- gauge_covar(ref_dist = "normal", sign_level = 0.01,
                        sign_level2 = 0.01, initial_est = "saturated",
                        iteration = 0, parameters = p, split = 0.33)
  expect_identical(avar7, covar7)
  avar8 <- gauge_avar(ref_dist = "normal", sign_level = 0.01,
                      initial_est = "saturated", iteration = 4,
                      parameters = p, split = 0.5)
  covar8 <- gauge_covar(ref_dist = "normal", sign_level = 0.01,
                        sign_level2 = 0.01, initial_est = "saturated",
                        iteration = 4, parameters = p, split = 0.5)
  expect_identical(avar8, covar8)
  avar9 <- gauge_avar(ref_dist = "normal", sign_level = 0.01,
                      initial_est = "saturated", iteration = "convergence",
                      parameters = p, split = 0.5)
  covar9 <- gauge_covar(ref_dist = "normal", sign_level = 0.01,
                        sign_level2 = 0.01, initial_est = "saturated",
                        iteration = "convergence", parameters = p, split = 0.5)
  expect_identical(avar9, covar9)
  # following setting is not completely identical, differ in digit 15
  avar10 <- gauge_avar(ref_dist = "normal", sign_level = 0.1,
                       initial_est = "saturated", iteration = "convergence",
                       parameters = p, split = 0.5)
  covar10 <- gauge_covar(ref_dist = "normal", sign_level = 0.1,
                         sign_level2 = 0.1, initial_est = "saturated",
                         iteration = "convergence", parameters = p, split = 0.5)
  expect_equal(avar10, covar10, tolerance = 0.000000000000001)

})