fixed wald test statistic calculation -- addresses #222

R/biasCorrectGEE.R

@@ -1,38 +1,38 @@
-#' Bias-correct the GEE sandwich variance-covariance matrix. 
-#' 
+#' Bias-correct the GEE sandwich variance-covariance matrix.
+#'
 #' @name biasCorrectGEE
 #' @author Jack Leary
-#' @description This functions implements several bias-correction methods for the GEE sandwich variance-covariance matrix; they are to be used when the number of subjects is small or the numer of timepoints per-subject is very large. 
-#' @importFrom stats fitted.values cor 
+#' @description This functions implements several bias-correction methods for the GEE sandwich variance-covariance matrix; they are to be used when the number of subjects is small or the numer of timepoints per-subject is very large.
+#' @importFrom stats fitted.values cor
 #' @importFrom dplyr with_groups summarise
 #' @importFrom MASS ginv
 #' @importFrom Matrix bdiag
 #' @param fitted.model The fitted model of class \code{geem} returned by \code{\link{marge2}}. Defaults to NULL.
 #' @param correction.method A string specifying the correction method to be used. Currently supported options are "df" and "kc". Defaults to "kc".
 #' @param id.vec A vector of subject IDs. Defaults to NULL.
-#' @param cor.structure A string specifying the correlation structure used in fitting the model. Defaults to "ar1".  
-#' @param verbose (Optional) A Boolean specifying whether or not verbose output should be printed to the console. Occasionally useful for debugging. Defaults to FALSE. 
-#' @return An object of class \code{matrix} containing the bias-corrected variance-covariance estimates. 
+#' @param cor.structure A string specifying the correlation structure used in fitting the model. Defaults to "ar1".
+#' @param verbose (Optional) A Boolean specifying whether or not verbose output should be printed to the console. Occasionally useful for debugging. Defaults to FALSE.
+#' @return An object of class \code{matrix} containing the bias-corrected variance-covariance estimates.
 #' @seealso \code{\link{waldTestGEE}}
 
 biasCorrectGEE <- function(fitted.model = NULL,
-                           correction.method = "kc", 
-                           id.vec = NULL, 
-                           cor.structure = "ar1", 
+                           correction.method = "kc",
+                           id.vec = NULL,
+                           cor.structure = "ar1",
                            verbose = FALSE) {
-  # check inputs 
+  # check inputs
   if (is.null(fitted.model)) { stop("Arguments to biasCorrectGEE() are missing.") }
   if (!inherits(fitted.model, "geem")) { stop("The fitted model must be of class 'geem'.") }
   correction.method <- tolower(correction.method)
   if (!correction.method %in% c("df", "kc")) { stop("Unsupported bias correction method in waldTestGEE().") }
-  if (correction.method == "kc" && is.null(id.vec)) { stop("Kauermann and Carroll requires the provision of a vector of subject IDs and a vector of raw gene counts.") }
-  cor.structure <- tolower(cor.structure) 
+  if (correction.method == "kc" && is.null(id.vec)) { stop("Kauermann and Carroll requires the provision of a vector of subject IDs.") }
+  cor.structure <- tolower(cor.structure)
   if (!cor.structure %in% c("exchangeable", "ar1")) { stop("Unrecognized correlation structure in biasCorrectGEE().") }
-  # extract sandwich variance-covariance matrix 
+  # extract sandwich variance-covariance matrix
   V_sandwich <- as.matrix(fitted.model$var)
   n_s <- length(fitted.model$clusz)
   if (correction.method == "df") {
-    # compute degrees-of-freedom bias-corrected variance-covariance matrix 
+    # compute degrees-of-freedom bias-corrected variance-covariance matrix
     p <- ncol(fitted.model$X)
     if (p >= n_s) {
       warning("Cannot perform DF bias correction on sandwich variance-covariance matrix as the number of subjects is less than or equal to the number of covariates.")
@@ -64,9 +64,9 @@ biasCorrectGEE <- function(fitted.model = NULL,
         }
         return(res)
       }
-      rho_estimates <- dplyr::with_groups(resid_df, 
-                                          subject, 
-                                          dplyr::summarise, 
+      rho_estimates <- dplyr::with_groups(resid_df,
+                                          subject,
+                                          dplyr::summarise,
                                           rho = rhoExchangeable(r_hat))
       rho_avg <- mean(rho_estimates$rho, na.rm = TRUE)
       if (verbose) {
@@ -78,15 +78,15 @@ biasCorrectGEE <- function(fitted.model = NULL,
         if (n_i < 2) {
           res <- NA_real_
         } else {
-          res <- stats::cor(residuals[-n_i], 
-                            residuals[-1], 
+          res <- stats::cor(residuals[-n_i],
+                            residuals[-1],
                             use = "complete.obs")
         }
         return(res)
       }
-      rho_estimates <- dplyr::with_groups(resid_df, 
-                                          subject, 
-                                          dplyr::summarise, 
+      rho_estimates <- dplyr::with_groups(resid_df,
+                                          subject,
+                                          dplyr::summarise,
                                           rho = rhoAR1(r_hat))
       rho_avg <- mean(rho_estimates$rho, na.rm = TRUE)
       if (verbose) {
@@ -126,8 +126,8 @@ biasCorrectGEE <- function(fitted.model = NULL,
       }, silent = TRUE)
       if (inherits(Var_Yi_inv, "try-error")) {
         if (verbose) {
-          warning(paste0("Covariance matrix inversion failed for subject ",  
-                         subjects[s], 
+          warning(paste0("Covariance matrix inversion failed for subject ",
+                         subjects[s],
                          ", applying diagonal matrix as fallback in biasCorrectGEE()."))
         }
         Var_Yi_inv <- diag(1, n_i, n_i)

R/waldTestGEE.R

@@ -9,7 +9,7 @@
 #' @param mod.0 The model under the null hypothesis. Must be of class \code{geem}. Defaults to NULL.
 #' @param bias.correct Boolean specifying whether a small-sample bias correction should be applied to the estimated sandwich variance-covariance matrix. Useful when the number of subjects is small. Defaults to FALSE.
 #' @param correction.method A string specifying the correction method to be used. Currently supported options are "df" and "kc". Defaults to "kc".
-#' @param id.vec A vector of subject IDs. Required when \code{correction.method} is "kc". Defaults to NULL. 
+#' @param id.vec A vector of subject IDs. Required when \code{correction.method} is "kc". Defaults to NULL.
 #' @param verbose (Optional) A Boolean specifying whether or not verbose output should be printed to the console. Occasionally useful for debugging. Defaults to FALSE.
 #' @return A list containing the Wald test statistic, a \emph{p}-value, and the degrees of freedom used in the test.
 #' @details
@@ -24,10 +24,11 @@
 #' @seealso \code{\link{modelLRT}}
 
 waldTestGEE <- function(mod.1 = NULL,
-                        mod.0 = NULL, 
-                        bias.correct = FALSE, 
-                        correction.method = "df", 
-                        id.vec = NULL, verbose = FALSE) {
+                        mod.0 = NULL,
+                        bias.correct = FALSE,
+                        correction.method = "kc",
+                        id.vec = NULL,
+                        verbose = FALSE) {
   # check inputs
   if (inherits(mod.1, "try-error") || inherits(mod.0, "try-error")) {
     res <- list(Wald_Stat = NA_real_,
@@ -38,9 +39,10 @@ waldTestGEE <- function(mod.1 = NULL,
   }
   correction.method <- tolower(correction.method)
   if (!correction.method %in% c("df", "kc")) { stop("Unsupported bias correction method in waldTestGEE().") }
-  if (bias.correct && correction.method == "kc" && is.null(id.vec)) { stop("The Kauermann and Carroll bias correction method requires a vector of subject IDs.") } 
+  if (bias.correct && correction.method == "kc" && is.null(id.vec)) { stop("The Kauermann and Carroll bias correction method requires a vector of subject IDs.") }
   mod.1 <- mod.1$final_mod
   if (is.null(mod.1) || is.null(mod.0) || !(inherits(mod.1, "geem") && inherits(mod.0, "geem"))) { stop("You must provide two geeM models to waldTestGee().") }
+  if (length(coef(mod.0)) != 1) { stop("Null GEE model must be intercept-only.") }
   if (length(coef(mod.1)) <= length(coef(mod.0))) {
     # can't calculate Wald statistic if both models are intercept-only
     res <- list(Wald_Stat = 0,
@@ -50,7 +52,7 @@ waldTestGEE <- function(mod.1 = NULL,
   } else {
     # compute test statistic & asymptotic p-value
     coef_alt_mod <- names(coef(mod.1))
-    coef_null_mod <- names(coef(mod.0))
+    coef_null_mod <- c("Intercept")
     coef_diff <- setdiff(coef_alt_mod, coef_null_mod)
     coef_idx <- rep(0, length(coef_diff))
     for (i in seq_len(length(coef_diff))) {
@@ -60,9 +62,9 @@ waldTestGEE <- function(mod.1 = NULL,
     vcov_mat <- as.matrix(mod.1$var)
     if (bias.correct) {
       vcov_mat <- biasCorrectGEE(mod.1,
-                                 correction.method = correction.method, 
-                                 id.vec = id.vec, 
-                                 cor.structure = mod.1$corr, 
+                                 correction.method = correction.method,
+                                 id.vec = id.vec,
+                                 cor.structure = mod.1$corr,
                                  verbose = verbose)
     }
     vcov_mat <- vcov_mat[coef_idx, coef_idx]