tests/testthat/test-rma-uni-methods.R

In wildmeta: Cluster Wild Bootstrapping for Meta-Analysis

suppressPackageStartupMessages(library(metafor))
suppressPackageStartupMessages(library(clubSandwich))

set.seed(20221002)
data("oswald2013", package = "robumeta")

# Create variables
oswald2013$yi <- atanh(oswald2013$R)
oswald2013$vi <- 1 / (oswald2013$N - 3)
oswald2013$esID <- 1:nrow(oswald2013)
oswald2013$wt <- 1 + rpois(nrow(oswald2013), lambda = 1)
table(oswald2013$wt)
oswald2013 <- oswald2013

# Use a subset of 30 studies
oswald2013_full <- oswald2013
study_sample <- sample(unique(oswald2013$Study), size = 30)
oswald2013 <- subset(oswald2013, Study %in% study_sample)


# Fit some models

mod_A <- rma.uni(yi = yi ~ 1, vi = vi, data = oswald2013)

mod_B <- rma.uni(yi = yi, mods = ~ 0 + Crit.Cat, vi = vi,
                 data = oswald2013)

mod_C <- rma.uni(yi ~ Crit.Cat + Crit.Domain + IAT.Focus + Scoring, vi = vi,
                 data = oswald2013)

mod_D <- rma.uni(yi = yi, mods = ~ Crit.Domain + IAT.Focus + Scoring, vi = vi,
                 data = oswald2013)

mod_F <- rma.uni(yi ~ Crit.Cat + Crit.Domain + IAT.Focus, vi = vi,
                 data = oswald2013)

Cmat_A <- constrain_equal("Crit.Cat", reg_ex = TRUE, coef(mod_B))
Cmat_B <- constrain_zero(7:10, coef(mod_C))
Cmat_D <- constrain_zero("Crit.Cat", reg_ex = TRUE, coef(mod_C))
Cmat_E <- constrain_zero("Crit.Domain", reg_ex = TRUE, coef(mod_C))
Cmat_F <- constrain_zero("Scoring", reg_ex = TRUE, coef(mod_C))

tol <- 1e-6

test_that("estimate_null() works for rma.uni objects.", {

  mod_A_con <- estimate_null(mod_B, C_mat = Cmat_A)
  expect_equal(get_fitted(mod_A), get_fitted(mod_A_con), tolerance = tol)
  expect_equal(get_res(mod_A), get_res(mod_A_con), tolerance = tol)

  mod_D_con <- estimate_null(mod_C, Cmat_D)
  expect_equal(get_fitted(mod_D), get_fitted(mod_D_con), tolerance = tol)
  expect_equal(get_res(mod_D), get_res(mod_D_con), tolerance = tol)

  mod_B_con <- estimate_null(mod_C, Cmat_B)
  expect_equal(get_fitted(mod_B), get_fitted(mod_B_con), tolerance = tol)
  expect_equal(get_res(mod_B), get_res(mod_B_con), tolerance = tol)

  oswald2013_mod <- oswald2013
  oswald2013_mod$X_null <- 2

  mod_G <- update(mod_C, data = oswald2013_mod)
  mod_F_Xnull <- estimate_null(mod_G, Cmat_F)
  expect_equal(get_fitted(mod_F), get_fitted(mod_F_Xnull), tolerance = tol)
  expect_equal(get_res(mod_F), get_res(mod_F_Xnull), tolerance = tol)

})

test_that("get_cluster() works for rma.uni objects.", {
  study_fac <- factor(1:nrow(oswald2013))
  expect_equal(get_cluster(mod_A), study_fac, tolerance = tol)
  expect_equal(get_cluster(mod_B), study_fac, tolerance = tol)
  expect_equal(get_cluster(mod_C), study_fac, tolerance = tol)
  expect_equal(get_cluster(mod_D), study_fac, tolerance = tol)
  expect_equal(get_cluster(mod_F), study_fac, tolerance = tol)
})

