From 4ab1f4dbb98a45a409b23cf7578842535a77d952 Mon Sep 17 00:00:00 2001
From: Shirley Mathur <shirmath556@gmail.com>
Date: Tue, 26 Aug 2025 16:09:00 -0700
Subject: [PATCH 01/17] implemented robust score test for uncorrelated linear
 reg

---
 R/fisher_info_contribution.R | 7 +++++++
 R/glm_test.R                 | 4 ++--
 R/score_contribution.R       | 7 +++++++
 3 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/R/fisher_info_contribution.R b/R/fisher_info_contribution.R
index 844f818..ffa2d64 100644
--- a/R/fisher_info_contribution.R
+++ b/R/fisher_info_contribution.R
@@ -18,5 +18,12 @@ fisher_info_contribution <- function(i, model_fits, yy, xx, family, link) { ###
     fisher_info_mat <- (1-model_fits[i])^(-1)*(model_fits[i])^(-1) * xx[i, ] %*% t(xx[i, ])
   }
   
+  if (family == "gaussian" & link == "identity") {
+    n <- length(yy)
+    p <- ncol(xx)
+    sigma2_tilde <- sum(yy - model_fits)/(n - p)
+    fisher_info_mat <- (1/sigma2_tilde)*xx[i, ] %*% t(xx[i, ])
+  }
+  
   return(fisher_info_mat)
 }
