tests/testthat/test-example_social_anxiety.R

In dmphillippo/multinma: Bayesian Network Meta-Analysis of Individual and Aggregate Data

# Generated by vignette example_social_anxiety.Rmd: do not edit by hand
# Instead edit example_social_anxiety.Rmd and then run precompile.R

skip_on_cran()



params <-
list(run_tests = FALSE, eval_multinomial = FALSE)

## ----code=readLines("children/knitr_setup.R"), include=FALSE------------------


## ----include=FALSE, setup-----------------------------------------------------
library(multinma)
library(dplyr)      # dplyr and tidyr for data manipulation
library(tidyr)
library(ggplot2)    # ggplot2 for plotting

## ----include=FALSE------------------------------------------------------------
options(mc.cores = parallel::detectCores())

## ----include=FALSE------------------------------------------------------------
nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), 
             "true" =, "warn" = 2, 
             parallel::detectCores())
options(mc.cores = nc)


## -----------------------------------------------------------------------------
head(social_anxiety)


## -----------------------------------------------------------------------------
sa_net <- set_agd_contrast(social_anxiety,
                           study = studyc, 
                           trt = trtc,
                           y = y, 
                           sample_size = 1,
                           se = se,
                           trt_class = classc,
                           trt_ref = "Waitlist")

sa_net


## -----------------------------------------------------------------------------
plot(sa_net, level = "class", weight_nodes = TRUE) + 
  theme(legend.position = "bottom", legend.box = "vertical")


## ----echo=TRUE, results='hide'------------------------------------------------
set.seed(951)
sa_fit_FE <- nma(sa_net,
                 trt_effects = "fixed",
                 prior_trt = normal(0, 100),
                 prior_het = half_normal(5),
)

sa_fit_RE <- nma(sa_net,
                 trt_effects = "random",
                 prior_trt = normal(0, 100),
                 prior_het = half_normal(5),
)


## -----------------------------------------------------------------------------
(sa_dic_FE <- dic(sa_fit_FE))
(sa_dic_RE <- dic(sa_fit_RE))


## ----echo=TRUE, results='hide'------------------------------------------------
sa_UME_RE <- nma(sa_net,
                 trt_effects = "random",
                 consistency = "ume",
                 prior_trt = normal(0, 100),
                 prior_het = half_normal(5))


## -----------------------------------------------------------------------------
(sa_dic_RE <- dic(sa_fit_RE))
(sa_dic_ume_RE <- dic(sa_UME_RE))

summary(sa_UME_RE, pars = "tau")
summary(sa_fit_RE, pars = "tau")


## -----------------------------------------------------------------------------
plot(sa_dic_RE, sa_dic_ume_RE, show_uncertainty = FALSE) +
  xlab("Residual deviance - No Class model") +
  ylab("Residual deviance - UME model")


## -----------------------------------------------------------------------------
as.data.frame(sa_dic_RE) %>%
  arrange(desc(resdev)) %>%
  head(5)


## ----echo=TRUE, results='hide', eval=!params$run_tests------------------------
# EMMELKAMP2006 <- data.frame(
#   Treatment_1 = c("CBT individual", "Waitlist", "Waitlist"),
#   Treatment_2 = c("Psychodynamic psychotherapy", "Psychodynamic psychotherapy", "CBT individual")
# )
# 
# ALDEN2011 <- data.frame(
#   Treatment_1 = c("Waitlist"),
#   Treatment_2 = c("CBT group")
# )
# 
# sa_fit_RE_nodesplit_EMMELKAMP <- nma(sa_net,
#                            consistency = "nodesplit",
#                            nodesplit = EMMELKAMP2006,
#                            trt_effects = "random",
#                            prior_trt = normal(0, 100),
#                            prior_het = half_normal(5),
# )
# 
# sa_fit_RE_nodesplit_ALDEN <- nma(sa_net,
#                              consistency = "nodesplit",
#                              nodesplit = ALDEN2011,
#                              trt_effects = "random",
#                              prior_trt = normal(0, 100),
#                              prior_het = half_normal(5),
# )


