Adds CR3 estimator to lm_robust and iv_robust

DESCRIPTION

@@ -20,7 +20,7 @@ Imports:
     Formula,
     rlang (>= 0.2.0)
 LinkingTo: Rcpp, RcppEigen
-RoxygenNote: 6.0.1
+RoxygenNote: 6.1.0
 LazyData: true
 Suggests:
     RcppEigen,

R/estimatr_iv_robust.R

@@ -22,7 +22,7 @@
 #' @param se_type The sort of standard error sought. If `clusters` is
 #' not specified the options are "HC0", "HC1" (or "stata", the equivalent),
 #'  "HC2" (default), "HC3", or
-#' "classical". If `clusters` is specified the options are "CR0", "CR2" (default), or "stata". Can also specify "none", which may speed up estimation of the coefficients.
+#' "classical". If `clusters` is specified the options are "CR0", "stata", "CR2" (default), or "CR3". Can also specify "none", which may speed up estimation of the coefficients.
 #' @param ci logical. Whether to compute and return p-values and confidence
 #' intervals, TRUE by default.
 #' @param alpha The significance level, 0.05 by default.

R/estimatr_lm_robust.R

@@ -19,7 +19,8 @@
 #' @param se_type The sort of standard error sought. If `clusters` is
 #' not specified the options are "HC0", "HC1" (or "stata", the equivalent),
 #'  "HC2" (default), "HC3", or
-#' "classical". If `clusters` is specified the options are "CR0", "CR2" (default), or "stata". Can also specify "none", which may speed up estimation of the coefficients.
+#' "classical". If `clusters` is specified the options are "CR0", "stata",
+#' "CR2" (default), or "CR3". Can also specify "none", which may speed up estimation of the coefficients.
 #' @param ci logical. Whether to compute and return p-values and confidence
 #' intervals, TRUE by default.
 #' @param alpha The significance level, 0.05 by default.

R/helper_lm_robust_fit.R