\ No newline at end of file
diff --git a/R/glm_test.R b/R/glm_test.R
index 8da6616..ac95798 100644
--- a/R/glm_test.R
+++ b/R/glm_test.R
@@ -28,8 +28,8 @@ 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")) {
-    stop(paste("This is only implemented this for Poisson family with log link and Binomial family with logit link.\n",
+  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",
                "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/score_contribution.R b/R/score_contribution.R
index 7084a78..372f9e2 100644
--- a/R/score_contribution.R
+++ b/R/score_contribution.R
@@ -19,6 +19,13 @@ score_contribution <- function(i, model_fits, yy, xx, family, link) { ### output
     #score_vec <- matrix(model_fits[i]*(1 - model_fits[i])^(-1)*(yy[i] - model_fits[i]) * xx[i,], ncol = 1)
     score_vec <- matrix((yy[i] - model_fits[i]) * xx[i, ], ncol = 1)
   }
+  
+  if (family == "gaussian" & link == "identity") {
+    n <- length(yy)
+    p <- ncol(xx)
+    sigma2_tilde <- sum(yy - model_fits)/(n - p)
+    score_vec <- matrix((1/sigma2_tilde)*(yy[i] - model_fits[i])*xx[i,], ncol = 1)
+  }
 
   return(score_vec)
 }
\ No newline at end of file

From ffd8e21190faaefc3e92549c50f50f8bd23b75fb Mon Sep 17 00:00:00 2001
From: Shirley Mathur <shirmath556@gmail.com>
Date: Wed, 27 Aug 2025 15:15:57 -0700
Subject: [PATCH 02/17] made code for computing robust score test statistics
 for correlated settings modular

---
 R/D_matrix_contribution.R | 27 +++++++++++++++++++++++
 R/S_matrix_contribution.R | 27 +++++++++++++++++++++++
 R/V_matrix_contribution.R | 45 +++++++++++++++++++++++++++++++++++++++
 R/robust_score_test.R     | 20 ++++++-----------
 4 files changed, 106 insertions(+), 13 deletions(-)
 create mode 100644 R/D_matrix_contribution.R
 create mode 100644 R/S_matrix_contribution.R
 create mode 100644 R/V_matrix_contribution.R

diff --git a/R/D_matrix_contribution.R b/R/D_matrix_contribution.R
new file mode 100644
index 0000000..34ab3db
--- /dev/null
+++ b/R/D_matrix_contribution.R
@@ -0,0 +1,27 @@
+#' Compute D matrix contribution for a given cluster to robust score statistic for glm robust score tests.
+#'
+#' @param indices Indices of observations in cluster.
+#' @param model_fits The fitted glm under the null hypothesis.
+#' @param yy Vector of observed responses.
+#' @param xx Design matrix for model.
+#' @param family The model family for the fitted glm.
+#' @param link The link function utilized in the fitted glm.
+#' 
+
+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,])
+  }
+  
+  if (family == "binomial" & link == "logit") {
+    D_i <- t(model_fits[indices]*(1-model_fits[indices])*xx[indices,])
+  }
+  
+  if (family == "gaussian" & link == "identity") {
+    D_i <- t(xx[indices,])
+  }
+  
+  return (D_i)
+}
diff --git a/R/S_matrix_contribution.R b/R/S_matrix_contribution.R
new file mode 100644
index 0000000..fe31518
--- /dev/null
+++ b/R/S_matrix_contribution.R
@@ -0,0 +1,27 @@
+#' Compute S matrix contribution for a given cluster to robust score statistic for glm robust score tests.
+#'
+#' @param indices Indices of observations in cluster.
+#' @param model_fits The fitted glm under the null hypothesis.
+#' @param yy Vector of observed responses.
+#' @param xx Design matrix for model.
+#' @param family The model family for the fitted glm.
+#' @param link The link function utilized in the fitted glm.
+#' 
+
+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)
+  }
+  
+  return (S_i)
+}
diff --git a/R/V_matrix_contribution.R b/R/V_matrix_contribution.R
new file mode 100644
index 0000000..a78a8d8
--- /dev/null
+++ b/R/V_matrix_contribution.R
@@ -0,0 +1,45 @@
+#' Compute V matrix contribution to robust score statistic for glm robust score tests.
+#'
+#' @param indices Indices of observations for cluster.
+#' @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 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_i <- matrix(NA, nrow = length(model_fits[indices]), ncol = length(model_fits[indices]))
+  n_i <- length(indices)
+  
+  if (family == "poisson" & link == "log") {
+    if (n_i > 1) {
+      V_i <- diag(sqrt(model_fits[indices])) %*% corr_mat %*% diag(sqrt(model_fits[indices]))
+    } else {
+      V_i <- model_fits[indices] * corr_mat
+    }
+  }
+  
+  if (family == "binomial" & link == "logit") {
+    if (n_i > 1) {
+      V_i <- diag(sqrt(model_fits[indices]*(1 - model_fits[indices]))) %*% corr_mat %*% diag(sqrt(model_fits[indices]*(1 - model_fits[indices])))
+    } else {
+      V_i <- model_fits[indices] * corr_mat
+    }
+  }
+  
+  if (family == "gaussian" & link == "identity") {
+    n <- length(indices)
+    p <- ncol(xx)
+    sigma2_tilde <- sum(yy - model_fits)/(n - p)
+    if (n_i > 1) {
+      V_i <- sqrt(diag(sigma2_tilde, n_i)) %*% corr_mat %*% sqrt(diag(sigma2_tilde, n_i))
+    } else {
+      V_i <- sigma2_tilde * corr_mat
+    }
+    
+  }
+  
+  return (V_i)
+}
diff --git a/R/robust_score_test.R b/R/robust_score_test.R
index 84a6d2d..899086a 100644
--- a/R/robust_score_test.R
+++ b/R/robust_score_test.R
@@ -67,25 +67,19 @@ robust_score_test <- function(glm_object, call_to_model, param = 1,
       xxi <- xx[indices, , drop = FALSE]
       model0_fits_i <- model0_fits[indices]
 
-      Di <- matrix(NA, nrow = pp, ncol = n_i)
-      for (j in 1:n_i) {
-        for (k in seq_len(pp)) {
-          Di[k, j] <- xxi[j, k] * model0_fits_i[j]
-        }
-      }
+      Di <- D_matrix_contribution(indices = indices, model_fits = model0_fits,
+                                  yy = yy, xx = xx, family = model1family, link = model1link)
 
       corr_matrix <- matrix(rep(glm_object$geese$alpha, n_i^2), nrow = n_i)
       diag(corr_matrix) <- 1
-      if (n_i > 1) {
-        Vi <- diag(sqrt(model0_fits_i)) %*% corr_matrix %*% diag(sqrt(model0_fits_i))
-      } else {
-        Vi <- sqrt(model0_fits_i) * corr_matrix * sqrt(model0_fits_i)
-      }
+      Vi <- V_matrix_contribution(indices = indices, model_fits = model0_fits, corr_mat = corr_matrix,
+                                  yy = yy, xx = xx, family = model1family, link = model1link)
 
-      Si <- yy[indices] - model0_fits_i
+      Si <- S_matrix_contribution(indices = indices, model_fits = model0_fits,
+                                  yy = yy, xx = xx, family = model1family, link = model1link)
 
       ## Recall solve(x1, x2) is the same as solve(x1) %*% x2
-      Umatrices[[ii]] <- Di %*% solve(Vi, matrix(Si, ncol = 1))
+      Umatrices[[ii]] <- Di %*% solve(Vi, Si)
       Amatrices[[ii]] <- Di %*% solve(Vi, t(Di))
 
     }

