Merge branch 'dev' of https://github.com/jr-leary7/scLANE into dev

R/fitGLMM.R

@@ -11,7 +11,7 @@
 #' @param Y A vector of raw single cell counts. Defaults to NULL.
 #' @param Y.offset (Optional) An offset to be included in the final model fit. Defaults to NULL.
 #' @param id.vec A vector of subject IDs. Defaults to NULL.
-#' @param adaptive Should basis functions be chosen adaptively? Defaults to FALSE.
+#' @param adaptive Should basis functions be chosen adaptively? Defaults to TRUE.
 #' @param approx.knot Should knot approximation be used in the calls to \code{\link{marge2}}? This speeds up computation somewhat. Defaults to TRUE.
 #' @param M.glm The number of possible basis functions to use in the calls to \code{\link{marge2}} when choosing basis functions adaptively.
 #' @param return.basis (Optional) Whether the basis model matrix (denoted \code{B_final}) should be returned as part of the \code{marge} model object. Defaults to FALSE.
@@ -35,7 +35,7 @@ fitGLMM <- function(X_pred = NULL,
                     Y = NULL,
                     Y.offset = NULL,
                     id.vec = NULL,
-                    adaptive = FALSE,
+                    adaptive = TRUE,
                     approx.knot = TRUE,
                     M.glm = 3,
                     return.basis = FALSE,

R/summarizeModel.R

@@ -5,9 +5,14 @@
 #' @author Rhonda Bacher
 #' @import magrittr
 #' @importFrom stats vcov
-#' @description This function summarizes the model for each gene and allows for quantiative interpretation of fitted gene dynamics.
+#' @importFrom dplyr mutate select
+#' @importFrom purrr map
+#' @importFrom tidyr pivot_longer
+#' @importFrom tidyselect everything
+#' @description This function summarizes the model for each gene and allows for quantitative interpretation of fitted gene dynamics.
 #' @param marge.model The fitted model output from \code{\link{marge2}} (this function is internally called by \code{\link{testDynamic}}). Defaults to NULL.
 #' @param pt The predictor matrix of pseudotime. Defaults to NULL.
+#' @param is.glmm Flag specifying whether or not data were generated using GLMM mode. Defaults to FALSE.
 #' @return A data.frame of the model coefficients, cutpoint intervals, and formatted equations.
 #' @seealso \code{\link{marge2}}
 #' @export
@@ -20,7 +25,9 @@
 #'                       Y.offset = cell_offset)
 #' model_summary <- summarizeModel(marge.model = marge_model, pt = sim_pseudotime)
 
-summarizeModel <- function(marge.model = NULL, pt = NULL) {
+summarizeModel <- function(marge.model = NULL, 
+                           pt = NULL, 
+                           is.glmm = FALSE) {
   # check inputs
   if (is.null(marge.model)) { stop("Please provide a non-NULL input argument for marge.model.") }
   if (is.null(pt)) { stop("Please provide a non-NULL input argument for pt.") }
@@ -35,73 +42,92 @@ summarizeModel <- function(marge.model = NULL, pt = NULL) {
     } else {
       coef_vcov <- stats::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)))
-
-    coef_df <- coef_df[-which(coef_df$coef_name == "Intercept"), ]
-
-    coef_df <- cbind(coef_df, extractBreakpoints(marge.model))
-
-    coef_df <- coef_df[order(coef_df$Breakpoint, coef_df$Direction), ]
-
+    if (is.glmm) {
+      coef_df_full <- stats::coef(marge.model$final_mod)[[1]][[1]] %>% 
+                      dplyr::mutate(subject = rownames(.), .before = 1)
+      coef_df_list <- split(coef_df_full, coef_df_full$subject) %>% 
+                      purrr::map(\(x) dplyr::select(x, -subject)) %>% 
+                      purrr::map(\(x) { colnames(x)[colnames(x) == "(Intercept)"] <- "Intercept"; x }) %>% 
+                      purrr::map(\(x) {
+                        tidyr::pivot_longer(x, 
+                                            cols = tidyselect::everything(), 
+                                            names_to = "coef_name", 
+                                            values_to = "coef_value") %>% 
+                        dplyr::mutate(coef_var = NA_real_) %>% 
+                        as.data.frame()
+                      })
+    } else {
+      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)))
+      coef_df_list <- list(coef_df)
+    }
+    coef_df_list <- purrr::map(coef_df_list, \(x) x[-which(x$coef_name == "Intercept"), ]) %>% 
+                    purrr::map(\(x) cbind(x, extractBreakpoints(marge.model))) %>% 
+                    purrr::map(\(x) x[order(x$Breakpoint, x$Direction), ])
+    # summarize hinge coefficients
     MIN <- min(pt[, 1])
     MAX <- max(pt[, 1])
