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
new file mode 100644
index 0000000..da2d63a
--- /dev/null
+++ b/NEWS.md
@@ -0,0 +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`. 
+
+## 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.
+
diff --git a/R/fit_mgx_model.R b/R/fit_mgx_model.R
index 7280289..1c95e44 100644
--- a/R/fit_mgx_model.R
+++ b/R/fit_mgx_model.R
@@ -21,6 +21,10 @@
 #' there are a small number of observations. Defaults to FALSE.
 #' @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)}.
+#'    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
@@ -49,9 +53,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 = FALSE)
+  control <- utils::modifyList(defaults, control)
+
   # first, check that everything is in enviro_df
   if (!(is.character(yy) &
         is.character(xstar) &
@@ -166,6 +175,25 @@ fit_mgx_model <- function(
     wts <- "wts"
   }
 
+  # center covariates if desired
+  if (control$center) {
+    # get variable names from the formula
+    vars <- setdiff(all.vars(update(formula, . ~ .))[-1], "predictor")
+
+    # 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]] <- my_df[[v]] - mean(my_df[[v]], na.rm = TRUE)
+      }
+    }
+
+    my_df <- data_centered
+  }
+
   my_df$offset <- log(my_df[[xx]])
   offset <- "offset"
 
@@ -177,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/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
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..4453ae4
--- /dev/null
+++ b/tests/testthat/test-fit_mgx_model_centering.R
@@ -0,0 +1,269 @@
+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,
+                           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)
+  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
+
+  # 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,
+                           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)
+
+  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,
+                                 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_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("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)
+
+})
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))
 })