## ----eval=!params$run_tests---------------------------------------------------
# summary(sa_fit_RE_nodesplit_ALDEN)
# summary(sa_fit_RE_nodesplit_EMMELKAMP)


## ----echo=TRUE, results='hide'------------------------------------------------
sa_fit_EXclass_RE <- nma(sa_net,
                         trt_effects = "random",
                         prior_trt = normal(0, 100),
                         prior_het = half_normal(5),
                         class_effects = "exchangeable",
                         prior_class_mean = normal(0, 10),
                         prior_class_sd = normal(0.33,0.1),
                         class_sd = 
                           list(`Exercise and SH no support` = 
                                  c("Exercise promotion", "Self-help no support"),
                                `SSRIs and NSSA` = 
                                  c("SSRI/SNRI", "NSSA"),
                                `Psychodynamic & Other psychological therapies` =
                                  c("Psychodynamic psychotherapy", "Other psychological therapies"))
                         )


## -----------------------------------------------------------------------------
(sa_dic_EXclass_RE <- dic(sa_fit_EXclass_RE))
(sa_dic_RE <- dic(sa_fit_RE))

summary(sa_fit_RE, pars = "tau")
summary(sa_fit_EXclass_RE, pars = "tau")


## -----------------------------------------------------------------------------
plot(sa_dic_EXclass_RE, sa_dic_RE, show_uncertainty = FALSE) +
  xlab("Residual deviance - Exchangeable Class model") +
  ylab("Residual deviance - No Class model")


## ----echo=TRUE, results='hide'------------------------------------------------
sa_fit_COclass_RE <- nma(sa_net,
                         trt_effects = "random",
                         prior_trt = normal(0, 100),
                         prior_het = half_normal(5),
                         class_effects = "common")

sa_fit_EXclass_FE <- nma(sa_net,
                         trt_effects = "fixed",
                         prior_trt = normal(0, 100),
                         prior_het = half_normal(5),
                         class_effects = "exchangeable",
                         prior_class_mean = normal(0, 10),
                         prior_class_sd = normal(0.33,0.1),
                         class_sd = 
                           list(`Exercise and SH no support` = 
                                  c("Exercise promotion", "Self-help no support"),
                                `SSRIs and NSSA` = 
                                  c("SSRI/SNRI", "NSSA"),
                                `Psychodynamic & Other psychological therapies` = 
                                  c("Psychodynamic psychotherapy", "Other psychological therapies"))
                         )


## -----------------------------------------------------------------------------
(sa_dic_COclass_RE <- dic(sa_fit_COclass_RE))
(sa_dic_EXclass_FE <- dic(sa_fit_EXclass_FE))
(sa_dic_EXclass_RE <- dic(sa_fit_EXclass_RE))

summary(sa_fit_COclass_RE, pars = "tau")
summary(sa_fit_EXclass_RE, pars = "tau")


## -----------------------------------------------------------------------------
plot(sa_dic_COclass_RE, sa_dic_EXclass_RE, show_uncertainty = FALSE) +
  xlab("Residual deviance - Common Class model") +
  ylab("Residual deviance - Exchangeable Class model")


## -----------------------------------------------------------------------------
plot(relative_effects(sa_fit_EXclass_RE), ref_line = 0)


## -----------------------------------------------------------------------------
plot(sa_fit_EXclass_RE,
     pars = "class_mean",
     ref_line = 0)


## ----fig.height=8-------------------------------------------------------------
# Relative treatment effects
trt_eff <- as_tibble(relative_effects(sa_fit_EXclass_RE)) %>% 
  # Add in class details
  mutate(Class = sa_net$classes[as.numeric(.trtb)],
         level = "treatment")

# Class effects
class_eff <- as_tibble(summary(sa_fit_EXclass_RE, pars = "class_mean")) %>% 
  # Extract class details
  mutate(Class = factor(gsub(".*\\[(.+)\\]", "\\1", parameter), levels = levels(sa_net$classes)),
         level = "class",
         .trtb = factor("Class Mean", levels = c(levels(sa_net$classes), "Class Mean")))