test_that("get_boot_F_f() works for rma.uni objects.", {

  y_boot <- oswald2013$yi

  Cmat_int <- constrain_zero(1, coefs = coef(mod_A))

  f_A <- get_boot_F_f(mod_A, C_mat = Cmat_int)
  expect_equal(
    f_A(y_boot = y_boot, cluster = get_cluster(mod_A)),
    Wald_test(mod_A, cluster = get_cluster(mod_A),
              constraints = Cmat_int, vcov = "CR0", test = "Naive-F")$Fstat, tolerance = tol
  )

  f_B <- get_boot_F_f(mod_B, C_mat = Cmat_A, type = "CR1", test = "HTZ")
  expect_equal(
    f_B(y_boot = y_boot, cluster = get_cluster(mod_B)),
    Wald_test(mod_B, cluster = get_cluster(mod_B),
              constraints = Cmat_A, vcov = "CR1", test = "HTZ")$Fstat, tolerance = tol
  )

  f_C1 <- get_boot_F_f(mod_C, C_mat = Cmat_B, type = "CR2", test = "EDT")
  expect_equal(
    f_C1(y_boot = y_boot, cluster = get_cluster(mod_C)),
    Wald_test(mod_C, cluster = get_cluster(mod_C),
              constraints = Cmat_B, vcov = "CR2", test = "EDT")$Fstat, tolerance = tol
  )

  f_C2 <- get_boot_F_f(mod_C, C_mat = Cmat_D, type = "CR3", test = "chi-sq")
  expect_equal(
    f_C2(y_boot = y_boot, cluster = get_cluster(mod_C)),
    Wald_test(mod_C, cluster = get_cluster(mod_C), constraints = Cmat_D, vcov = "CR3", test = "chi-sq")$Fstat, tolerance = tol
  )

  f_C3 <- get_boot_F_f(mod_C, C_mat = Cmat_E, type = "CR2", test = "Naive-Fp")
  expect_equal(
    f_C3(y_boot = y_boot, cluster = get_cluster(mod_C)),
    Wald_test(mod_C, cluster = get_cluster(mod_C), constraints = Cmat_E, vcov = "CR2", test = "Naive-Fp")$Fstat, tolerance = tol
  )

  f_C4 <- get_boot_F_f(mod_C, C_mat = Cmat_F, type = "CR0", test = "HTA")
  expect_equal(
    f_C4(y_boot = y_boot, cluster = get_cluster(mod_C)),
    Wald_test(mod_C, cluster = get_cluster(mod_C),
              constraints = Cmat_F, vcov = "CR0", test = "HTA")$Fstat, tolerance = tol
  )


  oswald2013$yi_boot <- oswald2013$yi + rnorm(nrow(oswald2013))

  mod_G <- rma.uni(yi = yi_boot, mods = ~ Crit.Cat + Crit.Domain + IAT.Focus, vi = vi,
                   data = oswald2013)
  mod_H <- rma.uni(yi = yi_boot ~ Crit.Cat + Crit.Domain + IAT.Focus, vi = vi,
                   data = oswald2013)
  expect_false(all(coef(mod_F) == coef(mod_G)))
  expect_equal(coef(mod_G), coef(mod_H), tolerance = tol)
  expect_false(all(y_boot == oswald2013$yi_boot))

  Cmat_G <- constrain_equal("Crit.Cat", reg_ex = TRUE, coefs = coef(mod_G))
  F_stat <- Wald_test(mod_F, cluster = get_cluster(mod_F),
                      constraints = Cmat_G, vcov = "CR0", test = "HTZ")$Fstat
  G_stat <- get_boot_F_f(mod_G, C_mat = Cmat_G, type = "CR0", test = "HTZ")(y_boot = y_boot, cluster = get_cluster(mod_G))
  H_stat <- get_boot_F_f(mod_H, C_mat = Cmat_G, type = "CR0", test = "HTZ")(y_boot = y_boot, cluster = get_cluster(mod_H))

  expect_equal(G_stat, F_stat, tolerance = tol)

  expect_equal(H_stat, F_stat, tolerance = tol)

  expect_false(G_stat == Wald_test(mod_G, cluster = get_cluster(mod_G), constraints = Cmat_G, vcov = "CR0", test = "HTZ")$Fstat)

})


