diff --git a/DESCRIPTION b/DESCRIPTION
index 5cbd8f8..16ded79 100644
--- a/DESCRIPTION
+++ b/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. 
diff --git a/NEWS.md b/NEWS.md
index 085dbd7..6ea9d0d 100644
--- a/NEWS.md
+++ b/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. 
diff --git a/R/backward_sel_WIC.R b/R/backward_sel_WIC.R
index ad11f1b..975051b 100644
--- a/R/backward_sel_WIC.R
+++ b/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 {
diff --git a/R/marge2.R b/R/marge2.R
index 5e3f74b..4e24a08 100644
--- a/R/marge2.R
+++ b/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") {
diff --git a/R/score_fun_gee.R b/R/score_fun_gee.R
index 0c29751..adde7df 100644
--- a/R/score_fun_gee.R
+++ b/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,
diff --git a/R/stat_out_score_gee_null.R b/R/stat_out_score_gee_null.R
index f684d39..4fea157 100644
--- a/R/stat_out_score_gee_null.R
+++ b/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
diff --git a/R/testDynamic.R b/R/testDynamic.R
index 735bd1d..9bbf015 100644
--- a/R/testDynamic.R
+++ b/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) {
diff --git a/man/marge2.Rd b/man/marge2.Rd
index f807e80..962edc1 100644
--- a/man/marge2.Rd
+++ b/man/marge2.Rd
@@ -17,7 +17,6 @@ marge2(
   approx.knot = TRUE,
   n.knot.max = 50,
   glm.backend = "MASS",
-  gee.scale.fix = FALSE,
   tols_score = 1e-05,
   minspan = NULL,
   return.basis = FALSE,
@@ -50,8 +49,6 @@ marge2(
 
 \item{glm.backend}{(Optional) Character specifying which GLM-fitting backend should be used. Must be one of "MASS" or "speedglm". Defaults to "MASS".}
 
-\item{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.}
-
 \item{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.}
 
 \item{minspan}{(Optional) A set minimum span value. Defaults to NULL.}
diff --git a/man/testDynamic.Rd b/man/testDynamic.Rd
index 91d89a2..34b34b8 100644
--- a/man/testDynamic.Rd
+++ b/man/testDynamic.Rd
@@ -12,7 +12,6 @@ testDynamic(
   is.gee = FALSE,
   cor.structure = "ar1",
   gee.bias.correction.method = NULL,
-  gee.scale.fix = FALSE,
   is.glmm = FALSE,
   glmm.adaptive = TRUE,
   id.vec = NULL,
@@ -38,8 +37,6 @@ testDynamic(
 
 \item{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.}
 
-\item{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.}
-
 \item{is.glmm}{Should a GLMM framework be used instead of the default GLM? Defaults to FALSE.}
 
 \item{glmm.adaptive}{(Optional) Should the basis functions for the GLMM be chosen adaptively? If not, uses 4 evenly spaced knots. Defaults to TRUE.}