Correcting statistical uncertainty when averages live on curved spaces

Directions on a sphere and covariance matrices have geometric structures that ordinary numerical regression does not capture. In a preprint, Ryunosuke Miyazaki of Hitotsubashi University and Chang Jun Im of Seoul National University develop a way to correct smoothing bias while respecting that geometry. Their proposed Riemannian robust bias correction, or RRBC, combines a corrected estimate with an uncertainty calculation that includes the randomness introduced by the correction itself.

Why smoothing can mislead

The target is a conditional Fréchet mean: the object that minimizes expected squared distance from the responses at a given predictor value. This extends the familiar conditional average to spaces where distance, rather than ordinary addition, defines the calculation. The paper considers finite-dimensional Riemannian manifolds, allowing both predictors and responses to occupy curved spaces.

Local regression borrows information from nearby predictor values. That smoothing reduces variability but can shift the estimate away from its target. The authors show that, at bandwidths of the order conventionally used to minimize mean squared error, this bias can be as large as sampling fluctuations. Confidence sets that ignore it can therefore miss their intended coverage even as samples grow.

One response is undersmoothing: narrowing the neighborhood until bias becomes negligible. But that leaves fewer nearby observations and increases variability.

Correcting the estimate—and its uncertainty

RRBC instead estimates the leading smoothing bias. It uses logarithm maps to express observations as vectors in tangent spaces, local linear regression for the initial estimate, and a local quadratic fit to estimate the correction. An exponential map moves the corrected estimate back onto the response manifold.

Crucially, estimating a correction creates additional uncertainty. The proposed leave-one-out covariance calculation tracks how the complete corrected estimate changes when each observation is removed. This accounts for variation in the correction and its relationship with the initial fit, rather than treating the correction as known.

Under the paper’s geometric, smoothness, moment and bandwidth assumptions, the authors establish large-sample validity for confidence sets at individual predictor values and simultaneous confidence tubes across evaluation regions. A Gaussian multiplier bootstrap calibrates the simultaneous procedure while preserving dependence across locations. These are asymptotic results, not guarantees of exact coverage in every finite dataset.

What the simulations show

The main simulations use six designs: two with spherical responses and four with positive-definite matrix predictors and responses. Each design has 1,000 replications at sample sizes of 250, 500, 1,000 and 2,000. Reported simultaneous coverage concerns finite evaluation grids.

In one spherical design with 250 observations, uncorrected nominal 95% pointwise coverage was 94.0% where smoothing displacement was close to zero, but only 29.2–47.0% at four other displayed locations. RRBC with leave-one-out covariance achieved 93.7–95.0% across those five points.

The matrix simulations demonstrate why correcting the center alone is insufficient. At 250 observations, bias correction paired with the uncorrected covariance achieved simultaneous coverage of just 14.7–17.4%, compared with 84.9–90.2% for RRBC with leave-one-out covariance.

Across all six designs at 2,000 observations, the latter procedure achieved 93.6–95.4% simultaneous coverage, versus 87.4–91.5% for undersmoothing. Its confidence sets were larger. Results also depended on implementation choices: replacing the uniform kernel with an Epanechnikov kernel reduced RRBC simultaneous Wald coverage in 23 of 24 design–sample-size settings.

Two applications

The cardiovascular analysis used 840 subjects from the Autonomic Aging database. Each response was a matrix summarizing low-frequency variation in systolic blood pressure and subsequent heartbeat intervals, including their correlation. The authors fitted separate relationships with body mass index for three age strata.

Their simultaneous tests rejected equality across age strata, with a Holm-adjusted p-value of 0.00015. They did not reject within-stratum constancy across BMI or BMI-invariance of pairwise separation between strata. Those non-rejections do not establish that the relationships are constant.

The ocean analysis used 1,107 El Niño and 1,879 La Niña drifter deployments, retaining at most one observation per deployment. Estimation incorporated spherical position and velocity direction; inference targeted headings at 85 tropical Pacific locations. Tests rejected equality between phase-specific heading fields and north–south mirror symmetry within each phase. These are phase-associated patterns in the analyzed data.

What remains unresolved

The authors explicitly distinguish first-order validity from optimal coverage accuracy. Their theory does not rank the finite-sample performance of the competing covariance estimators. They also identify joint selection of the initial-fit and bias-correction bandwidths as a topic needing further theory. The preprint’s contribution is therefore a framework supported by mathematical results, simulations and applications—not a universally optimal uncertainty procedure.