test_that("run_cwb options work for rma.uni objects.", {

  boot_yi_A <- run_cwb(mod_A, cluster = get_cluster(mod_A), R = 3,
                       auxiliary_dist = "Rademacher")
  fit_A <- get_fitted(mod_A)
  abs_res_A <- abs(get_res(mod_A))
  A_check <- sapply(boot_yi_A, function(x) all.equal(abs(x - fit_A), abs_res_A))
  expect_true(all(A_check))

  boot_yi_B <- run_cwb(mod_B, cluster = get_cluster(mod_B), R = 3,
                       auxiliary_dist = "Rademacher")
  fit_B <- get_fitted(mod_B)
  abs_res_B <- abs(get_res(mod_B))
  B_check <- sapply(boot_yi_B, function(x) all.equal(abs(x - fit_B), abs_res_B))
  expect_true(all(B_check))

  boot_yi_C <- run_cwb(mod_C, cluster = get_cluster(mod_C), R = 3,
                       auxiliary_dist = "Rademacher")
  fit_C <- get_fitted(mod_C)
  abs_res_C <- abs(get_res(mod_C))
  C_check <- sapply(boot_yi_C, function(x) all.equal(abs(x - fit_C), abs_res_C))
  expect_true(all(C_check))

})

compare_mod_results <- function(mod, scram, ord = 1:length(get_res(mod)), tol = 1e-6) {
  expect_equal(coef(mod), coef(scram), tolerance = tol)
  expect_equal(as.numeric(get_res(mod)[ord]), as.numeric(get_res(scram)), tolerance = tol)
  expect_equal(as.numeric(get_fitted(mod)[ord]), as.numeric(get_fitted(scram)), tolerance = tol)
}

test_that("Wald_test_cwb() results do not depend on sort order.", {

  skip_on_cran()

  tol <- 1e-5

  ord <- sample(1:nrow(oswald2013))
  oswald_scram <- oswald2013[ord,]

  scram_A <- rma.uni(yi = yi, vi = vi,
                     data = oswald_scram)
  compare_mod_results(mod_A, scram_A, ord)

  scram_B <- rma.uni(yi ~ 0 + Crit.Cat, vi = vi,
                     data = oswald_scram)
  compare_mod_results(mod_B, scram_B, ord)

  scram_C <- rma.uni(yi ~ Crit.Cat + Crit.Domain + IAT.Focus + Scoring, vi = vi,
                     data = oswald_scram)
  compare_mod_results(mod_C, scram_C, ord)

  orig_A <- Wald_test_cwb(mod_B, constraints = Cmat_A,
                          R = 4,
                          auxiliary_dist = "Rademacher",
                          adjust = "CR0",
                          type = "CR0",
                          test = "Naive-F",
                          seed = 1)

  scram_A <- Wald_test_cwb(scram_B, constraints = Cmat_A,
                           cluster = ord,
                           R = 4,
                           auxiliary_dist = "Rademacher",
                           adjust = "CR0",
                           type = "CR0",
                           test = "Naive-F",
                           seed = 1)

  expect_equal(attr(orig_A, "original"), attr(scram_A, "original"), tolerance = tol)
  expect_equal(attr(orig_A, "bootstraps"), attr(scram_A, "bootstraps"), tolerance = tol)


  orig_D <- Wald_test_cwb(mod_C, constraints = Cmat_D,
                          R = 3,
                          auxiliary_dist = "Mammen",
                          adjust = "CR1",
                          type = "CR1",
                          test = "Naive-Fp",
                          seed = 2)

  scram_D <- Wald_test_cwb(scram_C, constraints = Cmat_D,
                           cluster = ord,
                           R = 3,
                           auxiliary_dist = "Mammen",
                           adjust = "CR1",
                           type = "CR1",
                           test = "Naive-Fp",
                           seed = 2)

  expect_equal(attr(orig_D, "original"), attr(scram_D, "original"), tolerance = tol)
  expect_equal(attr(orig_D, "bootstraps"), attr(scram_D, "bootstraps"), tolerance = tol)

  orig_E <- Wald_test_cwb(mod_C, constraints = Cmat_E,
                          R = 5,
                          auxiliary_dist = "Webb six",
                          adjust = "CR2",
                          type = "CR2",
                          test = "HTZ",
                          seed = 3)

  scram_E <- Wald_test_cwb(scram_C, constraints = Cmat_E,
                           cluster = ord,
                           R = 5,
                           auxiliary_dist = "Webb six",
                           adjust = "CR2",
                           type = "CR2",
                           test = "HTZ",
                           seed = 3)

  expect_equal(attr(orig_E, "original"), attr(scram_E, "original"), tolerance = tol)
  expect_equal(attr(orig_E, "bootstraps"), attr(scram_E, "bootstraps"), tolerance = tol)

  orig_F <- Wald_test_cwb(mod_C, constraints = Cmat_F,
                          R = 8,
                          auxiliary_dist = "uniform",
                          adjust = "CR3",
                          type = "CR3",
                          test = "EDT",
                          seed = 4)

  scram_F <- Wald_test_cwb(scram_C, constraints = Cmat_F,
                           cluster = ord,
                           R = 8,
                           auxiliary_dist = "uniform",
                           adjust = "CR3",
                           type = "CR3",
                           test = "EDT",
                           seed = 4)

  expect_equal(attr(orig_F, "original"), attr(scram_F, "original"), tolerance = tol)
  expect_equal(attr(orig_F, "bootstraps"), attr(scram_F, "bootstraps"), tolerance = tol)

})


