Skip to content

Add more tests for formula w/ interactions, and for data w/ NA #79

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Mar 13, 2023
Merged

Add more tests for formula w/ interactions, and for data w/ NA #79

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
86943e8
add more tests re: interactions, including s(x) + s(y) + ti(x,y), and…
Mar 13, 2023
48c941e
add more tests on data with NA
Mar 13, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/pkgdown.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ articles:
installations: installations.html
voxel-wise_data: voxel-wise_data.html
walkthrough: walkthrough.html
last_built: 2023-03-10T22:37Z
last_built: 2023-03-13T20:36Z
urls:
reference: https://pennlinc.github.io/ModelArray/reference
article: https://pennlinc.github.io/ModelArray/articles
Expand Down
69 changes: 63 additions & 6 deletions tests/testthat/helper-expected_values.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -213,6 +213,12 @@ helper_generate_expect_lm <- function(fn.phenotypes,
data <- phenotypes
data$FD <- fd.simu

# phenotypes with NA:
phenotypes_wNA <- phenotypes
phenotypes_wNA$age[1] <- NA
data_wNA <- phenotypes_wNA
data_wNA$FD <- fd.simu

# different scenario:
thename <- "age"
formula <- FD ~ age
Expand All @@ -234,12 +240,11 @@ helper_generate_expect_lm <- function(fn.phenotypes,
dfout <- calcu_stat_lm(formula, data, idx.fixel.lm)
expected.results[[thename]] <- dfout

# phenotypes with NA:
phenotypes_wNA <- phenotypes
phenotypes_wNA$age[1] <- NA
data_wNA <- phenotypes_wNA
data_wNA$FD <- fd.simu

thename <- "age_phenotypeswNA"
formula <- FD ~ age
dfout <- calcu_stat_lm(formula, data_wNA, idx.fixel.lm)
expected.results[[thename]] <- dfout

thename <- "age_phenotypeswNA_na.action-na.omit"
formula <- FD ~ age
dfout <- calcu_stat_lm(formula, data_wNA, idx.fixel.lm,
Expand Down Expand Up @@ -298,6 +303,12 @@ helper_generate_expect_gam <- function(fn.phenotypes,

data <- phenotypes
data$FD <- fd.simu

# phenotypes with NA:
phenotypes_wNA <- phenotypes
phenotypes_wNA$age[1] <- NA
data_wNA <- phenotypes_wNA
data_wNA$FD <- fd.simu

# different scenario:
thename <- "s-age_sex"
Expand Down Expand Up @@ -339,6 +350,17 @@ helper_generate_expect_gam <- function(fn.phenotypes,
formula <- FD ~ s(age) + sex
dfout <- calcu_stat_gam(formula, data, idx.fixel.gam, method="REML") # MAKE SURE THIS IS IN! TODO
expected.results[[thename]] <- dfout

# with NA:
thename <- "s-age-wNA"
formula <- FD ~ s(age) + sex
dfout <- calcu_stat_gam(formula, data_wNA, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "s-age-wNA_na.action-na.omit"
formula <- FD ~ s(age) + sex
dfout <- calcu_stat_gam(formula, data_wNA, idx.fixel.gam, na.action = "na.omit")
expected.results[[thename]] <- dfout

thename <- "factorB_s-age_s-factorA"
formula <- FD ~ factorB + s(age) + s(factorA)
Expand Down Expand Up @@ -389,6 +411,11 @@ helper_generate_expect_gam <- function(fn.phenotypes,
formula <- FD ~ oSex + s(age,k=4, fx=TRUE) + s(age, by=oSex, fx=TRUE) + factorB
dfout <- calcu_stat_gam(formula, data, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "oSex_s-age-wNA-k-4-fx-T_s-age-byoSex-fx-T_factorB"
formula <- FD ~ oSex + s(age,k=4, fx=TRUE) + s(age, by=oSex, fx=TRUE) + factorB
dfout <- calcu_stat_gam(formula, data_wNA, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "oMultiLevels_s-age-k-4-fx-T_s-age-byoMultiLevels-fx-T_factorB"
formula <- FD ~ oMultiLevels + s(age,k=4, fx=TRUE) + s(age, by=oMultiLevels, fx=TRUE) + factorB
Expand All @@ -414,12 +441,42 @@ helper_generate_expect_gam <- function(fn.phenotypes,
formula <- FD ~ ti(age, fx=TRUE) + ti(factorB, fx=TRUE) + ti(age, factorB, fx=TRUE) + factorA
dfout <- calcu_stat_gam(formula, data, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "ti-age-wNA-fx-T_ti-factorB-fx-T_ti-age-factorB-fx-T_factorA"
formula <- FD ~ ti(age, fx=TRUE) + ti(factorB, fx=TRUE) + ti(age, factorB, fx=TRUE) + factorA
dfout <- calcu_stat_gam(formula, data_wNA, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "ti-age-fx-T_ti-factorB-fx-T_factorA"
formula <- FD ~ ti(age, fx=TRUE) + ti(factorB, fx=TRUE)+ factorA
dfout <- calcu_stat_gam(formula, data, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "ti-age-fx-T_ti-factorB-fx-F_ti-age-factorB-fx-T_factorA"
formula <- FD ~ ti(age, fx=TRUE) + ti(factorB, fx=FALSE) + ti(age, factorB, fx=TRUE) + factorA
dfout <- calcu_stat_gam(formula, data, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "ti-age-fx-T_ti-factorB-fx-F_factorA"
formula <- FD ~ ti(age, fx=TRUE) + ti(factorB, fx=FALSE) + factorA
dfout <- calcu_stat_gam(formula, data, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "s-age-fx-T_s-factorB-fx-T_ti-age-factorB-fx-T_factorA"
formula <- FD ~ s(age, fx=TRUE) + s(factorB, fx=TRUE) + ti(age, factorB, fx=TRUE) + factorA
dfout <- calcu_stat_gam(formula, data, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "s-age-wNA-fx-T_s-factorB-fx-T_ti-age-factorB-fx-T_factorA"
formula <- FD ~ s(age, fx=TRUE) + s(factorB, fx=TRUE) + ti(age, factorB, fx=TRUE) + factorA
dfout <- calcu_stat_gam(formula, data_wNA, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "s-age-fx-T_s-factorB-fx-T_factorA"
formula <- FD ~ s(age, fx=TRUE) + s(factorB, fx=TRUE)+ factorA
dfout <- calcu_stat_gam(formula, data, idx.fixel.gam)
expected.results[[thename]] <- dfout

thename <- "intercept" # not to replicate the process here.... otherwise dfout's class is "tbl_df"...
formula <- FD ~ 1
dfout <- calcu_stat_gam(formula, data, idx.fixel.gam)
Expand Down
83 changes: 83 additions & 0 deletions tests/testthat/test-ModelArray_gam.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -278,6 +278,49 @@ test_that("test that ModelArray.gam() works as expected", {
%>% isTRUE())

## different settings in mgcv::gam()'s additional arguments: test if the arguments have been passed into analyseOneElement.gam()
# handling inputs with NA:
phenotypes_wNA <- phenotypes
phenotypes_wNA$age[1] <- NA

# using default methods in `gam()` to handle NA:
mygam_agewNA <- ModelArray.gam(FD ~ s(age) + sex, data = modelarray, phenotypes = phenotypes_wNA, scalar = scalar_name, element.subset = element.subset,
n_cores = 2, pbar = FALSE)
compare_expected_results(mygam_agewNA, expected.results[["s-age-wNA"]])
expect_false(dplyr::all_equal(mygam_default %>% dplyr::select("s_age.statistic"),
mygam_agewNA %>% dplyr::select("s_age.statistic"))
%>% isTRUE())
# using default methods in `gam()` to handle NA + in a formula with interaction terms:
# WARNING: should not request effect size for `age` without removing NA observations from it!
# s(x, by=orderedFactor):
formula <- FD ~ oSex + s(age,k=4, fx=TRUE) + s(age, by=oSex, fx=TRUE) + factorB # ordered factor
mygam_agewNA_sby <- ModelArray.gam(formula = formula, data = modelarray, phenotypes = phenotypes_wNA, scalar = scalar_name, element.subset = element.subset,
n_cores = 2, pbar = FALSE)
compare_expected_results(mygam_agewNA_sby, expected.results[["oSex_s-age-wNA-k-4-fx-T_s-age-byoSex-fx-T_factorB"]])
# ti(x) + ti(z) + ti(x,z), where x has NA:
formula <- FD ~ ti(age, fx=TRUE) + ti(factorB, fx=TRUE) + ti(age, factorB, fx=TRUE) + factorA
mygam_agewNA_tiInteract <- ModelArray.gam(formula = formula, data = modelarray, phenotypes = phenotypes_wNA, scalar = scalar_name, element.subset = element.subset,
n_cores = 2, pbar = FALSE)
compare_expected_results(mygam_agewNA_tiInteract, expected.results[["ti-age-wNA-fx-T_ti-factorB-fx-T_ti-age-factorB-fx-T_factorA"]])
# s(x) + s(z) + ti(x,z), where x has NA:
formula <- FD ~ s(age, fx=TRUE) + s(factorB, fx=TRUE) + ti(age, factorB, fx=TRUE) + factorA
mygam_agewNA_s_tiInteract <- ModelArray.gam(formula = formula, data = modelarray, phenotypes = phenotypes_wNA, scalar = scalar_name, element.subset = element.subset,
n_cores = 2, pbar = FALSE)
compare_expected_results(mygam_agewNA_s_tiInteract, expected.results[["s-age-wNA-fx-T_s-factorB-fx-T_ti-age-factorB-fx-T_factorA"]])

# using "na.omit" to handle NA:
mygam_agewNA_na.omit <- ModelArray.gam(FD ~ s(age) + sex, data = modelarray, phenotypes = phenotypes_wNA, scalar = scalar_name, element.subset = element.subset,
na.action = "na.omit",
n_cores = 2, pbar = FALSE)
compare_expected_results(mygam_agewNA_na.omit, expected.results[["s-age-wNA_na.action-na.omit"]])
# default option to handle NA should be the same as using `na.omit`:
expect_true(dplyr::all_equal(mygam_agewNA %>% dplyr::select("s_age.statistic"),
mygam_agewNA_na.omit %>% dplyr::select("s_age.statistic"))
%>% isTRUE())
# using "na.fail" to handle NA:
expect_error(ModelArray.gam(FD ~ s(age) + sex, data = modelarray, phenotypes = phenotypes_wNA, scalar = scalar_name, element.subset = element.subset,
na.action = "na.fail",
n_cores = 2, pbar = FALSE))

# method:
mygam_methodREML <- ModelArray.gam(FD ~ s(age) + sex, data = modelarray, phenotypes = phenotypes, scalar = scalar_name, element.subset = element.subset,
method = "REML", # default = "GCV.Cp"
Expand Down Expand Up @@ -590,6 +633,7 @@ test_that("test that ModelArray.gam() works as expected", {
%in% colnames(mygam_sby_multiLevels))

## ti(x,z):
# ti(x) + ti(z) + ti(x,z):
formula <- FD ~ ti(age, fx=TRUE) + ti(factorB, fx=TRUE) + ti(age, factorB, fx=TRUE) + factorA
mygam_tiInteract <- ModelArray.gam(formula = formula, data = modelarray, phenotypes = phenotypes, scalar = scalar_name, element.subset = element.subset,
changed.rsq.term.index = c(3), var.model = c("dev.expl","adj.r.squared"),
Expand All @@ -608,6 +652,45 @@ test_that("test that ModelArray.gam() works as expected", {
expect_equal(mygam_tiInteract$model.adj.r.squared -mygam_tiInteract_red1$model.adj.r.squared,
mygam_tiInteract$ti_age_factorB.delta.adj.rsq)

# ti(x) + ti(z, fx=F) + ti(x,z):
formula <- FD ~ ti(age, fx=TRUE) + ti(factorB, fx=FALSE) + ti(age, factorB, fx=TRUE) + factorA
mygam_tifxF_tiInteract <- ModelArray.gam(formula = formula, data = modelarray, phenotypes = phenotypes, scalar = scalar_name, element.subset = element.subset,
changed.rsq.term.index = c(3), var.model = c("dev.expl","adj.r.squared"),
n_cores = 2, pbar = FALSE)
compare_expected_results(mygam_tifxF_tiInteract, expected.results[["ti-age-fx-T_ti-factorB-fx-F_ti-age-factorB-fx-T_factorA"]])
expect_true("ti_age.statistic" %in% colnames(mygam_tiInteract))
expect_true("ti_age_factorB.p.value" %in% colnames(mygam_tiInteract))
expect_true("ti_age_factorB.delta.adj.rsq" %in% colnames(mygam_tiInteract))
expect_true("ti_age_factorB.partial.rsq" %in% colnames(mygam_tiInteract))

# not recommend to request effect size for this formula
red.formula <- FD ~ ti(age, fx=TRUE) + ti(factorB, fx=FALSE)+ factorA
mygam_tifxF_tiInteract_red1 <- ModelArray.gam(formula = red.formula, data = modelarray, phenotypes = phenotypes, scalar = scalar_name, element.subset = element.subset,
var.model = c("adj.r.squared"),
n_cores = 2, pbar = FALSE)
compare_expected_results(mygam_tifxF_tiInteract_red1, expected.results[["ti-age-fx-T_ti-factorB-fx-F_factorA"]])
expect_equal(mygam_tifxF_tiInteract$model.adj.r.squared -mygam_tifxF_tiInteract_red1$model.adj.r.squared,
mygam_tifxF_tiInteract$ti_age_factorB.delta.adj.rsq)

# s(x) + s(z) + ti(x,z):
formula <- FD ~ s(age, fx=TRUE) + s(factorB, fx=TRUE) + ti(age, factorB, fx=TRUE) + factorA
mygam_s_tiInteract <- ModelArray.gam(formula = formula, data = modelarray, phenotypes = phenotypes, scalar = scalar_name, element.subset = element.subset,
changed.rsq.term.index = c(3), var.model = c("dev.expl","adj.r.squared"),
n_cores = 2, pbar = FALSE)
compare_expected_results(mygam_s_tiInteract, expected.results[["s-age-fx-T_s-factorB-fx-T_ti-age-factorB-fx-T_factorA"]])
expect_true("s_age.statistic" %in% colnames(mygam_s_tiInteract))
expect_true("ti_age_factorB.p.value" %in% colnames(mygam_s_tiInteract))
expect_true("ti_age_factorB.delta.adj.rsq" %in% colnames(mygam_s_tiInteract))
expect_true("ti_age_factorB.partial.rsq" %in% colnames(mygam_s_tiInteract))

red.formula <- FD ~ s(age, fx=TRUE) + s(factorB, fx=TRUE)+ factorA
mygam_s_tiInteract_red1 <- ModelArray.gam(formula = red.formula, data = modelarray, phenotypes = phenotypes, scalar = scalar_name, element.subset = element.subset,
var.model = c("adj.r.squared"),
n_cores = 2, pbar = FALSE)
compare_expected_results(mygam_s_tiInteract_red1, expected.results[["s-age-fx-T_s-factorB-fx-T_factorA"]])
expect_equal(mygam_s_tiInteract$model.adj.r.squared -mygam_s_tiInteract_red1$model.adj.r.squared,
mygam_s_tiInteract$ti_age_factorB.delta.adj.rsq)


# ## factorized, but not ordered: - NOT RECOMMEND
# formula <- FD ~ sexFactor + s(age) + s(age, by=sexFactor, fx=TRUE)
Expand Down
23 changes: 20 additions & 3 deletions tests/testthat/test-ModelArray_lm.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -277,9 +277,22 @@ test_that("ModelArray.lm() works as expected", {
expect_false(all(is.na(mylm_age_factorB$factorB.statistic)))

## Different optional arguments of lm: in order to test that the additional arguments have really been passed into the lm: #####
# test "na.action" with inputs with NA
# hanlding inputs with NA:
phenotypes_wNA <- phenotypes
phenotypes_wNA$age[1] <- NA

# test default `lm()` method for handling NA:
mylm_phenotypes_naActionDefault <- ModelArray.lm(FD ~ age, data = modelarray, phenotypes = phenotypes_wNA,
scalar = scalar_name, element.subset = element.subset, n_cores = 2, pbar=FALSE)
compare_expected_results(mylm_phenotypes_naActionDefault, expected.results[["age_phenotypeswNA"]])

# there should be differences in results, after changing one value to NA:
expect_false(all.equal(mylm %>% dplyr::select(age.estimate),
mylm_phenotypes_naActionDefault %>% dplyr::select(age.estimate))
%>% isTRUE())

# test different (not-default) `na.action` with inputs with NA:
# error if `na.action = "na.fail"`:
expect_error(ModelArray.lm(FD ~ age, data = modelarray, phenotypes = phenotypes_wNA, scalar = scalar_name, element.subset = element.subset,
var.terms = var.terms, var.model = var.model, n_cores = 1, pbar=FALSE,
na.action="na.fail")) # expect error of "missing values in object". If na.action was not passed into lm, there will not be error
Expand All @@ -295,7 +308,7 @@ test_that("ModelArray.lm() works as expected", {
var.terms = var.terms, var.model = var.model, n_cores = 2, pbar=TRUE,
na.action="na.fail")) # after updating ModelArray.lm with one row, specifying column names to keep, expect error (instead of warning for each core with error, expect_warning)

#mylm_phenotypes_naActionDefault <- ModelArray.lm(FD ~ age, data = modelarray, phenotypes = phenotypes_wNA, scalar = scalar_name, element.subset = element.subset, n_cores = 2, pbar=FALSE)
# okay if `na.action = "na.omit"`:
mylm_phenotypes_naActionOmit <- ModelArray.lm(FD ~ age, data = modelarray, phenotypes = phenotypes_wNA, scalar = scalar_name, element.subset = element.subset,
var.terms = var.terms, var.model = var.model, n_cores = 2, pbar=FALSE,
na.action="na.omit")
Expand All @@ -304,8 +317,12 @@ test_that("ModelArray.lm() works as expected", {
expect_false(all.equal(mylm %>% dplyr::select(age.estimate),
mylm_phenotypes_naActionOmit %>% dplyr::select(age.estimate))
%>% isTRUE())

# default option to handle NA should be the same as using `na.omit`:
expect_true(all.equal(mylm_phenotypes_naActionDefault %>% dplyr::select(age.estimate),
mylm_phenotypes_naActionOmit %>% dplyr::select(age.estimate))
%>% isTRUE())


# check if "weights" have been successfully passed into lm:
mylm_weights1 <- ModelArray.lm(FD ~ age, data = modelarray, phenotypes = phenotypes, scalar = scalar_name, element.subset = element.subset,
var.terms = var.terms, var.model = var.model,
Expand Down