spconform

# Install release version from CRAN:
install.packages("spconform")

# Or the development version from GitHub:
remotes::install_github("amjed-droid/spconform")

library(spconform)
library(sp)

data(meuse)
s <- as.matrix(meuse[, c("x", "y")])
y <- log(meuse$zinc)

# Any user-supplied point predictor
pred_fun <- function(s_train, y_train, s_new) {
  fit <- lm(y_train ~ s_train[, 1] + s_train[, 2] +
            I(s_train[, 1]^2) + I(s_train[, 2]^2))
  cbind(1, s_new[, 1], s_new[, 2],
        s_new[, 1]^2, s_new[, 2]^2) %*% coef(fit)
}

# 90% locally weighted conformal intervals
set.seed(123)
idx <- sample(nrow(s), floor(0.7 * nrow(s)))
out <- scp_geostatistical(s[idx, ], y[idx], s[-idx, ], pred_fun, alpha = 0.1)

print(out)
#> <spconform> geostatistical conformal prediction
#> Target coverage: 90.0%
#> Number of prediction points: 47

coverage_report(out, y[-idx])
#> $coverage
#> [1] 0.957
#> $mean_width
#> [1] 2.21

# spconform objects support standard S3 methods:
predict(out, interval = "prediction")
#>   fit   lwr   upr
#> 1 5.30  2.85  7.45
#> 2 5.12  2.70  7.32

# Extract residuals
residuals(out, y_true = y[-idx], type = "abs")

# Generate a comprehensive 4-panel diagnostic plot
diag <- diagnose(out, y_true = y[-idx], s_test = s[-idx])
plot(diag) # S3 plot method for visual inspection

# View S3 printed summary including Moran's I
print(diag)
#> === spconform Diagnostic Report ===
#> 
#> Marginal coverage:
#>   Empirical: 0.9574 (nominal: 0.9)
#>   Mean width: 2.2105 
#>   n = 47 , covered = 45 
#> 
#> Conditional coverage by spatial bin:
#>   Q1-1: 1.0000 (n=7, width=3.529)
#>   Q1-2: 1.0000 (n=6, width=2.975)
#>   Q4-4: 0.8889 (n=9, width=1.930)
#> 
#> Boundary effect:
#>   Near boundary:    0.9583 (n=24)
#>   Far from boundary: 0.9565 (n=23)

# Aggregate Meuse to a 6x6 grid (21 occupied cells)
xbreaks <- seq(min(meuse$x), max(meuse$x), length.out = 7)
ybreaks <- seq(min(meuse$y), max(meuse$y), length.out = 7)
meuse$cell_x  <- cut(meuse$x, xbreaks, include.lowest = TRUE, labels = FALSE)
meuse$cell_y  <- cut(meuse$y, ybreaks, include.lowest = TRUE, labels = FALSE)
meuse$cell_id <- (meuse$cell_y - 1) * 6 + meuse$cell_x

agg <- aggregate(log(zinc) ~ cell_id, data = meuse, FUN = mean)
names(agg) <- c("cell_id", "y")
cell_coords <- unique(meuse[, c("cell_id", "cell_x", "cell_y")])
agg <- merge(agg, cell_coords, by = "cell_id")
agg <- agg[order(agg$cell_id), ]

# Build adjacency matrix (Queen contiguity)
n_cells <- nrow(agg)
adj <- matrix(0, n_cells, n_cells)
for (i in 1:n_cells) {
  for (j in 1:n_cells) {
    if (i != j) {
      dx <- abs(agg$cell_x[i] - agg$cell_x[j])
      dy <- abs(agg$cell_y[i] - agg$cell_y[j])
      if (dx <= 1 && dy <= 1) adj[i, j] <- 1
    }
  }
}

# 80% neighbourhood-weighted conformal intervals
out_areal <- scp_areal(agg$y, adjacency = adj, alpha = 0.2, decay = 0.5)
coverage_report(out_areal, agg$y)
#> $coverage
#> [1] 0.81
#> $mean_width
#> [1] 1.79

library(mgcv)
library(bmstdr)

data("nysptime")
df <- nysptime[complete.cases(nysptime[, c("utmx", "utmy", "y8hrmax", "Day", "Month")]), ]
df$day_idx <- ifelse(df$Month == 7, df$Day, 31 + df$Day)

s <- as.matrix(df[, c("utmx", "utmy")])
t <- df$day_idx
s_3d <- cbind(s, t)
y <- df$y8hrmax

# Spatio-temporal GAM predictor
pred_fun_st <- function(s_train, y_train, s_new) {
  train_df <- data.frame(x = s_train[, 1], y = s_train[, 2],
                         day = s_train[, 3], z = y_train)
  fit <- gam(z ~ te(x, y, day, k = c(8, 8, 4)), data = train_df)
  new_df <- data.frame(x = s_new[, 1], y = s_new[, 2], day = s_new[, 3])
  as.numeric(predict(fit, newdata = new_df))
}

set.seed(123)
n <- nrow(s_3d)
train_idx <- sample(n, floor(0.7 * n))
out_st <- scp_geostatistical(
  s_train = s_3d[train_idx, ],
  y_train = y[train_idx],
  s0 = s_3d[-train_idx, ],
  pred_fun = pred_fun_st,
  t_train = t[train_idx],
  t0 = t[-train_idx],
  temporal_bandwidth = 5,
  alpha = 0.1,
  split = 0.5
)

coverage_report(out_st, y[-train_idx])
#> $coverage
#> [1] 0.909
#> $mean_width
#> [1] 41.7

citation("spconform")