fix: review findings for montecarlo/04_stata_coverage.R

.Rbuildignore

@@ -11,3 +11,6 @@
 ^src/.*\.so$
 ^src/.*\.o$
 ^src/.*\.dll$
+^scripts$
+^reviews$
+^montecarlo$

.Rhistory

@@ -1,3 +1,73 @@
+}
+})
+################################################################################
+# 3) Test bootstrap (boot=TRUE) with method='simple', check we have boot_out
+################################################################################
+test_that("r3d() boot=TRUE yields boot_out with confidence bands & tests (sharp)", {
+set.seed(106)
+test_data <- create_test_data(n = 40, sample_size = 20, fuzzy = FALSE)
+out_boot <- r3d(X = test_data$x, Y_list = test_data$y_list, method = "simple", p = 1,
+boot = TRUE, boot_reps = 10, test = "nullity", alpha = 0.2)
+expect_s3_class(out_boot, "r3d")
+expect_true("boot_out" %in% names(out_boot))
+expect_true("boot_out" %in% names(out_boot$bootstrap))
+expect_identical(out_boot$boot_out, out_boot$bootstrap$boot_out)
+bo <- out_boot$boot_out
+expect_true(all(c("cb_lower", "cb_upper", "test_stat", "test_crit_val", "p_value") %in% names(bo)))
+expect_equal(length(bo$cb_lower), length(out_boot$q_grid))
+expect_equal(length(bo$cb_upper), length(out_boot$q_grid))
+expect_equal(dim(bo$boot_taus), c(length(out_boot$q_grid), 10))
+expect_true(all(bo$cb_upper >= out_boot$tau, na.rm = TRUE))
+expect_true(all(bo$cb_lower <= out_boot$tau, na.rm = TRUE))
+})
+test_that("r3d() boot=TRUE yields boot_out with confidence bands & tests (fuzzy)", {
+set.seed(106)
+test_data <- create_test_data(n = 40, sample_size = 20, fuzzy = TRUE)
+out_boot_fuzzy <- r3d(X = test_data$x, Y_list = test_data$y_list, T = test_data$t, fuzzy = TRUE,
+method = "simple", p = 1, boot = TRUE, boot_reps = 10, test = "nullity", alpha = 0.2)
+expect_s3_class(out_boot_fuzzy, "r3d")
+expect_true("boot_out" %in% names(out_boot_fuzzy))
+expect_true("boot_out" %in% names(out_boot_fuzzy$bootstrap))
+expect_identical(out_boot_fuzzy$boot_out, out_boot_fuzzy$bootstrap$boot_out)
+bo <- out_boot_fuzzy$boot_out
+expect_true(all(c("cb_lower", "cb_upper", "test_stat", "test_crit_val", "p_value") %in% names(bo)))
+expect_equal(length(bo$cb_lower), length(out_boot_fuzzy$q_grid))
+expect_equal(length(bo$cb_upper), length(out_boot_fuzzy$q_grid))
+expect_equal(dim(bo$boot_taus), c(length(out_boot_fuzzy$q_grid), 10))
+expect_true(all(bo$cb_upper >= out_boot_fuzzy$tau, na.rm = TRUE))
+expect_true(all(bo$cb_lower <= out_boot_fuzzy$tau, na.rm = TRUE))
+})
+test_that("r3d() with boot=TRUE and test='homogeneity' works properly", {
+set.seed(1066)
+hetero_data <- create_hetero_effect_data(n = 60)
+q_grid <- c(0.1, 0.3, 0.5, 0.7, 0.9)
+out_test <- r3d(X = hetero_data$x, Y_list = hetero_data$y_list, cutoff = 0,
+method = "simple", p = 1, q_grid = q_grid,
+boot = TRUE, boot_reps = 20, test = "homogeneity")
+expect_true(!is.na(out_test$boot_out$test_stat))
+expect_true(!is.na(out_test$boot_out$p_value))
+expect_true(out_test$boot_out$test_stat > 0)
+})
+test_that("r3d_bootstrap() can be called directly with r3d object", {
+set.seed(107)
+test_data <- create_test_data(n = 40, sample_size = 20, fuzzy = FALSE)
+r3d_obj <- r3d(X = test_data$x, Y_list = test_data$y_list, method = "simple", p = 1, boot = FALSE)
+bo <- r3d_bootstrap(object = r3d_obj, X = test_data$x, Y_list = test_data$y_list,
+B = 5, alpha = 0.1, test = "homogeneity", cores = 1)
+expect_true(all(c("cb_lower", "cb_upper", "p_value") %in% names(bo)))
+expect_equal(length(bo$cb_lower), length(r3d_obj$q_grid))
+expect_true(!is.na(bo$p_value))
+})
+################################################################################
+# 4) Scenario-based test with random distributions
+################################################################################
+delta_mu <- 2
+delta_sigma <- 0.5
+base_mu_below <- 5
+slope_mu_below <- 0.5
+base_sigma_below <- 1
+slope_sigma_below <- 0.2
+c_cutoff <- 0
 simulate_data_scenario1 <- function(N, n_obs) {
 x <- runif(N, -1, 1)
 y <- vector("list", N)
@@ -403,110 +473,40 @@ expect_true(any(grepl("Gini coefficient above cutoff:", summary_output, fixed =
 expect_true(any(grepl("Gini coefficient below cutoff:", summary_output, fixed = TRUE)))
 expect_true(any(grepl("Difference in Gini coefficients:", summary_output, fixed = TRUE)))
 })
-?r3d_bwselect
-library(R3D)
-?r3d_bwselect
+devtools::check(cran = TRUE)
+devtools::document()
+devtools::document()
+devtools::build()
+devtools::check(cran = TRUE)
+devtools::document()
+devtools::build()
+devtools::check(cran = TRUE)
+devtools::document()
+devtools::build()
+devtools::check(cran = TRUE)
 devtools::document()
-devtools::install()
-devtools::install()
-devtools::install()
-devtools::install()
-#****************************************************************************************************************************************************
-# MASTER SCRIPT: DISTRIBUTIONAL REGRESSION DISCONTINUITY DESIGN
-# David Van Dijcke
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#### SET OVERALL PARAMETERS ####
-#### SET PATHS ####
-if (!require("here", character.only=T)) {install.packages("here", dependencies=TRUE)}; require("here")
-codeDir <- here::here()
-setwd(codeDir) # sets cd to program directory
-dir <- sub("/[^/]+$", "", codeDir)# get main directory
-dataIn <- file.path(dir, "data", "in")
-dataBy <- file.path(dir, "data", "by")
-dataOut <- file.path(dir, "data", "out")
-overleaf <- "/Users/davidvandijcke/Dropbox (University of Michigan)/Apps/Overleaf/r3d"
-tabs <- file.path(overleaf, "tabs")
-figs <- file.path(overleaf, "figs")
-#### USER-WRITTEN FUNCTIONS ####
-source("utils.R")
-?quantile
-#### robustness: use same-year income
-data_offset <- prepData(year_offset=0)
-?match.call
-devtools::install()
-devtools::install()
-devtools::install()
-devtools::install()
-devtools::install()
-remove.packages("modelsummary")
-source("r3d_simulation_scenarios.R")
-#****************************************************************************************************************************************************
-# MASTER SCRIPT: DISTRIBUTIONAL REGRESSION DISCONTINUITY DESIGN
-# David Van Dijcke
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#### SET OVERALL PARAMETERS ####
-#### SET PATHS ####
-if (!require("here", character.only=T)) {install.packages("here", dependencies=TRUE)}; require("here")
-codeDir <- here::here()
-setwd(codeDir) # sets cd to program directory
-dir <- sub("/[^/]+$", "", codeDir)# get main directory
-dataIn <- file.path(dir, "data", "in")
-dataBy <- file.path(dir, "data", "by")
-dataOut <- file.path(dir, "data", "out")
-overleaf <- "/Users/davidvandijcke/Dropbox (University of Michigan)/Apps/Overleaf/r3d"
-tabs <- file.path(overleaf, "tabs")
-figs <- file.path(overleaf, "figs")
-#### USER-WRITTEN FUNCTIONS ####
-source("utils.R")
-#****************************************************************************************************************************************************
-# MASTER SCRIPT: DISTRIBUTIONAL REGRESSION DISCONTINUITY DESIGN
-# David Van Dijcke
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#### SET OVERALL PARAMETERS ####
-#### SET PATHS ####
-if (!require("here", character.only=T)) {install.packages("here", dependencies=TRUE)}; require("here")
-codeDir <- here::here()
-setwd(codeDir) # sets cd to program directory
-dir <- sub("/[^/]+$", "", codeDir)# get main directory
-dataIn <- file.path(dir, "data", "in")
-dataBy <- file.path(dir, "data", "by")
-dataOut <- file.path(dir, "data", "out")
-overleaf <- "/Users/davidvandijcke/Dropbox (University of Michigan)/Apps/Overleaf/r3d"
-tabs <- file.path(overleaf, "tabs")
-figs <- file.path(overleaf, "figs")
-#### USER-WRITTEN FUNCTIONS ####
-source("utils.R")
-#****************************************************************************************************************************************************
-# MASTER SCRIPT: DISTRIBUTIONAL REGRESSION DISCONTINUITY DESIGN
-# David Van Dijcke
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#****************************************************************************************************************************************************
-#### SET OVERALL PARAMETERS ####
-#### SET PATHS ####
-if (!require("here", character.only=T)) {install.packages("here", dependencies=TRUE)}; require("here")
-codeDir <- here::here()
-setwd(codeDir) # sets cd to program directory
-dir <- sub("/[^/]+$", "", codeDir)# get main directory
-dataIn <- file.path(dir, "data", "in")
-dataBy <- file.path(dir, "data", "by")
-dataOut <- file.path(dir, "data", "out")
-overleaf <- "/Users/davidvandijcke/Dropbox (University of Michigan)/Apps/Overleaf/r3d"
-tabs <- file.path(overleaf, "tabs")
-figs <- file.path(overleaf, "figs")
-#### USER-WRITTEN FUNCTIONS ####
-source("utils.R")
+devtools::build()
+devtools::check(cran = TRUE)
+pkgdown::build_site()
+?pkgdown::builds_news
+?build_news
+?pkgdown::builds_news
+pkgdown::build_site()
+usethis::use_github_action("pkgdown")
+usethis::use_logo("docs/logo.png")
+pkgdown::build_site()
+usethis::use_logo("docs/logo.png")
+usethis::use_logo("docs/logo.png")
+pkgdown::build_site()
+usethis::use_logo("docs/logo.png")
+pkgdown::build_site()
+usethis::use_logo("docs/logo.png")
+usethis::use_logo("docs/logo.png")
+pkgdown::build_site()
+pkgdown::build_site()
+usethis::use_logo("docs/logo.png")
+usethis::use_logo("docs/logo.png")
+pkgdown::build_site()
+pkgdown::build_favicons()
+pkgdown::build_favicons(overwrite=TRUE)
+pkgdown::build_site()

.gitignore

@@ -16,3 +16,6 @@ montecarlo/test_*.log
 # Stata build artifacts
 stata_r3d/mata/*.log
 stata_r3d/mata/*.smcl
+
+# Review pipeline artifacts
+reviews/

montecarlo/04_stata_coverage.R

@@ -10,17 +10,17 @@ cat("=== R3D Monte Carlo: Stata Uniform Coverage ===\n\n")
 outdir <- file.path(getwd(), "output")
 stata_dir <- file.path(outdir, "results_Stata")
 
-NQ <- 9
 Q_GRID <- (1:9) / 10
+NQ <- length(Q_GRID)  # F4: derived from Q_GRID, not hardcoded independently
 N_SIMS <- 100
 
 # Pre-compute DGP2 true tau
 true_tau_dgp2 <- list()
+set.seed(99999)  # F1: seed set once before the loop, not reset each iteration
+large_n <- 1e6
+y_below <- rnorm(large_n, mean = 0, sd = 1) + 2 * rexp(large_n, rate = 1)  # F1: cached outside loop
 for (delta in c(0, 1, 2)) {
-  set.seed(99999)
-  large_n <- 1e6
   y_above <- rnorm(large_n, mean = delta, sd = 1) + 2 * rexp(large_n, rate = 1)
-  y_below <- rnorm(large_n, mean = 0, sd = 1) + 2 * rexp(large_n, rate = 1)
   true_tau_dgp2[[as.character(delta)]] <-
     quantile(y_above, probs = Q_GRID, names = FALSE) -
     quantile(y_below, probs = Q_GRID, names = FALSE)
@@ -35,6 +35,9 @@ cells <- expand.grid(
   stringsAsFactors = FALSE
 )
 
+# F6: check stata_dir exists before starting aggregation
+if (!dir.exists(stata_dir)) stop("Stata results directory not found: ", stata_dir)
+
 results <- list()
 
 for (i in seq_len(nrow(cells))) {
@@ -51,6 +54,7 @@ for (i in seq_len(nrow(cells))) {
   } else {
     true_tau <- true_tau_dgp2[[as.character(delta)]]
   }
+  stopifnot(length(true_tau) == NQ)  # F3: guard against silent recycling
 
   for (tag in c("sharedbw", "ownbw")) {
     # Collect results across sims
@@ -60,7 +64,7 @@ for (i in seq_len(nrow(cells))) {
     pval_null <- rep(NA, N_SIMS)
     pval_homo <- rep(NA, N_SIMS)
 
-    for (s in 1:N_SIMS) {
+    for (s in seq_len(N_SIMS)) {  # F5: seq_len avoids 1:0 pitfall
       f <- file.path(stata_dir, sprintf("mc_res_%s_%s_%d.csv", cell, tag, s))
       if (!file.exists(f)) next
       d <- read.csv(f)
@@ -77,11 +81,18 @@ for (i in seq_len(nrow(cells))) {
       }
     }
 
-    valid <- !is.na(tau_mat[, 1])
+    # F2: full-row validity across all three matrices
+    valid <- rowSums(is.na(tau_mat)) == 0 &
+             rowSums(is.na(cb_lo_mat)) == 0 &
+             rowSums(is.na(cb_hi_mat)) == 0
     n_valid <- sum(valid)
 
     if (n_valid == 0) next
 
+    # F7: separate validity mask for p-value-based rejection rates
+    pval_valid <- valid & !is.na(pval_null) & !is.na(pval_homo)
+    n_pval_valid <- sum(pval_valid)
+
     # Uniform coverage: true tau within [cb_lo, cb_hi] at ALL quantiles simultaneously
     covered_mat <- cb_lo_mat[valid, , drop = FALSE] <=
       matrix(true_tau, nrow = n_valid, ncol = NQ, byrow = TRUE) &
@@ -102,13 +113,13 @@ for (i in seq_len(nrow(cells))) {
                              matrix(true_tau, nrow = n_valid, ncol = NQ, byrow = TRUE))^2))
     mean_rmse <- mean(rmse)
 
-    # Rejection rates
-    rej_null <- mean(pval_null[valid] < 0.05, na.rm = TRUE)
-    rej_homo <- mean(pval_homo[valid] < 0.05, na.rm = TRUE)
+    # F7: Rejection rates use pval_valid mask (correct denominator)
+    rej_null <- mean(pval_null[pval_valid] < 0.05, na.rm = TRUE)
+    rej_homo <- mean(pval_homo[pval_valid] < 0.05, na.rm = TRUE)
 
     results[[length(results) + 1]] <- data.frame(
       cell = cell, dgp = dgp, method = method, n = n_obs, delta = delta,
-      bw_type = tag, n_valid = n_valid,
+      bw_type = tag, n_valid = n_valid, n_pval_valid = n_pval_valid,  # F7: report n_pval_valid
       uniform_coverage = uniform_coverage,
       min_pointwise_cov = min(pointwise_coverage),
       mean_pointwise_cov = mean(pointwise_coverage),
@@ -123,6 +134,9 @@ for (i in seq_len(nrow(cells))) {
 
 df <- do.call(rbind, results)
 
+# F6: fail fast with clear diagnostic if no results were collected
+if (is.null(df) || nrow(df) == 0) stop("No valid Stata results found in: ", stata_dir)
+
 # Print summary table
 cat(sprintf("%-45s  %4s  %6s  %6s  %6s  %6s  %6s\n",
             "Cell", "N", "UCov", "MnPCov", "Bias", "RMSE", "RejH0"))

stata_r3d/ado/r3d.ado

@@ -6,12 +6,19 @@ program define r3d, eclass
     version 18.0
 
     // Load Mata functions: try mlib first, then compile from source
+    local mata_loaded 0
     quietly capture mata: mata which r3d_compute_quantiles
-    if _rc != 0 {
+    if _rc == 0 {
+        local mata_loaded 1
+    }
+    if !`mata_loaded' {
         mata: mata mlib index
         quietly capture mata: mata which r3d_compute_quantiles
+        if _rc == 0 {
+            local mata_loaded 1
+        }
     }
-    if _rc != 0 {
+    if !`mata_loaded' {
         quietly capture do "`c(sysdir_plus)'r/r3d_mata.mata"
         if _rc != 0 {
             quietly capture do "mata/r3d_mata.mata"
@@ -152,14 +159,14 @@ program define r3d, eclass
     // Check if fuzzy RD
     local is_fuzzy 0
     if "`fuzzy'" != "" {
-        confirm numeric variable ``fuzzy''
+        confirm numeric variable `fuzzy'
         local is_fuzzy 1
         tempvar t_centered
         quietly gen double `t_centered' = `fuzzy' if `touse'
     }
 
     if "`weights'" != "" {
-        confirm numeric variable ``weights''
+        confirm numeric variable `weights'
     }
 
     // Compute empirical quantiles for each observation
@@ -396,8 +403,8 @@ program define r3d, eclass
                 replace tau = `tau'[1,`i'] in `i'
                 replace se = `se'[1,`i'] in `i'
                 if `bootstrap' > 0 {
-                    replace ci_lower = tau - invnormal(1-(1-`level_val'/100)/2)*se in `i'
-                    replace ci_upper = tau + invnormal(1-(1-`level_val'/100)/2)*se in `i'
+                    replace ci_lower = e(cb_lower)[1,`i'] in `i'
+                    replace ci_upper = e(cb_upper)[1,`i'] in `i'
                 }
             }
         }