test_that("Wald_test_cwb() works when rma.uni uses subset.", {

  oswald_sub <- subset(oswald2013_full, Crit.Cat == "Microbehavior")

  mod_full <- rma.uni(yi = yi, vi = vi,
                      mods = ~ 0 + IAT.Focus + Crit.Domain,
                      data = oswald_sub)

  mod_sub <- rma.uni(yi = yi, vi = vi,
                     mods = ~ 0 + IAT.Focus + Crit.Domain,
                     data = oswald2013_full,
                     subset = Crit.Cat == "Microbehavior")

  compare_mod_results(mod_full, mod_sub)

  test_full <- Wald_test_cwb(mod_full,
                             constraints = constrain_equal("IAT.Focus", reg_ex = TRUE),
                             R = 3,
                             auxiliary_dist = "Rademacher",
                             adjust = "CR0",
                             type = "CR0",
                             test = "Naive-F",
                             seed = 19)

  test_sub <- Wald_test_cwb(mod_sub,
                            constraints = constrain_equal("IAT.Focus", reg_ex = TRUE),
                            R = 3,
                            auxiliary_dist = "Rademacher",
                            adjust = "CR0",
                            type = "CR0",
                            test = "Naive-F",
                            seed = 19)

  expect_equal(attr(test_full, "original"), attr(test_sub, "original"), tolerance = tol)
  expect_equal(attr(test_full, "bootstraps"), attr(test_sub, "bootstraps"), tolerance = tol)

  test_sub2 <- Wald_test_cwb(mod_sub,
                            constraints = constrain_equal("IAT.Focus", reg_ex = TRUE),
                            R = 3,
                            auxiliary_dist = "Rademacher",
                            adjust = "CR0",
                            type = "CR0",
                            test = "Naive-F",
                            seed = 19)

  expect_equal(attr(test_sub2, "original"), attr(test_sub, "original"), tolerance = tol)
  expect_equal(attr(test_sub2, "bootstraps"), attr(test_sub, "bootstraps"), tolerance = tol)

})

