Merge pull request #101 from DoubleML/m-cleanup-tests

Refactoring, cleanup and extension of the unit test framework

Author: MalteKurz

R/double_ml.R

@@ -737,12 +737,12 @@ DoubleML = R6Class("DoubleML",
     #' a `matrix()` with adjusted p_values.
     p_adjust = function(method = "romano-wolf", return_matrix = TRUE) {
       if (all(is.na(self$coef))) {
-        stop("apply fit() before p_adust().")
+        stop("apply fit() before p_adjust().")
       }
 
       if (tolower(method) %in% c("rw", "romano-wolf")) {
         if (is.null(self$boot_t_stat) | all(is.na(self$coef))) {
-          stop("apply fit() & bootstrap() before p_adust().")
+          stop("apply fit() & bootstrap() before p_adjust().")
         }
         k = self$data$n_treat
         pinit = p_val_corrected = vector(mode = "numeric", length = k)

R/double_ml_data.R

@@ -245,12 +245,12 @@ DoubleMLData = R6Class("DoubleMLData",
         if (any(z_cols %in% d_cols)) {
           stop(paste(
             "At least one variable/column is set as treatment",
-            "variable ('d_cols') and instrumental variable in 'z_cols')."))
+            "variable ('d_cols') and instrumental variable in 'z_cols'."))
         }
         if (any(z_cols %in% x_cols)) {
           stop(paste(
             "At least one variable/column is set as covariate ('x_cols')",
-            "and instrumental variable in 'z_cols')."))
+            "and instrumental variable in 'z_cols'."))
         }
       }
     }

R/double_ml_iivm.R