# Combine and plot
bind_rows(trt_eff, class_eff) %>% 
  ggplot(aes(y = .trtb, 
             x = mean, xmin = `2.5%`, xmax = `97.5%`,
             colour = level, shape = level)) +
  geom_vline(xintercept = 0, colour = "grey60") +
  geom_pointrange() +
  facet_grid(rows = "Class", scales = "free", space = "free", labeller = label_wrap_gen(22)) +
  scale_shape_manual(values = c(15, 16), guide = guide_none()) +
  scale_colour_manual(values = c("#113259", "#55A480"), guide = guide_none()) +
  xlab("SMD") + ylab("") +
  theme_multinma() +
  theme(strip.text.y = element_text(angle = 0))


## -----------------------------------------------------------------------------
# Class means
EXclass_mean <- as.matrix(sa_fit_EXclass_RE, pars = "class_mean")

EXclass_mean <- cbind(`d[Reference]` = 0, EXclass_mean)

# Take ranks at each iteration
EXranks <- t(apply(EXclass_mean, 1, rank))

# Get median rank and 95% credible interval
EXresults <- t(apply(EXranks, 2, quantile, probs = c(0.025, 0.5, 0.975)))

# Convert to data frame
EXresults_df <- as.data.frame(EXresults)
EXresults_df$class <- rownames(EXresults_df)

EXresults_df$class <- factor(EXresults_df$class, 
                             levels = sort(unique(EXresults_df$class), decreasing = TRUE))

ggplot(EXresults_df, aes(x=class, y=`50%`, ymin=`2.5%`, ymax=`97.5%`)) +
geom_pointrange(size = 0.5) +
coord_flip() +
xlab("Class") + ylab("Posterior Ranks") +
theme_multinma()


## -----------------------------------------------------------------------------
EXranks_df <- as.data.frame(EXranks)

# Rank probabilities for class
rank_probs_EX <- apply(EXranks_df, 2, function(x) table(factor(x, levels = 1:ncol(EXranks_df))) / nrow(EXranks_df))

# Convert to data frame
rank_probs_df_EX <- as.data.frame(rank_probs_EX)

# Convert the data frame to a long format
rank_probs_long_EX <- rank_probs_df_EX %>%
  mutate(Rank = row_number()) %>%
  pivot_longer(
    cols = -Rank,
    names_to = "Class",
    values_to = "Probability"
  )

# Plot density
ggplot(rank_probs_long_EX, aes(x = Rank, y = Probability)) +
  geom_line() +
  facet_wrap(~ Class) +
  theme_multinma() +
  labs(x = "Rank",
       y = "Probability")


## ----smoking_tests, include=FALSE, eval=params$run_tests----------------------
#--- Test against TSD 4 results ---
library(testthat)
library(dplyr)

tol <- 0.05
tol_dic <- 0.1

# Relative effects
sa_EX_RE_releff <- as.data.frame(relative_effects(sa_fit_EXclass_RE))

