From d8d148aca1ef6986808a7e2315d9ed7bf355dda9 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Tue, 26 Aug 2025 11:23:36 -0700 Subject: [PATCH 1/8] updates to add control list --- R/fit_mgx_model.R | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/R/fit_mgx_model.R b/R/fit_mgx_model.R index 5c52d00..df58dfe 100644 --- a/R/fit_mgx_model.R +++ b/R/fit_mgx_model.R @@ -21,6 +21,10 @@ #' there is a small number of clusters. #' @param cluster_corr_coef when replicates are given, estimated value of the within-cluster correlation coefficient. This will only be used when gee estimation in `raoBust::gee_test` fails, and instead #' estimation is performed with a glm. This is set to NULL by default. +#' @param control a list of control parameters. +#' Defaults are \code{list(center = TRUE)}}. +#' Users can override some or all of these. +#' If \code{center = TRUE}, then covariates will be centered before being included in the model. #' #' @importFrom tibble tibble #' @importFrom dplyr bind_cols @@ -42,9 +46,14 @@ fit_mgx_model <- function( wts = NULL, replace_zeros = "minimum", use_jack_se = FALSE, - cluster_corr_coef = NULL + cluster_corr_coef = NULL, + control = list() ) { + # merge user-supplied control arguments with default + defaults <- list(center = TRUE) + control <- modifyList(defaults, control) + # Setting pseudocounts for xx and xstar if (replace_zeros == "minimum") { xx[xx == 0] <- min(xx[xx > 0]) # cheap pseudocount - take minimum @@ -100,6 +109,11 @@ fit_mgx_model <- function( } my_df$wts <- wts + # center covariates if desired + if (control$center) { + + } + ################################################ ## fit the model with Poisson regression ################################################ From c2ddea8a4db424a16ac7e6ed3f650900b842f5c8 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Tue, 26 Aug 2025 11:29:12 -0700 Subject: [PATCH 2/8] add NEWS.md file --- NEWS.md | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 NEWS.md diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..af9f493 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,9 @@ +# enviromtx 1.0.0 + +This is a major release that requires a new release of `raoBust` which performs many of the backend computation, changes the way data are input to `fit_mgx_model` and by default centers covariates. + +## Breaking changes + +* `raoBust` version 1.1.1 is required, which updates the computation of robust score tests that are used in `fit_mgx_model()` +* the inputs to `fit_mgx_model()` are updated so that all vectors need to be included within the `enviro_df` data frame, and cannot be taken directly from the global environment. +* a `control` argument is added for `fit_mgx_model()` with the option `center = TRUE` by default, which will center all covariates in the model. From f8b34074f546f316bb9558d684b6bfff63c97486 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Tue, 26 Aug 2025 13:49:14 -0700 Subject: [PATCH 3/8] start testing new centering --- R/fit_mgx_model.R | 20 ++++- tests/testthat/test-fit_mgx_model_centering.R | 84 +++++++++++++++++++ 2 files changed, 101 insertions(+), 3 deletions(-) create mode 100644 tests/testthat/test-fit_mgx_model_centering.R diff --git a/R/fit_mgx_model.R b/R/fit_mgx_model.R index fcc7bb9..e00b779 100644 --- a/R/fit_mgx_model.R +++ b/R/fit_mgx_model.R @@ -140,14 +140,28 @@ fit_mgx_model <- function( wts <- "wts" } - my_df$offset <- log(my_df[[xx]]) - offset <- "offset" - # center covariates if desired if (control$center) { + # get variable names from the formula + vars <- all.vars(update(formula, . ~ .))[-1] # drop response + + # copy data so we don't overwrite + data_centered <- my_df + + # center each variable in the formula + for (v in vars) { + is_binary <- length(unique(my_df[[v]])) < 3 + if (is.numeric(my_df[[v]]) & !is_binary) { + data_centered[[v]] <- scale(my_df[[v]], center = TRUE, scale = FALSE) + } + } + my_df <- data_centered } + my_df$offset <- log(my_df[[xx]]) + offset <- "offset" + ################################################ ## fit the model with Poisson regression ################################################ diff --git a/tests/testthat/test-fit_mgx_model_centering.R b/tests/testthat/test-fit_mgx_model_centering.R new file mode 100644 index 0000000..58bed20 --- /dev/null +++ b/tests/testthat/test-fit_mgx_model_centering.R @@ -0,0 +1,84 @@ +test_that("with small parameter values, centering does not affect results", { + + n <- 20 + + # no replicates + set.seed(4) + xx1 <- rpois(n, lambda=400) + xstar1 <- rpois(n, lambda=400) + beta0 <- 10 + beta1 <- 2 + beta2 <- 1 + beta3 <- -1 + xx_covariates1 <- rnorm(n) + xx_covariates2 <- rnorm(n) + yy1 <- rpois(n, xx1 * beta0 * (xstar1/xx1)^beta1 * exp(beta2 * xx_covariates1 + beta3 * xx_covariates2)) + output1 <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1, + xx_covariates1, + xx_covariates2), + formula= ~ xx_covariates1 + xx_covariates2, + replace_zeros=1) + output2 <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1, + xx_covariates1, + xx_covariates2), + formula= ~ xx_covariates1 + xx_covariates2, + replace_zeros=1, + control = list(center = FALSE)) + + expect_true(all.equal(output1, output2)) + + # with replicates + set.seed(5) + m <- 3 + n_clust <- 50 + sigma_b <- 0 + n <- n_clust * m + id <- rep(1:n_clust, each = m) + b <- rnorm(n_clust, mean = 0, sd = sigma_b) + + set.seed(3) + xx1 <- rpois(n, lambda=400) + xstar1 <- rpois(n, lambda=400) + beta0 <- 1 + beta1 <- 2 + beta2 <- 1 + beta3 <- -1 + xx_covariates1 <- rnorm(n) + xx_covariates2 <- rnorm(n) + eta <- xx1 * beta0 * (xstar1/xx1)^beta1 * exp(beta2 * xx_covariates1 + beta3 * xx_covariates2 + b[id]) + yy1 <- rnbinom(n, mu=eta, size=eta^2) + + my_df <- tibble(yy1, xx1, predictor = log((xstar1 + 1)/(xx1 + 1)), xx_covariates1, xx_covariates2, id) + + output3 <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1, xx_covariates1, xx_covariates2, id), + formula= ~ xx_covariates1 + xx_covariates2, + replicates="id", + replace_zeros=1) + output4 <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1, xx_covariates1, xx_covariates2, id), + formula= ~ xx_covariates1 + xx_covariates2, + replicates="id", + replace_zeros=1, + control = list(center = FALSE)) + expect_true(all.equal(output3, output4, tolerance = 1e-4)) + +}) + +test_that("with large parameter values, centering improves estimation", { + + skip("Skip longer sim study in automatic testing") + +}) + +test_that("with null parameter values, centering does not affect t1e control", { + + skip("Skip longer sim study in automatic testing") + +}) From eb010212188c606756eeb94ff4624a9ca1be6bdb Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Tue, 26 Aug 2025 16:18:15 -0700 Subject: [PATCH 4/8] keep investigating centering --- R/fit_mgx_model.R | 4 +- tests/testthat/test-fit_mgx_model_centering.R | 62 +++++++++++++++++++ 2 files changed, 64 insertions(+), 2 deletions(-) diff --git a/R/fit_mgx_model.R b/R/fit_mgx_model.R index a0ba40f..b647ebc 100644 --- a/R/fit_mgx_model.R +++ b/R/fit_mgx_model.R @@ -172,7 +172,7 @@ fit_mgx_model <- function( # center covariates if desired if (control$center) { # get variable names from the formula - vars <- all.vars(update(formula, . ~ .))[-1] # drop response + vars <- setdiff(all.vars(update(formula, . ~ .))[-1], "predictor") # copy data so we don't overwrite data_centered <- my_df @@ -181,7 +181,7 @@ fit_mgx_model <- function( for (v in vars) { is_binary <- length(unique(my_df[[v]])) < 3 if (is.numeric(my_df[[v]]) & !is_binary) { - data_centered[[v]] <- scale(my_df[[v]], center = TRUE, scale = FALSE) + data_centered[[v]] <- my_df[[v]] - mean(my_df[[v]], na.rm = TRUE) } } diff --git a/tests/testthat/test-fit_mgx_model_centering.R b/tests/testthat/test-fit_mgx_model_centering.R index 58bed20..b0763e5 100644 --- a/tests/testthat/test-fit_mgx_model_centering.R +++ b/tests/testthat/test-fit_mgx_model_centering.R @@ -75,6 +75,68 @@ test_that("with large parameter values, centering improves estimation", { skip("Skip longer sim study in automatic testing") + n <- 20 + + # large B values + set.seed(4) + xx1 <- rpois(n, lambda=400) + xstar1 <- rpois(n, lambda=400) + beta0 <- 10 + beta1 <- 20 + beta2 <- 30 + beta3 <- -1 + xx_covariates1 <- rnorm(n) + xx_covariates2 <- rnorm(n) + yy1 <- rpois(n, xx1 * beta0 * (xstar1/xx1)^beta1 * exp(beta2 * xx_covariates1 + beta3 * xx_covariates2)) + output1 <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1, + xx_covariates1, + xx_covariates2), + formula= ~ xx_covariates1 + xx_covariates2, + replace_zeros=1) + output2 <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1, + xx_covariates1, + xx_covariates2), + formula= ~ xx_covariates1 + xx_covariates2, + replace_zeros=1, + control = list(center = FALSE)) + + expect_true(all.equal(output1, output2)) + + + # large covariate values + n <- 20 + set.seed(3) + xx1 <- rpois(n, lambda=400) + xstar1 <- rpois(n, lambda=400) + beta0 <- 1.5 + beta1 <- -1 + beta2 <- 1 + beta3 <- 0 + xx_covariates1 <- rnorm(n) + 20 + xx_covariates2 <- rnorm(n) + yy1 <- rpois(n, xx1 * beta0 * (xstar1/xx1)^beta1 * exp(beta2 * xx_covariates1 + beta3 * xx_covariates2)) + output_center <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1, + xx_covariates1, + xx_covariates2), + formula= ~ xx_covariates1 + xx_covariates2, + replace_zeros=1) + output_nocenter <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1, + xx_covariates1, + xx_covariates2), + formula= ~ xx_covariates1 + xx_covariates2, + replace_zeros=1, + control = list(center = FALSE)) + + expect_true(all.equal(output1, output2)) + }) test_that("with null parameter values, centering does not affect t1e control", { From e96daf82b3a5245dfa583240d825edf49509458d Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Fri, 29 Aug 2025 15:55:57 -0700 Subject: [PATCH 5/8] add option to center and test this option, by default do not center --- NEWS.md | 6 +- R/fit_mgx_model.R | 4 +- tests/testthat/test-fit_mgx_model_centering.R | 143 ++++++++++++++++-- 3 files changed, 140 insertions(+), 13 deletions(-) diff --git a/NEWS.md b/NEWS.md index af9f493..b62add5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,9 +1,13 @@ # enviromtx 1.0.0 -This is a major release that requires a new release of `raoBust` which performs many of the backend computation, changes the way data are input to `fit_mgx_model` and by default centers covariates. +This is a major release that requires a new release of `raoBust` which performs many of the backend computation, changes the way data are input to `fit_mgx_model`. It also provides the option of centering covariates. ## Breaking changes * `raoBust` version 1.1.1 is required, which updates the computation of robust score tests that are used in `fit_mgx_model()` * the inputs to `fit_mgx_model()` are updated so that all vectors need to be included within the `enviro_df` data frame, and cannot be taken directly from the global environment. * a `control` argument is added for `fit_mgx_model()` with the option `center = TRUE` by default, which will center all covariates in the model. + +## Minor changes + +* `fit_mgx_model` now includes a `control` argument. One possible inclusion in the `control` argument is `center`. This is by default `FALSE`, but when it is `TRUE` all covariates will be centered before fitting the model. This should have very small if any effects on estimation, but could affect stability. diff --git a/R/fit_mgx_model.R b/R/fit_mgx_model.R index b647ebc..536cfaa 100644 --- a/R/fit_mgx_model.R +++ b/R/fit_mgx_model.R @@ -22,7 +22,7 @@ #' @param cluster_corr_coef When technical replicates are given, the estimated value of the within-cluster correlation coefficient. This will only be used when GEE estimation in `raoBust::gee_test` fails, and #' estimation is performed with a glm. Defaults to NULL by default (no robust score test is returned). #' @param control a list of control parameters. -#' Defaults are \code{list(center = TRUE)}}. +#' Defaults are \code{list(center = FALSE)}}. #' Users can override some or all of these. #' If \code{center = TRUE}, then covariates will be centered before being included in the model. #' @@ -58,7 +58,7 @@ fit_mgx_model <- function( ) { # merge user-supplied control arguments with default - defaults <- list(center = TRUE) + defaults <- list(center = FALSE) control <- modifyList(defaults, control) # first, check that everything is in enviro_df diff --git a/tests/testthat/test-fit_mgx_model_centering.R b/tests/testthat/test-fit_mgx_model_centering.R index b0763e5..4453ae4 100644 --- a/tests/testthat/test-fit_mgx_model_centering.R +++ b/tests/testthat/test-fit_mgx_model_centering.R @@ -59,20 +59,23 @@ test_that("with small parameter values, centering does not affect results", { xx = xx1, xx_covariates1, xx_covariates2, id), formula= ~ xx_covariates1 + xx_covariates2, replicates="id", - replace_zeros=1) + replace_zeros=1, + control = list(center = TRUE)) output4 <- fit_mgx_model(data.frame(yy = yy1, xstar = xstar1, xx = xx1, xx_covariates1, xx_covariates2, id), formula= ~ xx_covariates1 + xx_covariates2, replicates="id", - replace_zeros=1, - control = list(center = FALSE)) + replace_zeros=1) expect_true(all.equal(output3, output4, tolerance = 1e-4)) }) test_that("with large parameter values, centering improves estimation", { + # turns out this does not improve estimation, estimation is the same but less stable in some of these + # cases after centering + skip("Skip longer sim study in automatic testing") n <- 20 @@ -94,15 +97,15 @@ test_that("with large parameter values, centering improves estimation", { xx_covariates1, xx_covariates2), formula= ~ xx_covariates1 + xx_covariates2, - replace_zeros=1) + replace_zeros=1, + control = list(center = TRUE)) output2 <- fit_mgx_model(data.frame(yy = yy1, xstar = xstar1, xx = xx1, xx_covariates1, xx_covariates2), formula= ~ xx_covariates1 + xx_covariates2, - replace_zeros=1, - control = list(center = FALSE)) + replace_zeros=1) expect_true(all.equal(output1, output2)) @@ -125,22 +128,142 @@ test_that("with large parameter values, centering improves estimation", { xx_covariates1, xx_covariates2), formula= ~ xx_covariates1 + xx_covariates2, - replace_zeros=1) + replace_zeros=1, + control = list(center = TRUE)) output_nocenter <- fit_mgx_model(data.frame(yy = yy1, xstar = xstar1, xx = xx1, xx_covariates1, xx_covariates2), formula= ~ xx_covariates1 + xx_covariates2, + replace_zeros=1) + + expect_true(all.equal(output1, output2)) + +}) + +n <- 20 +test_that("reasonably accurate estimates", { + + # accuracy test 1 + set.seed(3) + xx1 <- rpois(n, lambda=400) + xstar1 <- rpois(n, lambda=400) + beta0 <- 100 + beta1 <- 20 + yy1 <- rpois(n, xx1 * beta0 * (xstar1/xx1)^beta1) + + output6 <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1), + replace_zeros=1, + control = list(center = TRUE)) + output6_noncenter <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1), replace_zeros=1, control = list(center = FALSE)) - expect_true(all.equal(output1, output2)) + expect_equal(output6, output6_noncenter) + + # accuracy test 2 + xx1 <- rpois(n, lambda=100) + xstar1 <- rpois(n, lambda=20) + beta0 <- 1 + beta1 <- -4 + yy1 <- rpois(n, xx1 * beta0 * (xstar1/xx1)^beta1) + yy1 + + output7 <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1), + replace_zeros=1, + control = list(center = TRUE)) + output7_nocenter <- fit_mgx_model(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1), + replace_zeros=1, + control = list(center = FALSE)) + + expect_true(all.equal(output7, output7_nocenter)) }) -test_that("with null parameter values, centering does not affect t1e control", { - skip("Skip longer sim study in automatic testing") +test_that("reasonably accurate estimates with covariates", { + + set.seed(4) + xx1 <- rpois(n, lambda=400) + xstar1 <- rpois(n, lambda=400) + beta0 <- 10 + beta1 <- 2 + beta2 <- 1 + beta3 <- -1 + xx_covariates1 <- rnorm(n) + xx_covariates2 <- rnorm(n) + yy1 <- rpois(n, xx1 * beta0 * (xstar1/xx1)^beta1 * exp(beta2 * xx_covariates1 + beta3 * xx_covariates2)) + + yy1 + output8 <- fit_mgx_model(cbind(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1), xx_covariates1, xx_covariates2), + formula= ~ xx_covariates1 + xx_covariates2, + replace_zeros=1, + control = list(center = TRUE)) + output8_nocenter <- fit_mgx_model(cbind(data.frame(yy = yy1, + xstar = xstar1, + xx = xx1), xx_covariates1, xx_covariates2), + formula= ~ xx_covariates1 + xx_covariates2, + replace_zeros=1, + control = list(center = FALSE)) + + expect_equal(output8, output8_nocenter) + +}) + +test_that("reasonably accurate estimates with covariates and correlation", { + + set.seed(5) + m <- 3 + n_clust <- 50 + sigma_b <- 0 + n <- n_clust * m + #### generate observations to be cluster correlated via random effect + id <- rep(1:n_clust, each = m) + b <- rnorm(n_clust, mean = 0, sd = sigma_b) + + set.seed(3) + xx1 <- rpois(n, lambda=400) + xstar1 <- rpois(n, lambda=400) + beta0 <- 1 + beta1 <- 2 + beta2 <- 5 + beta3 <- -5 + xx_covariates1 <- rnorm(n) + xx_covariates2 <- rnorm(n) + eta <- xx1 * beta0 * (xstar1/xx1)^beta1 * exp(beta2 * xx_covariates1 + beta3 * xx_covariates2 + b[id]) + # yy1 <- rpois(n, lambda=eta) + yy1 <- rnbinom(n, mu=eta, size=eta^2) + + my_df <- tibble(yy1, xx1, xstar1, xx_covariates1, xx_covariates2, id) + + output9 <- fit_mgx_model(yy = "yy1", + xstar = "xstar1", + xx = "xx1", + formula= ~ xx_covariates1 + xx_covariates2, + enviro_df=my_df, + replicates="id", + replace_zeros=1, + control = list(center = TRUE)) + output9_nocenter <- fit_mgx_model(yy = "yy1", + xstar = "xstar1", + xx = "xx1", + formula= ~ xx_covariates1 + xx_covariates2, + enviro_df=my_df, + replicates="id", + replace_zeros=1, + control = list(center = FALSE)) + + expect_equal(output9, output9_nocenter, tolerance = 1e-4) }) From bb492025a3ce2e646a40fafcab64d220524229b5 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Fri, 29 Aug 2025 16:06:09 -0700 Subject: [PATCH 6/8] update documentation --- R/fit_mgx_model.R | 4 ++-- man/fit_mgx_model.Rd | 8 +++++++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/R/fit_mgx_model.R b/R/fit_mgx_model.R index 536cfaa..e083daa 100644 --- a/R/fit_mgx_model.R +++ b/R/fit_mgx_model.R @@ -22,7 +22,7 @@ #' @param cluster_corr_coef When technical replicates are given, the estimated value of the within-cluster correlation coefficient. This will only be used when GEE estimation in `raoBust::gee_test` fails, and #' estimation is performed with a glm. Defaults to NULL by default (no robust score test is returned). #' @param control a list of control parameters. -#' Defaults are \code{list(center = FALSE)}}. +#' Defaults are \code{list(center = FALSE)}. #' Users can override some or all of these. #' If \code{center = TRUE}, then covariates will be centered before being included in the model. #' @@ -59,7 +59,7 @@ fit_mgx_model <- function( # merge user-supplied control arguments with default defaults <- list(center = FALSE) - control <- modifyList(defaults, control) + control <- utils::modifyList(defaults, control) # first, check that everything is in enviro_df if (!(is.character(yy) & diff --git a/man/fit_mgx_model.Rd b/man/fit_mgx_model.Rd index cd436b1..f67e735 100644 --- a/man/fit_mgx_model.Rd +++ b/man/fit_mgx_model.Rd @@ -14,7 +14,8 @@ fit_mgx_model( wts = NULL, replace_zeros = "minimum", use_jack_se = FALSE, - cluster_corr_coef = NULL + cluster_corr_coef = NULL, + control = list() ) } \arguments{ @@ -39,6 +40,11 @@ there are a small number of observations. Defaults to FALSE.} \item{cluster_corr_coef}{When technical replicates are given, the estimated value of the within-cluster correlation coefficient. This will only be used when GEE estimation in \code{raoBust::gee_test} fails, and estimation is performed with a glm. Defaults to NULL by default (no robust score test is returned).} + +\item{control}{a list of control parameters. +Defaults are \code{list(center = FALSE)}. +Users can override some or all of these. +If \code{center = TRUE}, then covariates will be centered before being included in the model.} } \description{ We would like to estimate the \eqn{\beta} parameters in the model From ece92734acf463b91d250bd38ca56717df2b27c9 Mon Sep 17 00:00:00 2001 From: amy Date: Thu, 4 Sep 2025 18:11:31 -0700 Subject: [PATCH 7/8] minor edits to NEWS --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index b62add5..2438cd3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,4 +10,4 @@ This is a major release that requires a new release of `raoBust` which performs ## Minor changes -* `fit_mgx_model` now includes a `control` argument. One possible inclusion in the `control` argument is `center`. This is by default `FALSE`, but when it is `TRUE` all covariates will be centered before fitting the model. This should have very small if any effects on estimation, but could affect stability. +* `fit_mgx_model` now includes a `control` argument. One possible inclusion in the `control` argument is `center`. This is by default `FALSE`. When `TRUE` all covariates will be centered before fitting the model. This should have negligible effects on estimation, but may improve stability in rare cases. From e3bab2f4c10f37277ba2bbaefe3700f123adc642 Mon Sep 17 00:00:00 2001 From: Sarah Teichman Date: Fri, 5 Sep 2025 09:30:03 -0700 Subject: [PATCH 8/8] fix error from missing data test, update version and news to reflect two different updates --- DESCRIPTION | 2 +- NEWS.md | 14 +++++++++----- R/fit_mgx_model.R | 7 +++++++ tests/testthat/test-missing-data.R | 7 +++++-- 4 files changed, 22 insertions(+), 8 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 514d0c6..0513d84 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: enviromtx Title: Predicting Species-Species Interactions from Environmental Metatranscriptomics Data -Version: 1.0.1 +Version: 1.1.0 Authors@R: c(person("Amy D", "Willis", email = "adwillis@uw.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-2802-4317")), person("Sarah", "Teichman", role = "aut"), diff --git a/NEWS.md b/NEWS.md index 2438cd3..da2d63a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,13 +1,17 @@ +# enviromtx 1.1.0 + +This is a minor release that provides the option of centering covariates. + +## Minor changes + +* `fit_mgx_model` now includes a `control` argument. One possible inclusion in the `control` argument is `center`. This is by default `FALSE`. When `TRUE` all covariates will be centered before fitting the model. This should have negligible effects on estimation, but may improve stability in rare cases. + # enviromtx 1.0.0 -This is a major release that requires a new release of `raoBust` which performs many of the backend computation, changes the way data are input to `fit_mgx_model`. It also provides the option of centering covariates. +This is a major release that requires a new release of `raoBust` which performs many of the backend computation, changes the way data are input to `fit_mgx_model`. ## Breaking changes * `raoBust` version 1.1.1 is required, which updates the computation of robust score tests that are used in `fit_mgx_model()` * the inputs to `fit_mgx_model()` are updated so that all vectors need to be included within the `enviro_df` data frame, and cannot be taken directly from the global environment. -* a `control` argument is added for `fit_mgx_model()` with the option `center = TRUE` by default, which will center all covariates in the model. - -## Minor changes -* `fit_mgx_model` now includes a `control` argument. One possible inclusion in the `control` argument is `center`. This is by default `FALSE`. When `TRUE` all covariates will be centered before fitting the model. This should have negligible effects on estimation, but may improve stability in rare cases. diff --git a/R/fit_mgx_model.R b/R/fit_mgx_model.R index 8ba28a7..1c95e44 100644 --- a/R/fit_mgx_model.R +++ b/R/fit_mgx_model.R @@ -205,6 +205,13 @@ fit_mgx_model <- function( if (n_incomplete > 0) { message(paste0(n_incomplete, " of your observations contain missing values for either `yy`, `xx`, `xstar`, or covariates within `formula`. These observations will be dropped from the analysis.")) } + # if after this cleaning there is still an NA, it cannot be wt, id, xx, xstar, predictor, of any + # covariate in the formula. Therefore it must be an irrelevant covariate, and we don't want raoBust + # to drop that sample, so we will remove this covariate + if (sum(is.na(my_df)) > 0) { + var_to_rm <- which(colSums(is.na(my_df)) > 0) + my_df <- my_df[, -var_to_rm] + } ################################################ ## fit the model with Poisson regression diff --git a/tests/testthat/test-missing-data.R b/tests/testthat/test-missing-data.R index debca76..ebcbb9d 100644 --- a/tests/testthat/test-missing-data.R +++ b/tests/testthat/test-missing-data.R @@ -82,7 +82,10 @@ test_that("when a value of yy, xx, xstar, or a relevant covariate is missing, th formula = ~ cov) expect_true(all.equal(glm_res, glm_res_oth)) - expect_silent(gee_res <- fit_mgx_model(my_df_irr_cov, formula = ~ cov, - replicates = "id")) + gee_res <- fit_mgx_model(my_df_irr_cov, formula = ~ cov, + replicates = "id") + gee_res_oth <- fit_mgx_model(subset(my_df_irr_cov, select = -irr_cov), + formula = ~ cov, replicates = "id") + expect_true(all.equal(gee_res, gee_res_oth)) })