@@ -72,9 +72,9 @@ lm_robust_fit <- function(y,
 
   se_type <- check_se_type(se_type, clustered)
 
-  if (weighted && se_type == "CR2" && fes) {
+  if (weighted && se_type %in% c("CR2", "CR3") && fes) {
     stop(
-      "Cannot use `fixed_effects` with weighted CR2 estimation at the moment. ",
+      "Cannot use `fixed_effects` with weighted CR2 or CR3 estimation at the moment. ",
       "Try setting `se_type` = \"stata\""
     )
   }
@@ -177,7 +177,7 @@ lm_robust_fit <- function(y,
 
       # For CR2 need X weighted by weights again
       # so that instead of having X * sqrt(W) we have X * W
-      if (se_type == "CR2") {
+      if (se_type %in% c("CR2", "CR3")) {
         data[["X"]] <- data[["weights"]] * data[["X"]]
         if (fes) {
           data[["femat"]] <- data[["weights"]] * data[["femat"]]
@@ -186,8 +186,6 @@ lm_robust_fit <- function(y,
 
     }
 
-
-
     # Also need second stage residuals for fstat
     if (iv_stage[[1]] == 2) {
       fit_vals[["fitted.values.iv"]] <- as.matrix(data[["X"]] %*% fit$beta_hat)
@@ -200,14 +198,14 @@ lm_robust_fit <- function(y,
     if (se_type != "none") {
 
       vcov_fit <- lm_variance(
-        X = if (se_type %in% c("HC2", "HC3", "CR2") && fes)
+        X = if (se_type %in% c("HC2", "HC3", "CR2", "CR3") && fes)
           cbind(data[["X"]], data[["femat"]])
           else data[["X"]],
-        Xunweighted = if (se_type %in% c("HC2", "HC3", "CR2") && fes && weighted)
+        Xunweighted = if (se_type %in% c("HC2", "HC3", "CR2", "CR3") && fes && weighted)
           cbind(data[["Xunweighted"]], data[["fematunweighted"]])
           else data[["Xunweighted"]],
         XtX_inv = fit$XtX_inv,
-        ei = if (se_type == "CR2" && weighted)
+        ei = if (se_type %in% c("CR2", "CR3") && weighted)
           fit_vals[["ei.unweighted"]]
           else fit_vals[["ei"]],
         weight_mean = data[["weight_mean"]],
@@ -388,7 +386,7 @@ lm_robust_fit <- function(y,
 check_se_type <- function(se_type, clustered) {
 
   # Allowable se_types with clustering
-  cl_se_types <- c("CR0", "CR2", "stata")
+  cl_se_types <- c("CR0", "CR2", "CR3", "CR1s", "stata")
   rob_se_types <- c("HC0", "HC1", "HC2", "HC3", "classical", "stata")
 
   # Parse cluster variable

R/helper_starprep.R

@@ -93,7 +93,7 @@ commarobust <- function(model,
     XtX_inv <- data[["weight_mean"]] * XtX_inv
 
     # Need unweighted resid and need to reweight X
-    if (se_type == "CR2") {
+    if (se_type %in% c("CR2", "CR3")) {
       eiunweighted <-
         as.matrix(data[["yunweighted"]] - data[["Xunweighted"]] %*% coefs)
       data[["X"]] <- data[["weights"]] * data[["X"]]
@@ -104,7 +104,7 @@ commarobust <- function(model,
     X = data[["X"]],
     Xunweighted = data[["Xunweighted"]],
     XtX_inv = XtX_inv,
-    ei = if (se_type == "CR2" && weighted) eiunweighted else ei,
+    ei = if (se_type %in% c("CR2", "CR3") && weighted) eiunweighted else ei,
     weight_mean = data[["weight_mean"]],
     cluster = data[["cluster"]],
     J = data[["J"]],
@@ -114,7 +114,6 @@ commarobust <- function(model,
     fe_rank = 0
   )
 
-
   ## Build return_list
   return_list <- list(
     coefficients = as.matrix(coef(model)),

src/lm_robust_helper.cpp

@@ -247,7 +247,7 @@ List lm_variance(Eigen::Map<Eigen::MatrixXd>& X,
   const int n(X.rows()), r(XtX_inv.cols()), ny(ei.cols());
   // Rcout << "X:" << std::endl << X << std::endl;
   int r_fe = r + fe_rank;
-  const bool clustered = ((se_type == "stata") || (se_type == "CR0") || (se_type == "CR2"));
+  const bool clustered = ((se_type == "stata") || (se_type == "CR0") || (se_type == "CR2") || (se_type == "CR3"));
   const int npars = r * ny;
   int sandwich_size = n;
   if (clustered) {
@@ -291,7 +291,7 @@ List lm_variance(Eigen::Map<Eigen::MatrixXd>& X,
     }
 
     Eigen::MatrixXd meatXtX_inv;
-    if ((se_type == "HC2") || (se_type == "HC3") || (se_type == "CR2")) {
+    if ((se_type == "HC2") || (se_type == "HC3") || (se_type == "CR2") || (se_type == "CR3")) {
       if (X.cols() > r) {
         meatXtX_inv = getMeatXtX(X, XtX_inv);
         r_fe = meatXtX_inv.cols();
@@ -324,7 +324,6 @@ List lm_variance(Eigen::Map<Eigen::MatrixXd>& X,
       }
       // Rcout << "temp_omega:" << std::endl << temp_omega << std::endl;
 
-
       for (int m = 0; m < ny; m++) {
         if (ny > 1) {
           // Preserve signs for off-diagonal vcov blocks in mlm
@@ -337,25 +336,28 @@ List lm_variance(Eigen::Map<Eigen::MatrixXd>& X,
     } else {
       // clustered
 
-      if (se_type == "CR2") {
+      if ((se_type == "CR2") || (se_type == "CR3")) {
         Xoriginal.resize(n, r);
         if (Xunweighted.isNotNull()) {
           Xoriginal = Rcpp::as<Eigen::Map<Eigen::MatrixXd> >(Xunweighted);
         } else {
           Xoriginal = X;
         }
 
-        H1s.resize(r_fe, r_fe*J);
-        H2s.resize(r_fe, r_fe*J);
-        H3s.resize(r_fe, r_fe*J);
-        P_diags.resize(r_fe, J);
+        if (se_type == "CR2") {
+          H1s.resize(r_fe, r_fe*J);
+          H2s.resize(r_fe, r_fe*J);
+          H3s.resize(r_fe, r_fe*J);
+          P_diags.resize(r_fe, J);
 
-        M_U_ct = meatXtX_inv.llt().matrixL();
-        MUWTWUM = meatXtX_inv * X.leftCols(r_fe).transpose() * X.leftCols(r_fe) * meatXtX_inv;
-        Omega_ct = MUWTWUM.llt().matrixL();
+          M_U_ct = meatXtX_inv.llt().matrixL();
+          MUWTWUM = meatXtX_inv * X.leftCols(r_fe).transpose() * X.leftCols(r_fe) * meatXtX_inv;
+          Omega_ct = MUWTWUM.llt().matrixL();
+        }
       }
 
       Eigen::Map<Eigen::ArrayXi> clusters = Rcpp::as<Eigen::Map<Eigen::ArrayXi> >(cluster);
+      // Rcout << "clusters: "<< std::endl << clusters << std::endl;
 
       double current_cluster = clusters(0);
       int clust_num = 0;
@@ -426,16 +428,34 @@ List lm_variance(Eigen::Map<Eigen::MatrixXd>& X,
               H2s.block(0, p_pos, r_fe, r_fe) = ME * X.block(start_pos, 0, len, r_fe) * M_U_ct;
               H3s.block(0, p_pos, r_fe, r_fe) = MEU * Omega_ct;
             }
+
+          } else if (se_type == "CR3") {
+
+            Eigen::CompleteOrthogonalDecomposition<Eigen::MatrixXd> At(
+                Eigen::MatrixXd::Identity(len, len) -
+                  Xoriginal.block(start_pos, 0, len, r_fe) *
+                  meatXtX_inv *
+                  X.block(start_pos, 0, len, r_fe).transpose()
+            );
+
+            if (Xunweighted.isNotNull()) {
+              // We require the transpose in both cases, but in the unweighted
+              // case `A` is symmetric and is not needed to be transposed first
+              // Solve against the identity matrix and then transpose in order
+              // to allow the more common unweighted case to be simpler
+              At_WX_inv =  At.solve(Eigen::MatrixXd::Identity(len, len)).transpose() * X.block(start_pos, 0, len, r_fe);
+            } else {
+              At_WX_inv =  At.solve(X.block(start_pos, 0, len, r_fe));
+            }
+
           }
 
           if (ny > 1) {
-
             // Stack residuals for this cluster from each model
-            // Rcout << "len: " << len << std::endl;
             Eigen::MatrixXd ei_block = ei.block(start_pos, 0, len, ny);
             Eigen::Map<const Eigen::MatrixXd> ei_long(ei_block.data(), 1, len*ny);
 
-            if (se_type == "CR2") {
+            if ((se_type == "CR2") || (se_type == "CR3")) {
               half_meat.block(clust_num, 0, 1, npars) =
                 ei_long *
                 Kr(Eigen::MatrixXd::Identity(ny, ny), At_WX_inv.leftCols(r));
@@ -447,7 +467,7 @@ List lm_variance(Eigen::Map<Eigen::MatrixXd>& X,
 
           } else {
 
-            if (se_type == "CR2") {
+            if ((se_type == "CR2") || (se_type == "CR3")) {
               half_meat.row(clust_num) =
                 ei.block(start_pos, 0, len, 1).transpose() *
                 At_WX_inv.leftCols(r);
@@ -456,7 +476,6 @@ List lm_variance(Eigen::Map<Eigen::MatrixXd>& X,
                 ei.block(start_pos, 0, len, 1).transpose() *
                 X.block(start_pos, 0, len, r);
             }
-
           }
           if (i < n) {
             current_cluster = clusters(i);

tests/testthat/helper-return-cleaners.R

@@ -0,0 +1,13 @@
+# This fn removes calls from function returns to make testing easier
+rmcall <- function(obj) {
+  if (!is.null(obj[["call"]])) {
+    obj[["call"]] <- NULL
+  }
+  return(obj)
+}
+
+# This fn casts tibbles as data.frames for equivalency tests
+# TODO implement this everywhere
+expect_equivalent_tbl <- function(obj1, obj2) {
+  expect_equivalent(as.data.frame(obj1), as.data.frame(obj2))
+}

tests/testthat/helper-se-types.R

@@ -0,0 +1,4 @@
+# se_types for various files
+
+se_types <- c("classical", "HC0", "HC1", "HC2", "HC3")
+cr_se_types <- c("CR0", "stata", "CR2", "CR3")