@@ -463,24 +463,35 @@ DoubleMLIIVM = R6Class("DoubleMLIIVM",
     },
     check_data = function(obj_dml_data) {
       one_treat = (obj_dml_data$n_treat == 1)
-      binary_treat = test_integerish(obj_dml_data$data[[obj_dml_data$d_cols]],
-        lower = 0, upper = 1)
-      if (!(one_treat & binary_treat)) {
-        stop(paste(
-          "Incompatible data.\n",
-          "To fit an IIVM model with DoubleML",
-          "exactly one binary variable with values 0 and 1",
-          "needs to be specified as treatment variable."))
+      err_msg = paste(
+        "Incompatible data.\n",
+        "To fit an IIVM model with DoubleML",
+        "exactly one binary variable with values 0 and 1",
+        "needs to be specified as treatment variable.")
+      if (one_treat) {
+        binary_treat = test_integerish(obj_dml_data$data[[obj_dml_data$d_cols]],
+                                       lower = 0, upper = 1)
+        if (!(one_treat & binary_treat)) {
+          stop(err_msg)
+        }
+      } else {
+        stop(err_msg)
       }
+
       one_instr = (obj_dml_data$n_instr == 1)
-      binary_instr = test_integerish(obj_dml_data$data[[obj_dml_data$z_cols]],
-        lower = 0, upper = 1)
-      if (!(one_instr & binary_instr)) {
-        stop(paste(
-          "Incompatible data.\n",
-          "To fit an IIVM model with DoubleML",
-          "exactly one binary variable with values 0 and 1",
-          "needs to be specified as instrumental variable."))
+      err_msg = paste(
+        "Incompatible data.\n",
+        "To fit an IIVM model with DoubleML",
+        "exactly one binary variable with values 0 and 1",
+        "needs to be specified as instrumental variable.")
+      if (one_instr) {
+        binary_instr = test_integerish(obj_dml_data$data[[obj_dml_data$z_cols]],
+                                       lower = 0, upper = 1)
+        if (!(one_instr & binary_instr)) {
+          stop(err_msg)
+        }
+      } else {
+        stop(err_msg)
       }
       return()
     }

R/double_ml_irm.R

@@ -351,9 +351,7 @@ DoubleMLIRM = R6Class("DoubleMLIRM",
         check_class(score, "function"))
       if (is.character(score)) {
         valid_score = c("ATE", "ATTE")
-        if (!(score %in% valid_score)) {
-          assertChoice(score, valid_score)
-        }
+        assertChoice(score, valid_score)
       }
       return()
     },
@@ -362,18 +360,23 @@ DoubleMLIRM = R6Class("DoubleMLIRM",
         stop(paste(
           "Incompatible data.\n", paste(obj_dml_data$z_cols, collapse = ", "),
           "has been set as instrumental variable(s).\n",
-          "To fit an interactive IV regression model use DoubleMLIIVM
-          instead of DoubleMLIRM."))
+          "To fit an interactive IV regression model use DoubleMLIIVM",
+          "instead of DoubleMLIRM."))
       }
       one_treat = (obj_dml_data$n_treat == 1)
-      binary_treat = test_integerish(obj_dml_data$data[[obj_dml_data$d_cols]],
-        lower = 0, upper = 1)
-      if (!(one_treat & binary_treat)) {
-        stop(paste(
-          "Incompatible data.\n",
-          "To fit an IRM model with DoubleML",
-          "exactly one binary variable with values 0 and 1",
-          "needs to be specified as treatment variable."))
+      err_msg = paste(
+        "Incompatible data.\n",
+        "To fit an IRM model with DoubleML",
+        "exactly one binary variable with values 0 and 1",
+        "needs to be specified as treatment variable.")
+      if (one_treat) {
+        binary_treat = test_integerish(obj_dml_data$data[[obj_dml_data$d_cols]],
+                                       lower = 0, upper = 1)
+        if (!(one_treat & binary_treat)) {
+          stop(err_msg)
+        }
+      } else {
+        stop(err_msg)
       }
       return()
     }

R/double_ml_pliv.R

@@ -348,8 +348,8 @@ DoubleMLPLIV = R6Class("DoubleMLPLIV",
           psi_b = psi_b)
       } else if (is.function(self$score)) {
         if (self$data$n_instr > 1) {
-          stop("Callable score not implemented for DoubleMLPLIV with
-               partialX=TRUE and partialZ=FALSE with several instruments")
+          stop(paste("Callable score not implemented for DoubleMLPLIV with",
+                     "partialX=TRUE and partialZ=FALSE with several instruments."))
         }
         psis = self$score(y, z, d, g_hat, m_hat, r_hat, smpls)
       }
@@ -415,8 +415,8 @@ DoubleMLPLIV = R6Class("DoubleMLPLIV",
           psi_a = psi_a,
           psi_b = psi_b)
       } else if (is.function(self$score)) {
-        stop("Callable score not implemented for DoubleMLPLIV
-             with partialX=TRUE and partialZ=TRUE.")
+        stop(paste("Callable score not implemented for DoubleMLPLIV",
+                   "with partialX=TRUE and partialZ=TRUE."))
         # res = self$score(y, d, g_hat, m_hat, m_hat_tilde)
       }
       res$preds = list(
@@ -454,8 +454,8 @@ DoubleMLPLIV = R6Class("DoubleMLPLIV",
         }
         res = list(psi_a = psi_a, psi_b = psi_b)
       } else if (is.function(self$score)) {
-        stop("Callable score not implemented for DoubleMLPLIV
-             with partialX=FALSE and partialZ=TRUE.")
+        stop(paste("Callable score not implemented for DoubleMLPLIV",
+                   "with partialX=FALSE and partialZ=TRUE."))
         # res = self$score(y, z, d, r_hat)
       }
       res$preds = list("ml_r" = r_hat)
@@ -674,6 +674,13 @@ DoubleMLPLIV = R6Class("DoubleMLPLIV",
       return()
     },
     check_data = function(obj_dml_data) {
+      if (obj_dml_data$n_instr == 0) {
+        stop(paste(
+          "Incompatible data.\n",
+          "At least one variable must be set as instrumental variable.\n",
+          "To fit a partially linear regression model without instrumental",
+          "variable(s) use DoubleMLPLR instead of DoubleMLPLIV."))
+      }
       return()
     }
   )

R/double_ml_plr.R

@@ -270,8 +270,8 @@ DoubleMLPLR = R6Class("DoubleMLPLR",
         stop(paste(
           "Incompatible data.\n", paste(obj_dml_data$z_cols, collapse = ", "),
           "has been set as instrumental variable(s).\n",
-          "To fit a partially linear IV regression model use
-          DoubleMLPLIV instead of DoubleMLPLR."))
+          "To fit a partially linear IV regression model use",
+          "DoubleMLPLIV instead of DoubleMLPLR."))
       }
       return()
     }

tests/testthat/helper-01-dgp.R

