tests/testthat/test-g_mlm.R

In scdhlm: Estimating Hierarchical Linear Models for Single-Case Designs

skip_if_not_installed("nlme")

# check if g_mlm() works

library(nlme)

test_that("g_mlm() is imported appropriately for Laski data.", {
  
  data(Laski)
  
  # simple model
  Laski_RML1 <- lme(fixed = outcome ~ treatment, 
                    random = ~ 1 | case, 
                    correlation = corAR1(0.01, ~ time | case), 
                    data = Laski)
  
  ## g_REML VS g_mlm
  Laski_g1_REML <- suppressWarnings(g_REML(Laski_RML1, p_const = c(0,1), r_const = c(1,0,1), returnModel = TRUE))
  Laski_g1_mlm <- g_mlm(Laski_RML1, p_const = c(0,1), r_const = c(1,0,1), returnModel = TRUE)
  
  expect_warning(g_REML(Laski_RML1, p_const = c(0,1), r_const = c(1,0,1), returnModel = TRUE))
  expect_equal(varcomp_vcov(Laski_RML1)[1, 1], 20214.7841577, tol = 1e-7)
  
  varcomp_ex <- extract_varcomp_lmeAR1(Laski_RML1) # extract varcomp using old function in scdhlm
  varcomp <- extract_varcomp(Laski_RML1) # extract varcomp using new function in lmeInfo
  expect_equal(varcomp_ex$sigma_sq, varcomp$sigma_sq)
  expect_equal(varcomp_ex$phi, varcomp$cor_params)
  expect_equal(varcomp_ex$Tau, varcomp$Tau$case)
  
  expect_output(summary(Laski_g1_REML))
  expect_output(print(Laski_g1_REML))
  expect_output(summary(Laski_g1_mlm))
  expect_output(print(Laski_g1_mlm))
  
  expect_equal(Laski_g1_mlm$p_beta, Laski_g1_REML$p_beta) # numerator of effect size
  expect_equal(Laski_g1_mlm$r_theta, Laski_g1_REML$r_theta) # squared denominator of effect size
  expect_equal(Laski_g1_mlm$delta_AB, Laski_g1_REML$delta_AB) # unadjusted (REML) effect size estimate
  expect_equal(as.numeric(Laski_g1_mlm$nu), Laski_g1_REML$nu) # degrees of freedom
  expect_equal(as.numeric(Laski_g1_mlm$kappa), Laski_g1_REML$kappa) # constant kappa
  expect_equal(as.numeric(Laski_g1_mlm$g_AB), Laski_g1_REML$g_AB) # corrected effect size estimate
  expect_equal(as.numeric(Laski_g1_mlm$SE_g_AB^2), Laski_g1_REML$V_g_AB) # Approximate variance estimate
  expect_equal(Laski_g1_mlm$theta$sigma_sq, Laski_g1_REML$sigma_sq) # Estimated level-1 variance
  expect_equal(Laski_g1_mlm$theta$cor_params, Laski_g1_REML$phi) # Estimated autocorrelation
  expect_equal(Laski_g1_mlm$theta$Tau$case, Laski_g1_REML$Tau) # Vector of level-2 variance components
  expect_equal(det(Laski_g1_mlm$info_inv), det(Laski_g1_REML$I_E_inv)) # Expected information matrix
  
  ## confidence interval
  expect_equal(CI_g(Laski_g1_REML), CI_g(Laski_g1_mlm))
  expect_equal(CI_g(Laski_g1_REML), lmeInfo:::CI_g.g_mlm(Laski_g1_mlm))
  expect_equal(CI_g(Laski_g1_REML, symmetric = FALSE), CI_g(Laski_g1_mlm, symmetric = FALSE))
  

  
  # complex model
  default_AB <- default_times(design = "MBP",
                              case = case, phase = treatment, session = time, 
                              data = Laski)
  D <- default_AB$B - default_AB$A
  
  dat <- preprocess_SCD(design = "MBP", 
                        case = case, phase = treatment, session = time, outcome = outcome, 
                        center = 4, data = Laski)
  Laski_RML2 <- suppressWarnings(lme(fixed = outcome ~ time + treatment + time_trt,
                                     random = ~ time | case, 
                                     correlation = corAR1(0.01, ~ time | case), 
                                     data = dat,
                                     control = lmeControl(msMaxIter = 50, apVar = FALSE, returnObject = TRUE)))
  
  ## g_REML VS g_mlm
  Laski_g2_REML <- suppressWarnings(g_REML(Laski_RML2, p_const = c(0,0,1,D), r_const = c(1,0,1,18,81), returnModel = TRUE))
  Laski_g2_mlm <- suppressWarnings(g_mlm(Laski_RML2, p_const = c(0,0,1,D), r_const = c(1,18,81,0,1), returnModel = TRUE))
  
  expect_equal(Laski_g2_mlm$p_beta, Laski_g2_REML$p_beta) # numerator of effect size
  expect_equal(Laski_g2_mlm$r_theta, Laski_g2_REML$r_theta) # squared denominator of effect size
  expect_equal(Laski_g2_mlm$delta_AB, Laski_g2_REML$delta_AB) # unadjusted (REML) effect size estimate
  expect_equal(as.numeric(Laski_g2_mlm$nu), Laski_g2_REML$nu) # degrees of freedom
  expect_equal(as.numeric(Laski_g2_mlm$kappa), Laski_g2_REML$kappa) # constant kappa
  expect_equal(as.numeric(Laski_g2_mlm$g_AB), Laski_g2_REML$g_AB) # corrected effect size estimate
  expect_equal(as.numeric(Laski_g2_mlm$SE_g_AB^2), Laski_g2_REML$V_g_AB) # Approximate variance estimate
  expect_equal(Laski_g2_mlm$theta$sigma_sq, Laski_g2_REML$sigma_sq) # Estimated level-1 variance
  expect_equal(Laski_g2_mlm$theta$cor_params, Laski_g2_REML$phi) # Estimated autocorrelation
  expect_equal(Laski_g2_mlm$theta$Tau$case, Laski_g2_REML$Tau) # Vector of level-2 variance components
  expect_equal(det(Laski_g2_mlm$info_inv), det(Laski_g2_REML$I_E_inv)) # Expected information matrix
  
  ## confidence interval
  expect_equal(CI_g(Laski_g2_REML), CI_g(Laski_g2_mlm))
  expect_equal(CI_g(Laski_g2_REML), lmeInfo:::CI_g.g_mlm(Laski_g2_mlm))
  expect_equal(CI_g(Laski_g2_REML, symmetric = FALSE), CI_g(Laski_g2_mlm, symmetric = FALSE))
  
})

