From 3247a66c8f91f06986e813ca665b66d9825319f1 Mon Sep 17 00:00:00 2001
From: Shirley Mathur <shirmath556@gmail.com>
Date: Wed, 10 Sep 2025 11:51:32 -0700
Subject: [PATCH 1/6] modified V_mat_contribution lin reg variance estimate

---
 R/V_matrix_contribution.R | 7 ++++---
 R/robust_score_test.R     | 2 +-
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/R/V_matrix_contribution.R b/R/V_matrix_contribution.R
index 6977d33..b35160c 100644
--- a/R/V_matrix_contribution.R
+++ b/R/V_matrix_contribution.R
@@ -4,12 +4,13 @@
 #' @param model_fits The fitted glm under the null hypothesis for the given indices.
 #' @param yy Vector of observed responses.
 #' @param xx Design matrix for model.
+#' @param m Number of parameters fixed under null hypothesis.
 #' @param corr_mat Working correlation matrix for model fit.
 #' @param family The model family for the fitted glm.
 #' @param link The link function utilized in the fitted glm.
 #' 
 
-V_matrix_contribution <- function(indices, model_fits, yy, xx, corr_mat, family, link) {
+V_matrix_contribution <- function(indices, model_fits, yy, xx, m = 1, corr_mat, family, link) {
   V_i <- matrix(NA, nrow = length(model_fits[indices]), ncol = length(model_fits[indices]))
   n_i <- length(indices)
   
@@ -31,8 +32,8 @@ V_matrix_contribution <- function(indices, model_fits, yy, xx, corr_mat, family,
   
   if (family == "gaussian" & link == "identity") {
     n <- length(yy)
-    p <- ncol(xx)
-    sigma2_tilde <- sum((yy - model_fits)^2)/(n - 1)
+    mod_df <- ncol(xx) - m
+    sigma2_tilde <- sum((yy - model_fits)^2)/(n - mod_df)
     if (n_i > 1) {
       V_i <- sqrt(diag(sigma2_tilde, n_i)) %*% corr_mat %*% sqrt(diag(sigma2_tilde, n_i))
     } else {
diff --git a/R/robust_score_test.R b/R/robust_score_test.R
index 899086a..02b30d4 100644
--- a/R/robust_score_test.R
+++ b/R/robust_score_test.R
@@ -73,7 +73,7 @@ robust_score_test <- function(glm_object, call_to_model, param = 1,
       corr_matrix <- matrix(rep(glm_object$geese$alpha, n_i^2), nrow = n_i)
       diag(corr_matrix) <- 1
       Vi <- V_matrix_contribution(indices = indices, model_fits = model0_fits, corr_mat = corr_matrix,
-                                  yy = yy, xx = xx, family = model1family, link = model1link)
+                                  yy = yy, xx = xx, m = pp0, family = model1family, link = model1link)
 
       Si <- S_matrix_contribution(indices = indices, model_fits = model0_fits,
                                   yy = yy, xx = xx, family = model1family, link = model1link)

From 015b4c7c21a28b0652a1c5e716a2dee04e10454f Mon Sep 17 00:00:00 2001
From: Shirley Mathur <shirmath556@gmail.com>
Date: Tue, 7 Oct 2025 18:56:41 -0700
Subject: [PATCH 2/6] implementing robust score test for linear reg with log
 link

---
 R/D_matrix_contribution.R    | 10 +++++++---
 R/S_matrix_contribution.R    | 17 +++--------------
 R/V_matrix_contribution.R    |  3 +--
 R/fisher_info_contribution.R | 14 ++++++++++++--
 R/gee_test.R                 |  4 ++--
 R/glm_test.R                 |  8 ++++----
 R/robust_score_test.R        | 34 +++++++++++++++++++++++-----------
 R/score_contribution.R       | 14 ++++++++++++--
 man/V_matrix_contribution.Rd | 13 ++++++++++++-
 9 files changed, 76 insertions(+), 41 deletions(-)

diff --git a/R/D_matrix_contribution.R b/R/D_matrix_contribution.R
index 34ab3db..a12e2ed 100644
--- a/R/D_matrix_contribution.R
+++ b/R/D_matrix_contribution.R
@@ -12,15 +12,19 @@ D_matrix_contribution <- function(indices, model_fits, yy, xx, family, link) {
   D_i <- matrix(NA, nrow = ncol(xx), ncol = length(model_fits))
   
   if (family == "poisson" & link == "log") {
-    D_i <- t(model_fits[indices]*xx[indices,])
+    D_i <- t(model_fits[indices]*xx[indices, , drop = FALSE])
   }
   
   if (family == "binomial" & link == "logit") {
-    D_i <- t(model_fits[indices]*(1-model_fits[indices])*xx[indices,])
+    D_i <- t(model_fits[indices]*(1-model_fits[indices])*xx[indices, , drop = FALSE])
   }
   
   if (family == "gaussian" & link == "identity") {
-    D_i <- t(xx[indices,])
+    D_i <- t(xx[indices, , drop = FALSE])
+  }
+  
+  if (family == "gaussian" & link == "log") {
+    D_i <- t(model_fits[indices]*xx[indices, , drop = FALSE])
   }
   
   return (D_i)
diff --git a/R/S_matrix_contribution.R b/R/S_matrix_contribution.R
index fe31518..c53b66a 100644
--- a/R/S_matrix_contribution.R
+++ b/R/S_matrix_contribution.R
@@ -9,19 +9,8 @@
 #' 
 
 S_matrix_contribution <- function(indices, model_fits, yy, xx, family, link) {
-  S_i <- matrix(NA, nrow = length(indices), ncol = 1)
-  
-  if (family == "poisson" & link == "log") {
-    S_i <- matrix(yy[indices] - model_fits[indices], ncol = 1)
-  }
-  
-  if (family == "binomial" & link == "logit") {
-    S_i <- matrix(yy[indices] - model_fits[indices], ncol = 1)
-  }
-  
-  if (family == "gaussian" & link == "identity") {
-    S_i <- matrix(yy[indices] - model_fits[indices], ncol = 1)
-  }
-  
+
+  S_i <- matrix(yy[indices] - model_fits[indices], ncol = 1)
+
   return (S_i)
 }
diff --git a/R/V_matrix_contribution.R b/R/V_matrix_contribution.R
index b35160c..e452588 100644
--- a/R/V_matrix_contribution.R
+++ b/R/V_matrix_contribution.R
@@ -30,7 +30,7 @@ V_matrix_contribution <- function(indices, model_fits, yy, xx, m = 1, corr_mat,
     }
   }
   
-  if (family == "gaussian" & link == "identity") {
+  if (family == "gaussian") {
     n <- length(yy)
     mod_df <- ncol(xx) - m
     sigma2_tilde <- sum((yy - model_fits)^2)/(n - mod_df)
@@ -39,7 +39,6 @@ V_matrix_contribution <- function(indices, model_fits, yy, xx, m = 1, corr_mat,
     } else {
       V_i <- sigma2_tilde * corr_mat
     }
-    
   }
   
   return (V_i)
diff --git a/R/fisher_info_contribution.R b/R/fisher_info_contribution.R
index ffa2d64..4a88202 100644
--- a/R/fisher_info_contribution.R
+++ b/R/fisher_info_contribution.R
@@ -4,11 +4,12 @@
 #' @param model_fits The fitted glm under the null hypothesis.
 #' @param yy Vector of observed responses.
 #' @param xx Design matrix for model.
+#' @param m Number of parameters fixed under null hypothesis.
 #' @param family The model family for the fitted glm.
 #' @param link The link function utilized in the fitted glm.
 #' 
 
-fisher_info_contribution <- function(i, model_fits, yy, xx, family, link) { ### output is p x p... Di^T Vi Di
+fisher_info_contribution <- function(i, model_fits, yy, xx, m = 1, family, link) { ### output is p x p... Di^T Vi Di
   fisher_info_mat <- matrix(NA, nrow = length(model_fits), ncol = 1)
   if (family == "poisson" & link == "log") {
     fisher_info_mat <- model_fits[i] * xx[i, ] %*% t(xx[i, ])
@@ -21,9 +22,18 @@ fisher_info_contribution <- function(i, model_fits, yy, xx, family, link) { ###
   if (family == "gaussian" & link == "identity") {
     n <- length(yy)
     p <- ncol(xx)
-    sigma2_tilde <- sum(yy - model_fits)/(n - p)
+    mod_df <- p-m
+    sigma2_tilde <- sum(yy - model_fits)/(n - mod_df)
     fisher_info_mat <- (1/sigma2_tilde)*xx[i, ] %*% t(xx[i, ])
   }
   
+  if (family == "gaussian" & link == "log") {
+    n <- length(yy)
+    p <- ncol(xx)
+    mod_df <- p-m
+    sigma2_tilde <- sum(yy - model_fits)/(n - mod_df)
+    fisher_info_mat <- ((model_fits[i]^2)/sigma2_tilde)*xx[i, ] %*% t(xx[i, ])
+  }
+  
   return(fisher_info_mat)
 }
\ No newline at end of file
diff --git a/R/gee_test.R b/R/gee_test.R
index f713ded..a058f0d 100644
--- a/R/gee_test.R
+++ b/R/gee_test.R
@@ -77,7 +77,7 @@ gee_test <- function(use_geeasy = TRUE, use_jack_se = FALSE, cluster_corr_coef =
     gee_fails <- TRUE
     new_cl <- cl
     new_cl[1] <- call("glm")
-    new_cl <- call_modify(new_cl, corstr = zap(), id = zap())
+    new_cl <- call_modify(new_cl, corstr = zap(), id = zap(), control = zap(), start = cl$control$init.beta)
     withCallingHandlers({
       glm_result <- try({eval(new_cl, envir = rlang::caller_env())})
     }, warning=function(w) {
@@ -92,7 +92,7 @@ gee_test <- function(use_geeasy = TRUE, use_jack_se = FALSE, cluster_corr_coef =
 
   gee_family <- gee_result$family$family
   gee_link <- gee_result$family$link
-  if ((gee_family != "poisson" | gee_link != "log") & (gee_family != "binomial" | gee_link != "logit") & (gee_family != "gaussian" | gee_link != "identity")) {
+  if ((gee_family != "poisson" | gee_link != "log") & (gee_family != "binomial" | gee_link != "logit") & (gee_family != "gaussian" | gee_link != "identity") &  (gee_family != "gaussian" | gee_link != "log")) {
     stop(paste("This is only implemented this for Poisson family with log link, Binomial family with logit link, and Gaussian family with identity link.\n",
                "You requested", gee_family, "family and", gee_link, "link \n",
                "Please open a GitHub issue if you're interested in other families."))
diff --git a/R/glm_test.R b/R/glm_test.R
index ac95798..aa67713 100644
--- a/R/glm_test.R
+++ b/R/glm_test.R
@@ -17,10 +17,9 @@ glm_test <- function(...) {
   # https://stackoverflow.com/questions/20680399/how-to-wrap-glm-in-a-function-pass-dotdotdot-directly-to-another-function-fail
   cl <- match.call()
   cl[1] <- call("glm")
-
   ### fit the glm, ignoring warnings about non integer inputs, as we take an estimating equations mindset
   withCallingHandlers({
-    glm_result <- eval(cl, envir = rlang::caller_env()) ### is the issue here??
+    glm_result <- eval(cl, envir = rlang::caller_env()) ### is the issue here?? 
   }, warning=function(w) {
     if (startsWith(conditionMessage(w), "non-integer x"))
       invokeRestart("muffleWarning")
@@ -28,8 +27,9 @@ glm_test <- function(...) {
 
   glm_family <- glm_result$family$family
   glm_link <- glm_result$family$link
-  if ((glm_family != "poisson" | glm_link != "log") & (glm_family != "binomial" | glm_link != "logit") & (glm_family != "gaussian" | glm_link != "identity")) {
-    stop(paste("This is only implemented this for Poisson family with log link, Binomial family with logit link, and Gaussian family with identity link.\n",
+  if ((glm_family != "poisson" | glm_link != "log") & (glm_family != "binomial" | glm_link != "logit") 
+      & (glm_family != "gaussian" | glm_link != "identity") & (glm_family != "gaussian" | glm_link != "log")) {
+    stop(paste("This is only implemented this for Poisson family with log link, Binomial family with logit link, and Gaussian family with identity or log link.\n",
                "You requested", glm_family, "family and", glm_link, "link \n",
                "Please open a GitHub issue if you're interested in other families."))
   }
diff --git a/R/robust_score_test.R b/R/robust_score_test.R
index 02b30d4..31d05b9 100644
--- a/R/robust_score_test.R
+++ b/R/robust_score_test.R
@@ -32,19 +32,30 @@ robust_score_test <- function(glm_object, call_to_model, param = 1,
   yy <- model1$y
   nn <- nrow(xx)
   pp <- ncol(xx)
-  xx0 <- xx[ , -param]
+  xx0 <- xx[ , -param, drop = FALSE]
   pp0 <- length(param)
 
 
 
   # stop("no")
   withCallingHandlers({
-    model0 <- glm.fit(x = xx0,
-                      y = yy,
-                      intercept = FALSE,
-                      offset = with(model1$data, eval(call_to_model$offset)),
-                      weights = with(model1$data, eval(call_to_model$weights)),
-                      family = eval(call_to_model$family))
+    if (model1family == "gaussian" & model1link == "log") {
+      model0 <- glm.fit(x = xx0,
+                        y = yy,
+                        intercept = FALSE,
+                        offset = with(model1$data, eval(call_to_model$offset)),
+                        weights = with(model1$data, eval(call_to_model$weights)),
+                        family = eval(call_to_model$family),
+                        start = rep(1, ncol(xx0)))
+    } else {
+      model0 <- glm.fit(x = xx0,
+                        y = yy,
+                        intercept = FALSE,
+                        offset = with(model1$data, eval(call_to_model$offset)),
+                        weights = with(model1$data, eval(call_to_model$weights)),
+                        family = eval(call_to_model$family))
+    }
+    
 
   }, warning = function(w) {
     if (startsWith(conditionMessage(w), "non-integer x"))
@@ -52,7 +63,7 @@ robust_score_test <- function(glm_object, call_to_model, param = 1,
   })
   model0_fits <- model0$fitted.values
 
-  if (is.factor(id)) {
+  if (length(id) == nn) {
 
     ## assume exchangeable
 
@@ -120,16 +131,17 @@ robust_score_test <- function(glm_object, call_to_model, param = 1,
                      solve(c_tilde %*% cov_U_hat %*% t(c_tilde), # (1 x p) x ((p x n) x (n x p)) x (p x 1)
                            (c_tilde %*% u_tilde_sum)))
 
-  } else if (is.na(id)) {
+  } else if (length(id) != nn) {
 
     u_tilde <- sapply(1:nn, score_contribution, model_fits = model0_fits,
                       family = model1family, link = model1link, xx = xx, yy = yy) ### p x n
 
     ## Reorder u tilde to match a
     u_tilde <- rbind(u_tilde[-param,], u_tilde[param,])
-
+  
     aa0 <- Reduce("+", sapply(1:nn, fisher_info_contribution, simplify=F,  model_fits = model0_fits,
-                              family = model1family, link = model1link, xx = xx, yy = yy)) ### p x p
+                              family = model1family, link = model1link, xx = xx, yy = yy, m = pp0)) ### p x p
+    
     aa0_11 <- aa0[setdiff(seq_len(pp), param), setdiff(seq_len(pp), param), drop = FALSE]
     aa0_21 <- aa0[param, setdiff(seq_len(pp), param), drop = FALSE]
     # c_tilde <- cbind(-aa0_21 %*% solve(aa0_11), diag(rep(1, pp0))) # 1 x p
diff --git a/R/score_contribution.R b/R/score_contribution.R
index 372f9e2..4294d50 100644
--- a/R/score_contribution.R
+++ b/R/score_contribution.R
@@ -4,11 +4,12 @@
 #' @param model_fits The fitted glm under the null hypothesis.
 #' @param yy Vector of observed responses.
 #' @param xx Design matrix for model.
+#' @param m Number of parameters fixed under null hypothesis.
 #' @param family The model family for the fitted glm.
 #' @param link The link function utilized in the fitted glm.
 #' 
 
-score_contribution <- function(i, model_fits, yy, xx, family, link) { ### output is p x 1
+score_contribution <- function(i, model_fits, yy, xx, m = 1, family, link) { ### output is p x 1
   score_mat <- matrix(NA, nrow = length(model_fits), ncol = 1)
   
   if (family == "poisson" & link == "log") {
@@ -23,9 +24,18 @@ score_contribution <- function(i, model_fits, yy, xx, family, link) { ### output
   if (family == "gaussian" & link == "identity") {
     n <- length(yy)
     p <- ncol(xx)
-    sigma2_tilde <- sum(yy - model_fits)/(n - p)
+    mod_df <- p-m
+    sigma2_tilde <- sum(yy - model_fits)/(n - mod_df)
     score_vec <- matrix((1/sigma2_tilde)*(yy[i] - model_fits[i])*xx[i,], ncol = 1)
   }
+  
+  if (family == "gaussian" & link == "log") {
+    n <- length(yy)
+    p <- ncol(xx)
+    mod_df <- p-m
+    sigma2_tilde <- sum(yy - model_fits)/(n - mod_df)
+    score_vec <- matrix((1/sigma2_tilde)*model_fits[i]*(yy[i] - model_fits[i])*xx[i,], ncol = 1)
+  }
 
   return(score_vec)
 }
\ No newline at end of file
diff --git a/man/V_matrix_contribution.Rd b/man/V_matrix_contribution.Rd
index d46754a..1e6237c 100644
--- a/man/V_matrix_contribution.Rd
+++ b/man/V_matrix_contribution.Rd
@@ -4,7 +4,16 @@
 \alias{V_matrix_contribution}
 \title{Compute V matrix contribution to robust score statistic for glm robust score tests.}
 \usage{
-V_matrix_contribution(indices, model_fits, yy, xx, corr_mat, family, link)
+V_matrix_contribution(
+  indices,
+  model_fits,
+  yy,
+  xx,
+  m = 1,
+  corr_mat,
+  family,
+  link
+)
 }
 \arguments{
 \item{indices}{Indices of observations for cluster.}
@@ -15,6 +24,8 @@ V_matrix_contribution(indices, model_fits, yy, xx, corr_mat, family, link)
 
 \item{xx}{Design matrix for model.}
 
+\item{m}{Number of parameters fixed under null hypothesis.}
+
 \item{corr_mat}{Working correlation matrix for model fit.}
 
 \item{family}{The model family for the fitted glm.}

From b9caf7986400934da40139ce0d37cbc4eea6895f Mon Sep 17 00:00:00 2001
From: Shirley Mathur <shirmath556@gmail.com>
Date: Fri, 10 Oct 2025 14:55:05 -0700
Subject: [PATCH 3/6] Revert "modified V_mat_contribution lin reg variance
 estimate"

This reverts commit 3247a66c8f91f06986e813ca665b66d9825319f1.
---
 R/V_matrix_contribution.R | 7 +++----
 R/robust_score_test.R     | 2 +-
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/R/V_matrix_contribution.R b/R/V_matrix_contribution.R
index e452588..98ba8fe 100644
--- a/R/V_matrix_contribution.R
+++ b/R/V_matrix_contribution.R
@@ -4,13 +4,12 @@
 #' @param model_fits The fitted glm under the null hypothesis for the given indices.
 #' @param yy Vector of observed responses.
 #' @param xx Design matrix for model.
-#' @param m Number of parameters fixed under null hypothesis.
 #' @param corr_mat Working correlation matrix for model fit.
 #' @param family The model family for the fitted glm.
 #' @param link The link function utilized in the fitted glm.
 #' 
 
-V_matrix_contribution <- function(indices, model_fits, yy, xx, m = 1, corr_mat, family, link) {
+V_matrix_contribution <- function(indices, model_fits, yy, xx, corr_mat, family, link) {
   V_i <- matrix(NA, nrow = length(model_fits[indices]), ncol = length(model_fits[indices]))
   n_i <- length(indices)
   
@@ -32,8 +31,8 @@ V_matrix_contribution <- function(indices, model_fits, yy, xx, m = 1, corr_mat,
   
   if (family == "gaussian") {
     n <- length(yy)
-    mod_df <- ncol(xx) - m
-    sigma2_tilde <- sum((yy - model_fits)^2)/(n - mod_df)
+    p <- ncol(xx)
+    sigma2_tilde <- sum((yy - model_fits)^2)/(n - 1)
     if (n_i > 1) {
       V_i <- sqrt(diag(sigma2_tilde, n_i)) %*% corr_mat %*% sqrt(diag(sigma2_tilde, n_i))
     } else {
diff --git a/R/robust_score_test.R b/R/robust_score_test.R
index 31d05b9..ec98498 100644
--- a/R/robust_score_test.R
+++ b/R/robust_score_test.R
@@ -84,7 +84,7 @@ robust_score_test <- function(glm_object, call_to_model, param = 1,
       corr_matrix <- matrix(rep(glm_object$geese$alpha, n_i^2), nrow = n_i)
       diag(corr_matrix) <- 1
       Vi <- V_matrix_contribution(indices = indices, model_fits = model0_fits, corr_mat = corr_matrix,
-                                  yy = yy, xx = xx, m = pp0, family = model1family, link = model1link)
+                                  yy = yy, xx = xx, family = model1family, link = model1link)
 
       Si <- S_matrix_contribution(indices = indices, model_fits = model0_fits,
                                   yy = yy, xx = xx, family = model1family, link = model1link)

From 68098b362a649e3765f6afb092b4f1c2c3c636f3 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Fri, 14 Nov 2025 16:59:28 -0800
Subject: [PATCH 4/6] Update test-basic.R

update tests for gaussian with log link (previously an error was expected for this, so fix that failing test since there is no longer an error now that this is implemented)
---
 tests/testthat/test-basic.R | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/tests/testthat/test-basic.R b/tests/testthat/test-basic.R
index 6900b15..49e4035 100644
--- a/tests/testthat/test-basic.R
+++ b/tests/testthat/test-basic.R
@@ -64,6 +64,9 @@ test_that("test gaussian example", {
   
   expect_type(glm_test(dist ~ speed, data = cars, family=gaussian(link = "identity"))$coef_tab,
               "double")
+
+  expect_type(glm_test(dist ~ speed, data = cars, family=gaussian(link = "log"))$coef_tab,
+              "double")
   
   cars$id <- c(rep(1:5, 9), 6:10)
   cars$bin <- rep(0:1, each = 25)
@@ -76,7 +79,7 @@ test_that("test gaussian example", {
 
 test_that("error with other model", {
   
-  expect_error(glm_test(dist ~ speed, data = cars, family=gaussian(link = "log")))
+  expect_error(glm_test(dist ~ speed, data = cars, family=gaussian(link = "logit")))
 
 })
 

From 51253b49d152d6cd9389eba20cee77e15819c95a Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Tue, 18 Nov 2025 09:40:15 -0800
Subject: [PATCH 5/6] update documentation for fisher_info and
 score_contribution

---
 man/fisher_info_contribution.Rd           |  4 +++-
 man/score_contribution.Rd                 |  4 +++-
 tests/testthat/test-multinom_pseudo_inv.R | 15 +++++++++++----
 3 files changed, 17 insertions(+), 6 deletions(-)

diff --git a/man/fisher_info_contribution.Rd b/man/fisher_info_contribution.Rd
index feeea10..2e84484 100644
--- a/man/fisher_info_contribution.Rd
+++ b/man/fisher_info_contribution.Rd
@@ -4,7 +4,7 @@
 \alias{fisher_info_contribution}
 \title{Compute fisher information contribution to robust score statistic for glm robust score tests.}
 \usage{
-fisher_info_contribution(i, model_fits, yy, xx, family, link)
+fisher_info_contribution(i, model_fits, yy, xx, m = 1, family, link)
 }
 \arguments{
 \item{i}{Index of observation.}
@@ -15,6 +15,8 @@ fisher_info_contribution(i, model_fits, yy, xx, family, link)
 
 \item{xx}{Design matrix for model.}
 
+\item{m}{Number of parameters fixed under null hypothesis.}
+
 \item{family}{The model family for the fitted glm.}
 
 \item{link}{The link function utilized in the fitted glm.}
diff --git a/man/score_contribution.Rd b/man/score_contribution.Rd
index 2d70ebd..da6ea1a 100644
--- a/man/score_contribution.Rd
+++ b/man/score_contribution.Rd
@@ -4,7 +4,7 @@
 \alias{score_contribution}
 \title{Compute score contribution to robust score statistic for glm robust score tests.}
 \usage{
-score_contribution(i, model_fits, yy, xx, family, link)
+score_contribution(i, model_fits, yy, xx, m = 1, family, link)
 }
 \arguments{
 \item{i}{Index of observation.}
@@ -15,6 +15,8 @@ score_contribution(i, model_fits, yy, xx, family, link)
 
 \item{xx}{Design matrix for model.}
 
+\item{m}{Number of parameters fixed under null hypothesis.}
+
 \item{family}{The model family for the fitted glm.}
 
 \item{link}{The link function utilized in the fitted glm.}
diff --git a/tests/testthat/test-multinom_pseudo_inv.R b/tests/testthat/test-multinom_pseudo_inv.R
index bcc2753..d0eecfe 100644
--- a/tests/testthat/test-multinom_pseudo_inv.R
+++ b/tests/testthat/test-multinom_pseudo_inv.R
@@ -32,10 +32,17 @@ test_that("the test stat with the pseudo inverse controls t1e when the inverse i
     df <- simulate_data_mult(nn = nn, strong=TRUE, ms = 10)
     df$Y[, 1] <- 0
     df$Y[, 3] <- 0
-    result <- suppressWarnings(multinom_test(df$X, df$Y, strong=TRUE))
-    result_alt <- suppressWarnings(multinom_test(df$X, df$Y, strong=TRUE, pseudo_inv = TRUE))
-    p[sim] <- result$p
-    p_alt[sim] <- result_alt$p
+    capture.output(
+      res <- suppressWarnings(multinom_test(df$X, df$Y, strong=TRUE)),
+      file = NULL
+    )
+    capture.output(
+      res_alt <- suppressWarnings(multinom_test(df$X, df$Y, strong=TRUE, pseudo_inv = TRUE)),
+      file = NULL
+    )
+    
+    if (!inherits(res, "try-error")) p[sim] <- res$p
+    if (!inherits(res_alt, "try-error")) p_alt[sim] <- res_alt$p
   }
 
   expect_true(sum(!is.na(p)) > 5)

From bc3387088ad4be1c2a75b42e90027f934b47de11 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Tue, 25 Nov 2025 12:15:29 -0800
Subject: [PATCH 6/6] update fisher_info_contribution to use sum of squared
 residuals for gaussian model

---
 R/fisher_info_contribution.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/fisher_info_contribution.R b/R/fisher_info_contribution.R
index 4a88202..c991619 100644
--- a/R/fisher_info_contribution.R
+++ b/R/fisher_info_contribution.R
@@ -23,7 +23,7 @@ fisher_info_contribution <- function(i, model_fits, yy, xx, m = 1, family, link)
     n <- length(yy)
     p <- ncol(xx)
     mod_df <- p-m
-    sigma2_tilde <- sum(yy - model_fits)/(n - mod_df)
+    sigma2_tilde <- sum((yy - model_fits)^2)/(n - mod_df)
     fisher_info_mat <- (1/sigma2_tilde)*xx[i, ] %*% t(xx[i, ])
   }
   
@@ -31,7 +31,7 @@ fisher_info_contribution <- function(i, model_fits, yy, xx, m = 1, family, link)
     n <- length(yy)
     p <- ncol(xx)
     mod_df <- p-m
-    sigma2_tilde <- sum(yy - model_fits)/(n - mod_df)
+    sigma2_tilde <- sum((yy - model_fits)^2)/(n - mod_df)
     fisher_info_mat <- ((model_fits[i]^2)/sigma2_tilde)*xx[i, ] %*% t(xx[i, ])
   }