From b87e353c9bdb20dee20a1ed4e776541fc5b8341f Mon Sep 17 00:00:00 2001 From: DGoettlich Date: Wed, 22 Apr 2026 01:55:24 +0200 Subject: [PATCH 1/2] [scpcR] added fe x iv x large-n branch; added back large-n fork --- DESCRIPTION | 2 +- R/core.R | 85 +++++- R/utils-data.R | 18 ++ R/utils-spatial.R | 292 ++++++++++++++++++++- man/scpc.Rd | 16 ++ tests/testthat/test-scpc-parity-stata.R | 78 ++++++ tests/testthat/test_scpc.R | 244 +++++++++++++++-- tests/testthat/test_scpc_fixest_iv_smoke.R | 50 ++++ tests/testthat/test_set_oms_wfin.R | 102 ++++++- 9 files changed, 834 insertions(+), 53 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index f336241..a0a0bf5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: scpcR Type: Package Title: Spatial Correlation-Robust Inference for Regression Coefficients -Version: 0.1.2 +Version: 0.1.3 Authors@R: c( person("David", "Boll", email = "david.boll@warwick.ac.uk", role = c("aut", "cre")), person("Daniel", "Goettlich", email = "daniel.goettlich@econ.uzh.ch", role = "aut") diff --git a/R/core.R b/R/core.R index a42d107..9e41916 100644 --- a/R/core.R +++ b/R/core.R @@ -27,6 +27,14 @@ #' @param avc Numeric; upper bound on the average pairwise correlation #' for which size is controlled. Must be in \code{(0.001, 0.99)}. #' Default is 0.03. +#' @param method Character string selecting the spatial algorithm. +#' \code{"auto"} chooses \code{"exact"} for smaller problems and +#' \code{"approx"} once the large-\eqn{n} threshold is reached. +#' \code{"exact"} always uses the full distance matrix. +#' \code{"approx"} uses the original large-\eqn{n} approximation branch. +#' @param large_n_seed Numeric; integer-valued seed used by the +#' large-\eqn{n} approximation branch. Ignored when +#' \code{method = "exact"}. Default is 1. #' @param uncond Logical; if \code{TRUE}, report unconditional critical #' values only (skip the conditional adjustment of Mueller and Watson #' 2023). Default is \code{FALSE}. @@ -46,6 +54,10 @@ #' \code{avc}.} #' \item{\code{cv}}{The unconditional 5\% critical value.} #' \item{\code{q}}{Number of spatial principal components selected.} +#' \item{\code{method}}{Spatial algorithm used: +#' \code{"exact"} or \code{"approx"}.} +#' \item{\code{large_n_seed}}{Seed used by the large-\eqn{n} +#' approximation branch.} #' \item{\code{call}}{The matched call.} #' } #' @@ -90,6 +102,8 @@ scpc <- function(model, cluster = NULL, ncoef = NULL, avc = 0.03, + method = "auto", + large_n_seed = 1, uncond = FALSE, cvs = FALSE) { @@ -99,6 +113,8 @@ scpc <- function(model, if (avc <= 0.001 || avc >= 0.99) { stop("`avc` must lie in (0.001, 0.99).") } + method <- .validate_scpc_method(method) + large_n_seed <- .validate_large_n_seed(large_n_seed) if ( !is.null(ncoef) && ( @@ -218,12 +234,19 @@ scpc <- function(model, if (ncol(coords) == 1) coords <- cbind(coords, 0) ## 2. Spatial kernel (unconditional) --------------------------------------- - D <- .getdistmat(coords, latlong) - spc <- .setOmsWfin(D, avc) + spc <- .setOmsWfin( + coords, + avc0 = avc, + latlong = latlong, + method = method, + large_n_seed = large_n_seed + ) Wfin <- spc$Wfin cvfin <- spc$cvfin Omsfin <- spc$Omsfin + perm <- spc$perm q <- ncol(Wfin) - 1 + large_n_random_state <- spc$random_state ## 3. Bread ---------------------------------------------------------------- bread_inv <- sandwich::bread(model) / n @@ -241,58 +264,87 @@ scpc <- function(model, for (j in seq_len(k_use)) { coef_j <- coef_vec[[j]] wj <- as.numeric(neff * bread_inv[j, ] %*% t(S)) + coef_j + wj_perm <- wj[perm] ## unconditional statistic ----------------------------------------------- - denom <- sqrt(sum((t(Wfin[, -1, drop = FALSE]) %*% wj)^2)) + denom <- sqrt(sum((t(Wfin[, -1, drop = FALSE]) %*% wj_perm)^2)) if (!is.finite(denom) || denom <= sqrt(.Machine$double.eps)) { stop("SCPC variance projection is degenerate for coefficient: ", coef_names[[j]]) } - tau_u <- as.numeric(sqrt(q) * crossprod(Wfin[, 1], wj) / denom) + tau_u <- as.numeric(sqrt(q) * crossprod(Wfin[, 1], wj_perm) / denom) SE <- as.numeric(denom / (sqrt(q) * sqrt(neff))) p_u <- .maxrp(Omsfin, q, abs(tau_u) / sqrt(q))$max ## conditional branch (skip when uncond = TRUE) -------------------------- if (!uncond) { + cl_vec_scpc <- NULL + if (!is.null(cluster)) { + cl_vec_scpc <- if (!identical(perm, seq_len(length(perm)))) { + factor(as.character(cl_vec), levels = levels(cl_vec)[perm]) + } else { + cl_vec + } + } + if (!is.null(cluster)) { xj_indiv <- as.numeric(neff * bread_inv[j, ] %*% t(model_mat_cond)) if (is_fixest_iv) { Wx <- .orthogonalize_W_cluster_iv( Wfin = Wfin, - cl_vec = cl_vec, + cl_vec = cl_vec_scpc, xj_indiv = xj_indiv, residualize = iv_residualizer ) } else { ## Clustered conditional: individual-level orthogonalization Wx <- .orthogonalize_W_cluster( - Wfin, cl_vec, xj_indiv, model_mat_cond, + Wfin, cl_vec_scpc, xj_indiv, model_mat_cond, include_intercept = cond_include_intercept ) } } else if (is_fixest_iv) { - xj <- as.numeric(neff * bread_inv[j, ] %*% t(model_mat_cond)) + xj_raw <- as.numeric(neff * bread_inv[j, ] %*% t(model_mat_cond)) + if (isTRUE(spc$large_n)) { + residualize <- .make_iv_residualizer( + iv_cond_design$X[perm, , drop = FALSE], + iv_cond_design$Z[perm, , drop = FALSE], + fixef_id = .permute_fixef_id(iv_cond_design$fixef_id, perm) + ) + xj <- xj_raw[perm] + } else { + residualize <- iv_residualizer + xj <- xj_raw + } xjs <- sign(xj) Wx <- .orthogonalize_W_iv( Wfin = Wfin, xj = xj, xjs = xjs, - residualize = iv_residualizer + residualize = residualize ) } else { ## Non-clustered conditional xj <- as.numeric(neff * bread_inv[j, ] %*% t(model_mat_cond)) + model_mat_cond_j <- model_mat_cond + if (isTRUE(spc$large_n)) { + xj <- xj[perm] + model_mat_cond_j <- model_mat_cond[perm, , drop = FALSE] + cond_fixef_id_j <- .permute_fixef_id(cond_fixef_id, perm) + } else { + cond_fixef_id_j <- cond_fixef_id + } xjs <- sign(xj) Wx <- .orthogonalize_W( - Wfin, xj, xjs, model_mat_cond, + Wfin, xj, xjs, model_mat_cond_j, include_intercept = cond_include_intercept, - fixef_id = cond_fixef_id + fixef_id = cond_fixef_id_j ) } if (!all(is.finite(Wx))) { stop("Conditional projection produced non-finite W values for coefficient: ", coef_names[[j]]) } - if (nrow(Wx) != nrow(D) || ncol(Wx) != ncol(Wfin)) { + if (nrow(Wx) != nrow(Wfin) || ncol(Wx) != ncol(Wfin)) { stop("Conditional projection produced W with incompatible dimensions for coefficient: ", coef_names[[j]]) } if ( @@ -302,7 +354,14 @@ scpc <- function(model, stop("Conditional projection collapsed all non-constant W columns for coefficient: ", coef_names[[j]]) } - Omsx <- .getOms(D, spc$c0, spc$cmax, Wx, 1.2) + omsx_res <- .get_conditional_Oms( + spc, + Wx, + latlong = latlong, + random_state = large_n_random_state + ) + Omsx <- omsx_res$Oms + large_n_random_state <- omsx_res$random_state p_c <- .maxrp(Omsx, q, abs(tau_u) / sqrt(q))$max cvx <- .getcv(Omsx, q, 0.05) p_final <- max(p_u, p_c) @@ -342,6 +401,8 @@ scpc <- function(model, c0 = spc$c0, cv = cvfin, q = q, + method = spc$method, + large_n_seed = large_n_seed, call = match.call() ), class = "scpc") } diff --git a/R/utils-data.R b/R/utils-data.R index d03dcb3..ccc3cf8 100644 --- a/R/utils-data.R +++ b/R/utils-data.R @@ -141,6 +141,24 @@ length(model$fixef_vars) > 0L } +.permute_fixef_id <- function(fixef_id, perm) { + ## the large-n branch reorders rows before it rebuilds the iv + ## residualizer. fixest keeps fixed-effect ids outside X and Z, so they + ## need the same row order as the permuted design matrices. + if (is.null(fixef_id)) { + return(NULL) + } + if (is.list(fixef_id)) { + out <- lapply(fixef_id, function(x) x[perm]) + names(out) <- names(fixef_id) + return(out) + } + if (is.null(dim(fixef_id))) { + return(fixef_id[perm]) + } + fixef_id[perm, , drop = FALSE] +} + .get_fixest_iv_design <- function(model) { if (!.is_fixest_iv_second_stage(model)) { stop("`.get_fixest_iv_design()` requires a fixest IV second-stage model.") diff --git a/R/utils-spatial.R b/R/utils-spatial.R index c407786..dad5af0 100644 --- a/R/utils-spatial.R +++ b/R/utils-spatial.R @@ -28,6 +28,26 @@ D } +.getdistvec <- function(S1, S2, latlong) { + ## Vector of paired distances between two coordinate matrices. + S1 <- as.matrix(S1) + S2 <- as.matrix(S2) + if (!identical(dim(S1), dim(S2))) { + stop("Internal error: paired distance inputs must have matching dimensions.") + } + if (latlong) { + lon1 <- S1[, 1] * pi / 180 + lat1 <- S1[, 2] * pi / 180 + lon2 <- S2[, 1] * pi / 180 + lat2 <- S2[, 2] * pi / 180 + dlon <- 0.5 * (lon1 - lon2) + dlat <- 0.5 * (lat1 - lat2) + asin(sqrt(sin(dlat)^2 + cos(lat1) * cos(lat2) * sin(dlon)^2)) / pi + } else { + sqrt(rowSums((S1 - S2)^2)) + } +} + .getc0fromavc <- function(dist, avc0) { ## Solve for exponential kernel parameter c0 such that ## mean(exp(-c0 * dist)) == avc0, using bisection. @@ -139,13 +159,210 @@ Oms } -.setOmsWfin <- function(distmat, avc0) { - n <- nrow(distmat) - distv <- .lvech(distmat) +.normalize_s <- function(S, latlong) { + ## This follows the original Princeton large-n branch closely so the + ## approximation uses the same row order as the reference code. + S <- as.matrix(S) + + if (!latlong) { + S <- sweep(S, 2, colMeans(S), `-`) + rot <- eigen(crossprod(S), symmetric = TRUE)$vectors + S <- S %*% rot + if (max(S[, 1]) != max(abs(S[, 1]))) { + S <- -S + } + perm <- do.call(order, as.data.frame(S)) + S <- S[perm, , drop = FALSE] + S <- sweep(S, 2, apply(S, 2, min), `-`) + smax <- max(S) + if (is.finite(smax) && smax > 0) { + S <- S / smax + } + } else { + S[, 1] <- S[, 1] - mean(S[, 1]) + S[, 1] <- ((S[, 1] + 180) %% 360) - 180 + perm <- order(S[, 2], S[, 1]) + S <- S[perm, , drop = FALSE] + } + + list(coords = S, perm = perm) +} + +.next_u <- function(random_t) { + random_t <- (64389 * random_t + 1) %% 2^32 + list(value = random_t / 2^32, state = random_t) +} + +.jumble_s <- function(S, m, random_t) { + n <- nrow(S) + for (i in seq_len(m)) { + nxt <- .next_u(random_t) + random_t <- nxt$state + j <- floor(nxt$value * n) + 1L + tmp <- S[j, ] + S[j, ] <- S[i, ] + S[i, ] <- tmp + } + list(coords = S, state = random_t) +} + +.ln_subset_evecs <- function(distmat, c0, qmax) { + Sig_d <- .demeanmat(exp(-c0 * distmat)) + n <- nrow(Sig_d) + if (requireNamespace("RSpectra", quietly = TRUE) && qmax < n - 1) { + RSpectra::eigs_sym(Sig_d, k = qmax, which = "LM")$vectors + } else { + eigen(Sig_d, symmetric = TRUE)$vectors[, seq_len(qmax), drop = FALSE] + } +} + +.lnset_wc0 <- function(S, avc0, qmax, minavc, latlong, + capN = 20L, m = 1000L, random_t = 1) { + n <- nrow(S) + m <- min(m, n) + ms <- vector("list", capN) + block_len <- m * (m - 1) / 2 + distvec <- numeric(capN * block_len) + + r <- S + for (i in seq_len(capN)) { + jumbled <- .jumble_s(r, m, random_t) + r <- jumbled$coords + random_t <- jumbled$state + ms[[i]] <- list( + coords = r[seq_len(m), , drop = FALSE], + distmat = .getdistmat(r[seq_len(m), , drop = FALSE], latlong) + ) + idx <- ((i - 1L) * block_len + 1L):(i * block_len) + distvec[idx] <- .lvech(ms[[i]]$distmat) + } + + c0 <- .getc0fromavc(distvec, avc0) + cmax <- .getc0fromavc(distvec, minavc) + + Wall <- matrix(0, n, capN * qmax) + for (i in seq_len(capN)) { + W0 <- .ln_subset_evecs(ms[[i]]$distmat, c0, qmax) + Wx <- matrix(0, n, qmax) + for (j in seq_len(m)) { + diff <- sweep(S, 2, ms[[i]]$coords[j, ], `-`) + v <- exp(-c0 * sqrt(rowSums(diff^2))) + Wx <- Wx + tcrossprod(v, W0[j, ]) + } + Wx <- sweep(Wx, 2, colMeans(Wx), `-`) + norms <- sqrt(colSums(Wx^2)) + norms[!is.finite(norms) | norms == 0] <- 1 + Wx <- sweep(Wx, 2, norms, `/`) + for (j in seq_len(qmax)) { + Wall[, (j - 1L) * capN + i] <- Wx[, j] + } + } + + W <- matrix(0, n, qmax) + for (i in seq_len(qmax)) { + Wx <- Wall[, seq_len(capN * i), drop = FALSE] + evec <- eigen(crossprod(Wx), symmetric = TRUE)$vectors[, 1, drop = FALSE] + W[, i] <- as.numeric(Wx %*% evec) + W[, i] <- W[, i] / sqrt(sum(W[, i]^2)) + Wall <- Wall - W[, i, drop = FALSE] %*% crossprod(W[, i], Wall) + } + + list( + W = cbind(rep(1 / sqrt(n), n), W), + c0 = c0, + cmax = cmax, + random_t = random_t + ) +} + +.raninds <- function(n, capM, random_t) { + v <- numeric(capM + 1L) + nxt <- .next_u(random_t) + random_t <- nxt$state + j <- floor(n * nxt$value) + + for (i in seq_len(capM + 1L)) { + v[i] <- j + 1L + nxt <- .next_u(random_t) + random_t <- nxt$state + j <- (j + 1 + floor(nxt$value * (n - 1))) %% n + } + + list(indices = as.integer(v), state = random_t) +} + +.lnget_Oms <- function(S, c0, cmax, W, cgridfac, + capM = 1000000L, random_t = 1, latlong = FALSE) { + nc <- .getnc(c0, cmax, cgridfac) + Oms <- vector("list", nc) + + n <- nrow(S) + inds_res <- .raninds(n, capM, random_t) + inds <- inds_res$indices + dist <- .getdistvec( + S[inds[seq_len(capM)], , drop = FALSE], + S[inds[2:(capM + 1L)], , drop = FALSE], + latlong + ) + W1 <- W[inds[seq_len(capM)], , drop = FALSE] + W2 <- W[inds[2:(capM + 1L)], , drop = FALSE] + + c <- c0 + for (i in seq_len(nc)) { + cd <- exp(-c * dist) + Oms[[i]] <- diag(ncol(W)) + + 0.5 * (n * (n - 1) / capM) * + (crossprod(W1, W2 * cd) + crossprod(W2, W1 * cd)) + c <- c * cgridfac + } + + list(Oms = Oms, state = inds_res$state) +} + +.validate_large_n_seed <- function(seed) { + if (!is.numeric(seed) || length(seed) != 1L || !is.finite(seed) || + seed < 0 || seed >= 2^32 || seed != floor(seed)) { + stop("`large_n_seed` must be a single integer-valued number in [0, 2^32).") + } + as.numeric(seed) +} + +.validate_scpc_method <- function(method) { + methods <- c("auto", "exact", "approx") + if (!is.character(method) || length(method) != 1L || !nzchar(method) || + !method %in% methods) { + stop("`method` must be one of \"auto\", \"exact\", or \"approx\".") + } + method +} + +.resolve_scpc_method <- function(n, method, large_n_threshold = 4500L) { + ## this is only the old exact-vs-approx routing rule pulled into one + ## place so the 4500 threshold can be tested without running scpc(). + method <- .validate_scpc_method(method) + if (identical(method, "auto")) { + if (n < large_n_threshold) "exact" else "approx" + } else { + method + } +} + +.setOmsWfin <- function(coords, avc0, latlong, method = "auto", large_n_seed = 1) { + ## keep the old large-n algorithm as-is, but let the routing and seed + ## checks happen up front so the split code stays easier to read. + method <- .validate_scpc_method(method) + large_n_seed <- .validate_large_n_seed(large_n_seed) + n <- nrow(coords) cgridfac <- 1.2 minavc <- 0.00001 + large_n_threshold <- 4500L + large_n_capN <- 20L + large_n_capM <- 1000000L + large_n_m <- 1000L + method_actual <- .resolve_scpc_method(n, method, large_n_threshold = large_n_threshold) + if (avc0 >= 0.05) { qmax <- 10 } else if (avc0 >= 0.01) { @@ -156,22 +373,81 @@ qmax <- 120 } - c0 <- .getc0fromavc(distv, avc0) - cmax <- .getc0fromavc(distv, minavc) + distmat <- NULL + coords_use <- coords + perm <- seq_len(n) + random_t <- NULL repeat { qmax <- min(qmax, n - 1) - W <- .getW(distmat, c0, qmax) - Oms <- .getOms(distmat, c0, cmax, W, cgridfac) + if (identical(method_actual, "exact")) { + distmat <- .getdistmat(coords, latlong) + distv <- .lvech(distmat) + c0 <- .getc0fromavc(distv, avc0) + cmax <- .getc0fromavc(distv, minavc) + W <- .getW(distmat, c0, qmax) + Oms <- .getOms(distmat, c0, cmax, W, cgridfac) + coords_use <- coords + perm <- seq_len(n) + random_t <- NULL + } else { + random_t <- large_n_seed + norm_s <- .normalize_s(coords, latlong) + coords_use <- norm_s$coords + perm <- norm_s$perm + ln_w <- .lnset_wc0( + coords_use, avc0, qmax, minavc, latlong, + capN = large_n_capN, m = large_n_m, random_t = random_t + ) + W <- ln_w$W + c0 <- ln_w$c0 + cmax <- ln_w$cmax + oms_res <- .lnget_Oms( + coords_use, c0, cmax, W, cgridfac, + capM = large_n_capM, random_t = ln_w$random_t, latlong = latlong + ) + Oms <- oms_res$Oms + random_t <- oms_res$state + } fin <- .setfinalW(Oms, W, qmax) if (fin$q < qmax || qmax == n - 1) break qmax <- round(qmax + qmax / 2) } + list( Wfin = fin$W, cvfin = fin$cv, Omsfin = Oms, c0 = c0, - cmax = cmax + cmax = cmax, + coords = coords_use, + perm = perm, + distmat = distmat, + method = method_actual, + large_n = identical(method_actual, "approx"), + random_state = random_t + ) +} + +.get_conditional_Oms <- function(spc, Wx, latlong, random_state) { + ## the old code did this inline in scpc(). this helper keeps the same + ## exact-vs-large-n logic in one place so core.R stays readable. + if (isTRUE(spc$large_n)) { + omsx_res <- .lnget_Oms( + spc$coords, + spc$c0, + spc$cmax, + Wx, + 1.2, + capM = 1000000L, + random_t = random_state, + latlong = latlong + ) + return(list(Oms = omsx_res$Oms, random_state = omsx_res$state)) + } + + list( + Oms = .getOms(spc$distmat, spc$c0, spc$cmax, Wx, 1.2), + random_state = random_state ) } diff --git a/man/scpc.Rd b/man/scpc.Rd index 88a4b53..6893dbb 100644 --- a/man/scpc.Rd +++ b/man/scpc.Rd @@ -17,6 +17,8 @@ scpc( cluster = NULL, ncoef = NULL, avc = 0.03, + method = "auto", + large_n_seed = 1, uncond = FALSE, cvs = FALSE ) @@ -58,6 +60,16 @@ in \code{data}. Coordinates must be constant within clusters.} for which size is controlled. Must be in \code{(0.001, 0.99)}. Default is 0.03.} +\item{method}{Character string selecting the spatial algorithm. +\code{"auto"} chooses \code{"exact"} for smaller problems and +\code{"approx"} once the large-\eqn{n} threshold is reached. +\code{"exact"} always uses the full distance matrix. +\code{"approx"} uses the original large-\eqn{n} approximation branch.} + +\item{large_n_seed}{Numeric; integer-valued seed used by the +large-\eqn{n} approximation branch. Ignored when +\code{method = "exact"}. Default is 1.} + \item{uncond}{Logical; if \code{TRUE}, report unconditional critical values only (skip the conditional adjustment of Mueller and Watson 2023). Default is \code{FALSE}.} @@ -92,6 +104,10 @@ An object of class \code{"scpc"} with components: \code{avc}.} \item{\code{cv}}{The unconditional 5\% critical value.} \item{\code{q}}{Number of spatial principal components selected.} + \item{\code{method}}{Spatial algorithm used: + \code{"exact"} or \code{"approx"}.} + \item{\code{large_n_seed}}{Seed used by the large-\eqn{n} + approximation branch.} \item{\code{call}}{The matched call.} } } diff --git a/tests/testthat/test-scpc-parity-stata.R b/tests/testthat/test-scpc-parity-stata.R index ae3c4d0..de9bb62 100644 --- a/tests/testthat/test-scpc-parity-stata.R +++ b/tests/testthat/test-scpc-parity-stata.R @@ -552,3 +552,81 @@ test_that("additional IV parity cases match Stata within tight tolerances", { } }) }) + +test_that("large-n synthetic scenarios match Stata in both unconditional and conditional modes", { + ctx <- skip_if_stata_scpc_unavailable() + skip_if( + !identical(Sys.getenv("SCPCR_RUN_STATA_HEAVY_TESTS", "false"), "true"), + "set SCPCR_RUN_STATA_HEAVY_TESTS=true to run large-n Stata parity tests" + ) + + set.seed(1) + n <- 4600 + dat <- data.frame( + y = stats::rnorm(n), + x = stats::rnorm(n), + lon = stats::runif(n, min = -125, max = -66), + lat = stats::runif(n, min = 25, max = 49) + ) + data_path <- tempfile("scpc_large_n_", fileext = ".csv") + utils::write.csv(dat, data_path, row.names = FALSE) + on.exit(unlink(data_path), add = TRUE) + + run_large_n_case <- function(id, uncond) { + stata_stats <- file.path(tempdir(), paste0("stata_stats_", id, ".csv")) + on.exit(unlink(stata_stats), add = TRUE) + + opts <- c("latlong", "avc(0.1)", "k(1)") + if (isTRUE(uncond)) { + opts <- c(opts, "uncond") + } + + stata_lines <- c( + "clear all", + "set more off", + sprintf("import delimited \"%s\", varnames(1) clear", data_path), + "regress y x, robust", + "gen s_1 = lat", + "gen s_2 = lon", + sprintf("scpc, %s", paste(opts, collapse = " ")) + ) + + stata_out <- run_stata_scpc_lines( + stata_lines, + write_cvs = FALSE, + stata_bin = ctx$stata_bin + ) + + fit_r <- stats::lm(y ~ x, data = dat) + out_r <- scpc( + model = fit_r, + data = dat, + lon = "lon", + lat = "lat", + ncoef = 2, + avc = 0.1, + method = "approx", + large_n_seed = 1, + uncond = uncond, + cvs = FALSE + ) + + expect_table_close( + stata_out$stats, + normalize_r_scpcstats(out_r), + vars = c("coef", "std_err", "t", "p", "ci_low", "ci_high"), + tolerance = c( + coef = 2e-4, + std_err = 2e-4, + t = 2e-4, + p = 2e-4, + ci_low = 2e-4, + ci_high = 2e-4 + ), + label = id + ) + } + + run_large_n_case("large_n_uncond", uncond = TRUE) + run_large_n_case("large_n_cond", uncond = FALSE) +}) diff --git a/tests/testthat/test_scpc.R b/tests/testthat/test_scpc.R index 19855d4..b03620d 100644 --- a/tests/testthat/test_scpc.R +++ b/tests/testthat/test_scpc.R @@ -1,19 +1,11 @@ -manual_fixest_iv_conditional_reference <- function(model, data, coord_cols, avc, coef_names) { +manual_fixest_iv_conditional_reference <- function(model, + data, + coord_cols, + avc, + coef_names, + method = "exact", + large_n_seed = 1) { design <- .get_fixest_iv_design(model) - residualize <- function(Y) { - X <- design$X - Z <- design$Z - Y <- as.matrix(Y) - if (!is.null(design$fixef_id)) { - X <- as.matrix(fixest::demean(X, f = design$fixef_id, nthreads = 1L)) - Z <- as.matrix(fixest::demean(Z, f = design$fixef_id, nthreads = 1L)) - Y <- as.matrix(fixest::demean(Y, f = design$fixef_id, nthreads = 1L)) - } - qrZ <- qr(Z) - PzX <- qr.fitted(qrZ, X) - PzY <- qr.fitted(qrZ, Y) - Y - X %*% solve(crossprod(X, PzX), crossprod(X, PzY)) - } n <- nrow(data) coef_vec <- stats::coef(model) @@ -23,13 +15,20 @@ manual_fixest_iv_conditional_reference <- function(model, data, coord_cols, avc, if (!is.null(design$fixef_id)) { model_mat <- as.matrix(fixest::demean(model_mat, f = design$fixef_id, nthreads = 1L)) } - D <- .getdistmat(as.matrix(data[, coord_cols, drop = FALSE]), latlong = FALSE) - spc <- .setOmsWfin(D, avc) + spc <- .setOmsWfin( + as.matrix(data[, coord_cols, drop = FALSE]), + avc0 = avc, + latlong = FALSE, + method = method, + large_n_seed = large_n_seed + ) Wfin <- spc$Wfin Omsfin <- spc$Omsfin + perm <- spc$perm q <- ncol(Wfin) - 1L levs <- c(0.32, 0.10, 0.05, 0.01) cvs_uncond <- vapply(levs, function(lv) .getcv(Omsfin, q, lv), 0.0) + large_n_random_state <- spc$random_state stats_out <- matrix( NA_real_, @@ -52,20 +51,51 @@ manual_fixest_iv_conditional_reference <- function(model, data, coord_cols, avc, } wj <- as.numeric(n * bread_inv[pos, , drop = TRUE] %*% t(S)) + coef_vec[[coef_name]] + wj_perm <- wj[perm] tau_u <- as.numeric( - sqrt(q) * crossprod(Wfin[, 1], wj) / - sqrt(sum((t(Wfin[, -1]) %*% wj)^2)) + sqrt(q) * crossprod(Wfin[, 1], wj_perm) / + sqrt(sum((t(Wfin[, -1]) %*% wj_perm)^2)) ) se <- as.numeric( - sqrt(sum((t(Wfin[, -1]) %*% wj)^2)) / (sqrt(q) * sqrt(n)) + sqrt(sum((t(Wfin[, -1]) %*% wj_perm)^2)) / (sqrt(q) * sqrt(n)) ) p_u <- .maxrp(Omsfin, q, abs(tau_u) / sqrt(q))$max - xj <- as.numeric(n * bread_inv[pos, , drop = TRUE] %*% t(model_mat)) + xj_raw <- as.numeric(n * bread_inv[pos, , drop = TRUE] %*% t(model_mat)) + if (spc$large_n) { + residualize <- .make_iv_residualizer( + design$X[perm, , drop = FALSE], + design$Z[perm, , drop = FALSE], + fixef_id = .permute_fixef_id(design$fixef_id, perm) + ) + xj <- xj_raw[perm] + } else { + residualize <- .make_iv_residualizer( + design$X, + design$Z, + fixef_id = design$fixef_id + ) + xj <- xj_raw + } xjs <- sign(xj) Wx <- .orthogonalize_W_iv(Wfin, xj, xjs, residualize = residualize) - Omsx <- .getOms(D, spc$c0, spc$cmax, Wx, 1.2) + if (spc$large_n) { + omsx_res <- .lnget_Oms( + spc$coords, + spc$c0, + spc$cmax, + Wx, + 1.2, + capM = 1000000L, + random_t = large_n_random_state, + latlong = FALSE + ) + Omsx <- omsx_res$Oms + large_n_random_state <- omsx_res$state + } else { + Omsx <- .getOms(spc$distmat, spc$c0, spc$cmax, Wx, 1.2) + } p_c <- .maxrp(Omsx, q, abs(tau_u) / sqrt(q))$max cvs_cond <- vapply(levs, function(lv) .getcv(Omsx, q, lv), 0.0) cv <- max(spc$cvfin, cvs_cond[[3]]) @@ -233,6 +263,117 @@ test_that("scpc validates cluster input and coordinate mode selection", { ) }) +test_that("scpc validates method and large_n_seed", { + dat <- make_scpc_data() + fit <- stats::lm(y ~ x, data = dat) + + expect_error( + scpc( + fit, + data = dat, + coords_euclidean = c("lon", "lat"), + avc = 0.1, + method = "bad", + uncond = TRUE + ), + "`method` must be one of \"auto\", \"exact\", or \"approx\"" + ) + expect_error( + scpc( + fit, + data = dat, + coords_euclidean = c("lon", "lat"), + avc = 0.1, + large_n_seed = 1.5, + uncond = TRUE + ), + "`large_n_seed` must be a single integer-valued number in \\[0, 2\\^32\\)" + ) +}) + +test_that("scpc method override can force exact or approximation branches", { + dat_small <- make_scpc_data(n = 40) + fit_small <- stats::lm(y ~ x, data = dat_small) + + out_auto_small <- scpc( + fit_small, + data = dat_small, + coords_euclidean = c("lon", "lat"), + ncoef = 1, + avc = 0.1, + method = "auto", + large_n_seed = 1, + uncond = TRUE + ) + out_exact_small <- scpc( + fit_small, + data = dat_small, + coords_euclidean = c("lon", "lat"), + ncoef = 1, + avc = 0.1, + method = "exact", + large_n_seed = 1, + uncond = TRUE + ) + out_approx_small_a <- scpc( + fit_small, + data = dat_small, + coords_euclidean = c("lon", "lat"), + ncoef = 1, + avc = 0.1, + method = "approx", + large_n_seed = 1, + uncond = TRUE + ) + out_approx_small_b <- scpc( + fit_small, + data = dat_small, + coords_euclidean = c("lon", "lat"), + ncoef = 1, + avc = 0.1, + method = "approx", + large_n_seed = 1, + uncond = TRUE + ) + + expect_identical(out_auto_small$method, "exact") + expect_identical(out_exact_small$method, "exact") + expect_identical(out_approx_small_a$method, "approx") + expect_equal(out_auto_small$scpcstats, out_exact_small$scpcstats, tolerance = 1e-10) + expect_equal(out_auto_small$c0, out_exact_small$c0, tolerance = 1e-10) + expect_equal(out_auto_small$cv, out_exact_small$cv, tolerance = 1e-10) + expect_equal(out_auto_small$q, out_exact_small$q) + expect_equal(out_approx_small_a$scpcstats, out_approx_small_b$scpcstats, tolerance = 1e-12) + expect_equal(out_approx_small_a$c0, out_approx_small_b$c0, tolerance = 1e-12) + expect_equal(out_approx_small_a$cv, out_approx_small_b$cv, tolerance = 1e-12) + expect_equal(out_approx_small_a$q, out_approx_small_b$q) + expect_true(all(is.finite(out_approx_small_a$scpcstats))) +}) + +test_that("scpc routes auto to the large-n path in a full run", { + skip_if( + !identical(Sys.getenv("SCPCR_RUN_HEAVY_TESTS", "false"), "true"), + "set SCPCR_RUN_HEAVY_TESTS=true to run the full large-n routing integration test" + ) + + dat_large <- make_scpc_data(n = 4500) + fit_large <- stats::lm(y ~ x, data = dat_large) + + out_auto_large <- scpc( + fit_large, + data = dat_large, + coords_euclidean = c("lon", "lat"), + ncoef = 1, + avc = 0.1, + method = "auto", + large_n_seed = 1, + uncond = TRUE + ) + + expect_identical(out_auto_large$method, "approx") + expect_true(all(is.finite(out_auto_large$scpcstats))) +}) + test_that("scpc does not accept the old k argument", { dat <- make_python_scpc_data() fit <- stats::lm(y ~ x, data = dat) @@ -407,6 +548,65 @@ test_that("conditional SCPC for fixest IV matches a direct 2SLS residualization }) }) +test_that("conditional SCPC for fixest IV matches a direct large-n reference", { + skip_if_not_installed("fixest") + + with_fixest_single_thread({ + set.seed(2031) + n_fe <- 30 + t_per_fe <- 6 + n <- n_fe * t_per_fe + fe <- rep(seq_len(n_fe), each = t_per_fe) + z <- stats::rnorm(n) + w <- stats::rnorm(n) + u <- stats::rnorm(n) + x <- 0.8 * z + 0.4 * w + 0.6 * u + stats::rnorm(n, sd = 0.2) + y <- 1 + 1.1 * x + 0.5 * w + stats::rnorm(n_fe)[fe] + u + dat <- data.frame( + y = y, + x = x, + z = z, + w = w, + fe = fe, + coord_x = runif(n), + coord_y = runif(n) + ) + + fit <- fixest::feols(y ~ w | fe | x ~ z, data = dat) + out <- scpc( + fit, + dat, + coords_euclidean = c("coord_x", "coord_y"), + ncoef = 2, + avc = 0.05, + method = "approx", + large_n_seed = 1, + uncond = FALSE, + cvs = TRUE + ) + ref <- manual_fixest_iv_conditional_reference( + fit, + dat, + coord_cols = c("coord_x", "coord_y"), + avc = 0.05, + coef_names = c("fit_x", "w"), + method = "approx", + large_n_seed = 1 + ) + + stats_out <- as.data.frame(out$scpcstats) + stats_out$term <- rownames(out$scpcstats) + stats_out <- stats_out[stats_out$term %in% rownames(ref$stats), , drop = FALSE] + stats_mat <- as.matrix(stats_out[, c("Coef", "Std_Err", "t", "P>|t|", "2.5 %", "97.5 %")]) + rownames(stats_mat) <- stats_out$term + + expect_identical(out$method, "approx") + expect_equal(out$large_n_seed, 1) + expect_equal(stats_mat[rownames(ref$stats), , drop = FALSE], ref$stats, tolerance = 1e-8) + expect_equal(out$scpccvs[rownames(ref$cvs), , drop = FALSE], ref$cvs, tolerance = 1e-8) + }) +}) + test_that("conditional SCPC for fixest FE matches lm with explicit FE dummies", { skip_if_not_installed("fixest") diff --git a/tests/testthat/test_scpc_fixest_iv_smoke.R b/tests/testthat/test_scpc_fixest_iv_smoke.R index 276b2ab..eb7f582 100644 --- a/tests/testthat/test_scpc_fixest_iv_smoke.R +++ b/tests/testthat/test_scpc_fixest_iv_smoke.R @@ -219,3 +219,53 @@ test_that("scpc handles clustered fixest IV models with absorbed fixed effects", expect_true(all(is.finite(out$scpccvs))) }) }) + +test_that("scpc handles clustered fixest IV models in the large-n approximation path", { + skip_if_not_installed("fixest") + + with_fixest_single_thread({ + set.seed(2006) + n_cluster <- 40 + cl_size <- 3 + n <- n_cluster * cl_size + cl <- rep(seq_len(n_cluster), each = cl_size) + z <- stats::rnorm(n) + w <- stats::rnorm(n) + u <- stats::rnorm(n) + x <- 0.8 * z + 0.2 * w + u + y <- 1 + 1.15 * x + 0.35 * w + u + lon_cl <- runif(n_cluster) + lat_cl <- runif(n_cluster) + dat <- data.frame( + y = y, + x = x, + w = w, + z = z, + cl = cl, + lon = lon_cl[cl], + lat = lat_cl[cl] + ) + + fit <- fixest::feols(y ~ w | x ~ z, data = dat) + out <- expect_no_warning( + scpc( + fit, + dat, + coords_euclidean = c("lon", "lat"), + cluster = "cl", + ncoef = 2, + avc = 0.1, + method = "approx", + large_n_seed = 1, + uncond = FALSE, + cvs = TRUE + ) + ) + + expect_s3_class(out, "scpc") + expect_identical(out$method, "approx") + expect_equal(out$large_n_seed, 1) + expect_true(all(is.finite(out$scpcstats))) + expect_true(all(is.finite(out$scpccvs))) + }) +}) diff --git a/tests/testthat/test_set_oms_wfin.R b/tests/testthat/test_set_oms_wfin.R index 3bea900..31f568f 100644 --- a/tests/testthat/test_set_oms_wfin.R +++ b/tests/testthat/test_set_oms_wfin.R @@ -1,20 +1,102 @@ -test_that("setOmsWfin returns a coherent spatial setup", { - distmat <- matrix( - c(0, 1, 2, - 1, 0, 1, - 2, 1, 0), - nrow = 3, - byrow = TRUE +test_that("setOmsWfin can force the exact path", { + coords <- cbind( + coord_x = c(0, 1, 0.5, 1.5, 2.0), + coord_y = c(0, 0, 1.0, 1.0, 1.5) ) - out <- .setOmsWfin(distmat, 0.1) + out <- .setOmsWfin( + coords, + avc0 = 0.1, + latlong = FALSE, + method = "exact", + large_n_seed = 1 + ) expect_type(out, "list") - expect_true(all(c("Wfin", "cvfin", "Omsfin", "c0", "cmax") %in% names(out))) - expect_equal(nrow(out$Wfin), nrow(distmat)) + expect_true(all( + c("Wfin", "cvfin", "Omsfin", "c0", "cmax", "coords", "perm", "distmat", "method", "large_n", "random_state") %in% + names(out) + )) + expect_identical(out$method, "exact") + expect_false(out$large_n) + expect_null(out$random_state) + expect_equal(out$perm, seq_len(nrow(coords))) + expect_equal(out$coords, coords) + expect_equal(out$distmat, .getdistmat(coords, latlong = FALSE), tolerance = 1e-12) + expect_equal(nrow(out$Wfin), nrow(coords)) expect_true(ncol(out$Wfin) >= 2) expect_gt(out$cvfin, 0) expect_gt(out$c0, 0) expect_gt(out$cmax, out$c0) expect_true(length(out$Omsfin) >= 2) }) + +test_that("setOmsWfin can force the large-n approximation path", { + coords <- cbind( + coord_x = seq(0, 1, length.out = 40), + coord_y = seq(1, 0, length.out = 40) + ) + + out <- .setOmsWfin( + coords, + avc0 = 0.1, + latlong = FALSE, + method = "approx", + large_n_seed = 1 + ) + + expect_identical(out$method, "approx") + expect_true(out$large_n) + expect_null(out$distmat) + expect_equal(sort(out$perm), seq_len(nrow(coords))) + expect_equal(nrow(out$coords), nrow(coords)) + expect_equal(nrow(out$Wfin), nrow(coords)) + expect_true(all(is.finite(out$coords))) + expect_true(all(is.finite(out$Wfin))) + expect_true(all(vapply(out$Omsfin, function(Om) all(is.finite(Om)), logical(1)))) + expect_true(is.numeric(out$random_state)) + expect_length(out$random_state, 1) + expect_true(is.finite(out$random_state)) +}) + +test_that("setOmsWfin auto chooses exact on small problems", { + coords <- cbind( + coord_x = c(0, 1, 0.5, 1.5, 2.0), + coord_y = c(0, 0, 1.0, 1.0, 1.5) + ) + + out <- .setOmsWfin( + coords, + avc0 = 0.1, + latlong = FALSE, + method = "auto", + large_n_seed = 1 + ) + + expect_identical(out$method, "exact") + expect_false(out$large_n) + expect_equal(out$perm, seq_len(nrow(coords))) +}) + +test_that("resolve_scpc_method routes auto at the large-n threshold", { + expect_identical(.resolve_scpc_method(4499L, "auto"), "exact") + expect_identical(.resolve_scpc_method(4500L, "auto"), "approx") + expect_identical(.resolve_scpc_method(100L, "exact"), "exact") + expect_identical(.resolve_scpc_method(100L, "approx"), "approx") +}) + +test_that("setOmsWfin validates method and large_n_seed", { + coords <- cbind( + coord_x = c(0, 1, 0.5, 1.5, 2.0), + coord_y = c(0, 0, 1.0, 1.0, 1.5) + ) + + expect_error( + .setOmsWfin(coords, avc0 = 0.1, latlong = FALSE, method = "bad", large_n_seed = 1), + "`method` must be one of \"auto\", \"exact\", or \"approx\"" + ) + expect_error( + .setOmsWfin(coords, avc0 = 0.1, latlong = FALSE, method = "exact", large_n_seed = 1.5), + "`large_n_seed` must be a single integer-valued number in \\[0, 2\\^32\\)" + ) +}) From c377421b1acbe46b0e9654a727cc3209d37e615d Mon Sep 17 00:00:00 2001 From: DGoettlich Date: Wed, 22 Apr 2026 01:56:04 +0200 Subject: [PATCH 2/2] [chores] corrected urls --- DESCRIPTION | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index a0a0bf5..b6e5ebe 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -11,8 +11,8 @@ Description: Implements spatial correlation-robust procedures for Mueller and Watson (2022) and Mueller and Watson (2023) . License: MIT + file LICENSE -URL: https://github.com/pdavidboll/scpcR -BugReports: https://github.com/pdavidboll/scpcR/issues +URL: https://github.com/spatial-spur/scpcR +BugReports: https://github.com/spatial-spur/scpcR/issues Encoding: UTF-8 LazyData: true RoxygenNote: 7.3.2