changed marge2() to properly use GEE mode in backward_sel_WIC()

R/backward_sel_WIC.R

@@ -7,7 +7,7 @@
 #' @importFrom gamlss gamlss
 #' @importFrom geeM geem
 #' @importFrom MASS negative.binomial
-#' @importFrom stats vcov
+#' @importFrom stats vcov coef
 #' @param Y The response variable. Defaults to NULL.
 #' @param B_new The model matrix. Defaults to NULL.
 #' @param is.gee Is the model a GEE? Defaults to FALSE.
@@ -47,7 +47,7 @@ backward_sel_WIC <- function(Y = NULL,
     vcov_mat <- try({ stats::vcov(fit, type = "all") }, silent = TRUE)
     if (inherits(vcov_mat, "try-error")) {
       covmat_unscaled <- chol2inv(fit$mu.qr$qr[1:(fit$mu.df - fit$mu.nl.df), 1:(fit$mu.df - fit$mu.nl.df), drop = FALSE])
-      wald_stat <- unname(coef(fit) / sqrt(diag(covmat_unscaled))^2)[-1]
+      wald_stat <- unname(stats::coef(fit) / sqrt(diag(covmat_unscaled))^2)[-1]
     } else {
       wald_stat <- unname((vcov_mat$coef / vcov_mat$se)[-c(1, length(vcov_mat$coef))]^2)
     }

R/marge2.R

@@ -84,12 +84,10 @@ marge2 <- function(X_pred = NULL,
   }
   q <- ncol(X_pred)  # Number of predictor variables
   B <- as.matrix(rep(1, NN))  # Start with the intercept model.
-  # estimate "known" theta for GEE models
-  if (is.gee || !is.glmm) {
-    theta_hat <- MASS::theta.mm(y = Y,
-                                mu = mean(Y),
-                                dfr = NN - 1)
-  }
+  # estimate "known" theta for NB distribution
+  theta_hat <- MASS::theta.mm(y = Y,
+                              mu = mean(Y),
+                              dfr = NN - 1)
 
   pen <- 2  # penalty for GCV criterion -- could also switch to log(N) later
   colnames(B) <- "Intercept"
@@ -777,7 +775,13 @@ marge2 <- function(X_pred = NULL,
   colnames(wic_mat_2)[(ncol_B  + 1)] <- "Forward pass model"
 
   wic_mat_2[1, (ncol_B + 1)] <- full.wic
-  wic1_2 <- backward_sel_WIC(Y = Y, B_new = B_new)
+  wic1_2 <- backward_sel_WIC(Y = Y, 
+                             B_new = B_new, 
+                             is.gee = is.gee, 
+                             id.vec = id.vec, 
+                             cor.structure = cor.structure, 
+                             theta.hat = theta_hat,
+                             sandwich.var = sandwich.var)
   wic_mat_2[2, 2:(length(wic1_2) + 1)] <- wic1_2
   WIC_2 <- sum(apply(wic_mat_2[seq(2), ], 1, min, na.rm = TRUE)) + 2 * ncol(B_new)
   WIC_vec_2 <- c(WIC_vec_2, WIC_2)
@@ -788,7 +792,12 @@ marge2 <- function(X_pred = NULL,
   cnames_2 <- c(cnames_2, list(colnames(B_new_2)))
   for (i in 2:(ncol_B - 1)) {
     wic1_2 <- backward_sel_WIC(Y = Y,
-                               B_new = B_new_2)
+                               B_new = B_new_2, 
+                               is.gee = is.gee, 
+                               id.vec = id.vec, 
+                               cor.structure = cor.structure, 
+                               theta.hat = theta_hat,
+                               sandwich.var = sandwich.var)
     if (i != (ncol_B - 1)) {
       wic_mat_2[(i + 1), colnames(B_new_2)[-1]] <- wic1_2
       WIC_2 <- sum(apply(wic_mat_2[seq(i + 1), ], 1, min, na.rm = TRUE)) + 2 * ncol(B_new_2)