# EX RE model parameters
test_ex_re <- tribble(
  ~parameter, ~mean, ~sd, ~`2.5%`, ~`50%`, ~`97.5%`,
  "d[Alprazolam]", -0.8314958, 0.28125207, -1.3755076, -0.8307162, -0.27711275,
  "d[CBT group]", -0.8142867, 0.09621211, -1.0065441, -0.8131570, -0.63148837,
  "d[CBT group + Fluoxetine]", -0.9211830, 0.19446183, -1.3156687, -0.9204697, -0.54765137,
  "d[CBT group + Moclobemide]", -1.2066027, 0.25000952, -1.6962705, -1.213031, -0.70271794,
  "d[CBT group + Phenelzine]", -1.6973883, 0.22347273, -2.1374866, -1.6956205, -1.26489286,
  "d[CBT group Enhanced]", -1.0473793, 0.21208307, -1.4793471, -1.0441864, -0.64037100,
  "d[CBT group Heimberg]", -0.7791976, 0.11644941, -1.0081920, -0.7790447, -0.55037257,
  "d[CBT individual]", -1.1730443, 0.15200737, -1.4664337, -1.1747683, -0.87424063,
  "d[CBT individual Heimberg]", -1.0210418, 0.19936795, -1.4083739, -1.0193569, -0.62629645,
  "d[Citalopram]", -0.8467242, 0.19163580, -1.2173637, -0.8569810, -0.43683285,
  "d[Clonazapam]", -1.0556030, 0.19780158, -1.4377446, -1.0581377, -0.67179776,
  "d[Cognitive therapy]", -1.5262932, 0.15985534, -1.8486260, -1.5249683, -1.21370893,
  "d[Cognitive therapy shortened sessions]", -0.9947664, 0.12715601, -1.2402119, -0.9965832, -0.73904361,
  "d[Escitalopram]", -0.8607573, 0.15488520, -1.1585547, -0.8654861, -0.54829102,
  "d[Exercise promotion]", -0.1805154, 0.40643886, -0.9872490, -0.1819663, 0.60912243,
  "d[Exposure in vivo]", -0.8190751, 0.13069686, -1.0765767, -0.8175697, -0.56103203,
  "d[Fluoxetine]", -0.8576191, 0.14304817, -1.1261241, -0.8605516, -0.57023950,
  "d[Fluvoxamine]", -0.9064442, 0.15317111, -1.2144374, -0.9042166, -0.60930856,
  "d[Gabapentin]", -0.8511784, 0.26831004, -1.3791457, -0.8477352, -0.32055234,
  "d[Interpersonal psychotherapy]", -0.4179080, 0.20660264, -0.8173726, -0.4137257, -0.01862694,
  "d[Levetiracetam]", -0.7959018, 0.32429393, -1.4507038, -0.7888249, -0.16125339,
  "d[Mindfulness]", -0.3321746, 0.21466120, -0.7685138, -0.3299451, 0.09519575,
  "d[Mirtazapine]", -0.7650427, 0.33297935, -1.4222047, -0.7676135, -0.11613886,
  "d[Moclobemide]", -0.7108082, 0.15076437, -1.0025218, -0.7098811, -0.41162089,
  "d[Paroxetine]", -0.9518310, 0.13803974, -1.2327792, -0.9530097, -0.68226256,
  "d[Paroxetine + Clonazapam]", -1.3265332, 0.29791193, -1.9052410, -1.3210094, -0.74984083,
  "d[Phenelzine]", -1.2715405, 0.15435898, -1.5720287, -1.2731792, -0.95776860,
  "d[Pill placebo]", -0.4342747, 0.12462281, -0.6844952, -0.4348165, -0.18907267,
  "d[Pregabalin]", -0.6868730, 0.18355833, -1.0606219, -0.6856703, -0.33094991,
  "d[Psychodynamic + Clonazepam]", -1.2161722, 0.28910209, -1.7912116, -1.2176852, -0.65840233,
  "d[Psychodynamic psychotherapy]", -0.6235235, 0.16406580, -0.9554527, -0.6227220, -0.29926675,
  "d[Psychological placebo]", -0.6467276, 0.15004879, -0.9335005, -0.6508577, -0.35282347,
  "d[Self-help book no support]", -0.7338743, 0.12963766, -0.9942650, -0.7345158, -0.48176798,
  "d[Self-help book with support]", -0.8347549, 0.17192892, -1.1769224, -0.8370861, -0.49302079,
  "d[Self-help internet no support]", -0.6028340, 0.14218310, -0.8756893, -0.6037919, -0.32938047,
  "d[Self-help internet with support]", -0.8435101, 0.08756184, -1.0161868, -0.8422352, -0.67559595,
  "d[Sertraline]", -0.8870449, 0.15228323, -1.1909221, -0.8867806, -0.57684247,
  "d[Social skills training]", -0.8263210, 0.27080686, -1.3600444, -0.8290302, -0.29391996,
  "d[Supportive therapy]", -0.2475428, 0.23241985, -0.6929505, -0.2487845, 0.21323466,
  "d[Venlafaxine]", -0.9027961, 0.14619543, -1.1883951, -0.9010945, -0.61143773)

