tests/testthat/test-wsub.R
In florale/multilevelcoda: Estimate Bayesian Multilevel Models for Compositional Data
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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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))))
})
florale/multilevelcoda documentation built on May 16, 2024, 4:25 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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 <- 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))))
})
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.