Two regressions of 0.11.0 in vcovSpHAC.felm(), found by a replication
pipeline within days of the release, are fixed. Results for fixest fits and
for felm fits with explicit unit/time are unchanged.
vcovSpHAC.felm() no longer infers unit and time from the absorbed
fixed effects. 0.11.0 used the first two absorbed effects as unit and
time whenever both arguments were omitted, which silently restricted
spatial pairs to within the second effect's groups (a region, an ethnic
group) and moved Conley standard errors by percentages for anyone relying
on the default; every earlier release, and the fixest method, treated
such a fit as one cross-sectional block. That behaviour is restored:
with both arguments omitted every pair within the cutoff enters and each
row is its own unit. Absorbed-effect names remain valid explicit
unit/time values. Regression test added. Reported from a replication
pipeline whose published values were computed with the pre-0.11.0
behaviour.vcovSpHAC.felm() no longer takes coordinate columns directly from a data
frame it found by name in the caller's or the formula's environment: when
data = is not supplied, rows are always recovered through the model-frame
alignment (expand.model.felm()), as in versions before 0.11.0. In 0.11.0 a
frame with the same name as the one in the felm call and the same row
count, for example a per-outcome frame reassigned in a loop or a pooled
frame in the global environment, could be used unaligned, moving Conley
standard errors silently, or produce NA coordinates that the new validation
then rejected. A user-supplied data = with as many rows as the fit is
still taken as aligned. Pass data = explicitly inside helper functions
and loops. Regression test added.vcovSpHAC() generic's help page gained a Value
section describing the returned matrix (CRAN request).Engine extraction, a Stata port sharing the same C++ engine, and a review round (six independent code reviews of the engine, the R layer, the Stata command, the plugin build, the reghdfe proposal, and the validation coverage) whose findings are fixed below.
felm fits got the wrong default ssc. lfe stores the
cluster-adjusted degrees of freedom (number of clusters minus one) in
df.residual when a fit carries a cluster specification, and the
small-sample scale n / df.residual inherited that value: on a 20-row
fixture the covariance was five times too large. The scale now uses the
regression's own residual degrees of freedom (N - p), so clustered and
unclustered fits of the same model give identical Conley covariances.vcovSpHAC.felm ignored the unit and time arguments and always
used the first two absorbed fixed effects. The names are now matched
against names(reg$fe) (with the model data as a fallback), a missing
name is an error, and the first-two-FE default applies only when nothing
was passed. A model with | region + unit + time used the wrong keys
before.felm fits); a time key that
does not parse is accepted for spatial blocking only and is an error when
lag_cutoff > 0.time without unit is honoured. The fixest method silently dropped
the time blocking when no unit was given; it now blocks by time with each
row as its own unit, matching the Stata command. lag_cutoff > 0 without
a unit is an error.fixest::obs(), so models extracted
from split = estimations and fits with several stacked selections align
their coordinates correctly. fixest_multi objects, lean = TRUE feols
fits, and fixed-effects-only models now fail with clear messages instead
of misleading ones.lat, lon, unit, time values, latitudes
outside [-90, 90], and longitudes outside [-180, 360] are rejected in R
before anything reaches C++; the engine repeats the checks (including
scores and cutoff) for its other front-ends. A NA time used to split
co-located observations into separate blocks.lat and latitude (or any two aliases) present in the data is now
an auto-detection error, as documented.data = either as the full original data (rows
addressed through fixest::obs()) or as a frame already aligned with the
fit, such as the fitted model frame after subset = or NA removal. A frame
with as many rows as the original data is always read as original data (so
subsets that only permute rows stay correct); any other row count is an
error.... are errors (typos such as psd_fxi = were
swallowed). maxobsmem warns once that it is ignored; neighbor = "band"
is documented as deprecated but still works.lfe fits made through the k-class path (any kclass = argument, even
kclass = 1) are rejected: lfe stores the raw endogenous regressors for
them instead of the projected 2SLS design, so the sandwich would be silently
wrong (a kclass = 1 fit gave 0.0028 where 2SLS gives 0.0097). Only OLS and
ordinary 2SLS are supported.ncores must be a single finite positive integer. The default remains all
logical cores, but getOption("fastconley.ncores") overrides it and a
true-ish _R_CHECK_LIMIT_CORES_ (CRAN's check convention) caps it at 2.Rcpp::checkUserInterrupt() in R,
SF_poll in the Stata plugin) and joins its workers before raising.std::thread
constructor used to terminate the process), ncores is capped at 1024,
and the sub-metre-cutoff cell cap no longer collapses all points into one
cell.dlat, dlon finite and positive, finite
lat0 in range): a zero longitude step used to spin forever.sin^2(theta/2) = |u_i - u_j|^2 / 4, distance 2R asin(sqrt a)), and the
spherical Bartlett weight uses the same asin form instead of acos(dot),
which is ill-conditioned at small angles. No per-pair sin/atan2 calls
remain, which speeds up the pairwise engine everywhere (about 2x on Linux,
6x on mingw-built Windows binaries, whose libm implements those functions
in x87 microcode). Results move at the 1e-14 relative level versus 0.10.0
on Bartlett haversine/spherical configurations.src/conley_core.h
(namespace conley, CONLEY_CORE_VERSION 0.11.1, bumped on every change
that alters results or entry points), shared with the Stata plugin. It
compiles without R against plain Armadillo under C++14
(tests/manual/core_standalone_check.cpp + .R verify bit-identity).std::thread pool inside
the engine with the same deterministic chunked reduction, so ncores
still never affects results. Depends: R (>= 4.0) declared for the
posix-threads Rtools toolchain on Windows.Three optimisation reviews (2026-09-05) led to preparation and engine changes that leave every result bit-identical (the 60-configuration bitwise battery, the standalone header check, and the plugin golden check all pass unchanged):
fepois/feglm, felm and fixest IV, the serial HAC, and pixel
aggregation; 1-vs-2-core identical() checks for every engine path; and
regression tests for each item above.tests/manual/bitwise-battery.R gained write-ref/check modes with a
committed reference (bitwise-battery-reference.txt), and
tests/manual/core_edge_probes.cpp exercises the edge cases in C++.R CMD check --as-cran workflow runs on Linux, macOS,
and Windows.method = "auto" as the default, chord versus great-circle wording, the
balanced-panel error (not a fallback), and the pixel-aggregation distance
bound (sqrt(2) * pixel for a pair, not pixel / 2).stata/ holds fastconley, a reghdfe-family command with a pure-Mata
engine and a compiled plugin built from src/conley_core.h for Linux,
Windows, and macOS (see stata/README.md and stata/CHANGELOG.md), and
stata/upstream/ two alternative reghdfe proposals: a generic
vce(external PROVIDER, ...) hook with fastconley as its first provider
(primary), and a self-contained native vce(conley ...) patch.CRAN preparation.
DESCRIPTION: package names quoted, Conley (1999) DOI added, copyright
holder role declared; LICENSE updated to match..Rbuildignore keeps development-only files (the agent instructions,
notes/, tests/manual/, stata/, .github/, and cran-comments.md)
out of the source tarball.felm fixture in the balanced-panel validation test no
longer emits "NaNs produced" warnings.GLM support and documented IV support.
vcovSpHAC.fixest now accepts GLM fits (any feglm family, including
fepois). The variance is the M-estimation sandwich H^{-1} B H^{-1},
built from the maximum-likelihood score matrix and inverse Hessian that
fixest stores on every (non-lean) fit — the same construction
fixest's own vcov_conley() uses for GLMs, verified against it at the
distance-formulation tolerance and against an exact same-distance
yardstick at ~1e-15. No estimation flag is needed (demeaned = TRUE is
only required for feols); weights, offsets, and the fixed-effect
profiling are already folded into the stored scores. Because the scores
ride through the existing engines unchanged, everything composes: pairwise
and grid/FFT engines, pixel aggregation, ssc, psd_fix, and — beyond
what fixest offers — the panel spatial + serial HAC via
lag_cutoff, now available for Poisson/GLM panels. lean = TRUE fits
(no stored scores) and femlm()/feNmlm() fits are rejected with clear
errors.
IV fits have in fact always produced the correct 2SLS Conley sandwich
through both methods — lfe and fixest store the projected
(second-stage) design in cX / X_demeaned and the structural residuals
in residuals, which is exactly what the score construction needs. This
is now documented and covered by the validation suite: felm IV and
feols IV agree with each other and with the exact yardstick at ~1e-15,
including weighted IV, multiple endogenous regressors, and IV panels with
serial HAC. See tests/manual/test-glm-iv-parity.R.
fixest feature parity: weighted fits, small-sample correction, PSD repair,
lat/lon auto-detection. Breaking: ssc and psd_fix default to TRUE
to match fixest's defaults out of the box; pass
ssc = FALSE, psd_fix = FALSE to reproduce earlier fastconley versions and
rbluhm/conley bit-for-bit.
vcovSpHAC now accepts weighted felm() and feols() fits. The meat
scores become s_i = w_i * e_i * x_i and the bread (X'WX)^{-1} — the
formula fixest's own weighted Conley vcov uses (verified exactly against
it with a self-pairs-only cutoff, rel. err ~1e-15). Weights enter only the
scores and the bread, so every engine — pairwise, grid/FFT, serial HAC,
pixel aggregation, balanced CSR reuse — works unchanged. Note lfe stores
sqrt(w) on the fit; vcovSpHAC squares it back.
Scales the variance matrix by n / (n - K), with K counting all
estimated parameters including absorbed fixed-effect levels (taken from
the fit's residual degrees of freedom). This is exactly fixest's default
Conley correction — its cluster adjustment (G.adj / cluster.adj) is a
no-op for Conley vcovs, so this one factor reproduces fixest defaults.
ssc = FALSE applies no correction, matching rbluhm/conley and previous
fastconley versions.
Conley spatial kernels do not guarantee a PSD variance matrix. With
psd_fix = TRUE (default, as in fixest) negative eigenvalues are
clamped to 1e-16 — the same semantics as fixest's vcov_fix — with a
warning when the fix noticeably changed the matrix (> 1e-8). With
psd_fix = FALSE the matrix is returned as computed and a warning is
emitted when it is noticeably non-PSD.
lat and lon now default to NULL and are auto-detected from the
data's column names (lat/latitude and lon/long/longitude/lng,
case-insensitive exact matches). A message reports the pick; ambiguous or
missing matches error with instructions.
ssc = FALSE, psd_fix = FALSE (the C++ engines and the default-off code
paths are untouched).fixest at self-pairs-only cutoff;
~1e-10 against a dense brute-force reference at 200 km for both
kernels; felm and fixest paths agree to ~1e-10.tests/testthat/test-weights-ssc.R (14 tests).Grid engine, part two: bartlett support (ring-FFT) and dateline wrap.
On a lattice the bartlett weight varies with the longitude offset, so the
per-ring-pair inner sum is a true 1D convolution rather than a boxcar.
FastGridMeat computes it via FFT (arma::fft): per ring pair, the
even-symmetric weight vector's (real) spectrum multiplies cached per-ring
score spectra, with one inverse FFT per target ring. Score spectra are
cached for a sliding latitude band plus a cutoff halo, and the reduction
is deterministically chunked — results remain bit-identical across
ncores.
Weights use the same per-distance arithmetic as the pairwise engine
(atan2 haversine, acos spherical, sqrt chord), and the same-cell
distance is hard-set to 0, so agreement with the pairwise engine is
~1e-15 (haversine) to ~1e-12 (spherical/chord — inherent conditioning of
acos/sqrt near zero distance, not algorithm error). The "auto" rule
uses an FFT-aware cost model for bartlett.
v0.7.0's grid engine clamped longitude windows at the lattice edges, so
on a raster spanning the full 360° circle it silently missed pairs that
are close "the short way" across the dateline (observed ~5e-3 relative
error on a global test raster). The engine now detects when the accept
window reaches across the dateline gap and switches both kernels to
circular windows (modular prefix-sum arcs for uniform; circular
convolution with period n_col_full for bartlett) — exact, validated
against the pairwise engine. When wrap would be needed but the lon step
does not tile 360° evenly (no consistent circular lattice exists),
method = "grid" stops with an informative error and method = "auto"
falls back to the pairwise engine. Non-wrapping rasters are unaffected:
results are bitwise identical to v0.7.0 (verified on the 30-config
battery).
Workstream C2: exact grid-native meat for raster data.
For the uniform kernel on a regular lat/lon lattice (raster cell centers,
gridded covariates), the within-cutoff accept set between two latitude
rings is a longitude-index interval, so the spatial meat reduces to
sliding-window sums over per-ring prefix sums — FastGridMeat. Cost is
O(n_ring * window * n_col * k), independent of the pair count, and the
accept threshold is the same dot-product constant the pairwise engine
uses, so the result is exact (agrees to FP summation order; no
approximation anywhere for natively gridded data).
dist_fns (monotone in the chord), multiple time
blocks (panels apply it per period), sparse occupancy, duplicate cells,
and is deterministic across ncores.method = "auto" (the new default) switches to the grid engine only when
a lattice is detected, the kernel is uniform, and a flop-balance estimate
says it wins; otherwise the pairwise engine runs as before. Scattered
(non-lattice) data is unaffected. method = "grid" errors informatively
when its requirements are not met.
The bartlett ring-FFT variant (exact per the C1 study) is planned as a follow-up; bartlett rasters currently stay on the pairwise engine.
Minor-backlog items M1-M4 from notes/OPTIMIZATION_PLAN.md.
vcovSpHAC.felm gains data =: when the passed frame has the same
row count as the fit (no NAs dropped, no subset), coordinates are taken
by direct column access with no model-frame re-evaluation; otherwise the
frame overrides the call-recovered data in the aligned model-frame path.ncores. 1M-row global cross-section at 16 threads:
0.46 -> 0.21 s (scaling 2.2x -> 4.8x).Phase 2 of notes/OPTIMIZATION_PLAN.md: memory diet, deterministic
reduction, and screen work. Results remain exact (same pairs, same
weights); summation order changed, so values differ from v0.5.0 by
<= ~5e-15 relative.
All meat accumulations use a deterministic chunked reduction (fixed-size
row/block chunks, partials summed in chunk order). ncores = 1 and
ncores = 16 produce bit-identical matrices — the old multicore
tolerance caveat is gone.
vcovSpHAC passes
the aggregated scores straight through; X/e arguments remain on the
internal wrappers for convenience. Measured peak RSS: 4M-row
cross-section 2321 -> 1558 MB (-33%); 40k x 40 balanced panel
1028 -> 748 MB (-27%). Same speed, identical checksums.csr_weight = c("double", "float"): opt-in float storage for the
balanced-path bartlett weights (8 -> 4 bytes per pair, <= ~6e-8
relative error per weight). Default stays exact double.a = (1 - dot)/2); the exact a-test remains the arbiter, so
results are bit-identical. CONUS 100k / 500 km bartlett/haversine:
15.1 -> 10.5 s single-threaded (1.4x).A unit-major T*k stacked score layout that streams the balanced CSR once instead of once per period was implemented and benchmarked. It lost to the period-major layout: spatial sorting makes neighbor gathers a sliding window of k-wide rows that stays L2-resident per period; any wider stacking pushes the window past L2 and thrashes the shared L3 at high thread counts (2.6 s vs 2.0 s at 16 cores). The period-major traversal stays, now deterministic and float-capable.
Phase 1 of notes/OPTIMIZATION_PLAN.md: 3D cell-grid neighbor search.
The spatial meat's candidate enumeration now defaults to a 3D cell grid.
Points are bucketed by their unit vectors into a cubic grid whose edge is
the unit-sphere chord equivalent of the cutoff; every supported distance is
monotone in the chord, so accepted pairs are never more than one cell apart
per axis — no pole or dateline special cases. Each row scans its own cell
plus five contiguous row ranges covering the 13 forward neighbor cells:
~3–4 candidates per accepted pair, independent of geographic extent, versus
the latitude band scan's 2*L_lon/(pi*r).
Both strategies call the identical per-pair accept test, so pair sets and
weights are exactly the same; results differ only by floating-point
summation order (observed <= 6e-15 relative across the validation matrix;
neighbor = "band" remains bitwise identical to v0.4.1). The band path is
kept for one release and will be removed in v0.6.0.
Measured single-threaded speedups (FastSpatialMeat, k = 10):
Phase 0 quick wins from notes/OPTIMIZATION_PLAN.md (Q1–Q6). Numerical
results are unchanged (verified bitwise against v0.4.0 on the balanced /
general / unbalanced × kernel × distance matrix at ncores = 1).
kernel = "uniform" (it was filled but never read), and column indices
are now 32-bit. Per stored pair: 16 bytes -> 4 (uniform) / 12 (bartlett).
A guard errors if a single period exceeds 2^32 - 1 units.a-test is unchanged).balanced_pnl = TRUE unit-set validation is now O(n) vectorized
instead of two grouping passes plus per-period string materialization.TimeDist internal export (unused since the serial
HAC moved to FastSerialHacPanel).~/.R/Makevars performance-build recipe.Development build for testing a faster spatial HAC path.
FastSpatialMeat() internal Rcpp routine.n x n spatial distance matrices in vcovSpHAC() with a screened spatial edge list.S = e * X, C_i = 0.5 S_i + sum_j w_ij S_j, and S'C + C'S.kernel = "bartlett" and kernel = "uniform".dist_fn choices: haversine, spherical, chord. The upstream flatearth option was dropped — its formula was not symmetric in (i, j) and the equirectangular approximation is no longer worth a separate code path now that spherical+uniform runs as a 3D dot-product threshold with no trig.(dist_fn, kernel) dispatch. The screen + distance + kernel-weight stack is specialised at compile time for each of the six (distance, kernel) combinations, removing per-pair runtime branches.(spherical, uniform) pair inclusion now reduces to a 3D dot product compared against cos(cutoff / R) — no trig in the inner loop. (chord, uniform) uses squared-Euclidean against (cutoff / R)². The bartlett variants only pay acos/sqrt on accepted pairs.s_i = e_i · X_i are stored row-major in a single flat std::vector<double> instead of being recomputed against column-major Armadillo memory on every pair.pixel argument for score pre-aggregation (R/vcovSpHAC.R). At pixel = 0 rows that share (lat, lon) within a time period are collapsed exactly; at pixel > 0 rows are first snapped to a uniform pixel-km grid (speed/accuracy trade-off).fixest::vcov_conley (kernel = "uniform", spherical, 500 km cutoff): fastconley is now 3-4× faster at pixel = 0 and 27-62× faster at pixel = 25 across n ∈ {5k … 100k}, and scales better than fixest past 4 threads.Rcpp::compileAttributes() is recommended after further edits.