From 569d6ffe44150f1b0593ac3eadbfdcaa92a7a509 Mon Sep 17 00:00:00 2001 From: Gregory Chen <65814078+hoppanda@users.noreply.github.com> Date: Wed, 26 Jun 2024 02:51:43 +0200 Subject: [PATCH] add in robincar vcov, its report and the tests (#16) * add in all draft functions * [skip style] [skip vbump] Restyle files * clean up design folder for merging * Update branch 9 with main (#23) * update data and fix design (#17) * update data and fix design * update treatment to factor * fix styler and lint * 7 counter factual prediction (#14) * add function prototype * add basic tests on functionality * add readme * add lm support * update tests * update prediction * simplify * fix typo * fix check issues * update vignettes skeleton * use dummy data in test * update print of cf predictions * fix checks * fix check issues * update bias args * [skip style] [skip vbump] Restyle files * add treatment to counter factual prediction * fix check * 18 sandwich variance (#19) * add vcov * fix check issues * correct issue --------- Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> * 20 add another developers design note (#22) * add in all draft functions * [skip style] [skip vbump] Restyle files * clean up design folder for merging * [skip style] [skip vbump] Restyle files * Update .gitignore --------- Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Liming <36079400+clarkliming@users.noreply.github.com> --------- Co-authored-by: Liming <36079400+clarkliming@users.noreply.github.com> Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> * add robincar vcov calc functions and report function * adding all the tests * [skip roxygen] [skip vbump] Roxygen Man Pages Auto Update * [skip style] [skip vbump] Restyle files * add Rds * delete unneeded RdS * add anhecova covariance (#24) * add anhecova covariance * fix error * use new implementation of anhecova variance (#25) * replace original function with faster version * update anhecova variance and remove report functions --------- Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Liming <36079400+clarkliming@users.noreply.github.com> Co-authored-by: 27856297+dependabot-preview[bot]@users.noreply.github.com <27856297+dependabot-preview[bot]@users.noreply.github.com> --- .gitignore | 3 + NAMESPACE | 1 + R/variance_anhecova.R | 81 +++++++++++++++++++ R/{variance.R => variance_hc.R} | 6 +- .../design_note.R} | 0 .../developer_note.R} | 1 - man/h_adjust_pi.Rd | 18 +++++ man/h_get_erb.Rd | 23 ++++++ man/vcovANHECOVA.Rd | 23 ++++++ tests/testthat/_snaps/variance.md | 30 +++++++ tests/testthat/test-variance.R | 16 ++++ 11 files changed, 199 insertions(+), 3 deletions(-) create mode 100644 R/variance_anhecova.R rename R/{variance.R => variance_hc.R} (71%) rename design/{package_structure/design_note_xgc_v1.R => code_replicate/design_note.R} (100%) rename design/{package_structure/developer_note_xgc_v1.R => code_replicate/developer_note.R} (99%) create mode 100644 man/h_adjust_pi.Rd create mode 100644 man/h_get_erb.Rd create mode 100644 man/vcovANHECOVA.Rd create mode 100644 tests/testthat/_snaps/variance.md create mode 100644 tests/testthat/test-variance.R diff --git a/.gitignore b/.gitignore index 7065c0b..d381fc6 100644 --- a/.gitignore +++ b/.gitignore @@ -31,3 +31,6 @@ coverage.* node_modules/ package-lock.json package.json +RobinCar2.Rcheck/ +RobinCar2*.tar.gz +RobinCar2*.tgz diff --git a/NAMESPACE b/NAMESPACE index 55f3734..2fce319 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/variance_anhecova.R b/R/variance_anhecova.R new file mode 100644 index 0000000..670f685 --- /dev/null +++ b/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 +} diff --git a/R/variance.R b/R/variance_hc.R similarity index 71% rename from R/variance.R rename to R/variance_hc.R index be5f633..7f0be20 100644 --- a/R/variance.R +++ b/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) diff --git a/design/package_structure/design_note_xgc_v1.R b/design/code_replicate/design_note.R similarity index 100% rename from design/package_structure/design_note_xgc_v1.R rename to design/code_replicate/design_note.R diff --git a/design/package_structure/developer_note_xgc_v1.R b/design/code_replicate/developer_note.R similarity index 99% rename from design/package_structure/developer_note_xgc_v1.R rename to design/code_replicate/developer_note.R index 697c938..f2db388 100644 --- a/design/package_structure/developer_note_xgc_v1.R +++ b/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) diff --git a/man/h_adjust_pi.Rd b/man/h_adjust_pi.Rd new file mode 100644 index 0000000..c32e099 --- /dev/null +++ b/man/h_adjust_pi.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance_anhecova.R +\name{h_adjust_pi} +\alias{h_adjust_pi} +\title{Obtain Adjustment for Proportion of Treatment Assignment} +\usage{ +h_adjust_pi(pi) +} +\arguments{ +\item{pi}{(\code{numeric}) vector of proportions.} +} +\value{ +Numeric matrix. +} +\description{ +Obtain Adjustment for Proportion of Treatment Assignment +} +\keyword{internal} diff --git a/man/h_get_erb.Rd b/man/h_get_erb.Rd new file mode 100644 index 0000000..7714c0b --- /dev/null +++ b/man/h_get_erb.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance_anhecova.R +\name{h_get_erb} +\alias{h_get_erb} +\title{Obtain Adjustment for Covariance Matrix} +\usage{ +h_get_erb(resi, group_idx, trt, pi, randomization) +} +\arguments{ +\item{resi}{(\code{numeric}) vector of residuals.} + +\item{group_idx}{(\code{list} of \code{integer}) index for each groups.} + +\item{trt}{(\code{factor}) of treatment assignment.} + +\item{pi}{(\code{numeric}) proportion of treatment assignment.} + +\item{randomization}{(\code{string}) name of the randomization schema.} +} +\description{ +Obtain Adjustment for Covariance Matrix +} +\keyword{internal} diff --git a/man/vcovANHECOVA.Rd b/man/vcovANHECOVA.Rd new file mode 100644 index 0000000..37b859d --- /dev/null +++ b/man/vcovANHECOVA.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variance_anhecova.R +\name{vcovANHECOVA} +\alias{vcovANHECOVA} +\title{ANHECOVA Covariance} +\usage{ +vcovANHECOVA(x, decompose = TRUE, randomization = "simple", ...) +} +\arguments{ +\item{x}{(\code{prediction_cf}) Counter-factual prediction.} + +\item{decompose}{(\code{flag}) whether to use decompose method to calculate the variance.} + +\item{randomization}{(\code{string}) randomization method.} + +\item{...}{Not used.} +} +\value{ +Named covariance matrix. +} +\description{ +ANHECOVA Covariance +} diff --git a/tests/testthat/_snaps/variance.md b/tests/testthat/_snaps/variance.md new file mode 100644 index 0000000..5c5c8c3 --- /dev/null +++ b/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 + diff --git a/tests/testthat/test-variance.R b/tests/testthat/test-variance.R new file mode 100644 index 0000000..217f294 --- /dev/null +++ b/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") + ) +})