test_that("g_mlm() is imported appropriately for Bryant 2018 data.", {
  
  data(Bryant2018)
  Bryant_RML1 <- lme(fixed = outcome ~ treatment,
                     random = ~ 1 | group/case,
                     correlation = corAR1(0.01, ~ session | group/case),
                     data = Bryant2018, 
                     na.action = na.omit)
  
  suppressWarnings(expect_error(g_REML(Bryant_RML1, p_const = c(0, 1), r_const = c(1, 0, 1, 1)))) # g_REML not available for 3-level data
  Bry_g1_mlm <- g_mlm(Bryant_RML1, p_const = c(0, 1), r_const = c(1, 1, 0, 1), infotype = "expected", returnModel = TRUE)
  
  expect_output(summary(Bry_g1_mlm))
  expect_output(print(Bry_g1_mlm))
  
})

test_that("Convergence indicator works correctly for g_mlm(), g_REML() and calc_BCSMD().", {
  
  data(Laski)
  
  # simple model
  Laski_BCSMD1 <- suppressWarnings(
    calc_BCSMD(design = "MBP", case = case, phase = treatment,
               session = time, outcome = outcome, center = 4,
               FE_base = 0, RE_base = 0, FE_trt = 0,
               data = Laski)
  )
  
  expect_identical(Laski_BCSMD1$Converged,"Yes")
  
  # complex model
  Laski_BCSMD2 <- suppressWarnings(
    calc_BCSMD(design = "MBP", case = case, phase = treatment,
               session = time, outcome = outcome, center = 4,
               FE_base = c(0,1), RE_base = c(0,1), FE_trt = c(0,1), RE_trt = c(0,1),
               data = Laski)
  )
  
  expect_identical(Laski_BCSMD2$Converged, "No")

})