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