diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index a3431f0..f982d5f 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 = "f2ee7950-3d1e-4553-9c61-4b23a73ed655")
         shell: Rscript {0}
diff --git a/DESCRIPTION b/DESCRIPTION
index f8d3d93..48386af 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: raoBust
 Title: Robust score tests for generalized linear models
-Version: 1.1.4.0
+Version: 1.2.0.0
 Authors@R: 
     c(person("Amy", "Willis", , "adwillis@uw.edu", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-2802-4317")), 
diff --git a/R/D_matrix_contribution.R b/R/D_matrix_contribution.R
new file mode 100644
index 0000000..e3ba4b6
--- /dev/null
+++ b/R/D_matrix_contribution.R
@@ -0,0 +1,26 @@
+#' 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) {
+  
+  if (family == "poisson" & link == "log") {
+    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, , drop = FALSE])
+  }
+  
+  if (family == "gaussian" & link == "identity") {
+    D_i <- t(xx[indices, , drop = FALSE])
+  }
+  
+  return (D_i)
+}
diff --git a/R/S_matrix_contribution.R b/R/S_matrix_contribution.R
new file mode 100644
index 0000000..26ed302
--- /dev/null
+++ b/R/S_matrix_contribution.R
@@ -0,0 +1,16 @@
+#' 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(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..781764c
--- /dev/null
+++ b/R/V_matrix_contribution.R
@@ -0,0 +1,49 @@
+#' 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 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.
+#' 
+
+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)
+  
+  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(yy)
+    m <- ncol(xx) - pp0
+    # 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 {
+      V_i <- sigma2_tilde * corr_mat
+    }
+    
+  }
+  
+  return (V_i)
+}
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/gee_test.R b/R/gee_test.R
index 7eacb69..7bd7576 100644
--- a/R/gee_test.R
+++ b/R/gee_test.R
@@ -107,8 +107,8 @@ gee_test <- function(..., use_geeasy = TRUE, use_jack_se = FALSE, cluster_corr_c
 
   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")) {
-    stop(paste("This is only implemented this for Poisson family with log link and Binomial family with logit link.\n",
+  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, 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 197fa84..3343ebc 100644
--- a/R/glm_test.R
+++ b/R/glm_test.R
@@ -32,8 +32,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/robust_score_test.R b/R/robust_score_test.R
index 424a2b7..cb48bf7 100644
--- a/R/robust_score_test.R
+++ b/R/robust_score_test.R
@@ -67,22 +67,16 @@ 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, pp0 = pp0, 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
       
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
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..490e927
--- /dev/null
+++ b/man/V_matrix_contribution.Rd
@@ -0,0 +1,37 @@
+% 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,
+  pp0 = 1,
+  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{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.}
+
+\item{link}{The link function utilized in the fitted glm.}
+}
+\description{
+Compute V matrix contribution to robust score statistic for glm robust score tests.
+}
diff --git a/tests/testthat/test-basic.R b/tests/testthat/test-basic.R
index bd4a3ba..6900b15 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")
@@ -40,6 +40,43 @@ test_that("offsets work", {
                        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")
+
+})
+
+test_that("test logistic example", {
+  
+  expect_type(glm_test((dist > 43) ~ speed, data = cars, family=binomial(link = "logit"))$coef_tab,
+              "double")
+  
+  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")
+  
+})
+
+test_that("test gaussian example", {
+  
+  expect_type(glm_test(dist ~ speed, data = cars, family=gaussian(link = "identity"))$coef_tab,
+              "double")
+  
+  cars$id <- c(rep(1:5, 9), 6:10)
+  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")
+  
+})
+
+test_that("error with other model", {
+  
+  expect_error(glm_test(dist ~ speed, data = cars, family=gaussian(link = "log")))
 
 })