From c6a48f7fd0fa01c6b3b4d43a4d6cf87e3bc67ab7 Mon Sep 17 00:00:00 2001
From: Shirley Mathur <shirmath556@gmail.com>
Date: Wed, 27 Aug 2025 15:48:49 -0700
Subject: [PATCH 03/17] fixed issue with linear reg correlated setting
 implementation

---
 R/V_matrix_contribution.R | 4 ++--
 R/gee_test.R              | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/R/V_matrix_contribution.R b/R/V_matrix_contribution.R
index a78a8d8..6977d33 100644
--- a/R/V_matrix_contribution.R
+++ b/R/V_matrix_contribution.R
@@ -30,9 +30,9 @@ V_matrix_contribution <- function(indices, model_fits, yy, xx, corr_mat, family,
   }
   
   if (family == "gaussian" & link == "identity") {
-    n <- length(indices)
+    n <- length(yy)
     p <- ncol(xx)
-    sigma2_tilde <- sum(yy - model_fits)/(n - p)
+    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/gee_test.R b/R/gee_test.R
index 4625ede..9ba151c 100644
--- a/R/gee_test.R
+++ b/R/gee_test.R
@@ -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")) {
+  if ((gee_family != "poisson" | gee_link != "log") & (gee_family != "binomial" | gee_link != "logit") & (gee_family != "gaussian" | gee_link != "identity")) {
     stop(paste("This is only implemented this for Poisson family with log link and Binomial family with logit link.\n",
                "You requested", gee_family, "family and", gee_link, "link \n",
                "Please open a GitHub issue if you're interested in other families."))

From 59224ccd304bc14c16a1d57b3260c889bb271785 Mon Sep 17 00:00:00 2001
From: Shirley Mathur <shirmath556@gmail.com>
Date: Thu, 28 Aug 2025 13:04:39 -0700
Subject: [PATCH 04/17] updated documentation for new functions

---
 man/D_matrix_contribution.Rd | 24 ++++++++++++++++++++++++
 man/S_matrix_contribution.Rd | 24 ++++++++++++++++++++++++
 man/V_matrix_contribution.Rd | 26 ++++++++++++++++++++++++++
 3 files changed, 74 insertions(+)
 create mode 100644 man/D_matrix_contribution.Rd
 create mode 100644 man/S_matrix_contribution.Rd
 create mode 100644 man/V_matrix_contribution.Rd

diff --git a/man/D_matrix_contribution.Rd b/man/D_matrix_contribution.Rd
new file mode 100644
index 0000000..a911961
--- /dev/null
+++ b/man/D_matrix_contribution.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/D_matrix_contribution.R
+\name{D_matrix_contribution}
+\alias{D_matrix_contribution}
+\title{Compute D matrix contribution for a given cluster to robust score statistic for glm robust score tests.}
+\usage{
+D_matrix_contribution(indices, model_fits, yy, xx, family, link)
+}
+\arguments{
+\item{indices}{Indices of observations in cluster.}
+
+\item{model_fits}{The fitted glm under the null hypothesis.}
+
+\item{yy}{Vector of observed responses.}
+
+\item{xx}{Design matrix for model.}
+
+\item{family}{The model family for the fitted glm.}
+
+\item{link}{The link function utilized in the fitted glm.}
+}
+\description{
+Compute D matrix contribution for a given cluster to robust score statistic for glm robust score tests.
+}
diff --git a/man/S_matrix_contribution.Rd b/man/S_matrix_contribution.Rd
new file mode 100644
index 0000000..4955a96
--- /dev/null
+++ b/man/S_matrix_contribution.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/S_matrix_contribution.R
+\name{S_matrix_contribution}
+\alias{S_matrix_contribution}
+\title{Compute S matrix contribution for a given cluster to robust score statistic for glm robust score tests.}
+\usage{
+S_matrix_contribution(indices, model_fits, yy, xx, family, link)
+}
+\arguments{
+\item{indices}{Indices of observations in cluster.}
+
+\item{model_fits}{The fitted glm under the null hypothesis.}
+
+\item{yy}{Vector of observed responses.}
+
+\item{xx}{Design matrix for model.}
+
+\item{family}{The model family for the fitted glm.}
+
+\item{link}{The link function utilized in the fitted glm.}
+}
+\description{
+Compute S matrix contribution for a given cluster to robust score statistic for glm robust score tests.
+}
diff --git a/man/V_matrix_contribution.Rd b/man/V_matrix_contribution.Rd
new file mode 100644
index 0000000..d46754a
--- /dev/null
+++ b/man/V_matrix_contribution.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/V_matrix_contribution.R
+\name{V_matrix_contribution}
+\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)
+}
+\arguments{
+\item{indices}{Indices of observations for cluster.}
+
+\item{model_fits}{The fitted glm under the null hypothesis for the given indices.}
+
+\item{yy}{Vector of observed responses.}
+
+\item{xx}{Design matrix for model.}
+
+\item{corr_mat}{Working correlation matrix for model fit.}
+
+\item{family}{The model family for the fitted glm.}
+
+\item{link}{The link function utilized in the fitted glm.}
+}
+\description{
+Compute V matrix contribution to robust score statistic for glm robust score tests.
+}

