Small updates to PR #40

.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}

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")), 

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)
+}

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)
+}

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)
+}

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)
 }

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."))
   }

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."))
   }

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
 

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)
 }

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")))
 
 })