add in robincar vcov, its report and the tests

.gitignore

@@ -31,3 +31,6 @@ coverage.*
 node_modules/
 package-lock.json
 package.json
+RobinCar2.Rcheck/
+RobinCar2*.tar.gz
+RobinCar2*.tgz

NAMESPACE

@@ -6,6 +6,7 @@ S3method(print,prediction_cf)
 S3method(vcovHC,prediction_cf)
 export(bias)
 export(predict_counterfactual)
+export(vcovANHECOVA)
 import(checkmate)
 importFrom(sandwich,vcovHC)
 importFrom(stats,coefficients)

R/variance_anhecova.R

@@ -0,0 +1,81 @@
+#' ANHECOVA Covariance
+#'
+#' @param x (`prediction_cf`) Counter-factual prediction.
+#' @param decompose (`flag`) whether to use decompose method to calculate the variance.
+#' @param randomization (`string`) randomization method.
+#' @param ... Not used.
+#'
+#' @return Named covariance matrix.
+#' @export
+vcovANHECOVA <- function(x, decompose = TRUE, randomization = "simple", ...) {
+  assert_class(x, "prediction_cf")
+  assert_flag(decompose)
+  assert_string(randomization)
+  resi <- attr(x, "residual")
+  est <- as.numeric(x)
+  preds <- attr(x, "predictions")
+  var_preds <- var(preds)
+  y <- attr(x, "response")
+  trt <- attr(x, "treatment")
+  pi_t <- as.numeric(table(trt) / length(trt))
+  trt_lvls <- levels(trt)
+  group_idx <- attr(x, "group_idx")
+
+  idx <- split(seq_len(length(trt)), trt)
+  cov_ymu <- vapply(idx, function(is) stats::cov(y[is], preds[is, ]), FUN.VALUE = rep(0, ncol(preds)))
+
+  if (decompose) {
+    vcov_sr <- (vapply(idx, function(is) stats::var(y[is]), FUN.VALUE = 0) + diag(var_preds) - 2 * diag(cov_ymu)) / pi_t
+  } else {
+    vcov_sr <- vapply(idx, function(is) stats::var(resi[is]), FUN.VALUE = 0) / pi_t
+  }
+
+  v <- diag(vcov_sr) + cov_ymu + t(cov_ymu) - var_preds
+  v <- v - h_get_erb(resi, group_idx, trt, pi_t, randomization)
+  ret <- v / length(resi)
+  dimnames(ret) <- list(trt_lvls, trt_lvls)
+  return(ret)
+}
+
+#' Obtain Adjustment for Proportion of Treatment Assignment
+#' @keywords internal
+#' @param pi (`numeric`) vector of proportions.
+#' @return Numeric matrix.
+h_adjust_pi <- function(pi) {
+  assert_numeric(pi)
+  diag(pi, nrow = length(pi)) - pi %*% t(pi)
+}
+
+#' Obtain Adjustment for Covariance Matrix
+#' @param resi (`numeric`) vector of residuals.
+#' @param group_idx (`list` of `integer`) index for each groups.
+#' @param trt (`factor`) of treatment assignment.
+#' @param pi (`numeric`) proportion of treatment assignment.
+#' @param randomization (`string`) name of the randomization schema.
+#' @keywords internal
+h_get_erb <- function(resi, group_idx, trt, pi, randomization) {
+  assert_list(group_idx, types = "integer")
+  if (length(group_idx) == 1L) {
+    return(0)
+  }
+  assert_string(randomization)
+  if (randomization %in% c("simple", "pocock-simon")) {
+    return(0)
+  }
+  assert_numeric(resi)
+  assert_factor(trt)
+  assert_numeric(pi)
+
+  omegaz_sr <- h_adjust_pi(pi)
+  n_tot <- length(resi)
+  resi_per_strata <- vapply(
+    group_idx,
+    function(ii) {
+      v <- vapply(split(resi[ii], trt[ii]), mean, FUN.VALUE = 0)
+      v * sqrt(length(ii) / n_tot)
+    },
+    FUN.VALUE = rep(0, length(levels(trt)))
+  )
+  rb_z <- resi_per_strata / as.numeric(pi)
+  tcrossprod(rb_z) * omegaz_sr
+}

R/variance.R → R/variance_hc.R

@@ -1,9 +1,11 @@
+#### Main functions ####
+
 #' @exportS3Method
-vcovHC.prediction_cf <- function(x, type = "HC3",...) {
+vcovHC.prediction_cf <- function(x, type = "HC3", ...) {
   fit <- attr(x, "fit")
   vc <- vcovHC(fit, type = type)
   mm <- attr(x, "model_matrix")
-  n <- length(md) / length(names(x))
+  n <- nrow(mm) / length(names(x))
   md <- family(fit)$mu.eta(attr(x, "predictions_linear")) / n
   z <- block_sum(as.vector(md) * mm, n)
   ret <- z %*% vc %*% t(z)

design/package_structure/developer_note_xgc_v1.R → design/code_replicate/developer_note.R

@@ -263,7 +263,6 @@ d <- sim_data(500)
 fit <- glm(y ~ trt:z1 + trt, family = binomial(link = "logit"), data = d)
 
 ## start of the pkg
-source(file.path(getwd(), "design", "package_structure", "all_draft_functions.R"))
 pred_counterfact <- get_countfact_pred(fit, "trt")
 fit.fvcov <- calculate_f_vcov(fit, "trt", pred_counterfact)
 report_fvcov(fit.fvcov, 3)

tests/testthat/_snaps/variance.md

@@ -0,0 +1,30 @@
+# vcovHC works
+
+    Code
+      vcovHC(pc)
+    Output
+                     pbo          trt1          trt2
+      pbo   1.106024e-03  4.523445e-07 -9.709004e-06
+      trt1  4.523445e-07  1.164889e-03 -7.709031e-07
+      trt2 -9.709004e-06 -7.709031e-07  1.170214e-03
+
+# vcovANHECOVA works
+
+    Code
+      vcovANHECOVA(pc)
+    Output
+                    pbo         trt1         trt2
+      pbo  1.128902e-03 1.856234e-05 1.333885e-05
+      trt1 1.856234e-05 1.184599e-03 2.178112e-05
+      trt2 1.333885e-05 2.178112e-05 1.157268e-03
+
+---
+
+    Code
+      vcovANHECOVA(pc, randomization = "permute_block")
+    Output
+                    pbo         trt1         trt2
+      pbo  1.128902e-03 1.856234e-05 1.333885e-05
+      trt1 1.856234e-05 1.184599e-03 2.178112e-05
+      trt2 1.333885e-05 2.178112e-05 1.157268e-03
+

tests/testthat/test-variance.R

@@ -0,0 +1,16 @@
+test_that("vcovHC works", {
+  pc <- predict_counterfactual(fit_binom, treatment ~ s1)
+  expect_snapshot(
+    vcovHC(pc)
+  )
+})
+
+test_that("vcovANHECOVA works", {
+  pc <- predict_counterfactual(fit_binom, treatment ~ s1)
+  expect_snapshot(
+    vcovANHECOVA(pc)
+  )
+  expect_snapshot(
+    vcovANHECOVA(pc, randomization = "permute_block")
+  )
+})