From a52b56a2625430017d3b91a3846ec01a0801b0a8 Mon Sep 17 00:00:00 2001
From: Shirley Mathur <shirmath556@gmail.com>
Date: Thu, 28 Aug 2025 15:53:14 -0700
Subject: [PATCH 05/17] updated error message to reflect linear regression as
 being implemented

---
 R/gee_test.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/gee_test.R b/R/gee_test.R
index 9ba151c..f713ded 100644
--- a/R/gee_test.R
+++ b/R/gee_test.R
@@ -93,7 +93,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")) {
-    stop(paste("This is only implemented this for Poisson family with log link and Binomial family with logit link.\n",
+    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."))
   }

From f714f88fa64b1d69e673549deb329738f4c43f31 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Mon, 1 Sep 2025 13:31:55 -0700
Subject: [PATCH 06/17] add tests of logistic and gaussian glms, update
 description to indicate minor addition (adding gaussian model option)

---
 DESCRIPTION                 |  2 +-
 tests/testthat/test-basic.R | 16 +++++++++++++++-
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 173aa41..4028635 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: raoBust
 Title: Robust score tests for generalized linear models
-Version: 1.1.1.0
+Version: 1.2.1.0
 Authors@R: 
     c(person("Amy", "Willis", , "adwillis@uw.edu", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-2802-4317")), 
diff --git a/tests/testthat/test-basic.R b/tests/testthat/test-basic.R
index bd4a3ba..c6dcb02 100644
--- a/tests/testthat/test-basic.R
+++ b/tests/testthat/test-basic.R
@@ -1,4 +1,4 @@
-test_that("offsets work", {
+test_that("test Poisson example", {
 
   expect_type(glm_test(dist ~ speed, data = cars, family=poisson(link="log"))$coef_tab,
               "double")
@@ -43,3 +43,17 @@ test_that("offsets work", {
 
 })
 
+test_that("test logistic example", {
+  
+  expect_type(glm_test((dist > 43) ~ speed, data = cars, family=binomial(link = "logit"))$coef_tab,
+              "double")
+  
+})
+
+test_that("test gaussian example", {
+  
+  expect_type(glm_test(dist ~ speed, data = cars, family=gaussian(link = "identity"))$coef_tab,
+              "double")
+  
+})
+

From fc7f85ba81f5c4b8567c4b538a7b19a265ed8ad9 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Mon, 1 Sep 2025 13:47:41 -0700
Subject: [PATCH 07/17] add codecov token

---
 .github/workflows/test-coverage.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index a3431f0..df45a95 100644
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -42,5 +42,5 @@ jobs:
         shell: Rscript {0}
 
       - name: Test coverage
-        run: covr::codecov()
+        run: covr::codecov(token = Sys.getenv("CODECOV_TOKEN"))
         shell: Rscript {0}

From 6362b4dc735fb56ba9264241abbca237dd5263f3 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Mon, 1 Sep 2025 13:52:20 -0700
Subject: [PATCH 08/17] try again to add token

---
 .github/workflows/test-coverage.yaml | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index df45a95..f710821 100644
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -42,5 +42,6 @@ jobs:
         shell: Rscript {0}
 
       - name: Test coverage
-        run: covr::codecov(token = Sys.getenv("CODECOV_TOKEN"))
-        shell: Rscript {0}
+        run: Rscript -e 'covr::codecov(token = Sys.getenv("CODECOV_TOKEN"))'
+        env:
+          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}

From 256d8269c96eb41e431b3fef54b1cab672b47469 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Mon, 1 Sep 2025 13:55:59 -0700
Subject: [PATCH 09/17] try one more time with codecov token

---
 .github/workflows/test-coverage.yaml | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index f710821..f982d5f 100644
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -42,6 +42,5 @@ jobs:
         shell: Rscript {0}
 
       - name: Test coverage
-        run: Rscript -e 'covr::codecov(token = Sys.getenv("CODECOV_TOKEN"))'
-        env:
-          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
+        run: covr::codecov(token = "f2ee7950-3d1e-4553-9c61-4b23a73ed655")
+        shell: Rscript {0}

From dd67f3752978db59eba54e9c268cb652bf4321d2 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Mon, 1 Sep 2025 14:16:18 -0700
Subject: [PATCH 10/17] tests for gee with logistic and gaussian models

---
 tests/testthat/test-basic.R | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/tests/testthat/test-basic.R b/tests/testthat/test-basic.R
index c6dcb02..fcb5f21 100644
--- a/tests/testthat/test-basic.R
+++ b/tests/testthat/test-basic.R
@@ -48,6 +48,11 @@ test_that("test logistic example", {
   expect_type(glm_test((dist > 43) ~ speed, data = cars, family=binomial(link = "logit"))$coef_tab,
               "double")
   
+  cars$id <- rep(1:5, 10)
+  expect_type(gee_test(formula = (dist > 43) ~ speed, data = cars, id = "id", 
+                       family=gaussian(link = "identity"))$coef_tab,
+              "list")
+  
 })
 
 test_that("test gaussian example", {
@@ -55,5 +60,10 @@ test_that("test gaussian example", {
   expect_type(glm_test(dist ~ speed, data = cars, family=gaussian(link = "identity"))$coef_tab,
               "double")
   
+  cars$id <- rep(1:5, 10)
+  expect_type(gee_test(formula = dist ~ speed, data = cars, id = "id", 
+                       family=gaussian(link = "identity"))$coef_tab,
+              "list")
+  
 })
 

From 019c48e61e9bfc3a67883f0287a91b255d8f962a Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Mon, 1 Sep 2025 14:17:32 -0700
Subject: [PATCH 11/17] check error when not included family and link combo

---
 tests/testthat/test-basic.R | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/tests/testthat/test-basic.R b/tests/testthat/test-basic.R
index fcb5f21..233d391 100644
--- a/tests/testthat/test-basic.R
+++ b/tests/testthat/test-basic.R
@@ -67,3 +67,9 @@ test_that("test gaussian example", {
   
 })
 
+test_that("error with other model", {
+  
+  expect_error(glm_test(dist ~ speed, data = cars, family=gaussian(link = "log")))
+
+})
+

From 24bdca465772833968228a386327c4bd4ba767fb Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Mon, 1 Sep 2025 17:14:03 -0700
Subject: [PATCH 12/17] simplify S, adjust D to work for clusters with size 1,
 increase code coverage for V, S, D matrix functions

---
 R/D_matrix_contribution.R   |  7 +++----
 R/S_matrix_contribution.R   | 13 +------------
 tests/testthat/test-basic.R |  9 +++++++--
 3 files changed, 11 insertions(+), 18 deletions(-)

diff --git a/R/D_matrix_contribution.R b/R/D_matrix_contribution.R
index 34ab3db..e3ba4b6 100644
--- a/R/D_matrix_contribution.R
+++ b/R/D_matrix_contribution.R
@@ -9,18 +9,17 @@
 #' 
 
 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])
   }
   
   return (D_i)
diff --git a/R/S_matrix_contribution.R b/R/S_matrix_contribution.R
index fe31518..26ed302 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/tests/testthat/test-basic.R b/tests/testthat/test-basic.R
index 233d391..b921015 100644
--- a/tests/testthat/test-basic.R
+++ b/tests/testthat/test-basic.R
@@ -40,6 +40,11 @@ test_that("test Poisson example", {
                        family=poisson(link="log"),
                        offset = log(rnorm(nrow(cars), 20)))$coef_tab,
               "double")
+  
+  cars$id <- c(rep(1:5, 9), 6:10)
+  expect_type(gee_test(formula = dist ~ speed, data = cars, id = "id", 
+                       family=poisson(link="log"))$coef_tab,
+              "list")
 
 })
 
@@ -48,7 +53,7 @@ test_that("test logistic example", {
   expect_type(glm_test((dist > 43) ~ speed, data = cars, family=binomial(link = "logit"))$coef_tab,
               "double")
   
-  cars$id <- rep(1:5, 10)
+  cars$id <- c(rep(1:5, 9), 6:10)
   expect_type(gee_test(formula = (dist > 43) ~ speed, data = cars, id = "id", 
                        family=gaussian(link = "identity"))$coef_tab,
               "list")
@@ -60,7 +65,7 @@ test_that("test gaussian example", {
   expect_type(glm_test(dist ~ speed, data = cars, family=gaussian(link = "identity"))$coef_tab,
               "double")
   
-  cars$id <- rep(1:5, 10)
+  cars$id <- c(rep(1:5, 9), 6:10)
   expect_type(gee_test(formula = dist ~ speed, data = cars, id = "id", 
                        family=gaussian(link = "identity"))$coef_tab,
               "list")

From f776f5d846bc72151d7854caf30ba5d1e29d5f09 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Wed, 10 Sep 2025 13:50:14 -0700
Subject: [PATCH 13/17] updating variance estimation for gaussian model, n-m
 instead of n-1

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

diff --git a/R/V_matrix_contribution.R b/R/V_matrix_contribution.R
index 6977d33..6105ead 100644
--- a/R/V_matrix_contribution.R
+++ b/R/V_matrix_contribution.R
@@ -9,7 +9,7 @@
 #' @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, pp0 = 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 +31,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)
+    m <- ncol(xx) - pp0
+    sigma2_tilde <- sum((yy - model_fits)^2)/(n - m)
     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..0a1c69d 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, pp0 = pp0, family = model1family, link = model1link)
 
       Si <- S_matrix_contribution(indices = indices, model_fits = model0_fits,
                                   yy = yy, xx = xx, family = model1family, link = model1link)

From 0e0804f54870181cd12bd9da8269e06c0e3f35d1 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Wed, 10 Sep 2025 13:51:05 -0700
Subject: [PATCH 14/17] updating documentation accordingly

---
 man/V_matrix_contribution.Rd | 11 ++++++++++-
 man/gee_test.Rd              |  2 +-
 2 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/man/V_matrix_contribution.Rd b/man/V_matrix_contribution.Rd
index d46754a..b67a888 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,
+  pp0 = 1,
+  corr_mat,
+  family,
+  link
+)
 }
 \arguments{
 \item{indices}{Indices of observations for cluster.}
diff --git a/man/gee_test.Rd b/man/gee_test.Rd
index ebf3487..dfb5812 100644
--- a/man/gee_test.Rd
+++ b/man/gee_test.Rd
@@ -22,7 +22,7 @@ instead be performed with a GLM.}
 
 \item{skip_gee}{When TRUE doesn't try to optimize with a GEE (just uses a GLM). This should only be used internally for testing.}
 
-\item{...}{Arguments that you would pass to a regular \code{geepack::geeglm} call. Any observations with \code{NA} values in the data (response or covariates) will be dropped.}
+\item{...}{Arguments that you would pass to a regular \code{geepack::geeglm} call. Any observations with \code{NA} values in the data (response or covariates or id) will be dropped.}
 }
 \description{
 Generalized Estimating Equations under technical replication with robust and non-robust Wald and Rao (score) tests

From ab5103acecd8896e8cbd7c81911d7c9fcda80f4d Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Wed, 10 Sep 2025 14:01:53 -0700
Subject: [PATCH 15/17] fix relevant documentation

---
 R/V_matrix_contribution.R    | 1 +
 man/V_matrix_contribution.Rd | 2 ++
 2 files changed, 3 insertions(+)

diff --git a/R/V_matrix_contribution.R b/R/V_matrix_contribution.R
index 6105ead..29243e6 100644
--- a/R/V_matrix_contribution.R
+++ b/R/V_matrix_contribution.R
@@ -4,6 +4,7 @@
 #' @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 pp0 Number of fixed parameters 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.
diff --git a/man/V_matrix_contribution.Rd b/man/V_matrix_contribution.Rd
index b67a888..490e927 100644
--- a/man/V_matrix_contribution.Rd
+++ b/man/V_matrix_contribution.Rd
@@ -24,6 +24,8 @@ V_matrix_contribution(
 
 \item{xx}{Design matrix for model.}
 
+\item{pp0}{Number of fixed parameters under null hypothesis.}
+
 \item{corr_mat}{Working correlation matrix for model fit.}
 
 \item{family}{The model family for the fitted glm.}

From 16553781277cad65ab03e4b636b6cabce0c6d10c Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Fri, 12 Sep 2025 15:03:44 -0700
Subject: [PATCH 16/17] update denominator for sigma^2 for gaussian case in
 V_matrix_contribution, because it cancels out and is therefore not needed

---
 R/V_matrix_contribution.R   | 5 ++++-
 tests/testthat/test-basic.R | 4 +++-
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/R/V_matrix_contribution.R b/R/V_matrix_contribution.R
index 29243e6..781764c 100644
--- a/R/V_matrix_contribution.R
+++ b/R/V_matrix_contribution.R
@@ -33,7 +33,10 @@ V_matrix_contribution <- function(indices, model_fits, yy, xx, pp0 = 1, corr_mat
   if (family == "gaussian" & link == "identity") {
     n <- length(yy)
     m <- ncol(xx) - pp0
-    sigma2_tilde <- sum((yy - model_fits)^2)/(n - m)
+    # note, while sigma^2_tilde should in theory have a denominator, this would cancel out
+    # in the test statistic, so we set it to 1
+    # sigma2_tilde <- sum((yy - model_fits)^2)/(n - m)
+    sigma2_tilde <- sum((yy - model_fits)^2)
     if (n_i > 1) {
       V_i <- sqrt(diag(sigma2_tilde, n_i)) %*% corr_mat %*% sqrt(diag(sigma2_tilde, n_i))
     } else {
diff --git a/tests/testthat/test-basic.R b/tests/testthat/test-basic.R
index b921015..6900b15 100644
--- a/tests/testthat/test-basic.R
+++ b/tests/testthat/test-basic.R
@@ -66,7 +66,9 @@ test_that("test gaussian example", {
               "double")
   
   cars$id <- c(rep(1:5, 9), 6:10)
-  expect_type(gee_test(formula = dist ~ speed, data = cars, id = "id", 
+  cars$bin <- rep(0:1, each = 25)
+  cars$cat <- rep(c("A", "B", "C", "D", "E"), 10)
+  expect_type(gee_test(formula = dist ~ speed + bin + cat, data = cars, id = "id", 
                        family=gaussian(link = "identity"))$coef_tab,
               "list")
   

From a8b57319f4a8727db6275620717190c7274a2c98 Mon Sep 17 00:00:00 2001
From: Sarah Teichman <teichs@uw.edu>
Date: Mon, 27 Oct 2025 13:04:51 -0700
Subject: [PATCH 17/17] fix gee_test.Rd file

---
 man/gee_test.Rd | 2 --
 1 file changed, 2 deletions(-)

diff --git a/man/gee_test.Rd b/man/gee_test.Rd
index c565709..27c3d24 100644
--- a/man/gee_test.Rd
+++ b/man/gee_test.Rd
@@ -23,8 +23,6 @@ gee_test(
 instead be performed with a GLM.}
 
 \item{skip_gee}{When TRUE doesn't try to optimize with a GEE (just uses a GLM). This should only be used internally for testing.}
-
-\item{...}{Arguments that you would pass to a regular \code{geepack::geeglm} call. Any observations with \code{NA} values in the data (response or covariates or id) will be dropped.}
 }
 \description{
 Generalized Estimating Equations under technical replication with robust and non-robust Wald and Rao (score) tests