Merge pull request #268 from jr-leary7/dev

Dev

NEWS.md

@@ -1,9 +1,11 @@
 # Changes in v0.8.7 
 
-+ Switched GEE fitting back to use `scale.fix = FALSE` and a fixed value for the Negative-binomial overdispersion parameter as it improves model fits. 
++ Switched GEE fitting back to use `scale.fix = FALSE` and substituted a fixed value for the Negative-binomial overdispersion parameter (instead of estimating via method-of-moments) as it improves model fits. 
 + Added option to use a Lagrange Multiplier (Score) test for GEE mode instead of the default Wald test. The relevant argument is `gee.test` in `testDynamic()`.
 + Updated documentation and some tests. 
 + Added column called `Null_Fit_Notes` to output from `getResultsDE()` to describe when and how null models fail. This doesn't happen frequently, but it's good info to have when it does. 
++ Expanded test suite to include C++ functions. 
++ Fixed incorrectly-implemented Lagrange Multiplier test after initial code was seen to be incorrect. 
 
 # Changes in v0.8.6
 

R/biasCorrectGEE.R

@@ -136,12 +136,13 @@ biasCorrectGEE <- function(fitted.model = NULL,
     }
     W <- as.matrix(Matrix::bdiag(cov_matrices))
     X <- fitted.model$X
-    XWX <- t(X) %*% W %*% X
+    X_t <- t(X)
+    XWX <- X_t %*% W %*% X
     XWX_inv <- try({ eigenMapMatrixInvert(XWX, n_cores = 1L) }, silent = TRUE)
     if (inherits(XWX_inv, "try-error")) {
       XWX_inv <- eigenMapPseudoInverse(XWX, n_cores = 1L)
     }
-    H <- X %*% XWX_inv %*% t(X) %*% W
+    H <- X %*% XWX_inv %*% X_t %*% W
     tr_H <- sum(diag(H))
     if (verbose) {
       message(paste0("Trace of projection matrix H estimated at: ", round(tr_H, 5)))

R/fitGLMM.R

@@ -14,9 +14,9 @@
 #' @param id.vec A vector of subject IDs. Defaults to NULL.
 #' @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 M.glm The number of possible basis functions to use in the calls to \code{\link{marge2}} when choosing basis functions adaptively. Defaults to 3. 
 #' @param reg.penalty (Optional) String specifying the penalty type to be used when fitting a regularized negative-binomial model to select optimal basis functions. Defaults to "snet". 
-#' @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.
+#' @param return.basis (Optional) Whether the basis model matrix (denoted \code{basis_mtx}) should be returned as part of the \code{marge} model object. Defaults to FALSE.
 #' @param return.GCV (Optional) Whether the final GCV value should be returned as part of the \code{marge} model object. Defaults to FALSE.
 #' @param verbose (Optional) Should intermediate output be printed to the console? Defaults to FALSE.
 #' @return An object of class \code{marge} containing the fitted model & other optional quantities of interest (basis function matrix, GCV, etc.).

R/getResultsDE.R

@@ -40,7 +40,7 @@ getResultsDE <- function(test.dyn.res = NULL,
                                      function(x) {
                                        purrr::map_dfr(x,
                                                       function(y) {
-                                                        as.data.frame(rbind(y[seq(14)])) %>%
+                                                        as.data.frame(rbind(y[seq(15)])) %>%
                                                           dplyr::mutate(dplyr::across(tidyselect::everything(), \(z) unname(unlist(z))))
                                                       })
                                      }) %>%

R/testDynamic.R

@@ -324,6 +324,7 @@ testDynamic <- function(expr.mat = NULL,
                                                        ifelse(gee.test == "wald", "Wald", "Score"), 
                                                        "LRT"),
                                Test_Stat_Note = NA_character_,
+                               Degrees_Freedom = NA_real_, 
                                P_Val = NA_real_,
                                LogLik_MARGE = marge_sumy$ll_marge,
                                LogLik_Null = null_sumy$null_ll,
@@ -366,6 +367,7 @@ testDynamic <- function(expr.mat = NULL,
                                            ifelse(gee.test == "wald", test_res$Wald_Stat, test_res$Score_Stat), 
                                            test_res$LRT_Stat)
      lineage_list[[j]]$Test_Stat_Note <- test_res$Notes
+     lineage_list[[j]]$Degrees_Freedom <- test_res$DF 
      lineage_list[[j]]$P_Val <- test_res$P_Val
     }
     names(lineage_list) <- paste0("Lineage_", LETTERS[seq(n_lineages)])
