Merge pull request #259 from jr-leary7/dev

changed scale.fix to TRUE throughout -- related to #241

Author: jr-leary7

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: scLANE
 Type: Package
 Title: Model Gene Expression Dynamics with Spline-Based NB GLMs, GEEs, & GLMMs
-Version: 0.8.5
+Version: 0.8.6
 Authors@R: c(person(given = c("Jack", "R."), family = "Leary", email = "j.leary@ufl.edu", role = c("aut", "cre"), comment = c(ORCID = "0009-0004-8821-3269")), 
              person(given = "Rhonda", family = "Bacher", email = "rbacher@ufl.edu", role = c("ctb", "fnd"), comment = c(ORCID = "0000-0001-5787-476X")))
 Description: Our scLANE model uses truncated power basis spline models to build flexible, interpretable models of single cell gene expression over pseudotime or latent time. 

NEWS.md

@@ -1,3 +1,7 @@
+# Changes in v0.8.6
+
++ Changed GEE fitting to use `scale.fix = TRUE` throughout the package, as it appears to be faster and more statistically efficient based on simulated data benchmarking. 
+
 # Changes in v0.8.5
 
 + Minor bug fixes. 

R/backward_sel_WIC.R

@@ -37,6 +37,7 @@ backward_sel_WIC <- function(Y = NULL,
                       id = id.vec,
                       corstr = cor.structure,
                       family = MASS::negative.binomial(theta.hat, link = "log"),
+                      scale.fix = TRUE,
                       sandwich = sandwich.var)
     wald_stat <- unname(summary(fit)$wald.test[-1])^2
   } else {

R/marge2.R

@@ -26,7 +26,6 @@
 #' @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 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.
@@ -66,7 +65,6 @@ marge2 <- function(X_pred = NULL,
                    approx.knot = TRUE,
                    n.knot.max = 50,
                    glm.backend = "MASS",
-                   gee.scale.fix = FALSE,
                    tols_score = 1e-5,
                    minspan = NULL,
                    return.basis = FALSE,
@@ -838,7 +836,7 @@ marge2 <- function(X_pred = NULL,
                               id = id.vec,
                               family = MASS::negative.binomial(theta_hat, link = log),
                               corstr = cor.structure,
-                              scale.fix = gee.scale.fix,
+                              scale.fix = TRUE,
                               sandwich = sandwich.var)
     } else {
       if (glm.backend == "MASS") {

R/score_fun_gee.R

@@ -99,7 +99,7 @@ score_fun_gee <- function(Y = NULL,
     Sigma <- Sigma22 - temp_prod_1 - temp_prod_2 + temp_prod_3
     Sigma_inv <- try({ eigenMapMatrixInvert(Sigma, n_cores = 1L) }, silent = TRUE)
     if (inherits(Sigma_inv, "try-error")) {
-      Sigma_inv <- MASS::ginv(Sigma)
+      Sigma_inv <- eigenMapPseudoInverse(Sigma, n_cores = 1L)
     }
     temp_prod <- eigenMapMatMult(A = t(B.est),
                                  B = Sigma_inv,

R/stat_out_score_gee_null.R

@@ -36,8 +36,8 @@ stat_out_score_gee_null <- function(Y = NULL,
                      data = NULL,
                      corstr = cor.structure,
                      family = MASS::negative.binomial(theta.hat),
-                     scale.fix = FALSE,
-                     sandwich = FALSE, 
+                     scale.fix = TRUE,
+                     sandwich = FALSE,
                      maxit = 10)
   alpha_est <- ests$alpha
   sigma_est <- ests$phi

R/testDynamic.R

@@ -25,7 +25,6 @@
 #' @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 (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.
@@ -64,7 +63,6 @@ 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,
@@ -188,7 +186,6 @@ 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),
@@ -244,7 +241,7 @@ testDynamic <- function(expr.mat = NULL,
                      data = null_mod_df,
                      family = MASS::negative.binomial(theta_hat),
                      corstr = cor.structure,
-                     scale.fix = gee.scale.fix,
+                     scale.fix = TRUE,
                      sandwich = ifelse(is.null(gee.bias.correction.method), FALSE, TRUE))
         }, silent = TRUE)
       } else if (is.glmm) {