@@ -0,0 +1,60 @@
+se_repeated = function(se_s, coefficients, theta_s) {
+    se = sqrt(stats::median(se_s^2 + (theta_s - coefficients)^2))
+  return(se)
+}
+
+
+sample_splitting = function(k, data) {
+
+  resampling = mlr3::ResamplingCV$new()
+  resampling$param_set$values$folds = k
+
+  dummy_task = mlr3::Task$new("dummy_resampling", "regr", data)
+  resampling = resampling$instantiate(dummy_task)
+
+  n_iters = resampling$iters
+  train_ids = lapply(1:n_iters, function(x) resampling$train_set(x))
+  test_ids = lapply(1:n_iters, function(x) resampling$test_set(x))
+
+  return(list(train_ids = train_ids, test_ids = test_ids))
+}
+
+
+draw_bootstrap_weights = function(bootstrap, n_rep_boot, n_obs) {
+  if (bootstrap == "Bayes") {
+    weights = stats::rexp(n_rep_boot * n_obs, rate = 1) - 1
+  } else if (bootstrap == "normal") {
+    weights = stats::rnorm(n_rep_boot * n_obs)
+  } else if (bootstrap == "wild") {
+    weights = stats::rnorm(n_rep_boot * n_obs) / sqrt(2) + (stats::rnorm(n_rep_boot * n_obs)^2 - 1) / 2
+  } else {
+    stop("invalid boot method")
+  }
+  weights = matrix(weights, nrow = n_rep_boot, ncol = n_obs, byrow = TRUE)
+  
+  return(weights)
+}
+
+
+functional_bootstrap = function(theta, se, psi, psi_a, k, smpls,
+                                n_rep_boot, weights) {
+  score = psi
+  J = mean(psi_a)
+  boot_coef = matrix(NA, nrow = 1, ncol = n_rep_boot)
+  boot_t_stat = matrix(NA, nrow = 1, ncol = n_rep_boot)
+  for (i in seq(n_rep_boot)) {
+    boot_coef[1, i] = mean(weights[i, ] * 1 / J * score)
+    boot_t_stat[1, i] = boot_coef[1, i] / se
+  }
+
+  res = list(boot_coef = boot_coef, boot_t_stat = boot_t_stat)
+  return(res)
+}
+
+trim_vec = function(values, trimming_threshold) {
+  if (trimming_threshold > 0) {
+    values[values < trimming_threshold] = trimming_threshold
+    values[values > 1 - trimming_threshold] = 1 - trimming_threshold
+  }
+  return(values)
+}

tests/testthat/helper-01-helper_functions.R

@@ -132,3 +132,41 @@ dgp1_toeplitz = function(n, p, betamax = 4, decay = 0.99, threshold = 0, noiseva
 
   return(data)
 }
+
+make_data_pliv_partialZ = function(n_obs, alpha = 1, dim_x = 5, dim_z = 150) {
+  sigma_e_u = matrix(c(1, 0.6, 0.6, 1), ncol = 2)
+  mu_e_u = rep(0, 2)
+  e_u = mvtnorm::rmvnorm(n = n_obs, mean = mu_e_u, sigma = sigma_e_u)
+  epsilon = e_u[, 1]
+  u = e_u[, 2]
+  
+  sigma_x = toeplitz(0.5^(0:(dim_x - 1)))
+  mu_x = rep(0, dim_x)
+  x = mvtnorm::rmvnorm(n = n_obs, mean = mu_x, sigma = sigma_x)
+  
+  I_z = diag(x = 1, ncol = dim_z, nrow = dim_z)
+  mu_xi = rep(0, dim_z)
+  xi = mvtnorm::rmvnorm(n = n_obs, mean = mu_xi, sigma = 0.25 * I_z)
+  
+  beta = 1 / (1:dim_x)^2
+  gamma = beta
+  delta = 1 / (1:dim_z)^2
+  
+  zeros = matrix(0, nrow = dim_x, ncol = (dim_z - dim_x))
+  I_x = diag(x = 1, ncol = dim_x, nrow = dim_x)
+  Pi = cbind(I_x, zeros)
+  
+  z = x %*% Pi + xi
+  d = x %*% gamma + z %*% delta + u
+  y = alpha * d + x %*% beta + epsilon
+  
+  
+  colnames(x) = paste0("X", 1:dim_x)
+  colnames(z) = paste0("Z", 1:dim_z)
+  colnames(y) = "y"
+  colnames(d) = "d"
+  
+  data = data.frame(x, y, d, z)
+  return(data)
+}
+