@@ -393,6 +395,7 @@ testDynamic <- function(expr.mat = NULL,
                                      ifelse(gee.test == "wald", "Wald", "Score"), 
                                      "LRT"),
              Test_Stat_Note = NA_character_,
+             Degrees_Freedom = NA_real_, 
              P_Val = NA_real_,
              LogLik_MARGE = NA_real_,
              LogLik_Null = NA_real_,

tests/testthat/test_scLANE.R

@@ -25,10 +25,18 @@ null_stat_gee <- stat_out_score_gee_null(Y = Y_exp,
                                          theta.hat = 1)
 tp1_res <- tp1(x = rnorm(30), t = 0)
 tp2_res <- tp2(x = rnorm(30), t = 0)
+# generate PD matrix for use in testing C++ matrix functions 
+n <- 25
+A <- matrix(rnorm(n^2), nrow = n, ncol = n)
+A <- t(A) %*% A
 
 # generate scLANE results w/ all three modes
 withr::with_output_sink(tempfile(), {
-  # choose candidate genes
+  # C++ matrix operations 
+  B <- eigenMapMatMult(A, A)
+  C <- eigenMapMatrixInvert(A)
+  D <- eigenMapPseudoInverse(A)
+  # candidate gene selection
   candidate_genes <- chooseCandidateGenes(sim_data_seu,
                                           group.by.subject = TRUE,
                                           id.vec = sim_data_seu$subject,
@@ -113,6 +121,7 @@ withr::with_output_sink(tempfile(), {
                           is.gee = TRUE,
                           id.vec = sim_data$subject,
                           cor.structure = "ar1",
+                          sandwich.var = TRUE, 
                           return.basis = TRUE,
                           return.GCV = TRUE,
                           return.WIC = TRUE)
@@ -142,6 +151,10 @@ withr::with_output_sink(tempfile(), {
                              alt.df = as.data.frame(marge_mod_GEE_offset$basis_mtx), 
                              null.df = data.frame(Y = counts_test[, 3]), 
                              id.vec = sim_data$subject)
+  # bias-correct GEE sandwich variance-covariance matrix 
+  V_kc <- biasCorrectGEE(marge_mod_GEE$final_mod, 
+                         correction.method = "kc", 
+                         id.vec = sim_data$subject)
   # run GLMM model -- no offset
   glmm_mod <- fitGLMM(pt_test,
                       Y = counts_test[, 4],
@@ -282,6 +295,18 @@ withr::with_output_sink(tempfile(), {
 })
 
 # run tests
+test_that("internal C++ functions", {
+  expect_type(B, "double")
+  expect_equal(ncol(B), 25)
+  expect_equal(nrow(B), 25)
+  expect_type(C, "double")
+  expect_equal(ncol(C), 25)
+  expect_equal(nrow(C), 25)
+  expect_type(D, "double")
+  expect_equal(ncol(D), 25)
+  expect_equal(nrow(D), 25)
+})
+
 test_that("internal marge functions", {
   expect_type(min_span_res, "double")
   expect_type(max_span_res, "double")
@@ -386,6 +411,7 @@ test_that("marge2() output -- GEE backend", {
   expect_equal(marge_mod_GEE_offset$model_type, "GEE")
   expect_true(marge_mod_GEE$final_mod$converged)
   expect_true(marge_mod_GEE_offset$final_mod$converged)
+  expect_type(V_kc, "double")
 })
 
 test_that("Statistical testing output", {