diff --git a/NAMESPACE b/NAMESPACE
index 95f1dcb..c022fc3 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,6 +7,7 @@ export(enrichDynamicGenes)
 export(extractBreakpoints)
 export(fitGLMM)
 export(getFittedValues)
+export(getKnotDist)
 export(getResultsDE)
 export(marge2)
 export(modelLRT)
@@ -86,10 +87,10 @@ importFrom(glmmTMB,nbinom2)
 importFrom(parallel,clusterEvalQ)
 importFrom(parallel,clusterExport)
 importFrom(parallel,clusterSetRNGStream)
-importFrom(parallel,detectCores)
 importFrom(parallel,makeCluster)
 importFrom(parallel,stopCluster)
 importFrom(purrr,discard)
+importFrom(purrr,imap)
 importFrom(purrr,map)
 importFrom(purrr,map2)
 importFrom(purrr,map_chr)
diff --git a/R/summarizeModel.R b/R/summarizeModel.R
index e54e994..5282b21 100644
--- a/R/summarizeModel.R
+++ b/R/summarizeModel.R
@@ -24,7 +24,11 @@ summarizeModel <- function(marge.model = NULL, pt = NULL) {
                      Trend.Segment = NA_real_)
   } else {
     # extract model equation, slopes, & covariances
-    coef_vcov <- vcov(marge.model$final_mod)
+    if (marge.model$model_type == "GEE") {
+      coef_vcov <- as.matrix(marge.model$final_mod$var)
+    } else {
+      coef_vcov <- vcov(marge.model$final_mod)
+    }
     coef_df <- data.frame(coef_name = names(coef(marge.model$final_mod)),
                           coef_value = unname(coef(marge.model$final_mod)),
                           coef_var = unname(diag(coef_vcov)))
diff --git a/R/testDynamic.R b/R/testDynamic.R
index 06e9217..1111db5 100644
--- a/R/testDynamic.R
+++ b/R/testDynamic.R
@@ -8,7 +8,7 @@
 #' @importFrom bigstatsr as_FBM
 #' @importFrom foreach foreach %dopar% registerDoSEQ
 #' @importFrom doParallel registerDoParallel
-#' @importFrom parallel makeCluster detectCores stopCluster clusterEvalQ clusterExport clusterSetRNGStream
+#' @importFrom parallel makeCluster stopCluster clusterEvalQ clusterExport clusterSetRNGStream
 #' @importFrom withr with_output_sink
 #' @importFrom MASS glm.nb negative.binomial theta.mm
 #' @importFrom dplyr rename mutate relocate
@@ -93,6 +93,9 @@ testDynamic <- function(expr.mat = NULL,
   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)
+  }
   if (inherits(expr.mat, "SingleCellExperiment")) {
     expr.mat <- BiocGenerics::counts(expr.mat)[genes, ]
     expr.mat <- as.matrix(expr.mat)
@@ -102,15 +105,12 @@ testDynamic <- function(expr.mat = NULL,
                                      assay = Seurat::DefaultAssay(expr.mat))
     expr.mat <- as.matrix(expr.mat[genes, ])
   } else if (inherits(expr.mat, "dgCMatrix")) {
-    expr.mat <- as.matrix(expr.mat)
+    expr.mat <- as.matrix(expr.mat[genes, ])
+  } else {
+    expr.mat <- expr.mat[genes, ]
   }
   if (!(inherits(expr.mat, "matrix") || inherits(expr.mat, "array"))) { stop("Input expr.mat must be coerceable to a matrix of integer counts.") }
   expr.mat <- t(expr.mat)  # transpose to cell x gene matrix
-  if (is.null(genes)) {
-    genes <- colnames(expr.mat)
-  } else {
-    expr.mat <- expr.mat[, genes]
-  }
 
   # extract pseudotime dataframe if input is results from Slingshot
   if (inherits(pt, "SlingshotDataSet")) {
@@ -165,8 +165,7 @@ testDynamic <- function(expr.mat = NULL,
   package_list <- c("glm2", "scLANE", "MASS",  "bigstatsr", "broom.mixed", "dplyr", "stats")
   if (is.gee) {
     package_list <- c(package_list, "geeM")
-  }
-  if (is.glmm) {
+  } else if (is.glmm) {
     package_list <- c(package_list, "glmmTMB")
   }
 
diff --git a/tests/testthat/test_scLANE.R b/tests/testthat/test_scLANE.R
index 29bd595..f07554b 100644
--- a/tests/testthat/test_scLANE.R
+++ b/tests/testthat/test_scLANE.R
@@ -34,7 +34,7 @@ withr::with_output_sink(tempfile(), {
                                 size.factor.offset = cell_offset,
                                 n.cores = 2,
                                 track.time = TRUE)
-  gee_gene_stats <- testDynamic(sim_data,
+  gee_gene_stats <- testDynamic(expr.mat = sim_data,
                                 pt = pt_test,
                                 genes = genes_to_test,
                                 n.potential.basis.fns = 5,
@@ -53,7 +53,7 @@ withr::with_output_sink(tempfile(), {
                                  glmm.adaptive = TRUE,
                                  id.vec = sim_data$subject,
                                  n.cores = 2,
-                                 track.time = TRUE)
+                                 track.time = FALSE)
   # get results tables overall
   glm_test_results <- getResultsDE(glm_gene_stats)
   gee_test_results <- getResultsDE(gee_gene_stats)
@@ -191,6 +191,7 @@ withr::with_output_sink(tempfile(), {
   gsea_res <- enrichDynamicGenes(glm_test_results, species = "hsapiens")
   coef_summary_glm <- summarizeModel(marge_mod_offset, pt = pt_test)
   coef_summary_gee <- summarizeModel(marge_mod_GEE_offset, pt = pt_test)
+  knot_df <- getKnotDist(glm_gene_stats)
 })
 
 # run tests
@@ -359,3 +360,8 @@ test_that("summarizeModels() output", {
   expect_type(coef_summary_glm$Slope.Segment, "double")
   expect_type(coef_summary_gee$Slope.Segment, "double")
 })
+
+test_that("getKnotDist() output", {
+  expect_s3_class(knot_df, "data.frame")
+  expect_equal(ncol(knot_df), 3)
+})