-
-    coef_ranges <- mapply(function(x, y) {
-      if (x == "Right") {
-        c(y, MAX)
-      } else if (x == "Left") {
-        c(MIN, y)
-      }
-    }, coef_df$Direction, coef_df$Breakpoint)
-
-    num_segments <- length(unique(coef_df$Breakpoint)) + 1
-    mod_seg <- list()
-    for (i in seq(num_segments)) {
-      if (i == 1) {
-        mod_seg[[i]] <- c(MIN, coef_df$Breakpoint[i])
-      } else if (i == num_segments) {
-        mod_seg[[i]] <- c(coef_df$Breakpoint[i - 1], MAX)
-      } else {
-        mod_seg[[i]] <- c(coef_df$Breakpoint[i - 1], coef_df$Breakpoint[i])
-      }
-    }
-    coef_ranges <- t(coef_ranges)
-    mod_seg_overlaps <- lapply(mod_seg, function(x) {
-      overlap_ind <- c()
-      for (i in seq(nrow(coef_ranges))) {
-        if (x[1] < coef_ranges[i, 2] && x[2] > coef_ranges[i, 1]) {
-          overlap_ind <- c(overlap_ind, i)
+    mod_summ_list <- purrr::map(coef_df_list, \(z) {
+      coef_ranges <- mapply(function(x, y) {
+        if (x == "Right") {
+          c(y, MAX)
+        } else if (x == "Left") {
+          c(MIN, y)
         }
-      }
-      return(overlap_ind)
-    })
-
-    seg_slopes <- lapply(mod_seg_overlaps, function(x) {
-      if (length(x) == 1) {
-        if (coef_df$Direction[x] == "Right") {
-          temp_slp <- coef_df$coef_value[x]
-        } else if (coef_df$Direction[x] == "Left") {
-          temp_slp <- -1 * coef_df$coef_value[x]
+      }, z$Direction, z$Breakpoint)
+      num_segments <- length(unique(z$Breakpoint)) + 1
+      mod_seg <- list()
+      for (i in seq(num_segments)) {
+        if (i == 1) {
+          mod_seg[[i]] <- c(MIN, z$Breakpoint[i])
+        } else if (i == num_segments) {
+          mod_seg[[i]] <- c(z$Breakpoint[i - 1], MAX)
+        } else {
+          mod_seg[[i]] <- c(z$Breakpoint[i - 1], z$Breakpoint[i])
         }
-      } else {
-        to_rev <- which(coef_df$Direction[x] == "Left")
-        coef_df$coef_value[x][to_rev] <- -1 * coef_df$coef_value[x][to_rev]
-        temp_slp <- sum(coef_df$coef_value[x])
       }
-      return(temp_slp)
+      coef_ranges <- t(coef_ranges)
+      mod_seg_overlaps <- lapply(mod_seg, function(x) {
+        overlap_ind <- c()
+        for (i in seq(nrow(coef_ranges))) {
+          if (x[1] < coef_ranges[i, 2] && x[2] > coef_ranges[i, 1]) {
+            overlap_ind <- c(overlap_ind, i)
+          }
+        }
+        return(overlap_ind)
+      })
+
+      seg_slopes <- lapply(mod_seg_overlaps, function(x) {
+        if (length(x) == 1) {
+          if (z$Direction[x] == "Right") {
+            temp_slp <- z$coef_value[x]
+          } else if (z$Direction[x] == "Left") {
+            temp_slp <- -1 * z$coef_value[x]
+          }
+        } else {
+          to_rev <- which(z$Direction[x] == "Left")
+          z$coef_value[x][to_rev] <- -1 * z$coef_value[x][to_rev]
+          temp_slp <- sum(z$coef_value[x])
+        }
+        return(temp_slp)
+      })
+      # combine results
+      seg_slopes <- do.call(c, seg_slopes)
+      # create discretized trend summary
+      seg_trends <- ifelse(seg_slopes > 0, 1, -1)
+      seg_trends[seg_slopes == 0] <- 0
+      # prepare results
+      mod_summ <- list(Breakpoint = unique(z$Breakpoint),
+                       Slope.Segment = seg_slopes,
+                       Trend.Segment = seg_trends)
+      return(mod_summ)
     })
-    # combine results
-    seg_slopes <- do.call(c, seg_slopes)
-    # create discretized trend summary
-    seg_trends <- ifelse(seg_slopes > 0, 1, -1)
-    seg_trends[seg_slopes == 0] <- 0
-    # prepare results
-    mod_summ <- list(Breakpoint = unique(coef_df$Breakpoint),
-                     Slope.Segment = seg_slopes,
-                     Trend.Segment = seg_trends)
-
+    if (!is.glmm) {
+      mod_summ <- mod_summ_list[[1]]
+    } else {
+      mod_summ <- mod_summ_list
+    }
   }
   return(mod_summ)
 }

R/testDynamic.R

@@ -15,6 +15,7 @@
 #' @importFrom MASS glm.nb negative.binomial theta.mm
 #' @importFrom dplyr rename mutate relocate
 #' @importFrom broom.mixed tidy
+#' @importFrom purrr imap reduce
 #' @importFrom stats predict logLik deviance offset
 #' @importFrom geeM geem
 #' @importFrom glmmTMB glmmTMB nbinom2
@@ -287,16 +288,24 @@ testDynamic <- function(expr.mat = NULL,
                                       Gene = genes[i],
                                       Lineage = LETTERS[j],
                                       .before = 1)
-     # solve values of slopes across pseudotime intervals -- TODO: add support for GLMM backend
-     if (!is.glmm) {
-       marge_dynamic_df <- summarizeModel(marge.model = marge_mod,
-                                                   pt = pt[lineage_cells, j, drop = FALSE])
+     # solve values of slopes across pseudotime intervals
+     marge_dynamic_df <- summarizeModel(marge.model = marge_mod, 
+                                        pt = pt[lineage_cells, j, drop = FALSE], 
+                                        is.glmm = is.glmm)
+     if (is.glmm) {
+       marge_dynamic_df <- purrr::imap(marge_dynamic_df, \(x, y) {
+         data.frame(t(unlist(x))) %>%
+         dplyr::mutate(Subject = y, 
+                       Gene = genes[i],
+                       Lineage = LETTERS[j],
+                       .before = 1)
+       })
+       marge_dynamic_df <- purrr::reduce(marge_dynamic_df, rbind)
+     } else {
        marge_dynamic_df <- data.frame(t(unlist(marge_dynamic_df))) %>%
                            dplyr::mutate(Gene = genes[i],
                                          Lineage = LETTERS[j],
                                          .before = 1)
-     } else {
-       marge_dynamic_df <- NULL
      }
 
      # format results list