added gee.scale.fix argument to testDynamic() and marge2() -- related…

… to #255

R/marge2.R

@@ -25,7 +25,8 @@
 #' @param sandwich.var (Optional) Should the sandwich variance estimator be used instead of the model-based estimator? Default to FALSE.
 #' @param approx.knot (Optional) Should the set of candidate knots be subsampled in order to speed up computation? This has little effect on the final fit, but can improve computation time somewhat. Defaults to TRUE.
 #' @param n.knot.max (Optional) The maximum number of candidate knots to consider. Uses random sampling (don't worry, a random seed is set internally) to select this number of unique values from the reduced set of all candidate knots. Defaults to 50.
-#' @param glm.backend (Optional) Character specifying which GLM-fitting backend should be used. Must be one of "MASS" or "speedglm". Defaults to "MASS". 
+#' @param glm.backend (Optional) Character specifying which GLM-fitting backend should be used. Must be one of "MASS" or "speedglm". Defaults to "MASS".
+#' @param gee.scale.fix (Optional) Boolean specifying whether the dispersion should be estimated from the data or held fixed at 1 when fitting in GEE mode. Defaults to FALSE.
 #' @param tols_score (Optional) The set tolerance for monitoring the convergence for the difference in score statistics between the parent and candidate model (this is the lack-of-fit criterion used for MARGE). Defaults to 0.00001.
 #' @param minspan (Optional) A set minimum span value. Defaults to NULL.
 #' @param return.basis (Optional) Whether the basis model matrix should be returned as part of the \code{marge} model object. Defaults to FALSE.
@@ -64,7 +65,8 @@ marge2 <- function(X_pred = NULL,
                    sandwich.var = FALSE,
                    approx.knot = TRUE,
                    n.knot.max = 50,
-                   glm.backend = "MASS", 
+                   glm.backend = "MASS",
+                   gee.scale.fix = FALSE,
                    tols_score = 1e-5,
                    minspan = NULL,
                    return.basis = FALSE,
@@ -836,7 +838,7 @@ marge2 <- function(X_pred = NULL,
                               id = id.vec,
                               family = MASS::negative.binomial(theta_hat, link = log),
                               corstr = cor.structure,
-                              scale.fix = FALSE,
+                              scale.fix = gee.scale.fix,
                               sandwich = sandwich.var)
     } else {
       if (glm.backend == "MASS") {
@@ -849,12 +851,12 @@ marge2 <- function(X_pred = NULL,
                                   model = FALSE)
         final_mod <- stripGLM(glm.obj = final_mod)
       } else if (glm.backend == "speedglm") {
-        final_mod <- speedglm::speedglm(model_formula, 
-                                        data = model_df, 
-                                        family = MASS::negative.binomial(theta_hat, link = "log"), 
-                                        trace = FALSE, 
-                                        model = FALSE, 
-                                        y = FALSE, 
+        final_mod <- speedglm::speedglm(model_formula,
+                                        data = model_df,
+                                        family = MASS::negative.binomial(theta_hat, link = "log"),
+                                        trace = FALSE,
+                                        model = FALSE,
+                                        y = FALSE,
                                         fitted = TRUE)
       }
     }

R/testDynamic.R

@@ -24,7 +24,8 @@
 #' @param size.factor.offset (Optional) An offset to be included in the final model fit. Can be generated easily with \code{\link{createCellOffset}}. Defaults to NULL.
 #' @param is.gee Should a GEE framework be used instead of the default GLM? Defaults to FALSE.
 #' @param cor.structure If the GEE framework is used, specifies the desired working correlation structure. Must be one of "ar1", "independence", or "exchangeable". Defaults to "ar1".
-#' @param gee.bias.correction.method Specify which small-sample bias correction to be used on the sandwich variance-covariance matrix prior to test statistic estimation. Options are "kc" and "df". Defaults to NULL, indicating the use of the model-based variance.
+#' @param gee.bias.correction.method (Optional) Specify which small-sample bias correction to be used on the sandwich variance-covariance matrix prior to test statistic estimation. Options are "kc" and "df". Defaults to NULL, indicating the use of the model-based variance.
+#' @param gee.scale.fix (Optional) Boolean specifying whether the dispersion should be estimated from the data or held fixed at 1 when fitting in GEE mode. Defaults to FALSE.
 #' @param is.glmm Should a GLMM framework be used instead of the default GLM? Defaults to FALSE.
 #' @param id.vec If a GEE or GLMM framework is being used, a vector of subject IDs to use as input to \code{\link[geeM]{geem}} or \code{\link[glmmTMB]{glmmTMB}}. Defaults to NULL.
 #' @param glmm.adaptive (Optional) Should the basis functions for the GLMM be chosen adaptively? If not, uses 4 evenly spaced knots. Defaults to TRUE.
@@ -63,6 +64,7 @@ testDynamic <- function(expr.mat = NULL,
                         is.gee = FALSE,
                         cor.structure = "ar1",
                         gee.bias.correction.method = NULL,
+                        gee.scale.fix = FALSE,
                         is.glmm = FALSE,
                         glmm.adaptive = TRUE,
                         id.vec = NULL,
@@ -73,7 +75,7 @@ testDynamic <- function(expr.mat = NULL,
                         random.seed = 312) {
   # check inputs
   if (is.null(expr.mat) || is.null(pt)) { stop("You forgot some inputs to testDynamic().") }
-  
+
   # get raw counts from SingleCellExperiment or Seurat object & transpose to cell x gene dense matrix
   if (is.null(genes)) {
     genes <- rownames(expr.mat)
@@ -144,7 +146,7 @@ testDynamic <- function(expr.mat = NULL,
 
   # build list of objects to prevent from being sent to parallel workers
   necessary_vars <- c("expr.mat", "genes", "pt", "n.potential.basis.fns", "approx.knot", "is.glmm", "gee.bias.correction.method",
-                      "verbose", "n_lineages", "id.vec", "cor.structure", "is.gee", "glmm.adaptive", "size.factor.offset")
+                      "verbose", "n_lineages", "id.vec", "cor.structure", "is.gee", "gee.scale.fix", "glmm.adaptive", "size.factor.offset")
   if (any(ls(envir = .GlobalEnv) %in% necessary_vars)) {
     no_export <- c(ls(envir = .GlobalEnv)[-which(ls(envir = .GlobalEnv) %in% necessary_vars)],
                    ls()[-which(ls() %in% necessary_vars)])
@@ -186,11 +188,12 @@ testDynamic <- function(expr.mat = NULL,
                  Y = expr.mat[lineage_cells, i],
                  Y.offset = size.factor.offset[lineage_cells],
                  is.gee = is.gee,
+                 gee.scale.fix = gee.scale.fix,
                  id.vec = id.vec[lineage_cells],
                  cor.structure = cor.structure,
                  sandwich.var = ifelse(is.null(gee.bias.correction.method), FALSE, TRUE),
                  M = n.potential.basis.fns,
-                 approx.knot = approx.knot, 
+                 approx.knot = approx.knot,
                  return.basis = TRUE)
         }, silent = TRUE)
       } else if (is.glmm) {
@@ -241,7 +244,7 @@ testDynamic <- function(expr.mat = NULL,
                      data = null_mod_df,
                      family = MASS::negative.binomial(theta_hat),
                      corstr = cor.structure,
-                     scale.fix = FALSE,
+                     scale.fix = gee.scale.fix,
                      sandwich = ifelse(is.null(gee.bias.correction.method), FALSE, TRUE))
         }, silent = TRUE)
       } else if (is.glmm) {