test_that("FE relative effects", {
  expect_equivalent(sa_EX_RE_releff$mean, test_ex_re$mean, tolerance = tol)
  expect_equivalent(sa_EX_RE_releff$sd, test_ex_re$sd, tolerance = tol)
  expect_equivalent(sa_EX_RE_releff$`2.5%`, test_ex_re$`2.5%`, tolerance = tol)
  expect_equivalent(sa_EX_RE_releff$`50%`, test_ex_re$`50%`, tolerance = tol)
  expect_equivalent(sa_EX_RE_releff$`97.5%`, test_ex_re$`97.5%`, tolerance = tol)
})

# DIC checks
test_that("FE DIC", {
  expect_equivalent(sa_dic_FE$resdev, 288.3, tolerance = tol_dic)
  expect_equivalent(sa_dic_FE$pd, 40.1, tolerance = tol_dic)
  expect_equivalent(sa_dic_FE$dic, 328.4, tolerance = tol_dic)
})

test_that("RE DIC", {
  expect_equivalent(sa_dic_RE$resdev, 162.6, tolerance = tol_dic)
  expect_equivalent(sa_dic_RE$pd, 94.8, tolerance = tol_dic)
  expect_equivalent(sa_dic_RE$dic, 257.4, tolerance = tol_dic)
})

test_that("UME RE DIC", {
  expect_equivalent(sa_dic_ume_RE$resdev, 160.7, tolerance = tol_dic)
  expect_equivalent(sa_dic_ume_RE$pd, 108.3, tolerance = tol_dic)
  expect_equivalent(sa_dic_ume_RE$dic, 269, tolerance = tol_dic)
})

test_that("EX RE DIC", {
  expect_equivalent(sa_dic_EXclass_RE$resdev, 162.5, tolerance = tol_dic)
  expect_equivalent(sa_dic_EXclass_RE$pd, 87.6, tolerance = tol_dic)
  expect_equivalent(sa_dic_EXclass_RE$dic, 250.1, tolerance = tol_dic)
})

test_that("CO RE DIC", {
  expect_equivalent(sa_dic_COclass_RE$resdev, 158.4, tolerance = tol_dic)
  expect_equivalent(sa_dic_COclass_RE$pd, 93.1, tolerance = tol_dic)
  expect_equivalent(sa_dic_COclass_RE$dic, 251.5, tolerance = tol_dic)
})

test_that("EX FE DIC", {
  expect_equivalent(sa_dic_EXclass_FE$resdev, 285, tolerance = tol_dic)
  expect_equivalent(sa_dic_EXclass_FE$pd, 34.5, tolerance = tol_dic)
  expect_equivalent(sa_dic_EXclass_FE$dic, 319.4, tolerance = tol_dic)
})

#--- Tau checks ---
CO_RE_tau <- summary(sa_fit_COclass_RE, pars = "tau")
EX_RE_tau <- summary(sa_fit_EXclass_RE, pars = "tau")
RE_tau <- summary(sa_fit_RE, pars = "tau")
UME_RE_tau <- summary(sa_UME_RE, pars = "tau")

test_that("CO RE Tau", {
  expect_equivalent(CO_RE_tau$summary$mean, 0.252, tolerance = tol)
  expect_equivalent(CO_RE_tau$summary$sd, 0.0291, tolerance = tol)
  expect_equivalent(CO_RE_tau$summary$`2.5%`, 0.198, tolerance = tol)
  expect_equivalent(CO_RE_tau$summary$`50%`, 0.251, tolerance = tol)
  expect_equivalent(CO_RE_tau$summary$`97.5%`, 0.312, tolerance = tol)
})

test_that("EX RE Tau", {
  expect_equivalent(EX_RE_tau$summary$mean, 0.196, tolerance = tol)
  expect_equivalent(EX_RE_tau$summary$sd, 0.0298, tolerance = tol)
  expect_equivalent(EX_RE_tau$summary$`2.5%`, 0.142, tolerance = tol)
  expect_equivalent(EX_RE_tau$summary$`50%`, 0.195, tolerance = tol)
  expect_equivalent(EX_RE_tau$summary$`97.5%`, 0.257, tolerance = tol)
})

