Getting Help and Troubleshooting

Q1: My empirical coverage is slightly below \(1 - \alpha\) (e.g. 88.5% instead of 90%). Is my model invalid?

Not necessarily. Conformal prediction guarantees valid marginal coverage in expectation: \[P\left(Y_{n+1} \in \widehat{C}(S_{n+1})\right) \ge 1 - \alpha\] For any finite test evaluation set of size \(m\), the observed empirical coverage follows a Binomial distribution \(\text{Binom}(m, 1 - \alpha)\) with standard error: \[\text{SE} = \sqrt{\frac{\alpha(1 - \alpha)}{m}}\] For example, with \(m = 100\) and nominal coverage \(90\%\) (\(\alpha = 0.10\)), the standard error is \(\sqrt{0.1 \times 0.9 / 100} = 0.03\) (\(3\%\)). An observed coverage between \(84\%\) and \(96\%\) falls well within normal \(95\%\) binomial sampling variation. You can compute the exact binomial confidence interval via:

report <- coverage_report(out, y_test)
binom.test(round(report$coverage * length(y_test)), length(y_test), p = 1 - out$alpha)

Q2: Should I use scp_geostatistical() or scp_areal()?

• Use scp_geostatistical() when observations have continuous coordinates \((x, y)\) or \((\text{lon}, \text{lat})\), such as weather stations, air quality sensors, or soil core samples.
• Use scp_areal() when data are aggregated over discrete spatial units (counties, municipalities, census tracts, or lattice grid cells) where relationships are defined by an adjacency matrix \(W\).

Q3: How do spatial bandwidth and \(k\)-NN affect interval widths?

• Bandwidth (\(h\)): In scp_geostatistical(), weights decay as \(w_i = \exp(-\|s_0 - s_i\|^2 / (2h^2))\).
  • A smaller \(h\) restricts calibration to immediately proximate neighbors, producing localized intervals that adapt to local heteroscedasticity, but may be noisy if local data are sparse.
  • A larger \(h\) borrows strength from more calibration points, yielding smoother, more stable intervals that approach standard split conformal prediction.
  • The default (bandwidth = NULL) employs the median pairwise distance heuristic, which provides a robust balance for isotropic spatial domains.
• \(k\)-Nearest Neighbors (k_neighbors): In spconform 0.1.2+, for large calibration datasets (\(N > 1000\)), setting e.g. k_neighbors = 50 restricts kernel weighting to the 50 closest calibration stations, accelerating computation up to 50x while preserving asymptotic validity.

Q4: Does spconform support spatio-temporal data?

Yes. Supply t_train and t0 alongside coordinates in scp_geostatistical():

out <- scp_geostatistical(
  s_train = s_train, y_train = y_train, s0 = s_test,
  pred_fun = my_st_predictor,
  t_train = t_train, t0 = t_test,
  temporal_bandwidth = 7,  # e.g. 7 days
  alpha = 0.1
)

1. Marginal Coverage & Interval Width

• coverage: Fraction of true values falling within \([L(s_0), U(s_0)]\). Should be \(\ge 1 - \alpha\).
• mean_width vs median_width: If the mean width is substantially larger than the median, spatial outliers or high-uncertainty boundary clusters are inflating interval widths.

2. Winkler Interval Score (WIS)

A strictly proper scoring rule balancing interval sharpness and miscoverage penalties: \[\text{WIS}_\alpha(l, u, y) = (u - l) + \frac{2}{\alpha}(l - y)\mathbf{1}(y < l) + \frac{2}{\alpha}(y - u)\mathbf{1}(y > u)\] Lower scores indicate superior predictive performance (narrower intervals with minimal violations).

3. Spatial Quadrant Stratification

Partitions the 2D spatial bounding box into \(n \times n\) quadrants (default \(4 \times 4 = 16\)). * Inspect whether any specific quadrant suffers severe localized under-coverage (e.g. coverage \(< 75\%\)). * Local deficits often highlight unobserved spatial covariates or sharp topographical transitions.

4. Convex Hull Boundary Proximity

Examines whether intervals near the edge of the spatial domain experience degraded coverage due to spatial extrapolation. diagnose() partitions test points into Interior vs. Boundary units.

5. Spatial Autocorrelation (Moran’s \(I\))

Computes Moran’s \(I\) on absolute conformal residuals: * \(I \approx 0\) (\(p > 0.05\)): Nonconformity scores are spatially independent; local exchangeability holds well. * \(I > 0\) (\(p < 0.05\)): Prediction errors cluster spatially, indicating unmodeled spatial autocorrelation in the underlying point predictor.

Recipe 1: Random Forest (ranger)

library(ranger)

pred_rf <- function(s_train, y_train, s_new) {
  df_tr <- data.frame(x = s_train[, 1], y = s_train[, 2], target = y_train)
  df_te <- data.frame(x = s_new[, 1],   y = s_new[, 2])
  
  fit <- ranger(target ~ x + y, data = df_tr, num.trees = 100)
  predict(fit, data = df_te)$predictions
}

Recipe 2: Generalized Additive Model (mgcv::gam)

library(mgcv)

pred_gam <- function(s_train, y_train, s_new) {
  df_tr <- data.frame(x = s_train[, 1], y = s_train[, 2], target = y_train)
  df_te <- data.frame(x = s_new[, 1],   y = s_new[, 2])
  
  fit <- gam(target ~ s(x, y, k = 20), data = df_tr)
  as.numeric(predict(fit, newdata = df_te))
}

Recipe 3: Ordinary Spatial Kriging (gstat)

library(sp)
library(gstat)

pred_kriging <- function(s_train, y_train, s_new) {
  sp_tr <- SpatialPointsDataFrame(coords = s_train, data = data.frame(y = y_train))
  sp_te <- SpatialPoints(coords = s_new)
  
  vgm_mod <- fit.variogram(variogram(y ~ 1, sp_tr), vgm("Exp"))
  krig_out <- krige(y ~ 1, sp_tr, sp_te, model = vgm_mod, debug.level = 0)
  krig_out$var1.pred
}

1. Coordinate & Dimension Mismatches

2. Missing & Non-Finite Values (NA, NaN, Inf)

3. Custom Predictor (pred_fun) Return Types & Single-Row Drops

4. Coordinate Reference System (CRS) & Distance Scale Discrepancies

5. Areal Adjacency Matrix Structural Errors

6. Calibration Sample Size vs. Nominal Miscoverage (\(\alpha\))

7. Evaluation & Diagnostic Alignment

8. Graphics Layout & Multi-Panel Plot Devices