test_that("Wald_test_cwb() works when rma.uni has missing observations.", {

  set.seed(20230212)
  oswald_missX <- oswald2013_full
  oswald_missX$IAT.Focus[sample(1:nrow(oswald2013_full), size = 40)] <- NA
  oswald_sub <- subset(oswald_missX, Crit.Cat == "Microbehavior")
  oswald_complete <- subset(oswald_missX, Crit.Cat == "Microbehavior" & !is.na(IAT.Focus))

  mod_complete <- rma.uni(yi = yi, vi = vi,
                          mods = ~ 0 + IAT.Focus + Crit.Domain,
                          data = oswald_complete)

  expect_warning(
    mod_sub <- rma.uni(yi = yi, vi = vi,
                       mods = ~ 0 + IAT.Focus + Crit.Domain,
                       data = oswald_sub)
  )

  expect_warning(
    mod_missX <- rma.uni(yi = yi, vi = vi,
                         mods = ~ 0 + IAT.Focus + Crit.Domain,
                         data = oswald_missX,
                         subset = Crit.Cat == "Microbehavior")
  )

  compare_mod_results(mod_complete, mod_sub)
  compare_mod_results(mod_complete, mod_missX)

  test_complete <- Wald_test_cwb(mod_sub,
                                 constraints = constrain_equal("IAT.Focus", reg_ex = TRUE),
                                 R = 5,
                                 auxiliary_dist = "Rademacher",
                                 adjust = "CR0",
                                 type = "CR0",
                                 test = "Naive-F",
                                 seed = 19)

  test_sub <- Wald_test_cwb(mod_sub,
                            constraints = constrain_equal("IAT.Focus", reg_ex = TRUE),
                            R = 5,
                            auxiliary_dist = "Rademacher",
                            adjust = "CR0",
                            type = "CR0",
                            test = "Naive-F",
                            seed = 19)

  expect_equal(attr(test_complete, "original"), attr(test_sub, "original"), tolerance = tol)
  expect_equal(attr(test_complete, "bootstraps"), attr(test_sub, "bootstraps"), tolerance = tol)

  test_missX <- Wald_test_cwb(mod_missX,
                              constraints = constrain_equal("IAT.Focus", reg_ex = TRUE),
                              R = 5,
                              auxiliary_dist = "Rademacher",
                              adjust = "CR0",
                              type = "CR0",
                              test = "Naive-F",
                              seed = 19)

  expect_equal(attr(test_complete, "original"), attr(test_missX, "original"), tolerance = tol)
  expect_equal(attr(test_complete, "bootstraps"), attr(test_missX, "bootstraps"), tolerance = tol)

  test_missX2 <- Wald_test_cwb(mod_missX,
                               constraints = constrain_equal("IAT.Focus", reg_ex = TRUE),
                               cluster = oswald_missX$esID,
                               R = 5,
                               auxiliary_dist = "Rademacher",
                               adjust = "CR0",
                               type = "CR0",
                               test = "Naive-F",
                               seed = 19)

  expect_equal(attr(test_missX2, "original"), attr(test_missX, "original"), tolerance = tol)
  expect_equal(attr(test_missX2, "bootstraps"), attr(test_missX, "bootstraps"), tolerance = tol)

  test_sub2 <- Wald_test_cwb(mod_sub,
                             constraints = constrain_equal("IAT.Focus", reg_ex = TRUE),
                             cluster = oswald_sub$esID,
                             R = 5,
                             auxiliary_dist = "Rademacher",
                             adjust = "CR0",
                             type = "CR0",
                             test = "Naive-F",
                             seed = 19)

  expect_equal(attr(test_sub2, "original"), attr(test_sub, "original"), tolerance = tol)
  expect_equal(attr(test_sub2, "bootstraps"), attr(test_sub, "bootstraps"), tolerance = tol)


})



test_that("Wald_test_cwb() works with user-weighted rma.uni models.", {

  mod_wt1 <- rma.uni(yi ~ 0 + Crit.Cat + Crit.Domain + IAT.Focus + Scoring,
                     vi = vi, weights = wt,
                     data = oswald2013)


  test_wt1 <- Wald_test_cwb(mod_wt1,
                            constraints = constrain_equal("Crit.Cat", reg_ex = TRUE),
                            R = 3,
                            seed = 19)

  expect_s3_class(test_wt1, "Wald_test_wildmeta")
  expect_true(!is.na(test_wt1$p_val))

  mod_wt2 <- rma.uni(yi ~ 0 + Crit.Cat + Crit.Domain + IAT.Focus + Scoring,
                     vi = vi, weights = wt,
                     data = oswald2013)

  test_wt2 <- Wald_test_cwb(mod_wt2,
                            constraints = constrain_equal("Crit.Cat", reg_ex = TRUE),
                            R = 3,
                            seed = 19)

  expect_s3_class(test_wt2, "Wald_test_wildmeta")
  expect_true(!is.na(test_wt2$p_val))

  expect_equal(test_wt1, test_wt2, tolerance = tol)

})