test_that("RE Tau", {
  expect_equivalent(RE_tau$summary$mean, 0.207, tolerance = tol)
  expect_equivalent(RE_tau$summary$sd, 0.0308, tolerance = tol)
  expect_equivalent(RE_tau$summary$`2.5%`, 0.149, tolerance = tol)
  expect_equivalent(RE_tau$summary$`50%`, 0.207, tolerance = tol)
  expect_equivalent(RE_tau$summary$`97.5%`, 0.269, tolerance = tol)
})

test_that("UME RE Tau", {
  expect_equivalent(EX_RE_tau$summary$mean, 0.221, tolerance = tol)
  expect_equivalent(EX_RE_tau$summary$sd, 0.0333, tolerance = tol)
  expect_equivalent(EX_RE_tau$summary$`2.5%`, 0.158, tolerance = tol)
  expect_equivalent(EX_RE_tau$summary$`50%`, 0.221, tolerance = tol)
  expect_equivalent(EX_RE_tau$summary$`97.5%`, 0.289, tolerance = tol)
})

# nodesplit_ALDEN <- summary(sa_fit_RE_nodesplit_ALDEN)
# nodesplit_EMMELKAMP <- summary(sa_fit_RE_nodesplit_EMMELKAMP)
# 
# ALDEN_test <- tribble(
#   ~parameter, ~mean, ~sd, ~`2.5%`, ~`50%`, ~`97.5%`,
#   "d_dir[CBT group vs. Waitlist]", -0.90, 0.13, -1.16, -0.90, -0.64,
#   "d_ind[CBT group vs. Waitlist]", -0.73, 0.13, -0.99, -0.73, -0.46,
#   "omega", -0.17, 0.17, -0.51, -0.17, 0.18,
#   "tau", 0.21, 0.03, 0.15, 0.21, 0.27
# )
# 
# 
# EMMELKAMP_test1 <- tribble(
#   ~parameter, ~mean, ~sd, ~`2.5%`, ~`25%`, ~`50%`, ~`75%`, ~`97.5%`, ~Bulk_ESS, ~Tail_ESS, ~Rhat,
#   "d_dir", 0.71, 0.33, 0.06, 0.49, 0.71, 0.94, 1.35, 5354, 2759, 1.00,
#   "d_ind", 0.45, 0.26, -0.07, 0.27, 0.45, 0.63, 0.96, 4159, 3058, 1.00,
#   "omega", 0.26, 0.39, -0.51, 0.00, 0.25, 0.53, 1.02, 4849, 3147, 1.00,
#   "tau", 0.21, 0.03, 0.15, 0.19, 0.21, 0.23, 0.27, 1291, 2516, 1.01
# )
# 
# EMMELKAMP_test2 <- tribble(
#   ~parameter, ~mean, ~sd, ~`2.5%`, ~`25%`, ~`50%`, ~`75%`, ~`97.5%`, ~Bulk_ESS, ~Tail_ESS, ~Rhat,
#   "d_dir", -0.61, 0.19, -0.96, -0.73, -0.61, -0.48, -0.22, 5103, 3262, 1.00,
#   "d_ind", -0.68, 0.42, -1.49, -0.96, -0.68, -0.40, 0.15, 4157, 3138, 1.00,
#   "omega", 0.08, 0.46, -0.84, -0.23, 0.07, 0.39, 0.98, 4293, 2877, 1.00,
#   "tau", 0.21, 0.03, 0.15, 0.19, 0.21, 0.23, 0.27, 959, 1921, 1.00
# )
# 
# EMMELKAMP_test3 <- tribble(
#   ~parameter, ~mean, ~sd, ~`2.5%`, ~`25%`, ~`50%`, ~`75%`, ~`97.5%`, ~Bulk_ESS, ~Tail_ESS, ~Rhat,
#   "d_dir", -0.84, 0.25, -1.32, -1.00, -0.84, -0.67, -0.36, 4676, 3392, 1,
#   "d_ind", -1.46, 0.24, -1.93, -1.62, -1.46, -1.30, -0.99, 1466, 2413, 1,
#   "omega", 0.62, 0.34, -0.04, 0.39, 0.62, 0.85, 1.31, 1999, 2802, 1,
#   "tau", 0.20, 0.03, 0.15, 0.18, 0.20, 0.22, 0.26, 1238, 1834, 1
# )
# 
# test_that("Node-splitting estimates_ALDEN", {
#   expect_equal(ALDEN_test$mean, as.numeric(nodesplit_ALDEN$summary[[1]]$summary$mean), tolerance = tol)
#   expect_equal(ALDEN_test$sd, as.numeric(nodesplit_ALDEN$summary[[1]]$summary$sd), tolerance = tol)
#   expect_equal(ALDEN_test$`2.5%`, as.numeric(nodesplit_ALDEN$summary[[1]]$summary$`2.5%`), tolerance = tol)
#   expect_equal(ALDEN_test$`50%`, as.numeric(nodesplit_ALDEN$summary[[1]]$summary$`50%`), tolerance = tol)
#   expect_equal(ALDEN_test$`97.5%`, as.numeric(nodesplit_ALDEN$summary[[1]]$summary$`97.5%`), tolerance = tol)
# })
# 
# test_that("Node-splitting estimates_EMMEL1", {
#   expect_equal(EMMELKAMP_test1$mean, as.numeric(nodesplit_EMMELKAMP$summary[[1]]$summary$mean), tolerance = tol)
#   expect_equal(EMMELKAMP_test1$sd, as.numeric(nodesplit_EMMELKAMP$summary[[1]]$summary$sd), tolerance = tol)
#   expect_equal(EMMELKAMP_test1$`2.5%`, as.numeric(nodesplit_EMMELKAMP$summary[[1]]$summary$`2.5%`), tolerance = tol)
#   expect_equal(EMMELKAMP_test1$`50%`, as.numeric(nodesplit_EMMELKAMP$summary[[1]]$summary$`50%`), tolerance = tol)
#   expect_equal(EMMELKAMP_test1$`97.5%`, as.numeric(nodesplit_EMMELKAMP$summary[[1]]$summary$`97.5%`), tolerance = tol)
# })
# 
# test_that("Node-splitting estimates_EMMEL2", {
#   expect_equal(EMMELKAMP_test2$mean, as.numeric(nodesplit_EMMELKAMP$summary[[2]]$summary$mean), tolerance = tol)
#   expect_equal(EMMELKAMP_test2$sd, as.numeric(nodesplit_EMMELKAMP$summary[[2]]$summary$sd), tolerance = tol)
#   expect_equal(EMMELKAMP_test2$`2.5%`, as.numeric(nodesplit_EMMELKAMP$summary[[2]]$summary$`2.5%`), tolerance = tol)
#   expect_equal(EMMELKAMP_test2$`50%`, as.numeric(nodesplit_EMMELKAMP$summary[[2]]$summary$`50%`), tolerance = tol)
#   expect_equal(EMMELKAMP_test2$`97.5%`, as.numeric(nodesplit_EMMELKAMP$summary[[2]]$summary$`97.5%`), tolerance = tol)
# })
# 
# test_that("Node-splitting estimates_EMMEL3", {
#   expect_equal(EMMELKAMP_test3$mean, as.numeric(nodesplit_EMMELKAMP$summary[[3]]$summary$mean), tolerance = tol)
#   expect_equal(EMMELKAMP_test3$sd, as.numeric(nodesplit_EMMELKAMP$summary[[3]]$summary$sd), tolerance = tol)
#   expect_equal(EMMELKAMP_test3$`2.5%`, as.numeric(nodesplit_EMMELKAMP$summary[[3]]$summary$`2.5%`), tolerance = tol)
#   expect_equal(EMMELKAMP_test3$`50%`, as.numeric(nodesplit_EMMELKAMP$summary[[3]]$summary$`50%`), tolerance = tol)
#   expect_equal(EMMELKAMP_test3$`97.5%`, as.numeric(nodesplit_EMMELKAMP$summary[[3]]$summary$`97.5%`), tolerance = tol)
# })


# Force clean up
rm(list = ls())
gc()