# skip_on_cran() # # if (!requireNamespace("cmdstanr", quietly = TRUE)) { # backend <- "rstan" # ## if using rstan backend, models can crash on Windows and MAC OS # ## so skip if on windows and cannot use cmdstanr # skip_on_os("windows") # } else { # if (isFALSE(is.null(cmdstanr::cmdstan_version(error_on_NA = FALSE)))) { # backend <- "cmdstanr" # } # } # # # Packages # library(testthat) # library(data.table) # library(multilevelcoda) # library(extraoperators) # library(brms) # library(lme4) # # # model # #--------------------------------------------------------------------------------------------------- # data(mcompd) # data(sbp) # data(psub) # # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID", total = 1440) # # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + bilr2 + bilr3 + bilr4 + # wilr1 + wilr2 + wilr3 + wilr4 + Female + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # foreach::registerDoSEQ() # # x <- wsub(object = m, basesub = psub, delta = 2) # # Testing #--------------------------------------------------------------------------------------------------- # test_that("wsub errors for invalid input", { # # ## missing object # expect_error(x <- wsub(basesub = psub, delta = 2)) # # ## missing basesub # expect_error(x <- wsub(object = m, delta = 2)) # # ## not brmcoda model # m1 <- lmer(Stress ~ 1 + (1 | ID), data = mcompd) # expect_error(x <- wsub(object = m1, basesub = psub, delta = 2)) # # ## invalid delta # expect_error(x <- wsub(object = m, basesub = psub, delta = -10)) # expect_error(x <- wsub(object = m, basesub = psub, delta = 1:10)) # # ## missing delta # expect_error(x <- substitution(object = m1, basesub = psub)) # # ## basesub does not have the same components as parts in cilr # ps <- build.basesub(c("WAKE", "MVPA", "LPA", "SB")) # expect_error(x <- wsub(object = m, basesub = ps, delta = 2)) # # ## basesub does have the same names as parts in cilr # ps <- build.basesub(parts = c("Sleep", "WAKE", "MVPA", "LPA", "SB")) # expect_error(x <- wsub(object = m, basesub = ps, delta = 2)) # # }) # test_that("wsub works as expected for adjusted/unadjusted model", { # # ## reference grid is provided for unadjusted model # suppressWarnings( # m2 <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + bilr2 + bilr3 + bilr4 + # wilr1 + wilr2 + wilr3 + wilr4 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # rg <- data.table(Age = 1) # expect_warning(x <- wsub(object = m2, basesub = psub, delta = 2, regrid = rg)) # # ## incorect reference grid 1 # rg <- data.table(Age = 1) # expect_error(x <- wsub(object = m, basesub = psub, delta = 2, regrid = rg)) # # ## reference grid has matching names with ILRs # rg <- data.table(bilr1 = 1) # expect_error(x <- wsub(object = m, basesub = psub, delta = 2, regrid = rg)) # # ## incorect reference grid 2 # rg <- data.table(bilr1 = 1, Age = 1) # expect_error(x <- wsub(object = m, basesub = psub, delta = 2, regrid = rg)) # # # delta out of range # expect_error(x <- wsub(object = m, basesub = psub, delta = 1000)) # # ## function knows to use correct user's specified reference grid # rg <- data.table(Female = 1) # x3 <- wsub(object = m, basesub = psub, delta = 2, regrid = rg) # expect_true(all(x3$TST$Female == 1)) # expect_true(all(x3$WAKE$Female == 1)) # expect_true(all(x3$MVPA$Female == 1)) # expect_true(all(x3$LPA$Female == 1)) # expect_true(all(x$SB$Female == 1)) # # expect_true(all(x3$TST$Female != 0)) # expect_true(all(x3$WAKE$Female != 0)) # expect_true(all(x3$MVPA$Female != 0)) # expect_true(all(x3$LPA$Female != 0)) # expect_true(all(x3$SB$Female != 0)) # # ## model with unspecified reference grid works as expected # expect_equal(x$TST$Female, NULL) # expect_equal(x$WAKE$Female, NULL) # expect_equal(x$MVPA$Female, NULL) # expect_equal(x$LPA$Female, NULL) # expect_equal(x$SB$Female, NULL) # # ## model with unspecified reference grid works as expected # expect_true("Female" %nin% colnames(x$TST)) # expect_true("Female" %nin% colnames(x$WAKE)) # expect_true("Female" %nin% colnames(x$MVPA)) # expect_true("Female" %nin% colnames(x$LPA)) # expect_true("Female" %nin% colnames(x$SB)) # # ## average across reference grid as default # x4 <- wsub(object = m, basesub = psub, delta = 2, summary = TRUE) # x5 <- wsub(object = m, basesub = psub, delta = 2) # expect_equal(x4, x5) # # ## keep prediction at each level of refrence grid # cilr <- complr(data = mcompd[ID %in% c(1:5, 185:190), .SD[1:3], by = ID], sbp = sbp, # parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID", total = 1440) # # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + bilr2 + bilr3 + bilr4 + # wilr1 + wilr2 + wilr3 + wilr4 + Female + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # # x6 <- wsub(object = m, basesub = psub, delta = 2, summary = FALSE) # # expect_equal(nrow(x6$TST), nrow(x5$TST) * 2) # expect_equal(nrow(x6$WAKE), nrow(x5$WAKE) * 2) # expect_equal(nrow(x6$MVPA), nrow(x5$MVPA) * 2) # expect_equal(nrow(x6$LPA), nrow(x5$LPA) * 2) # expect_equal(nrow(x6$SB), nrow(x5$SB) * 2) # # expect_true("Female" %in% colnames(x6$TST)) # expect_true("Female" %in% colnames(x6$WAKE)) # expect_true("Female" %in% colnames(x6$MVPA)) # expect_true("Female" %in% colnames(x6$LPA)) # expect_true("Female" %in% colnames(x6$SB)) # # expect_true(all(x6$TST$Female %in% c(0, 1))) # expect_true(all(x6$WAKE$Female %in% c(0, 1))) # expect_true(all(x6$MVPA$Female %in% c(0, 1))) # expect_true(all(x6$LPA$Female %in% c(0, 1))) # expect_true(all(x6$SB$Female %in% c(0, 1))) # # }) # # test_that("wsub checks for user-specified reference composition", { # # # incorrect length # ref1 <- c(400, 60, 500, 60) # expect_error(wsub(object = m, basesub = psub, recomp = ref1, delta = 2)) # # # incorrect class # ref2 <- c("400", "100", "500", "200", "200") # expect_error(wsub(object = m, basesub = psub, recomp = ref2, delta = 2)) # # # incorrect class # ref3 <- c(400, 100, 500, 200, 200) # expect_error(x <- wsub(object = m, basesub = psub, recomp = ref3, delta = 2)) # # # values outside of possible range # ref4 <- c(100, 100, 900, 100, 240) # expect_error(x <- wsub(object = m, basesub = psub, recomp = ref4, delta = 2)) # # # include 0 # ref5 <- c(100, 200, 900, 0, 240) # expect_error(x <- wsub(object = m, basesub = psub, recomp = ref5, delta = 2)) # # }) # # test_that("wsub outputs what expected", { # # ## types # expect_type(x, "list") # expect_equal(length(x), length(m$complr$parts)) # expect_s3_class(x$TST, "data.table") # expect_s3_class(x$WAKE, "data.table") # expect_s3_class(x$MVPA, "data.table") # expect_s3_class(x$LPA, "data.table") # expect_s3_class(x$SB, "data.table") # # expect_type(x$TST$Mean, "double") # expect_type(x$TST$CI_low, "double") # expect_type(x$TST$CI_high, "double") # expect_type(x$TST$Delta, "double") # expect_type(x$TST$From, "character") # expect_type(x$TST$To, "character") # # expect_type(x$WAKE$Mean, "double") # expect_type(x$WAKE$CI_low, "double") # expect_type(x$WAKE$CI_high, "double") # expect_type(x$WAKE$Delta, "double") # expect_type(x$WAKE$From, "character") # expect_type(x$WAKE$To, "character") # # expect_type(x$MVPA$Mean, "double") # expect_type(x$MVPA$CI_low, "double") # expect_type(x$MVPA$CI_high, "double") # expect_type(x$MVPA$Delta, "double") # expect_type(x$MVPA$From, "character") # expect_type(x$MVPA$To, "character") # # expect_type(x$LPA$Mean, "double") # expect_type(x$LPA$CI_low, "double") # expect_type(x$LPA$CI_high, "double") # expect_type(x$LPA$Delta, "double") # expect_type(x$LPA$From, "character") # expect_type(x$LPA$To, "character") # # expect_type(x$SB$Mean, "double") # expect_type(x$SB$CI_low, "double") # expect_type(x$SB$CI_high, "double") # expect_type(x$SB$Delta, "double") # expect_type(x$SB$From, "character") # expect_type(x$SB$To, "character") # # expect_true(ncol(x$TST) >= 8) # expect_true(ncol(x$WAKE) >= 8) # expect_true(ncol(x$MVPA) >= 8) # expect_true(ncol(x$LPA) >= 8) # expect_true(ncol(x$SB) >= 8) # # expect_true(all(x$TST$To == "TST")) # expect_true(all(x$WAKE$To == "WAKE")) # expect_true(all(x$MVPA$To == "MVPA")) # expect_true(all(x$LPA$To == "LPA")) # expect_true(all(x$SB$To == "SB")) # # expect_true(all(x$TST$Level == "within")) # expect_true(all(x$WAKE$Level == "within")) # expect_true(all(x$MVPA$Level == "within")) # expect_true(all(x$LPA$Level == "within")) # expect_true(all(x$SB$Level == "within")) # # }) # # test_that("wsub gives results in sensible range", { # # ## difference in outcome # expect_true(x$TST$Mean %ae% "[-0.5, 0) | (0, 0.5]") # expect_true(x$WAKE$Mean %ae% "[-0.5, 0) | (0, 0.5]") # expect_true(x$MVPA$Mean %ae% "[-0.5, 0) | (0, 0.5]") # expect_true(x$LPA$Mean %ae% "[-0.5, 0) | (0, 0.5]") # expect_true(x$SB$Mean %ae% "[-0.5, 0) | (0, 0.5]") # # expect_true(x$TST$CI_low %ae% "[-1, 0) | (0, 1]") # expect_true(x$WAKE$CI_low %ae% "[-1, 0) | (0, 1]") # expect_true(x$MVPA$CI_low %ae% "[-1, 0) | (0, 1]") # expect_true(x$LPA$CI_low %ae% "[-1, 0) | (0, 1]") # expect_true(x$SB$CI_low %ae% "[-1, 0) | (0, 1]") # # expect_true(x$TST$CI_high %ae% "[-1, 0) | (0, 1]") # expect_true(x$WAKE$CI_high %ae% "[-1, 0) | (0, 1]") # expect_true(x$MVPA$CI_high %ae% "[-1, 0) | (0, 1]") # expect_true(x$LPA$CI_high %ae% "[-1, 0) | (0, 1]") # expect_true(x$SB$CI_high %ae% "[-1, 0) | (0, 1]") # # }) # # test_that("wsub gives results in expected direction and magnitude", { # # ## values are opposite sign for opposite substitution # for (i in seq_along(x)) { # expect_true(all(x[[i]][, sign(Mean[sign(Delta) == 1]) # %a!=% sign(Mean[sign(Delta) == -1]), by = From]$V1)) # } # # ## results for 1 min have smaller magnitude than 2 mins # for (i in seq_along(x)) { # expect_true(all(x[[i]][, abs(Mean[abs(Delta) == 1]) # < abs(Mean[abs(Delta) == 2])])) # } # }) # #--------------------------------------------------------------------------------------------------- # # Check the results from substitution model align with brm model # # using results from pairwise substitution # ## Estimates should be in the direction between pairwise coordinates and pairwise substitution # ## CIs should indicate consistent significance between pairwise coordinates and substitution # # sbp <- matrix(c( # 1, 1, -1,-1, -1, # 1, -1, 0, 0, 0, # 0, 0, -1, -1, 1, # 0, 0, 1, -1, 0), ncol = 5, byrow = TRUE) # # cilr <- complr(data = mcompd, sbp = sbp, # parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID", total = 1440) # # suppressWarnings(m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + bilr2 + bilr3 + bilr4 + # wilr1 + wilr2 + wilr3 + wilr4 + Female + (1 | ID), # seed = 123)) # # a <- wsub(object = m, basesub = psub, delta = 2, summary = TRUE) # # test_that("wsub's estimates matches with brm model's (TST vs WAKE and MVPA vs LPA)", { # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[7, 1] > 0))) { # wilr2 = more TST less WAKE # expect_true(all(a$TST[From == "WAKE" & Delta > 1]$Mean > 0)) # more TST less WAKE # expect_true(all(a$WAKE[From == "TST" & Delta > 1]$Mean < 0)) # more WAKE less TST # # } else { # expect_true(all(a$TST[From == "WAKE" & Delta > 1]$Mean < 0)) # expect_true(all(a$WAKE[From == "TST" & Delta > 1]$Mean > 0)) # # } # # if (isTRUE(suppressWarnings(summary(m$model)$fixed[9, 1] > 0))) { # wilr4 = more MVPA less LPA # expect_true(all(a$MVPA[From == "LPA" & Delta > 1]$Mean > 0)) # more MVPA less LPA # expect_true(all(a$LPA[From == "MVPA" & Delta > 1]$Mean < 0)) # more LPA less MVPA # # } else { # expect_true(all(a$MVPA[From == "LPA" & Delta > 1]$Mean < 0)) # expect_true(all(a$LPA[From == "MVPA" & Delta > 1]$Mean > 0)) # } # # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[7, 3], summary(m$model)$fixed[7, 4])) # == (0 %agele% c(a$TST[From == "WAKE" & Delta == 1]$CI_low, # a$TST[From == "WAKE" & Delta == 1]$CI_high)))) # # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[9, 3], summary(m$model)$fixed[9, 4])) # == (0 %agele% c(a$MVPA[From == "LPA" & Delta == 1]$CI_low, # a$MVPA[From == "LPA" & Delta == 1]$CI_high)))) # # }) # # sbp <- matrix(c( # 1, -1, 1,-1, -1, # 1, 0, -1, 0, 0, # 0, 1, 0, -1, -1, # 0, 0, 0, 1, -1), ncol = 5, byrow = TRUE) # # cilr <- complr(data = mcompd, sbp = sbp, # parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID", total = 1440) # # suppressWarnings(m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + bilr2 + bilr3 + bilr4 + # wilr1 + wilr2 + wilr3 + wilr4 + Female + (1 | ID), # seed = 123)) # # s <- wsub(object = m, basesub = psub, delta = 2, summary = TRUE) # # test_that("wsub's results matches with brm model (TST vs MVPA and LPA vs SB)", { # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[7, 1] > 0))) { # wilr2 = more TST less MVPA # expect_true(all(s$TST[From == "MVPA" & Delta > 1]$Mean > 0)) # more TST less MVPA # expect_true(all(s$MVPA[From == "TST" & Delta > 1]$Mean < 0)) # more MVPA less TST # } else { # expect_true(all(s$TST[From == "MVPA" & Delta > 1]$Mean < 0)) # expect_true(all(s$MVPA[From == "TST" & Delta > 1]$Mean > 0)) # } # # if (isTRUE(suppressWarnings(summary(m$model)$fixed[9, 1] > 0))) { # wilr4 = more LPA less SB # expect_true(all(s$LPA[From == "SB" & Delta > 1]$Mean > 0)) # more LPA less SB # expect_true(all(s$SB[From == "LPA" & Delta > 1]$Mean < 0)) # more SB less LPA # } else { # expect_true(all(s$LPA[From == "SB" & Delta > 1]$Mean < 0)) # expect_true(all(s$SB[From == "LPA" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[7, 3], summary(m$model)$fixed[7, 4])) # == (0 %agele% c(s$TST[From == "MVPA" & Delta == 1]$CI_low, # s$TST[From == "MVPA" & Delta == 1]$CI_high)))) # # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[9, 3], summary(m$model)$fixed[9, 4])) # == (0 %agele% c(s$LPA[From == "SB" & Delta == 1]$CI_low, # s$LPA[From == "SB" & Delta == 1]$CI_high)))) # # }) # # sbp <- matrix(c( # 1, -1, -1, 1, -1, # 1, 0, 0, -1, 0, # 0, -1, 1, 0, -1, # 0, 1, 0, 0, -1), ncol = 5, byrow = TRUE) # # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID", total = 1440) # # suppressWarnings(m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + bilr2 + bilr3 + bilr4 + # wilr1 + wilr2 + wilr3 + wilr4 + Female + (1 | ID), # chain = 1, iter = 500, seed = 123)) # # d <- wsub(object = m, basesub = psub, delta = 2) # # test_that("wsub's results matches with brm model (TST vs LPA and WAKE vs SB)", { # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[7, 1] > 0))) { # bilr2 = more TST less LPA # expect_true(all(d$TST[From == "LPA" & Delta > 1]$Mean > 0)) # more TST less LPA # expect_true(all(d$LPA[From == "TST" & Delta > 1]$Mean < 0)) # more LPA less TST # } else { # expect_true(all(d$TST[From == "LPA" & Delta > 1]$Mean < 0)) # expect_true(all(d$LPA[From == "TST" & Delta > 1]$Mean > 0)) # } # # if (isTRUE(suppressWarnings(summary(m$model)$fixed[9, 1] > 0))) { # bilr4 = more WAKE less SB # expect_true(all(d$WAKE[From == "SB" & Delta > 1]$Mean > 0)) # more WAKE less SB # expect_true(all(d$SB[From == "WAKE" & Delta > 1]$Mean < 0)) # more SB less WAKE # } else { # expect_true(all(d$WAKE[From == "SB" & Delta > 1]$Mean < 0)) # expect_true(all(d$SB[From == "WAKE" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[7, 3], summary(m$model)$fixed[7, 4])) # == (0 %agele% c(d$TST[From == "LPA" & Delta == 1]$CI_low, # d$TST[From == "LPA" & Delta == 1]$CI_high)))) # # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[9, 3], summary(m$model)$fixed[9, 4])) # == (0 %agele% c(d$WAKE[From == "SB" & Delta == 1]$CI_low, # d$WAKE[From == "SB" & Delta == 1]$CI_high)))) # # }) # # sbp <- matrix(c( # 1, -1, -1, -1, 1, # 1, 0, 0, 0, -1, # 0, -1, -1, 1, 0, # 0, 1, -1, 0, 0), ncol = 5, byrow = TRUE) # # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID", total = 1440) # # suppressWarnings(m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + bilr2 + bilr3 + bilr4 + # wilr1 + wilr2 + wilr3 + wilr4 + Female + (1 | ID), # chain = 1, iter = 500, seed = 123)) # # f <- wsub(object = m, basesub = psub, delta = 2) # # test_that("wsub's results matches with brm model (TST vs SB and WAKE vs MVPA)", { # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[7, 1] > 0))) { # bilr2 = more TST less SB # expect_true(all(f$TST[From == "SB" & Delta > 1]$Mean > 0)) # more TST less SB # expect_true(all(f$SB[From == "TST" & Delta > 1]$Mean < 0)) # more SB less TST # } else { # expect_true(all(f$TST[From == "SB" & Delta > 1]$Mean < 0)) # expect_true(all(f$SB[From == "TST" & Delta > 1]$Mean > 0)) # } # # if (isTRUE(suppressWarnings(summary(m$model)$fixed[9, 1] > 0))) { # bilr4 = more WAKE less MVPA # expect_true(all(f$WAKE[From == "MVPA" & Delta > 1]$Mean > 0)) # more WAKE less MVPA # expect_true(all(f$MVPA[From == "WAKE" & Delta > 1]$Mean < 0)) # more MVPA less WAKE # } else { # expect_true(all(f$WAKE[From == "MVPA" & Delta > 1]$Mean < 0)) # expect_true(all(f$MVPA[From == "WAKE" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[7, 3], summary(m$model)$fixed[7, 4])) # == (0 %agele% c(f$TST[From == "SB" & Delta == 1]$CI_low, # f$TST[From == "SB" & Delta == 1]$CI_high)))) # # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[9, 3], summary(m$model)$fixed[9, 4])) # == (0 %agele% c(f$WAKE[From == "MVPA" & Delta == 1]$CI_low, # f$WAKE[From == "MVPA" & Delta == 1]$CI_high)))) # # }) # # sbp <- matrix(c( # -1, -1, 1, -1, 1, # 0, 0, 1, 0, -1, # 1, -1, 0, -1, 0, # 0, 1, 0, -1, 0), ncol = 5, byrow = TRUE) # # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("TST", "WAKE", "MVPA", "LPA", "SB"), idvar = "ID", total = 1440) # # suppressWarnings(m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + bilr2 + bilr3 + bilr4 + # wilr1 + wilr2 + wilr3 + wilr4 + Female + (1 | ID), # seed = 123)) # # g <- wsub(object = m, basesub = psub, delta = 2) # # test_that("wsub's results matches with brm model (MVPA vs SB) and (WAKE vs LPA)", { # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[7, 1] > 0))) { # wilr2 = more MVPA less SB # expect_true(all(g$MVPA[From == "SB" & Delta > 1]$Mean > 0)) # more MVPA less SB # expect_true(all(g$SB[From == "MVPA" & Delta > 1]$Mean < 0)) # more SB less MVPA # } else { # expect_true(all(g$MVPA[From == "SB" & Delta > 1]$Mean < 0)) # expect_true(all(g$SB[From == "MVPA" & Delta > 1]$Mean > 0)) # } # # if (isTRUE(suppressWarnings(summary(m$model)$fixed[9, 1] > 0))) { # wilr4 = more WAKE less LPA # expect_true(all(g$WAKE[From == "LPA" & Delta > 1]$Mean > 0)) # more WAKE less LPA # expect_true(all(g$LPA[From == "WAKE" & Delta > 1]$Mean < 0)) # more LPA less WAKE # } else { # expect_true(all(g$WAKE[From == "LPA" & Delta > 1]$Mean < 0)) # expect_true(all(g$LPA[From == "WAKE" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[7, 3], summary(m$model)$fixed[7, 4])) # == (0 %agele% c(g$MVPA[From == "SB" & Delta == 1]$CI_low, # g$MVPA[From == "SB" & Delta == 1]$CI_high)))) # # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[9, 3], summary(m$model)$fixed[9, 4])) # == (0 %agele% c(g$WAKE[From == "LPA" & Delta == 1]$CI_low, # g$WAKE[From == "LPA" & Delta == 1]$CI_high)))) # # }) # # #--------------------------------------------------------------------------------------------------- # # Test 2-component composition for consistency between brm model and substitution model # ## TST vs WAKE # test_that("wsub's results matches with brm model for 2-component composition (TST vs WAKE)", { # # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("TST", "WAKE"), idvar = "ID", total = 1440) # psub <- build.basesub(c("TST", "WAKE")) # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # a <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(a$TST[From == "WAKE" & Delta > 1]$Mean > 0)) # expect_true(all(a$WAKE[From == "TST" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(a$TST[From == "WAKE" & Delta > 1]$Mean < 0)) # expect_true(all(a$WAKE[From == "TST" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(a$TST[From == "WAKE" & Delta == 1]$CI_low, # a$TST[From == "WAKE" & Delta == 1]$CI_high)))) # # }) # # ## TST vs MVPA # test_that("wsub's results matches with brm model for 2-component composition (TST vs MVPA)", { # # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("TST", "MVPA"), idvar = "ID", total = 1440) # psub <- build.basesub(c("TST", "MVPA")) # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # b <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(b$TST[From == "MVPA" & Delta > 1]$Mean > 0)) # expect_true(all(b$MVPA[From == "TST" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(b$TST[From == "MVPA" & Delta > 1]$Mean < 0)) # expect_true(all(b$MVPA[From == "TST" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(b$TST[From == "MVPA" & Delta == 1]$CI_low, # b$TST[From == "MVPA" & Delta == 1]$CI_high)))) # # }) # # ## TST vs LPA # test_that("wsub's results matches with brm model for 2-component composition (TST vs LPA)", { # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("TST", "LPA"), idvar = "ID", total = 1440) # psub <- build.basesub(c("TST", "LPA")) # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # c <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(c$TST[From == "LPA" & Delta > 1]$Mean > 0)) # expect_true(all(c$LPA[From == "TST" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(c$TST[From == "LPA" & Delta > 1]$Mean < 0)) # expect_true(all(c$LPA[From == "TST" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(c$TST[From == "LPA" & Delta == 1]$CI_low, # c$TST[From == "LPA" & Delta == 1]$CI_high)))) # # }) # # ## TST vs SB # test_that("wsub's results matches with brm model for 2-component composition (TST vs SB)", { # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("TST", "SB"), idvar = "ID", total = 1440) # psub <- build.basesub(c("TST", "SB")) # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # d <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(d$TST[From == "SB" & Delta > 1]$Mean > 0)) # expect_true(all(d$SB[From == "TST" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(d$TST[From == "SB" & Delta > 1]$Mean < 0)) # expect_true(all(d$SB[From == "TST" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(d$TST[From == "SB" & Delta == 1]$CI_low, # d$TST[From == "SB" & Delta == 1]$CI_high)))) # # }) # # ## WAKE vs MVPA # test_that("wsub's results matches with brm model for 2-component composition (WAKE vs MVPA)", { # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("WAKE", "MVPA"), idvar = "ID", total = 1440) # psub <- build.basesub(c("WAKE", "MVPA")) # # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # e <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(e$WAKE[From == "MVPA" & Delta > 1]$Mean > 0)) # expect_true(all(e$MVPA[From == "WAKE" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(e$WAKE[From == "MVPA" & Delta > 1]$Mean < 0)) # expect_true(all(e$MVPA[From == "WAKE" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(e$WAKE[From == "MVPA" & Delta == 1]$CI_low, # e$WAKE[From == "MVPA" & Delta == 1]$CI_high)))) # # }) # # ## WAKE vs LPA # test_that("wsub's results matches with brm model for 2-component composition (WAKE vs LPA)", { # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("WAKE", "LPA"), idvar = "ID", total = 1440) # psub <- build.basesub(c("WAKE", "LPA")) # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # f <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(f$WAKE[From == "LPA" & Delta > 1]$Mean > 0)) # expect_true(all(f$LPA[From == "WAKE" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(f$WAKE[From == "LPA" & Delta > 1]$Mean < 0)) # expect_true(all(f$LPA[From == "WAKE" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(f$WAKE[From == "LPA" & Delta == 1]$CI_low, # f$WAKE[From == "LPA" & Delta == 1]$CI_high)))) # # }) # # ## WAKE vs SB # test_that("wsub's results matches with brm model for 2-component composition (WAKE vs SB)", { # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("WAKE", "SB"), idvar = "ID", total = 1440) # psub <- build.basesub(c("WAKE", "SB")) # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # g <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(g$WAKE[From == "SB" & Delta > 1]$Mean > 0)) # expect_true(all(g$SB[From == "WAKE" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(g$WAKE[From == "SB" & Delta > 1]$Mean < 0)) # expect_true(all(g$SB[From == "WAKE" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(g$WAKE[From == "SB" & Delta == 1]$CI_low, # g$WAKE[From == "SB" & Delta == 1]$CI_high)))) # # }) # # ## MVPA vs LPA # test_that("wsub's results matches with brm model for 2-component composition (MVPA vs LPA)", { # # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("MVPA", "LPA"), idvar = "ID", total = 1440) # psub <- build.basesub(c("MVPA", "LPA")) # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # h <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(h$MVPA[From == "LPA" & Delta > 1]$Mean > 0)) # expect_true(all(h$LPA[From == "MVPA" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(h$MVPA[From == "LPA" & Delta > 1]$Mean < 0)) # expect_true(all(h$LPA[From == "MVPA" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(h$MVPA[From == "LPA" & Delta == 1]$CI_low, # h$MVPA[From == "LPA" & Delta == 1]$CI_high)))) # # }) # # ## MVPA vs SB # test_that("wsub's results matches with brm model for 2-component composition (MVPA vs SB)", { # # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("MVPA", "SB"), idvar = "ID", total = 1440) # psub <- build.basesub(c("MVPA", "SB")) # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # i <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(i$MVPA[From == "SB" & Delta > 1]$Mean > 0)) # expect_true(all(i$SB[From == "MVPA" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(i$MVPA[From == "SB" & Delta > 1]$Mean < 0)) # expect_true(all(i$SB[From == "MVPA" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(i$MVPA[From == "SB" & Delta == 1]$CI_low, # i$MVPA[From == "SB" & Delta == 1]$CI_high)))) # # }) # # ## LPA vs SB # test_that("wsub's results matches with brm model for 2-component composition (LPA vs SB)", { # # sbp <- as.matrix(data.table(1, -1)) # cilr <- complr(data = mcompd[ID %in% 1:10, .SD[1:3], by = ID], sbp = sbp, # parts = c("LPA", "SB"), idvar = "ID", total = 1440) # psub <- build.basesub(c("LPA", "SB")) # suppressWarnings( # m <- brmcoda(complr = cilr, # formula = Stress ~ bilr1 + wilr1 + (1 | ID), # chain = 1, iter = 500, seed = 123, # backend = backend)) # j <- wsub(object = m, basesub = psub, delta = 1:2) # # ## Estimates # if (isTRUE(suppressWarnings(summary(m$model)$fixed[3, 1] > 0))) { # expect_true(all(j$LPA[From == "SB" & Delta > 1]$Mean > 0)) # expect_true(all(j$SB[From == "LPA" & Delta > 1]$Mean < 0)) # } else { # expect_true(all(j$LPA[From == "SB" & Delta > 1]$Mean < 0)) # expect_true(all(j$SB[From == "LPA" & Delta > 1]$Mean > 0)) # } # # # CIs # suppressWarnings(expect_true( # (0 %gele% c(summary(m$model)$fixed[3, 3], summary(m$model)$fixed[3, 4])) # == (0 %agele% c(j$LPA[From == "SB" & Delta == 1]$CI_low, # j$LPA[From == "SB" & Delta == 1]$CI_high)))) # # })