Package {soilcnn}

ALE: how the prediction changes along each variable's range.

Description

Accumulated local effects (Apley & Zhu 2020). A variable's range is cut into bins by its quantiles at the points; within each bin, every point is moved to the bin's lower edge and to its upper edge, and the difference of the two predictions is the local effect. Accumulated across the bins and centred, it is the curve. A point is moved only across its own bin, so the model is asked about values close to the ones it saw – which is why ALE is read with correlated predictors, where a partial dependence plot would ask about combinations that do not occur.

Usage

ale_effect(variables = NULL, bins = 20L)

Arguments

Details

A point's whole patch is moved, every pixel and every window by the same amount: the neighbourhood keeps its texture and changes its level. Moving the centre pixel alone would build a spike no landscape has. A categorical – a one-hot set, or a binary map alone – has no range: its effect is the prediction with the whole patch made each class, less the prediction as it is. That is a partial dependence, which asks about a class also where it does not occur; the number of points in each class says where it does.

Value

An importance_spec, for dsm_importance(). Its importance is the spread of the curve over the points – the standard deviation of the effect at their values; flat is no effect (Greenwell et al. 2018 for the idea, on partial dependence).

Examples

ale_effect(bins = 10)

The DI of new rows of RAW predictor values, QC'd and scaled as the reference was.

Description

The DI of new rows of RAW predictor values, QC'd and scaled as the reference was.

Usage

aoa_di(aref, values)

Arguments

Value

The DI of each row; NA where a channel is missing after QC.

Examples

ex <- example_landscape()
store <- dsm_prepare(ex$profiles, target = "soc_stock", raster_dir = ex$raster_dir,
                     windows = c(3, 7), out_dir = file.path(tempdir(), "landscape"),
                     percentage = "^clay_pct$", transform = "log1p", n_cores = 1,
                     overwrite = TRUE, verbose = FALSE)
data <- dsm_load(store, windows = integer(0), verbose = FALSE)
predictors <- as.character(data$store$predictors)
qc <- utils::read.csv2(file.path(data$patch_dir, "qc_table.csv"))
qc <- qc[match(predictors, qc$predictor), ]
# a fitted model's scaling is its predictor_scaling.csv; here, the profiles' own
x <- as.matrix(data$points[, predictors])
scaling <- data.frame(predictor = predictors, center = colMeans(x),
                      scale = apply(x, 2, sd))
plan <- spatial_folds(data$store$meta, k = 3, block_size = 0.05)
aref <- aoa_reference(data$points, predictors, qc, scaling, plan)
aoa_di(aref, data$points[1:5, predictors])

The dissimilarity reference of a fitted model.

Description

The centre pixel of every point that trained or validated in the tuning plan, QC'd and scaled with the MODEL's scaling (stage 07's construction). Each point's fold is the one it was held out in, so its cross-validated DI is the distance to the nearest point OUTSIDE its fold: the dissimilarity the model that predicted it during cross-validation actually faced. That is the DI a calibration residual comes with. A point that only ever trained – the training side of a holdout – has fold NA and no cross-validated DI: it is in the reference, a neighbour to every fold, and out of the threshold.

Usage

aoa_reference(points, predictors, qc_table, scaling, plan, weights = NULL)

Arguments

Value

An aoa_reference: the DI reference, the AOA threshold, and each used point's fold and cross-validated DI (NA for a point never held out).

Examples

ex <- example_landscape()
store <- dsm_prepare(ex$profiles, target = "soc_stock", raster_dir = ex$raster_dir,
                     windows = c(3, 7), out_dir = file.path(tempdir(), "landscape"),
                     percentage = "^clay_pct$", transform = "log1p", n_cores = 1,
                     overwrite = TRUE, verbose = FALSE)
data <- dsm_load(store, windows = integer(0), verbose = FALSE)
predictors <- as.character(data$store$predictors)
qc <- utils::read.csv2(file.path(data$patch_dir, "qc_table.csv"))
qc <- qc[match(predictors, qc$predictor), ]
# a fitted model's scaling is its predictor_scaling.csv; here, the profiles' own
x <- as.matrix(data$points[, predictors])
scaling <- data.frame(predictor = predictors, center = colMeans(x),
                      scale = apply(x, 2, sd))
plan <- spatial_folds(data$store$meta, k = 3, block_size = 0.05)
aref <- aoa_reference(data$points, predictors, qc, scaling, plan)
as.numeric(aref$threshold)     # beyond it, the cross-validated error does not apply
head(aref$cv)

The DI threshold that separates the area of applicability.

Description

Derived from the CROSS-VALIDATED training data: for each training row, the distance to the nearest training row in a DIFFERENT fold. That set is the dissimilarity the model actually coped with during cross-validation, so the outlier-removed maximum of it is the largest dissimilarity for which the reported error has been demonstrated.

Usage

aoa_threshold(ref, folds, k_iqr = 1.5)

Arguments

Details

A row never held out – the training side of a single split – has fold NA: it is a neighbour to every fold, and has no cross-validated DI of its own, since no model predicted it.

Value

The threshold, with the cross-validated DI attached as "cv_di" (NA for the rows never held out).

Examples

set.seed(1)
x_train <- matrix(rnorm(200), 100, 2)    # 100 training rows of 2 scaled predictors
ref <- di_reference(x_train)
th <- aoa_threshold(ref, folds = rep(1:5, 20))
as.numeric(th)                 # the threshold
summary(attr(th, "cv_di"))     # the cross-validated DI it was taken from

A metric, or two models' difference in it, with an interval from whole blocks.

Description

The points are drawn by blocks: each draw takes as many blocks as there are, with replacement, and every point of a block drawn comes with it. With against, the statistic is the difference pred less against, at the same points. Two weightings answer two questions: "profile", every point one vote; "block", every block one vote, its points weighing 1/n each. They differ in the estimate as well as in the interval, and where they disagree the advantage rests on the densely sampled blocks.

Usage

block_bootstrap(
  obs,
  pred,
  blocks,
  against = NULL,
  metric = c("mae", "rmse", "ccc", "bias"),
  weights = c("profile", "block"),
  n_boot = 2000L,
  conf = 0.95,
  seed = 42L
)

Arguments

Value

A block_bootstrap: a tibble, one row per metric and weighting, with estimate (pred's), against (the other model's), difference, ci_low and ci_high (of the difference when against is given, of the estimate otherwise), n_points, n_blocks, n_effective (for "profile": the number of independent points that would give the same interval) and share_better (the share of points, or of blocks, where pred errs less).

Examples

ex <- example_landscape()
obs <- ex$profiles$soc_stock
set.seed(1)
pred_a <- obs * exp(rnorm(length(obs), 0, 0.20))   # two models' predictions
pred_b <- obs * exp(rnorm(length(obs), 0, 0.25))
blk <- equal_area_blocks(ex$profiles$x, ex$profiles$y, size_km = 2)
block_bootstrap(obs, pred_a, blk, against = pred_b, metric = "mae", n_boot = 500)

The interval as the blocks grow.

Description

block_bootstrap() for blocks of growing side. Blocks smaller than the distance over which the values resemble each other are themselves dependent, and their interval is still too narrow; where the width stops growing, they no longer are.

Usage

block_bootstrap_by_size(
  x,
  y,
  obs,
  pred,
  against = NULL,
  sizes_km = c(25, 50, 100, 200, 400),
  metric = "mae",
  weights = "profile",
  n_boot = 2000L,
  conf = 0.95,
  seed = 42L,
  coords = c("lonlat", "metres"),
  centre = NULL
)

Arguments

Value

A block_bootstrap_by_size: a tibble, one row per size, with size_km, n_blocks, the estimate (or the difference), ci_low, ci_high, width and n_effective.

Examples

ex <- example_landscape()
obs <- ex$profiles$soc_stock
set.seed(1)
pred_a <- obs * exp(rnorm(length(obs), 0, 0.20))   # two models' predictions
pred_b <- obs * exp(rnorm(length(obs), 0, 0.25))
s <- block_bootstrap_by_size(ex$profiles$x, ex$profiles$y, obs, pred_a, against = pred_b,
                             sizes_km = c(1, 2, 5, 10), n_boot = 500)
s
plot(s)

Compute all regression metrics for one obs/pred pair.

Description

Compute all regression metrics for one obs/pred pair.

Usage

calc_metrics(obs, pred)

Arguments

Value

A one-row tibble with columns: n, ccc, r2, mae, nse, rmse, rpd, mqi, bias (signed, native units) and bias_pct (relative to mean(obs)).

Examples

obs  <- c(12, 25, 31, 48, 60)
pred <- c(15, 22, 35, 40, 66)
calc_metrics(obs, pred)

Every caret method that could be borrowed here.

Description

Every caret method that could be borrowed here.

Usage

caret_available(pattern = NULL)

Arguments

Value

A tibble of method, label and the parameters it tunes.

Examples

caret_available("^rf$|^ranger$")

Turn a caret method into a model_spec.

Description

Turn a caret method into a model_spec.

Usage

caret_spec(method, name = method, search = c("grid", "random"), ...)

Arguments

Value

A model_spec with input == "table".

Examples

caret_spec("rf")

Lin's Concordance Correlation Coefficient.

Description

Lin's Concordance Correlation Coefficient.

Usage

ccc(obs, pred)

Arguments

Value

A single numeric, or NA when fewer than two finite pairs remain.

Examples

obs <- c(10, 20, 30, 40)
ccc(obs, obs)        # 1: perfect agreement
ccc(obs, obs + 5)    # below 1: an offset is disagreement too, though r is 1

Check a plan's folds are well formed, and describe them.

Description

Run this before any training: a broken plan caught here costs a second, caught later it costs the whole run. Checks that, within each fold, train and validation are disjoint; that test (when present) appears in neither; and that across the folds every pooled row is validated exactly once – the property that makes the k metrics a partition of the pool rather than an arbitrary set of overlapping subsets.

Usage

check_fold_plan(plan, meta = NULL, group = "auto")

Arguments

Value

tibble, one row per fold.

Examples

ex <- example_landscape()
meta <- data.frame(sample_id = seq_len(nrow(ex$profiles)), x = ex$profiles$x,
                   y = ex$profiles$y)
plan <- spatial_folds(meta, k = 3, test_frac = 0.2, block_size = 0.05)
check_fold_plan(plan, meta)

Several importances, side by side.

Description

Aligns importances computed by different methods, variants, rows or runs – permutation over all rows against within blocks, the test set against the folds, the spatial design against the random one – by what each perturbed (the variable, for a permutation), and says how far their rankings agree. Where they disagree is the reading: a variable that keeps its importance within blocks discriminates locally; one that loses it carries a regional gradient.

Usage

compare_importance(..., n = 20L)

Arguments

Value

An importance_comparison: table (one row per target: each importance, its rank, and its share of that importance's largest), agreement (Spearman's correlation between the rankings, over the variables they share), labels, and for two SHAP importances of the same points points_agreement (per variable, the correlation of their values point by point and their mean difference against the second's size).

Examples

run <- example_run()    # a small fitted run, made once a session
perm <- dsm_importance(run$final, run$data, permutation_importance(draws = 2),
                       verbose = FALSE)
shap <- dsm_importance(run$final, run$data, shap_importance(samples = 20, background = 20),
                       verbose = FALSE)
compare_importance(permutation = perm, shap = shap)

Calibrate a conformal interval from held-out residuals.

Description

Calibrate a conformal interval from held-out residuals.

Usage

conformal_calibrate(obs, pred, alpha = 0.1, difficulty = NULL, group = NULL)

Arguments

Value

An object of class "conformal_cal".

Examples

set.seed(1)
obs  <- rlnorm(300, 3, 0.4)
pred <- obs * exp(rnorm(300, 0, 0.25))
cal  <- conformal_calibrate(obs[1:150], pred[1:150], alpha = 0.1)
cal
# the same residuals, every one of 30 groups weighing alike
conformal_calibrate(obs[1:150], pred[1:150], alpha = 0.1, group = rep(1:30, 5))

Calibrate on the folds of a resampling run, then report coverage.

Description

The natural calibration set in this framework is the VALIDATION side of each fold: it is held out of training, it is separated spatially with a buffer, and it already exists. This walks the per-unit predictions a run wrote, calibrates on one fold and measures coverage on the others, so the reported PICP is always out-of-calibration.

Usage

conformal_cv(pred_obs, alpha = 0.1, difficulty = NULL, group = NULL)

Arguments

Value

A picp_report over the pooled out-of-calibration points.

Examples

set.seed(1)
pred_obs <- data.frame(fold = rep(1:5, each = 60), obs = rlnorm(300, 3, 0.4))
pred_obs$pred <- pred_obs$obs * exp(rnorm(300, 0, 0.25))
conformal_cv(pred_obs, alpha = 0.1)

Turn predictions into intervals.

Description

Turn predictions into intervals.

Usage

conformal_interval(cal, pred, difficulty = NULL, lower_limit = -Inf)

Arguments

Value

A tibble: pred, lower, upper and width, one row per prediction.

Examples

set.seed(1)
obs  <- rlnorm(300, 3, 0.4)
pred <- obs * exp(rnorm(300, 0, 0.25))
cal  <- conformal_calibrate(obs[1:150], pred[1:150], alpha = 0.1)
iv <- conformal_interval(cal, pred[151:300], lower_limit = 0)
head(iv)

Fit the scale of an interval: |residual| on covariates that say how wrong.

Description

The scale of conformal_scaled_calibrate(), fitted on its own: a + b1 x1 + b2 x2 + ..., by least squares on |obs - pred|, with a floor. Fitted on points the quantile is NOT taken on – a calibration set's quantile, scaled by a fit to the cross-validated residuals, keeps the whole calibration set for q; fitted and calibrated on the same points, the fit would absorb the residuals the quantile measures.

Usage

conformal_scale_fit(obs, pred, covariates, floor_frac = 0.05)

Arguments

Value

A conformal_scale: the coefficients, the floor, and how well the fit explained |residual| (r2_fit).

Examples

set.seed(1)
level <- runif(300, 1, 4)
di    <- runif(300)                      # a dissimilarity index
obs   <- exp(level + rnorm(300, 0, 0.1 + 0.3 * di))
pred  <- exp(level)
conformal_scale_fit(obs, pred, data.frame(level = pred, di = di))

Calibrate an interval whose width is a fitted scale.

Description

Calibrate an interval whose width is a fitted scale.

Usage

conformal_scaled_calibrate(
  obs,
  pred,
  covariates,
  alpha = 0.1,
  fit_frac = 0.5,
  floor_frac = 0.05,
  seed = 42L,
  scale = NULL,
  group = NULL
)

Arguments

Value

A conformal_scaled (also a conformal_cal).

Examples

set.seed(1)
level <- runif(300, 1, 4)
di    <- runif(300)                      # a dissimilarity index
obs   <- exp(level + rnorm(300, 0, 0.1 + 0.3 * di))
pred  <- exp(level)
covariates <- data.frame(level = pred, di = di)
cal <- conformal_scaled_calibrate(obs[1:200], pred[1:200], covariates[1:200, ],
                                  alpha = 0.1)
cal
# the scale fitted on other points, and all 100 of these calibrating q
sc <- conformal_scale_fit(obs[1:200], pred[1:200], covariates[1:200, ])
conformal_scaled_calibrate(obs[201:300], pred[201:300], covariates[201:300, ],
                           scale = sc)

Turn predictions into intervals with a fitted scale.

Description

Turn predictions into intervals with a fitted scale.

Usage

conformal_scaled_interval(cal, pred, covariates, lower_limit = -Inf)

Arguments

Value

A tibble: pred, lower, upper and width, one row per prediction.

Examples

set.seed(1)
level <- runif(300, 1, 4)
di    <- runif(300)                      # a dissimilarity index
obs   <- exp(level + rnorm(300, 0, 0.1 + 0.3 * di))
pred  <- exp(level)
covariates <- data.frame(level = pred, di = di)
cal <- conformal_scaled_calibrate(obs[1:200], pred[1:200], covariates[1:200, ],
                                  alpha = 0.1)
iv <- conformal_scaled_interval(cal, pred[201:300], covariates[201:300, ],
                                lower_limit = 0)
picp_report(obs[201:300], iv$lower, iv$upper)

Context importance: how far from the point, and at which scale, the model reads.

Description

The patch is a neighbourhood. This permutes part of it – the point's own pixel, a ring of pixels at one distance from it, every ring but the centre, or one window's whole input – and measures the loss of skill, as permutation_importance() does for a variable. The windows are cut around the same point, so a ring is the same ground in every window, and it takes its values from one donor in all of them.

Usage

context_importance(
  by = c("ring", "window"),
  bands = NULL,
  per_variable = FALSE,
  draws = 5L,
  metric = c("ccc", "rmse", "rmse_transform"),
  seed = 42L
)

Arguments

Details

A large importance of the context says the network reads the neighbourhood. It does not say the neighbourhood adds what the centre lacks: at 250 m most covariates are smooth, a pixel three cells away is close to a copy of the centre, and a network can lean on the rim and learn nothing the centre would not have told it. Ring d holds 8d pixels, so the table gives the cost per pixel beside the total: a wider ring costs more for its area alone.

Value

An importance_spec, for dsm_importance().

Examples

context_importance(by = "ring")

Calibrate a CV+ interval from the folds of a cross-validation.

Description

Calibrate a CV+ interval from the folds of a cross-validation.

Usage

cv_plus_calibrate(obs, pred, fold, alpha = 0.1, difficulty = NULL)

Arguments

Value

A cv_plus_cal: each fold's sorted scores, the rank, the folds.

Examples

set.seed(1)
x <- runif(200)
fold <- rep(1:5, 40)
y <- 10 + 5 * x + rnorm(200)
# each fold's model: a line fitted without that fold
fits <- lapply(1:5, function(k) lm(y ~ x, data = data.frame(x, y)[fold != k, ]))
oof  <- vapply(seq_along(y), function(i)
  unname(predict(fits[[fold[i]]], data.frame(x = x[i]))), numeric(1))
cal <- cv_plus_calibrate(y, oof, fold, alpha = 0.1)
cal
# the five models at three new points, one column per fold
new <- data.frame(x = c(0.1, 0.5, 0.9))
cv_plus_interval(cal, sapply(fits, predict, newdata = new))

Turn the fold models' predictions into CV+ intervals.

Description

Turn the fold models' predictions into CV+ intervals.

Usage

cv_plus_interval(cal, fold_pred, difficulty = NULL, lower_limit = -Inf)

Arguments

Value

A tibble: lower, upper and width, one row per new point.

Examples

set.seed(1)
fold <- rep(1:4, 25)
oof  <- rnorm(100, 20, 3)              # each point's out-of-fold prediction
obs  <- oof + rnorm(100, 0, 2)
cal  <- cv_plus_calibrate(obs, oof, fold, alpha = 0.1)
# the four fold models at two new points: they disagree a little
cv_plus_interval(cal, rbind(c(19.5, 20.2, 20.0, 20.4), c(30.1, 29.0, 29.8, 30.6)),
                 lower_limit = 0)

Cross-validated residuals from a tuning run, for calibration.

Description

Every point appears once as validation, pooled over folds and averaged over seeds, so the residual distribution covers the whole study area instead of one fold's worth of it.

Usage

cv_residuals(run_dir, config_id, role = "validation")

Arguments

Value

A tibble with sample_id, obs, pred (the seed ensemble's median), and n_seeds; or NULL when the run wrote no usable predictions.

Examples

run <- example_run()    # a small fitted run, made once a session
res <- cv_residuals(run$fit$run_dir, run$final$selected_config_ids)
head(res)
conformal_calibrate(res$obs, res$pred, alpha = 0.1)

Cross-validated residuals of a configuration, from any tuning run.

Description

Cross-validated residuals of a configuration, from any tuning run.

Usage

cv_residuals_for_config(run_dir, cfg_row, required = TRUE)

Arguments

Value

What cv_residuals() returns, with the run's own config_id attached as attribute "config_id".

Examples

run <- example_run()    # a small fitted run, made once a session
grid <- readRDS(file.path(run$fit$run_dir, "tune_grid.rds"))
cfg  <- grid[grid$config_id == run$final$selected_config_ids, ]
# found by its hyperparameters, so the same call reads any other tuning run
head(cv_residuals_for_config(run$fit$run_dir, cfg))

Summarise a training set for dissimilarity-index computation.

Description

Computed once and reused for every prediction chunk: the training data does not change, and the pairwise mean is the expensive part.

Usage

di_reference(x_train, weights = NULL, max_pairs = 2e+06, seed = 42L)

Arguments

Value

An object of class "di_reference".

Examples

set.seed(1)
x_train <- matrix(rnorm(200), 100, 2)    # 100 training rows of 2 scaled predictors
ref <- di_reference(x_train)
ref

The dissimilarity index of new data.

Description

The dissimilarity index of new data.

Usage

dissimilarity_index(ref, x, chunk = 2000L)

Arguments

Value

Numeric vector, one DI per row. 0 means identical to a training point; 1 means as far from the training data as two training points are from each other on average.

Examples

set.seed(1)
x_train <- matrix(rnorm(200), 100, 2)    # 100 training rows of 2 scaled predictors
ref <- di_reference(x_train)
dissimilarity_index(ref, rbind(c(0, 0), c(4, 4)))   # near the data, and far from it

Refit the configuration a tuning run selects, under N seeds.

Description

Refit the configuration a tuning run selects, under N seeds.

Usage

dsm_final(
  tuning,
  data = NULL,
  config = "auto",
  rule = c("one_se", "rank1"),
  metric = "val_ccc",
  seeds = 10L,
  validation_frac = 0.15,
  predpoints = NULL,
  training = list(),
  transform = NULL,
  n_cores = NULL,
  threads_per_unit = 5L,
  max_ram_gb = NULL,
  conformal_alpha = c(0.1, 0.05),
  intervals = c("cv", "split", "cv_plus"),
  strata = NULL,
  output_dir = NULL,
  run_id = NULL,
  resume = TRUE,
  verbose = TRUE
)

Arguments

Value

A dsm_final, printed with the report.

Examples

run <- example_run()    # a small fitted run, made once a session
final <- dsm_final(run$fit, seeds = 2, n_cores = 1, threads_per_unit = 1,
                   training = list(n_epochs = 10, patience = 5), output_dir = tempdir(),
                   run_id = "final_example", verbose = FALSE)
final

Variable importance of a fitted model, by the method you choose.

Description

The method is an argument, as the validation design is in dsm_train(): permutation_importance() measures how much the model relies on each variable, over all rows or within blocks or classes; context_importance() how far from the point, and through which window, it reads; shap_importance() how each variable pushes each prediction, up or down; sage_importance() how much of the model's skill each variable carries, shared fairly among variables that carry the same information; refit_importance() how much skill a network trained without the variable loses; ale_effect() how the prediction changes along each variable's range. Each variable is a channel, or the channels of one categorical (see importance_groups()), or a group of yours.

Usage

dsm_importance(
  final,
  data,
  method = permutation_importance(),
  rows = c("test", "folds"),
  groups = "auto",
  config = NULL,
  threads = NULL,
  batch_size = 512L,
  verbose = TRUE,
  at = NULL,
  rasters = NULL,
  seeds = NULL,
  chunk_points = 5000L
)

Arguments

Details

Before anything is perturbed, each model's own score is computed again and must reproduce the one its run wrote; if it does not, nothing is measured.

Value

A dsm_importance: table (one row per target – a variable, a band, a variable at a band or a window – ranked), by_model (one row per target and model), baseline (each model's unperturbed score), groups, targets (what each target is), gate (by window, for a model with a gate: the gate point by point), units (the models scored) and the settings. For SHAP: table (mean |SHAP|, its share and direction per variable), points (every point's values), patch and pixels (where in the patch they sit) and completeness (how well they add up). For ALE: table (the spread, trend and extremes of each effect), curves (each continuous variable's curve, with its spread between models) and classes (each categorical's class effects). For SAGE: table (the loss each variable takes away, and its share), by_model and losses (each model's loss, and the mean prediction's, which SAGE splits). For a refit: table, by_model, raw (each unit's scores), check (the seeds trained again with nothing left out, against the run), noise (what a refit moves a score by with nothing to lose) and refit_dir.

Examples

run <- example_run()    # a small fitted run, made once a session
imp <- dsm_importance(run$final, run$data, permutation_importance(draws = 2),
                      verbose = FALSE)
imp
plot(imp)

Load everything a run needs, and check that it fits together.

Description

Replaces the six-step preamble every pipeline script used to repeat: open the store, read the points, read the predictor types, align the points to the store, read the raster resolution, and verify the store can serve this configuration.

Usage

dsm_load(
  patch_dir,
  points = NULL,
  type_table = NULL,
  windows = NULL,
  cell_size = NULL,
  raster_table = NULL,
  target_col = NULL,
  verbose = TRUE
)

Arguments

Value

A dsm_data object. ⁠$transform⁠ is the target transform the store was built under – name, forward and inverse – or NULL when the store did not record one.

Examples

ex <- example_landscape()
store <- dsm_prepare(ex$profiles, target = "soc_stock", raster_dir = ex$raster_dir,
                     windows = c(3, 7), out_dir = file.path(tempdir(), "landscape"),
                     percentage = "^clay_pct$", transform = "log1p", n_cores = 1,
                     overwrite = TRUE, verbose = FALSE)
dsm_load(store, windows = integer(0))   # the table alone: loading windows needs torch

Predict a map from a final model.

Description

Predict a map from a final model.

Usage

dsm_predict(
  final,
  data,
  rasters = NULL,
  qc_table = NULL,
  extent = NULL,
  config = NULL,
  calibration = NULL,
  aoa_weights = NULL,
  alpha = 0.1,
  intervals = NULL,
  weighting = c("point", "group"),
  clamp = NULL,
  bands = NULL,
  engine = c("auto", "patch"),
  output_dir = NULL,
  run_id = NULL,
  resume = TRUE,
  n_cores = NULL,
  threads_per_worker = 5L,
  max_ram_gb = NULL,
  unit_rows = NULL,
  step_rows = NULL,
  chunk_cols = 2048L,
  probe = TRUE,
  verbose = TRUE
)

Arguments

Value

A dsm_prediction, printed.

Examples

run <- example_run()    # a small fitted run, made once a session
map <- dsm_predict(run$final, run$data, output_dir = tempdir(), run_id = "map_example",
                   n_cores = 1, threads_per_worker = 1, verbose = FALSE)
map

Build a patch store from a point table and a folder of aligned rasters.

Description

Build a patch store from a point table and a folder of aligned rasters.

Usage

dsm_prepare(
  points,
  target,
  raster_dir,
  windows,
  out_dir = NULL,
  store_dir = NULL,
  metadata_dir = NULL,
  points_file = NULL,
  profile_id = "profile_id",
  coords = c("x", "y"),
  crs = NULL,
  percentage = NULL,
  dummy = "auto",
  drop = NULL,
  na_below = NULL,
  percentage_limits = c(0, 100),
  transform = c("none", "log1p"),
  target_min = NULL,
  subsample = NULL,
  target_label = NULL,
  target_unit = NULL,
  raster_pattern = "\\.tif$",
  chunk_nrows = 1000L,
  n_cores = NULL,
  max_ram_gb = NULL,
  read_gap = 256L,
  read_max_cols = 4096L,
  overwrite = FALSE,
  verbose = TRUE
)

Arguments

Value

A dsm_store, which dsm_load() accepts directly.

Examples

ex <- example_landscape()
store <- dsm_prepare(ex$profiles, target = "soc_stock", raster_dir = ex$raster_dir,
                     windows = c(3, 7), out_dir = file.path(tempdir(), "landscape"),
                     percentage = "^clay_pct$", transform = "log1p", n_cores = 1,
                     overwrite = TRUE, verbose = FALSE)
store

The declaration of a final model that already exists.

Description

dsm_final() writes final_report.md as it fits. A final model fitted by stage 04, before dsm_final() existed, has the seeds and the files but no declaration – and retraining ten seeds to obtain one would also change them (T1: the thread count changes the numbers). This writes it from what is on disk, re-assembling the seeds the way P3 proved dsm_final() does, and changes nothing that is there: it ADDS final_report.md and selected_hyperparameters.csv to the run directory.

Usage

dsm_report_final(
  run_dir,
  tuning_dir,
  conformal_alpha = c(0.1, 0.05),
  verbose = TRUE
)

Arguments

Value

A dsm_final describing the run, printed with the declaration.

Examples

run <- example_run()    # a small fitted run, made once a session
# the declaration of a final run, written from what is on disk
dsm_report_final(run$final$run_dir, run$fit$run_dir)

Fit and tune a model over a resampling plan.

Description

Fit and tune a model over a resampling plan.

Usage

dsm_train(
  data,
  model = "cnn",
  resampling = spatial_cv(),
  tune_grid = NULL,
  tune_length = 20L,
  n_seeds = 3L,
  transform = NULL,
  clamp = c(0, Inf),
  features = c("centre", "window_mean"),
  output_dir,
  run_id = format(Sys.time(), "%Y%m%d_%H%M%S"),
  base_seed = 42L,
  device = NULL,
  n_cores = NULL,
  in_session = !is.null(device),
  threads_per_unit = 5L,
  max_ram_gb = NULL,
  test_ids = NULL,
  calibration_ids = NULL,
  resume = TRUE,
  evaluate_test = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Value

The runner's result, plus the plan and the data it used.

Examples

ex <- example_landscape()
store <- dsm_prepare(ex$profiles, target = "soc_stock", raster_dir = ex$raster_dir,
                     windows = c(3, 7), out_dir = file.path(tempdir(), "landscape"),
                     percentage = "^clay_pct$", transform = "log1p", n_cores = 1,
                     overwrite = TRUE, verbose = FALSE)
data <- dsm_load(store, verbose = FALSE)
fit <- dsm_train(data, resampling = spatial_cv(k = 2, block_size = 0.05, test_frac = 0.2),
                 tune_length = 2, n_seeds = 1, output_dir = tempdir(),
                 in_session = TRUE, n_cores = 2, n_epochs = 5, verbose = FALSE)
fit

Square blocks of equal area, in kilometres.

Description

Each point falls in a square cell of size_km on a Lambert azimuthal equal-area projection centred on the points (on the sphere of WGS84's authalic radius), so every block covers the same ground wherever it lies – which a grid in degrees does not: a degree of longitude is 111 km at the equator and 71 km at 50 degrees.

Usage

equal_area_blocks(x, y, size_km, coords = c("lonlat", "metres"), centre = NULL)

Arguments

Value

A character vector, one block per point (the cell's column and row, "i_j"), with the size and the projection as attributes.

Examples

ex <- example_landscape()
blk <- equal_area_blocks(ex$profiles$x, ex$profiles$y, size_km = 2)
length(unique(blk))    # blocks holding a profile
table(table(blk))      # how many blocks hold 1, 2, ... profiles

A small synthetic landscape, for the examples and the tests.

Description

80 x 80 cells of 0.0025 degrees (about 250 m) in coordinates of south-eastern Brazil – the place is borrowed, the landscape is made up – with eight predictors and 160 soil profiles, 90 of them in six survey clusters. The soil organic carbon stock (t/ha, 0-30 cm) was drawn from vegetation, clay, temperature, geology and the topographic position of the cell in its 7 x 7 neighbourhood: a network that reads the neighbourhood has something to find that the centre pixel does not show.

Usage

example_landscape()

Value

A list: raster_dir, the folder of the eight GeoTIFFs (elevation, temperature, precipitation, ndvi, clay_pct, and geology as three 0/1 dummies); profiles, a data frame of profile_id, x, y, survey and soc_stock; and truth, the formula the stock was drawn from.

Examples

ex <- example_landscape()
head(ex$profiles)
list.files(ex$raster_dir)

A small fitted run on the example landscape, made once a session.

Description

The store, a tuning run and a final model on example_landscape(), in dir: what the examples of the functions that need a fitted model start from. Made at the first call – a minute or two on one core – and kept for the session; later calls return it at once. It is an example, not a model: 160 profiles, one small configuration, 30 epochs.

Usage

example_run(dir = file.path(tempdir(), "soilcnn_example"), verbose = FALSE)

Arguments

Value

A list: data (from dsm_load()), fit (from dsm_train()), final (from dsm_final()) and dir.

Examples

run <- example_run()
run$final

How much does each fold still leak?

Description

Reports, per fold, the same quantities the pipeline check reports for the fixed split: the share of validation points sharing a raster cell with a training point, and sharing patch pixels at each window. Deliberately the SAME criterion, so a number here can be compared with a number there rather than being a second, incompatible notion of leakage.

Usage

fold_leakage_report(plan, meta, cell_size, windows = c(3L, 9L, 15L))

Arguments

Value

A tibble with a fold column: the overlap shares per role and window, fold by fold.

Examples

ex <- example_landscape()
meta <- data.frame(sample_id = seq_len(nrow(ex$profiles)), x = ex$profiles$x,
                   y = ex$profiles$y)
blocks   <- spatial_folds(meta, k = 3, block_size = 0.05)
buffered <- spatial_folds(meta, k = 3, block_size = 0.05, buffer = 7 * 0.0025)
fold_leakage_report(blocks, meta, cell_size = 0.0025, windows = c(3, 7))
fold_leakage_report(buffered, meta, cell_size = 0.0025, windows = c(3, 7))

Record which config was chosen, and when.

Description

Writes the selection into the run directory so a later test-set report can prove the choice preceded it. Call it at the moment the choice is made.

Usage

freeze_selection(
  run_dir,
  config_id,
  rule = "one_se",
  metric = "val_ccc",
  note = NA_character_
)

Arguments

Value

The selection record, invisibly.

Examples

run_dir <- file.path(tempdir(), "tuning_example")
freeze_selection(run_dir, "cfg_002", note = "one_se on the cross-validation")
# the same choice again is accepted; another is refused
try(freeze_selection(run_dir, "cfg_005"))

Fetch a registered model.

Description

Fetch a registered model.

Usage

get_model(name)

Arguments

Value

The model_spec registered under name.

Examples

get_model("rf")

A single train/validation/test split of a point table.

Description

The fold constructor behind holdout_cv(). Whole groups are drawn, never single rows (see group), so the rows of one profile never straddle the split.

Usage

holdout(
  meta,
  validation_frac = 0.15,
  test_frac = 0.15,
  test_ids = NULL,
  seed = 42L,
  group = "auto",
  calibration_frac = 0,
  calibration_ids = NULL
)

Arguments

Value

A fold_plan with one fold.

Examples

ex <- example_landscape()
meta <- data.frame(sample_id = seq_len(nrow(ex$profiles)), x = ex$profiles$x,
                   y = ex$profiles$y)
holdout(meta, validation_frac = 0.2, test_frac = 0.2)
# with a calibration set for the intervals beside the test set
holdout(meta, validation_frac = 0.2, test_frac = 0.2, calibration_frac = 0.15)

A single train/validation/test split.

Description

A single train/validation/test split.

Usage

holdout_cv(
  validation_frac = 0.15,
  test_frac = 0.15,
  group = "auto",
  seed = 42L,
  calibration_frac = 0
)

Arguments

Value

A resample_spec, as spatial_cv() returns.

Examples

holdout_cv(validation_frac = 0.2)

Which channels the importance treats as one variable.

Description

A categorical raster arrives as several 0/1 channels, one per class. Permuting one of them alone makes points that belong to two classes or to none, so the channels of one categorical are one variable. The package is not told which raster a dummy came from, so it infers the sets: dummies that share a name prefix (cut at "_") and are never 1 together at the points. Every other channel is a variable of its own. The sets found are listed; check them, and give your own grouping when they are not what the rasters were.

Usage

importance_groups(data, groups = "auto")

Arguments

Details

A grouping of your own can be anything that makes sense to perturb as one: the classes of a categorical named outright, or themes (climate, relief, vegetation...), which answers how much the model relies on each theme as a block – the way to read many correlated variables.

Value

A tibble in channel order: channel, variable, and rule ("alone", "one-hot set" or "yours").

Examples

ex <- example_landscape()
store <- dsm_prepare(ex$profiles, target = "soc_stock", raster_dir = ex$raster_dir,
                     windows = c(3, 7), out_dir = file.path(tempdir(), "landscape"),
                     percentage = "^clay_pct$", transform = "log1p", n_cores = 1,
                     overwrite = TRUE, verbose = FALSE)
data <- dsm_load(store, windows = integer(0), verbose = FALSE)
importance_groups(data)          # the geology dummies, found as one variable
themes <- data.frame(channel = c("temperature", "precipitation"), variable = "climate")
importance_groups(data, themes)  # and two channels made one by you

Maps of SHAP values, from an importance computed at points of the map.

Description

Lays the values of dsm_importance()⁠(at = )⁠ back on a grid: each cell the mean of the points in it. With the points' own grid (from importance_points()) every point is its own cell; with a coarser resolution, each cell averages the points it holds.

Usage

importance_map(x, resolution = NULL, output_dir = NULL)

Arguments

Value

An importance_map, a list: shap (one layer per variable: the mean SHAP value of the points in each cell, in the network's units), dominant (in each cell, the variable with the largest mean |SHAP|, as the number in legend), legend, prediction (the mean prediction, native units) and files. The three rasters are packed with terra::wrap(), so the map survives base::saveRDS(); terra::unwrap() gives each as a SpatRaster. plot() draws it.

Examples

run <- example_run()    # a small fitted run, made once a session
pts <- importance_points(run$final, run$data, every = 4)
imp <- dsm_importance(run$final, run$data, shap_importance(samples = 10, background = 20),
                      at = pts, seeds = 42, verbose = FALSE)
m <- importance_map(imp)
plot(m)

Points of the map to explain: a regular grid of the model's raster cells.

Description

The cells of the rasters the model was fitted on, every every rows and columns, inside extent, where the first channel has a value. With every = 1 every cell – the map at full resolution, for a tile; with 40 at 250 m, one point every 10 km, for a region. The points carry their grid, so importance_map() can lay the values back on it.

Usage

importance_points(final, data, extent = NULL, every = 1L, rasters = NULL)

Arguments

Value

A tibble of x, y (cell centres), with the grid as an attribute.

Examples

run <- example_run()    # a small fitted run, made once a session
pts <- importance_points(run$final, run$data, every = 4)   # every 4th cell
nrow(pts)

One weight per channel, for the area of applicability, from an importance.

Description

The dissimilarity index measures how far a pixel is from the training data, one axis per channel. Unweighted, a channel the model ignores counts as much as the one it leans on, and a pixel unlike the training data in an ignored channel falls outside the area of applicability for nothing. Meyer & Pebesma (2021) weight each axis by the predictor's importance; this gives each channel its variable's importance – every channel of a one-hot set, or of a group of yours, its variable's – and a negative importance, noise around zero, weight zero.

Usage

importance_weights(x)

Arguments

Value

A named numeric vector, one weight per channel in the store's order, for dsm_predict(aoa_weights = ) or aoa_reference(weights = ).

Examples

run <- example_run()    # a small fitted run, made once a session
imp <- dsm_importance(run$final, run$data, shap_importance(samples = 20, background = 20),
                      verbose = FALSE)
importance_weights(imp)   # for dsm_predict(aoa_weights = ) or aoa_reference(weights = )

Is each row inside the area of applicability?

Description

Is each row inside the area of applicability?

Usage

inside_aoa(di, threshold)

Arguments

Value

Logical vector. TRUE means the cross-validated error applies here.

Examples

set.seed(1)
x_train <- matrix(rnorm(200), 100, 2)    # 100 training rows of 2 scaled predictors
ref <- di_reference(x_train)
threshold <- aoa_threshold(ref, folds = rep(1:5, 20))
inside_aoa(dissimilarity_index(ref, rbind(c(0, 0), c(4, 4))), threshold)

Folds matched to where the map will be predicted (kNNDM).

Description

The better-founded alternative to spatial_cv(). Blocks need a block size and a buffer width, and nothing in the data says whether the chosen ones were right. kNNDM instead shapes the folds so the distance from a validation point to its nearest training point is distributed like the distance from a PREDICTION pixel to its nearest training point.

Usage

knndm_cv(
  k = 5L,
  predpoints = NULL,
  hold_out_test = FALSE,
  crs = 4326,
  project_to = "+proj=moll +lon_0=0 +datum=WGS84 +units=m",
  seed = 42L,
  hold_out_calibration = FALSE,
  ...
)

Arguments

Details

It therefore needs predpoints: a sample of where the map will be drawn. There is no default, because a default would turn the method into an expensive random split. prediction_sample(raster) produces one.

When the samples are well spread over the prediction area it converges by itself to ordinary random k-fold – it does not impose separation that prediction will not face. That is the property blocks cannot have.

Needs the CAST and sf packages. See R/knndm.R for the projection question, which is not optional on lon/lat data.

Value

A resample_spec, as spatial_cv() returns.

Examples

ex <- example_landscape()
elevation <- terra::rast(file.path(ex$raster_dir, "elevation.tif"))
knndm_cv(k = 5, predpoints = prediction_sample(elevation, size = 500))

Fold plan from k-fold Nearest Neighbour Distance Matching.

Description

Fold plan from k-fold Nearest Neighbour Distance Matching.

Usage

knndm_folds(
  meta,
  k = 5L,
  predpoints = NULL,
  test_ids = NULL,
  hold_out_test = FALSE,
  crs = 4326,
  project_to = "+proj=moll +lon_0=0 +datum=WGS84 +units=m",
  seed = 42L,
  calibration_ids = NULL,
  hold_out_calibration = FALSE,
  ...
)

Arguments

Value

A fold_plan.

Examples

ex <- example_landscape()
meta <- data.frame(sample_id = seq_len(nrow(ex$profiles)), x = ex$profiles$x,
                   y = ex$profiles$y)
elevation <- terra::rast(file.path(ex$raster_dir, "elevation.tif"))
knndm_folds(meta, k = 3, predpoints = prediction_sample(elevation, size = 500))

The most recent run under base, by time, and only if it finished.

Description

WHY THIS IS NOT sort(dirs, decreasing = TRUE)[1].

Usage

latest_run_dir(
  base,
  prefix,
  require_file = NULL,
  require_pattern = NULL,
  label = "run",
  on_none = c("stop", "null")
)

Arguments

Details

That expression means "last alphabetically", which equals "most recent" only while every run id is a timestamp sharing one prefix. Runs given names broke it silently: soc_0_5cm_design_spatial sorts ahead of soc_0_5cm_20260916_232318 because 'd' > '2', so "latest" started resolving to a run that had died on its first unit, and stage 04 would have refit the final model against it without a word.

It also never asked whether the run FINISHED. Stage 04 creates its output directory before its own validations run, so a failure leaves a ⁠final_<timestamp>⁠ behind that "latest" would then deploy – which is what B2's review found from the other direction.

Value

The run id (basename), or stop() / NULL when nothing qualifies.

Examples

base <- file.path(tempdir(), "runs")
for (id in c("run_a", "run_b")) {
  dir.create(file.path(base, id, "comparison"), recursive = TRUE, showWarnings = FALSE)
  writeLines("done", file.path(base, id, "comparison", "comparison_ranked.csv"))
}
Sys.setFileTime(file.path(base, "run_a", "comparison", "comparison_ranked.csv"),
                Sys.time() - 3600)
dir.create(file.path(base, "run_c"), showWarnings = FALSE)   # started, never finished
latest_run_dir(base, "run_", require_file = "comparison/comparison_ranked.csv")

Every registered model, as a table.

Description

Every registered model, as a table.

Usage

list_models()

Value

A tibble: name, input ("table" or "patches"), tunable (whether it draws its own grid) and description.

Examples

list_models()

Build a manual grid from explicit lists of values per parameter.

Description

Build a manual grid from explicit lists of values per parameter.

Usage

make_manual_tune_grid(...)

Arguments

Value

A tibble, one configuration per row, with a config_id column.

Examples

make_manual_tune_grid(
  embedding_dim = c(256L, 384L),
  gate_type     = c("vector_featurewise", "no_gate_concat"),
  base_lr       = c(0.001, 0.0005)
)

Generate a random hyperparameter grid (like caret's tuneLength).

Description

Generate a random hyperparameter grid (like caret's tuneLength).

Usage

make_tune_grid(
  tune_length = 20L,
  seed = NULL,
  fixed = list(),
  windows = NULL,
  n_train = NULL
)

Arguments

Value

A tibble with one row per configuration. window_sizes and conv_channels are stored as list-columns.

Examples

make_tune_grid(tune_length = 3, seed = 1, windows = c(3, 7))

Describe a model the framework can fit.

Description

Describe a model the framework can fit.

Usage

model_spec(
  name,
  input,
  fit,
  predict,
  default_grid = NULL,
  count_params = NULL,
  description = "",
  fit_args = NULL
)

Arguments

Value

A model_spec.

Examples

mean_only <- model_spec(
  "mean_only", input = "table",
  fit = function(x, y, cfg, ...) list(mu = mean(y)),
  predict = function(object, x, ...) rep(object$mu, nrow(x)),
  description = "the training mean: a floor every model must beat")
mean_only

What each part of the patch is worth to one trained unit of a tuning run.

Description

Hides a part of every patch – the centre pixel, each ring around it, the whole context – and measures what the unit loses for it. The caller names the unit and nothing more: the checkpoint is found, the architecture rebuilt from the grid and the fold's cache from the plan.

Usage

occlusion_report(
  run_dir,
  data,
  config_id,
  fold = 1L,
  seed_i = 1L,
  role = "validation",
  transform = NULL,
  device,
  ...
)

Arguments

Value

A spatial_occlusion: table, one row per part hidden, with the pixels hidden, the CCC and MAE without them and the change in CCC, and baseline_ccc, the unit's own. Printed with what it does and does not mean.

Examples

run <- example_run()    # a small fitted run, made once a session
occlusion_report(run$fit$run_dir, run$data, config_id = run$final$selected_config_ids,
                 device = torch::torch_device("cpu"))

Pick the simplest config within one standard error of the best.

Description

Pick the simplest config within one standard error of the best.

Usage

one_se(by_config, metric = "val_ccc", complexity = "n_params", maximise = NULL)

Arguments

Value

One row of by_config, with within_one_se (how many configs were tied) and simpler_than_best (whether the rule actually moved the choice) attached – a selection rule that silently returns the same answer as the default should say so.

Examples

# A tuning table: three configurations, each on three folds under two seeds
set.seed(1)
cfgs <- data.frame(config_id = c("cfg_001", "cfg_002", "cfg_003"),
                   skill = c(0.615, 0.62, 0.55), n_params = c(2e5, 2e6, 5e4))
cmp <- merge(expand.grid(config_id = cfgs$config_id, fold = 1:3, seed = 1:2,
                         stringsAsFactors = FALSE), cfgs)
cmp$status  <- "success"
cmp$val_ccc <- cmp$skill + rnorm(nrow(cmp), 0, 0.02)
pick <- one_se(summarise_resamples(cmp, metrics = "val_ccc"), metric = "val_ccc")
pick$config_id

Paired comparison between two model families.

Description

WHY PAIRED, AND WHY IT MATTERS HERE.

Usage

paired_family_test(
  a,
  b,
  metric = "val_ccc",
  config_a = NULL,
  config_b = NULL,
  label_a = "a",
  label_b = "b",
  conf = 0.95
)

Arguments

Details

Every family in stage 03b runs on the SAME fold plan with the SAME seeds, by construction. So each family's units come in matched pairs: rf_context on fold 2 seed 3 and the CNN on fold 2 seed 3 saw the identical training rows.

Comparing the two MEANS and their separate spreads throws that away. Most of the spread across units is the fold – some folds are simply harder, and both families suffer on them together. A paired comparison subtracts that shared difficulty out before asking whether anything is left.

The first 03b run made exactly this mistake in the reader's favour: it reported +0.0091 CCC against a seed noise floor of 0.0383 and concluded "smaller than the noise". The conclusion happened to hold, but the reasoning did not – the noise floor is the spread of ONE family across seeds, which is not the standard error of a DIFFERENCE, and comparing a difference to it can hide a real effect as easily as it can invent one.

What is reported is a mean difference with its own standard error and a paired t interval. Not a p-value on its own: with 9 pairs, "not significant" is mostly a statement about 9, and the interval says how large an effect is still compatible with the data – which is the actual question when the answer is "the convolution buys nothing".

Value

An object of class "paired_comparison".

Examples

# A tuning table: three configurations, each on three folds under two seeds
set.seed(1)
cfgs <- data.frame(config_id = c("cfg_001", "cfg_002", "cfg_003"),
                   skill = c(0.615, 0.62, 0.55), n_params = c(2e5, 2e6, 5e4))
cmp <- merge(expand.grid(config_id = cfgs$config_id, fold = 1:3, seed = 1:2,
                         stringsAsFactors = FALSE), cfgs)
cmp$status  <- "success"
cmp$val_ccc <- cmp$skill + rnorm(nrow(cmp), 0, 0.02)
# a second family, trained on the same folds under the same seeds
forest <- cmp
forest$config_id <- sub("cfg", "rf", forest$config_id)
forest$val_ccc <- forest$val_ccc - 0.03 + rnorm(nrow(forest), 0, 0.01)
paired_family_test(cmp, forest, metric = "val_ccc", config_a = "cfg_002",
                   config_b = "rf_002", label_a = "cnn", label_b = "forest")

Permutation importance: how much a fitted model relies on each variable.

Description

Each point takes one variable's whole patch – every pixel, every window – from another point; the model predicts again, and the loss of skill is the variable's importance to that model. Nothing is retrained.

Usage

permutation_importance(
  draws = 5L,
  within = NULL,
  fill = c("permute", "mean"),
  metric = c("ccc", "rmse", "rmse_transform"),
  seed = 42L
)

Arguments

Details

A variable whose permutation costs nothing is one the model does not use. It may still carry information other variables also carry: permutation measures reliance, not necessity, and a model that can find the same thing elsewhere leans on whichever it found first.

The mean of a z-scored channel is 0 and of a dummy its class frequency, so fill = "mean" gives every point the average landscape, flat. A large drop under it with a small one under permutation says the model reacts to the input being unusual rather than to its content.

Value

An importance_spec, for dsm_importance().

Examples

permutation_importance()
permutation_importance(draws = 3, within = 0.05)   # donors from the same 0.05 block

Coverage, and whether it holds where it is needed.

Description

Coverage, and whether it holds where it is needed.

Usage

picp_report(obs, lower, upper, group = NULL, alpha = 0.1)

Arguments

Value

An object of class "picp_report".

Examples

set.seed(1)
obs  <- rlnorm(300, 3, 0.4)
pred <- obs * exp(rnorm(300, 0, 0.25))
cal  <- conformal_calibrate(obs[1:150], pred[1:150], alpha = 0.1)
iv <- conformal_interval(cal, pred[151:300], lower_limit = 0)
# one width for every level: right on average, and by level?
picp_report(obs[151:300], iv$lower, iv$upper, alpha = 0.1,
            group = ifelse(pred[151:300] > median(pred), "high", "low"))

Plot the interval as the blocks grow.

Description

Plot the interval as the blocks grow.

Usage

## S3 method for class 'block_bootstrap_by_size'
plot(x, ...)

Arguments

Value

x, invisibly.

Plot a variable importance.

Description

The figure that answers the importance's own question: bars for a permutation, the cost per pixel by ring (or per variable, or per window) for a context, the summary plot for SHAP – a point per profile at its SHAP value, coloured by the variable's value there – the curves for ALE, for SAGE the loss each variable takes away, and for a refit the skill a network trained without it loses.

Usage

## S3 method for class 'dsm_importance'
plot(x, n = 20L, ...)

Arguments

Value

x, invisibly.

Plot importances side by side.

Description

The first two importances' shares of their largest, one against the other, each variable a point on a 1:1 line; and, for two SHAP importances of the same points, how well they agree variable by variable.

Usage

## S3 method for class 'importance_comparison'
plot(x, ...)

Arguments

Value

x, invisibly.

Plot maps of SHAP values.

Description

One panel per variable on a diverging scale centred on zero – where the variable raises the prediction and where it lowers it – then the dominant variable and the prediction.

Usage

## S3 method for class 'importance_map'
plot(x, n = 12L, ...)

Arguments

Value

x, invisibly.

Plot a spatial correlogram.

Description

Plot a spatial correlogram.

Usage

## S3 method for class 'spatial_correlogram'
plot(x, ...)

Arguments

Value

x, invisibly.

A regular sample of the prediction area, as kNNDM's predpoints.

Description

kNNDM needs the places the map will cover, not the places it was trained on. A few thousand regularly spaced pixels describe that well enough – the method compares DISTRIBUTIONS of distances, and a regular sample estimates one cheaply. A random sample would work too and is noisier for the same size.

Usage

prediction_sample(raster, size = 5000L)

Arguments

Value

A data frame with x and y, in the raster's own CRS.

Examples

ex <- example_landscape()
elevation <- terra::rast(file.path(ex$raster_dir, "elevation.tif"))
head(prediction_sample(elevation, size = 200))

Print a block_bootstrap

Description

Print a block_bootstrap

Usage

## S3 method for class 'block_bootstrap'
print(x, ...)

Arguments

Value

x, invisibly.

Print a caret_fitted

Description

Print a caret_fitted

Usage

## S3 method for class 'caret_fitted'
print(x, ...)

Arguments

Value

x, invisibly.

Print a check_ledger

Description

Print a check_ledger

Usage

## S3 method for class 'check_ledger'
print(x, ...)

Arguments

Value

x, invisibly.

Print a conformal_cal

Description

Print a conformal_cal

Usage

## S3 method for class 'conformal_cal'
print(x, ...)

Arguments

Value

x, invisibly.

Print a conformal_scale

Description

Print a conformal_scale

Usage

## S3 method for class 'conformal_scale'
print(x, ...)

Arguments

Value

x, invisibly.

Print a conformal_scaled

Description

Print a conformal_scaled

Usage

## S3 method for class 'conformal_scaled'
print(x, ...)

Arguments

Value

x, invisibly.

Print a cv_plus_cal

Description

Print a cv_plus_cal

Usage

## S3 method for class 'cv_plus_cal'
print(x, ...)

Arguments

Value

x, invisibly.

Print a di_reference

Description

Print a di_reference

Usage

## S3 method for class 'di_reference'
print(x, ...)

Arguments

Value

x, invisibly.

Print a dsm_data

Description

Print a dsm_data

Usage

## S3 method for class 'dsm_data'
print(x, ...)

Arguments

Value

x, invisibly.

Print a dsm_final

Description

Print a dsm_final

Usage

## S3 method for class 'dsm_final'
print(x, ...)

Arguments

Value

x, invisibly.

Print a dsm_fit

Description

Print a dsm_fit

Usage

## S3 method for class 'dsm_fit'
print(x, ...)

Arguments

Value

x, invisibly.

Print a variable importance.

Description

Print a variable importance.

Usage

## S3 method for class 'dsm_importance'
print(x, n = 20L, ...)

Arguments

Value

x, invisibly.

Print a dsm_prediction

Description

Print a dsm_prediction

Usage

## S3 method for class 'dsm_prediction'
print(x, ...)

Arguments

Value

x, invisibly.

Print a dsm_store

Description

Print a dsm_store

Usage

## S3 method for class 'dsm_store'
print(x, ...)

Arguments

Value

x, invisibly.

Print a fold plan, with its per-fold sizes.

Description

Print a fold plan, with its per-fold sizes.

Usage

## S3 method for class 'fold_plan'
print(x, ...)

Arguments

Value

x, invisibly.

Print a comparison of importances.

Description

Print a comparison of importances.

Usage

## S3 method for class 'importance_comparison'
print(x, ...)

Arguments

Value

x, invisibly.

Print an importance method.

Description

Print an importance method.

Usage

## S3 method for class 'importance_spec'
print(x, ...)

Arguments

Value

x, invisibly.

Print a mlp_fitted

Description

Print a mlp_fitted

Usage

## S3 method for class 'mlp_fitted'
print(x, ...)

Arguments

Value

x, invisibly.

Print a model_spec

Description

Print a model_spec

Usage

## S3 method for class 'model_spec'
print(x, ...)

Arguments

Value

x, invisibly.

Print a paired_comparison

Description

Print a paired_comparison

Usage

## S3 method for class 'paired_comparison'
print(x, ...)

Arguments

Value

x, invisibly.

Print a picp_report

Description

Print a picp_report

Usage

## S3 method for class 'picp_report'
print(x, ...)

Arguments

Value

x, invisibly.

Print a resample_spec

Description

Print a resample_spec

Usage

## S3 method for class 'resample_spec'
print(x, ...)

Arguments

Value

x, invisibly.

Print a rf_fitted

Description

Print a rf_fitted

Usage

## S3 method for class 'rf_fitted'
print(x, ...)

Arguments

Value

x, invisibly.

Print a smearing_cal

Description

Print a smearing_cal

Usage

## S3 method for class 'smearing_cal'
print(x, ...)

Arguments

Value

x, invisibly.

Print a smearing_check

Description

Print a smearing_check

Usage

## S3 method for class 'smearing_check'
print(x, ...)

Arguments

Value

x, invisibly.

Print a spatial_correlogram

Description

Print a spatial_correlogram

Usage

## S3 method for class 'spatial_correlogram'
print(x, ...)

Arguments

Value

x, invisibly.

Say what the occlusion found, including when it found nothing.

Description

Say what the occlusion found, including when it found nothing.

Usage

## S3 method for class 'spatial_occlusion'
print(x, ...)

Arguments

Value

x, invisibly.

Say what the grid did on the test set, and what that does and does not mean.

Description

Say what the grid did on the test set, and what that does and does not mean.

Usage

## S3 method for class 'test_optimism'
print(x, ...)

Arguments

Value

x, invisibly.

Report what an AOA covers.

Description

Prints the threshold, the share of cells inside it and the quantiles of the DI, and says so plainly when less than half the map is inside.

Usage

print_aoa(di, threshold, label = "prediction area")

Arguments

Value

The logical mask of the cells inside (inside_aoa()), invisibly.

Examples

set.seed(1)
x_train <- matrix(rnorm(200), 100, 2)    # 100 training rows of 2 scaled predictors
ref <- di_reference(x_train)
threshold <- aoa_threshold(ref, folds = rep(1:5, 20))
print_aoa(dissimilarity_index(ref, matrix(rnorm(400, sd = 1.5), 200, 2)), threshold)

Print a noise-floor report in the terms it should be read in.

Description

Print a noise-floor report in the terms it should be read in.

Usage

print_noise_floor(nf, digits = 4L)

Arguments

Value

nf, invisibly; NULL when no floor could be estimated.

Examples

# A tuning table: three configurations, each on three folds under two seeds
set.seed(1)
cfgs <- data.frame(config_id = c("cfg_001", "cfg_002", "cfg_003"),
                   skill = c(0.615, 0.62, 0.55), n_params = c(2e5, 2e6, 5e4))
cmp <- merge(expand.grid(config_id = cfgs$config_id, fold = 1:3, seed = 1:2,
                         stringsAsFactors = FALSE), cfgs)
cmp$status  <- "success"
cmp$val_ccc <- cmp$skill + rnorm(nrow(cmp), 0, 0.02)
print_noise_floor(seed_noise_floor(cmp, metric = "val_ccc"))

Say what one_se() did, including when it did nothing.

Description

Say what one_se() did, including when it did nothing.

Usage

print_one_se(pick, metric = "val_ccc", digits = 4L)

Arguments

Value

pick, invisibly.

Examples

# A tuning table: three configurations, each on three folds under two seeds
set.seed(1)
cfgs <- data.frame(config_id = c("cfg_001", "cfg_002", "cfg_003"),
                   skill = c(0.615, 0.62, 0.55), n_params = c(2e5, 2e6, 5e4))
cmp <- merge(expand.grid(config_id = cfgs$config_id, fold = 1:3, seed = 1:2,
                         stringsAsFactors = FALSE), cfgs)
cmp$status  <- "success"
cmp$val_ccc <- cmp$skill + rnorm(nrow(cmp), 0, 0.02)
print_one_se(one_se(summarise_resamples(cmp, metrics = "val_ccc"), metric = "val_ccc"))

Random k-fold. Ignores geography by construction.

Description

Right when the rows really are independent, and the cleanest way to MEASURE what geography is worth: run it against spatial_cv() on the same points and the gap is the spatial optimism.

Usage

random_cv(
  k = 5L,
  test_frac = 0.15,
  group = "auto",
  seed = 42L,
  calibration_frac = 0
)

Arguments

Value

A resample_spec, as spatial_cv() returns.

Examples

random_cv(k = 10)

Random k-fold over the train+validation pool.

Description

Ignores geography by construction. On spatially clustered data that makes the validation metric optimistic – the pipeline check measures it: a random split put ~27% of validation in the SAME 250 m pixel as a training point. Offered anyway, for two reasons: it is the right answer when the rows really are independent, and running it against spatial_folds() on the same data is the cleanest way to MEASURE how much geography is worth here.

Usage

random_folds(
  meta,
  k = 5L,
  test_frac = 0,
  test_ids = NULL,
  seed = 42L,
  group = "auto",
  calibration_frac = 0,
  calibration_ids = NULL
)

Arguments

Value

A fold_plan with k folds.

Examples

ex <- example_landscape()
meta <- data.frame(sample_id = seq_len(nrow(ex$profiles)), x = ex$profiles$x,
                   y = ex$profiles$y)
random_folds(meta, k = 3, test_frac = 0.2)

Refit importance: what the model loses when trained without a variable.

Description

Leave one covariate out (LOCO; Lei et al. 2018, Williamson et al. 2023): each variable's channels are set to their training mean – every pixel, every window, every row – and the final run's seeds are trained again, everything else as the run trained them: its split, scaling, schedule, thread count, and the seed. The test skill lost is the variable's importance. A permutation asks how much the fitted network relies on a variable; this asks whether the variable is needed at all: variables that carry the same information each cost little, because a network trained without one finds it in the other.

Usage

refit_importance(
  metric = c("ccc", "rmse", "rmse_transform"),
  check_seeds = 1L,
  run_id = NULL,
  resume = TRUE,
  output_dir = NULL,
  n_cores = NULL,
  max_ram_gb = NULL
)

Arguments

Details

Before any variable is left out, the first check_seeds seeds are trained again with nothing left out, and must give the run's own predictions – or the refit is not the run's training, and it stops. It trains (variables x seeds + check_seeds) networks, side by side as dsm_final() does: give dsm_importance() themes as groups, and a few of the run's seeds.

Value

An importance_spec, for dsm_importance().

Examples

refit_importance(check_seeds = 1)

Leave-region-out, on a grouping that already exists (biome, catchment, ...).

Description

Leave-region-out, on a grouping that already exists (biome, catchment, ...).

Usage

region_cv(group, k = NULL, test_frac = 0.15, seed = 42L, calibration_frac = 0)

Arguments

Value

A resample_spec, as spatial_cv() returns.

Examples

ex <- example_landscape()
region_cv(group = ex$profiles$survey)

Fold by a grouping column: no group is split across folds.

Description

For when the unit that must not leak is named rather than geometric – a survey campaign, a country, a laboratory, a soil map unit. With k = number of groups this is leave-one-group-out.

Usage

region_folds(
  meta,
  group,
  k = NULL,
  test_frac = 0,
  test_ids = NULL,
  seed = 42L,
  calibration_frac = 0,
  calibration_ids = NULL
)

Arguments

Value

A fold_plan.

Examples

ex <- example_landscape()
meta <- data.frame(sample_id = seq_len(nrow(ex$profiles)), x = ex$profiles$x,
                   y = ex$profiles$y)
region_folds(meta, group = ex$profiles$survey, k = 3)

Add a model to the registry.

Description

Add a model to the registry.

Usage

register_model(spec, overwrite = FALSE)

Arguments

Value

spec, invisibly.

Examples

mean_only <- model_spec(
  "mean_only", input = "table",
  fit = function(x, y, cfg, ...) list(mu = mean(y)),
  predict = function(object, x, ...) rep(object$mu, nrow(x)),
  description = "the training mean: a floor every model must beat")
register_model(mean_only, overwrite = TRUE)
list_models()      # dsm_train(data, model = "mean_only", ...) now fits it

Turn a spec into a fold plan against real points.

Description

Exported because a plan is worth looking at before spending a night on it: ⁠plan <- resolve_resampling(spatial_cv(k = 5), data); print(plan)⁠.

Usage

resolve_resampling(
  spec,
  data,
  test_ids = NULL,
  windows = NULL,
  verbose = TRUE,
  calibration_ids = NULL
)

Arguments

Value

A fold_plan, already checked with check_fold_plan().

Examples

ex <- example_landscape()
store <- dsm_prepare(ex$profiles, target = "soc_stock", raster_dir = ex$raster_dir,
                     windows = c(3, 7), out_dir = file.path(tempdir(), "landscape"),
                     percentage = "^clay_pct$", transform = "log1p", n_cores = 1,
                     overwrite = TRUE, verbose = FALSE)
data <- dsm_load(store, windows = integer(0), verbose = FALSE)
resolve_resampling(spatial_cv(k = 3, block_size = 0.05, test_frac = 0.2), data,
                   windows = c(3, 7))

SAGE: how much of the model's skill each variable carries, shared fairly.

Description

Shapley Additive Global importancE (Covert, Lundberg & Lee 2020): the Shapley values of the loss the model explains. The value of a set of variables is the loss of the model's predictions when only they are known – the others taken from background points of the model's own training data, the whole patch, the same point in every window, as in shap_importance()("kernel"). A variable's SAGE value is the loss it takes away; variables that carry the same information split it, where a permutation credits the shared part to neither and a refit without one gives each nothing (the other stands in). The values add up exactly to the loss of the mean prediction – no variable known – less the model's own.

Usage

sage_importance(
  loss = c("mse", "mae"),
  background = 16L,
  max_points = 200L,
  permutations = 64L,
  seed = 42L
)

Arguments

Details

Exact over every coalition up to 14 variables, by sampled permutations up to 40: meant for a few variables – give groups themes.

Value

An importance_spec, for dsm_importance().

Examples

sage_importance(background = 8, max_points = 100)

Score every trained unit of a tuning run on the held-out test set.

Description

Score every trained unit of a tuning run on the held-out test set.

Usage

score_test_grid(
  run_dir,
  data,
  transform = NULL,
  device,
  config_ids = NULL,
  allow_unfrozen = FALSE
)

Arguments

Value

An object of class "test_optimism".

Examples

run <- example_run()    # a small fitted run, made once a session
# every configuration on the test set, once the choice is frozen
score_test_grid(run$fit$run_dir, run$data, device = torch::torch_device("cpu"))

How much does a metric move when ONLY the seed changes?

Description

The noise floor of the whole experiment. Computed within (config, fold), so nothing but the random draw differs between the runs being compared; then summarised across all configs.

Usage

seed_noise_floor(comparison, metric = "val_ccc")

Arguments

Details

Any difference between two configs smaller than this is not evidence. That sentence is the entire reason the function exists, and it is why oneSE-style selection is worth having: the best mean is frequently the luckiest draw.

Value

A list: metric; by_config, one row per (config, fold) trained under more than one seed, with the mean, sd and range over the seeds; n_comparable, how many such rows; median_sd, the floor; and max_range, the widest spread seen.

Examples

# A tuning table: three configurations, each on three folds under two seeds
set.seed(1)
cfgs <- data.frame(config_id = c("cfg_001", "cfg_002", "cfg_003"),
                   skill = c(0.615, 0.62, 0.55), n_params = c(2e5, 2e6, 5e4))
cmp <- merge(expand.grid(config_id = cfgs$config_id, fold = 1:3, seed = 1:2,
                         stringsAsFactors = FALSE), cfgs)
cmp$status  <- "success"
cmp$val_ccc <- cmp$skill + rnorm(nrow(cmp), 0, 0.02)
nf <- seed_noise_floor(cmp, metric = "val_ccc")
nf$median_sd

The configuration a final run deployed.

Description

Read from the list of configurations in the order they were chosen – not from the grid's order, in which a runner-up can come first.

Usage

selected_config_id(summary, label = "this final run")

Arguments

Value

the config id, or stop() – "auto" must never travel on unresolved.

Examples

summary <- list(selected_config_ids = c("cfg_004", "cfg_011"))
selected_config_id(summary)   # the first chosen; the grid's order can differ

Set torch's thread pools.

Description

Split out of setup_torch_device() so that a caller who already holds a device can change the threads without building a second one – dsm_train(device = d, n_cores = 8) is that caller. The threads belong to the R session, not to a device: whatever set them last is what every later torch call gets.

Usage

set_torch_threads(n_threads = NULL)

Arguments

Details

The interop pool can be sized once per session, before torch first runs work in parallel; later calls leave it as it is and say so once.

Value

The number of threads, invisibly.

Examples

before <- Sys.getenv(c("OMP_NUM_THREADS", "MKL_NUM_THREADS"), unset = NA)
n <- torch::torch_get_num_threads()
set_torch_threads(2)
torch::torch_get_num_threads()
# the session as it was: torch's threads, and the two variables set here
torch::torch_set_num_threads(n)
Sys.unsetenv(names(before)[is.na(before)])
if (any(!is.na(before))) do.call(Sys.setenv, as.list(before[!is.na(before)]))

Configure torch threads and select compute device.

Description

Configure torch threads and select compute device.

Usage

setup_torch_device(n_threads = NULL, use_cuda = TRUE)

Arguments

Value

A torch_device object.

Examples

before <- Sys.getenv(c("OMP_NUM_THREADS", "MKL_NUM_THREADS"), unset = NA)
n <- torch::torch_get_num_threads()
device <- setup_torch_device(n_threads = 2, use_cuda = FALSE)
device
# the session as it was: torch's threads, and the two variables set here
torch::torch_set_num_threads(n)
Sys.unsetenv(names(before)[is.na(before)])
if (any(!is.na(before))) do.call(Sys.setenv, as.list(before[!is.na(before)]))

SHAP importance: how each variable pushes each prediction, up or down.

Description

Shapley values (Lundberg & Lee 2017) of the network's predictions, estimated from its gradients. Each point's prediction, less a reference prediction, is shared among the variables, with a sign: this one raised it, that one lowered it. Averaged in absolute value over the points they are a global importance – of how much each variable MOVES the predictions, which is not how much it improves them (permutation_importance() measures that); kept point by point (⁠$points⁠) they are the map of what drives the prediction where (Padarian et al. 2020; Wadoux et al. 2023 for soil carbon).

Usage

shap_importance(
  estimator = c("expected_gradients", "integrated_gradients", "kernel"),
  samples = 100L,
  background = NULL,
  steps = 50L,
  permutations = 64L,
  max_points = NULL,
  seed = 42L
)

Arguments

Details

The values are in the units the network predicts in: for a log1p target, a SHAP value of +0.1 raises the predicted value by about 10%. Each point's values add up to its prediction less the background's mean prediction (expected gradients) or less the prediction at the baseline (integrated gradients); dsm_importance() checks that they do and stops if they do not. SHAP by deep-network rules (DeepSHAP) is not offered: it needs a propagation rule written for every layer of this architecture.

Value

An importance_spec, for dsm_importance().

Examples

shap_importance()
shap_importance("kernel", background = 8, max_points = 50)

The conditional MEAN surface, from transform-space predictions.

Description

The conditional MEAN surface, from transform-space predictions.

Usage

smear(
  pred_transform,
  cal,
  lower_limit = 0,
  method = c("global", "level", "total")
)

Arguments

Value

Numeric, in native units: an estimate of E[y | x] rather than of its median.

Examples

set.seed(1)
f <- runif(500, 2, 5)                    # predictions, in log1p space
z <- f + rnorm(500, 0, 0.4)              # what was observed, in log1p space
cal <- smearing_factor(z, f)
expm1(c(2, 3, 4))                          # the median surface
smear(c(2, 3, 4), cal)                     # the mean surface, by Duan's factor
smear(c(2, 3, 4), cal, method = "level")   # by a factor that follows the level

The mean surface of a median map already made.

Description

The conditional mean is ⁠exp(z) S(z) - 1⁠ of the median z in transform space, so a mean map by any of smear()'s factors comes from the median map alone, cell by cell, without running the network again.

Usage

smear_map(
  median,
  cal,
  method = c("global", "level", "total"),
  lower_limit = 0,
  filename = ""
)

Arguments

Value

A SpatRaster: the conditional mean, cell by cell.

Examples

set.seed(1)
f <- runif(500, 2, 5)                    # predictions, in log1p space
z <- f + rnorm(500, 0, 0.4)              # what was observed, in log1p space
cal <- smearing_factor(z, f)
median_map <- terra::rast(nrows = 10, ncols = 10, vals = expm1(runif(100, 2, 5)))
smear_map(median_map, cal)

The smearing factors on held-out points: which one unbiases the mean surface.

Description

Each factor of smear() applied to held-out predictions – points none of them was calibrated on, a test set – against what was observed there: the bias of the total, which a total taken from the map inherits, and the bias within quintiles of the prediction, where a factor by level should help. The median surface, with no factor, is the first row. On the global SOC model, calibrated on three folds, the global factor won (+2.7% against -3.9% by level and -5.8% for a total), because those folds saw half the points the deployed model saw; this measures it again for a run of yours.

Usage

smearing_check(
  cal,
  obs_transform,
  pred_transform,
  methods = c("global", "level", "total"),
  lower_limit = 0
)

Arguments

Value

A smearing_check: a tibble, one row per surface, with bias_pct (its total against the observed total, in %), mae, rmse, and bias_pct_q1 to bias_pct_q5 (within quintiles of the prediction, low to high).

Examples

set.seed(1)
f <- runif(1000, 2, 5)
z <- f + rnorm(1000, 0, 0.2 + 0.1 * (f - 2))   # the spread grows with the level
cal <- smearing_factor(z[1:500], f[1:500])
smearing_check(cal, z[501:1000], f[501:1000])  # on points it was not calibrated on

The smearing factor, from held-out residuals in transform space.

Description

The smearing factor, from held-out residuals in transform space.

Usage

smearing_factor(
  obs_transform,
  pred_transform,
  transform = c("log1p", "log"),
  bins = 5L
)

Arguments

Value

An object of class "smearing_cal".

Examples

set.seed(1)
f <- runif(500, 2, 5)                    # predictions, in log1p space
z <- f + rnorm(500, 0, 0.4)              # what was observed, in log1p space
cal <- smearing_factor(z, f)
cal

The smearing factor from a tuning run's out-of-fold predictions.

Description

ONE RESIDUAL PER POINT, not one per (point, seed).

Usage

smearing_from_run(run_dir, config_id, role = "validation", transform = "log1p")

Arguments

Details

The deployed prediction is the ensemble MEDIAN over seeds, so the residual being calibrated has to be the ensemble's and not one member's. This is the same reasoning cv_residuals() applies for the conformal interval, and for a while this function claimed to use "the same out-of-fold residuals" while quietly reading the per-seed rows instead: 9,276 rows for 3,092 points.

Two consequences, and the second is the one that matters:

• exp() is convex, so mean(exp(e)) grows with var(e). A single member's residual is noisier than the ensemble's (sd 0.6698 against 0.6507 here), which inflated S by 0.99% – almost exactly the 1.00% that the variance difference alone predicts.

• n read 9,276 when there were 3,092 exchangeable units, overstating the evidence threefold in the one line a reader would use to judge it.

The 1% is NOT why this was changed. The target is bracketed between 1.259 and 1.389 by two held-out measurements of the deployed ensemble itself, so no choice inside that range is defensible as more accurate. It was changed because the claim in the docs was false and the printed n was wrong.

Value

A smearing_cal, or NULL when the run wrote no usable predictions.

Examples

run <- example_run()    # a small fitted run, made once a session
smearing_from_run(run$fit$run_dir, run$final$selected_config_ids)

Moran's I by distance class: how far apart two points stop resembling each other.

Description

For each class of distance, Moran's I of value over the pairs of points that far apart: the mean product of their deviations from the mean, over the variance. Near zero, points that far apart are as unlike as any two; the distance where it falls to zero is the scale blocks should exceed. For two models compared at the same points, value is the difference of their errors – the quantity whose mean is compared.

Usage

spatial_correlogram(
  x,
  y,
  value,
  breaks_km = c(0, 5, 10, 25, 50, 100, 200, 400, 800),
  coords = c("lonlat", "metres")
)

Arguments

Value

A spatial_correlogram: a tibble with from_km, to_km, n_pairs, moran_i, and expected, Moran's I with no spatial structure (-1/(n-1)).

Examples

ex <- example_landscape()
obs <- ex$profiles$soc_stock
set.seed(1)
pred_a <- obs * exp(rnorm(length(obs), 0, 0.20))   # two models' predictions
pred_b <- obs * exp(rnorm(length(obs), 0, 0.25))
d <- abs(pred_a - obs) - abs(pred_b - obs)    # what a comparison of the two compares
cg <- spatial_correlogram(ex$profiles$x, ex$profiles$y, d,
                          breaks_km = c(0, 1, 2, 5, 10, 20))
cg
plot(cg)

Spatially blocked k-fold: whole blocks of ground go to one fold.

Description

Spatially blocked k-fold: whole blocks of ground go to one fold.

Usage

spatial_cv(
  k = 5L,
  block_size = "auto",
  buffer = "auto",
  test_frac = 0.15,
  max_share = 0.1,
  buffer_metric = c("chebyshev", "euclidean"),
  seed = 42L,
  calibration_frac = 0
)

Arguments

Value

A resample_spec. resolve_resampling() turns it into folds against the points, and dsm_train() does that itself.

Examples

spatial_cv(k = 5, block_size = "auto", buffer = "auto")
# 15% of the points held out to calibrate split conformal intervals
spatial_cv(k = 5, calibration_frac = 0.15)

Spatially blocked folds of a point table.

Description

The fold constructor behind spatial_cv(): the points are binned into square blocks, whole blocks go to one fold, and the buffer then drops the training points too close to a validation or test point.

Usage

spatial_folds(
  meta,
  k = 5L,
  test_frac = 0,
  block_size = NULL,
  buffer = NULL,
  buffer_metric = c("chebyshev", "euclidean"),
  test_ids = NULL,
  seed = 42L,
  blocks_per_fold = 10L,
  calibration_frac = 0,
  calibration_ids = NULL
)

Arguments

Value

A fold_plan with k folds.

Examples

ex <- example_landscape()
meta <- data.frame(sample_id = seq_len(nrow(ex$profiles)), x = ex$profiles$x,
                   y = ex$profiles$y)
# blocks of 0.05 degrees, and a buffer of the widest window (7 cells of 0.0025)
spatial_folds(meta, k = 3, test_frac = 0.2, block_size = 0.05, buffer = 7 * 0.0025)

Aggregate a unit-level comparison table into one row per config.

Description

Aggregate a unit-level comparison table into one row per config.

Usage

summarise_resamples(
  comparison,
  metrics = c("val_ccc", "val_mae", "val_rmse", "val_r2", "val_mqi", "test_ccc",
    "test_mae"),
  by = "config_id"
)

Arguments

Value

One row per config: n, mean, sd and se of each metric, ranked by the first metric. Failed units are excluded from the statistics but counted in n_failed, because a config that crashes 2 runs in 3 is not the same as one that completed all three.

Examples

# A tuning table: three configurations, each on three folds under two seeds
set.seed(1)
cfgs <- data.frame(config_id = c("cfg_001", "cfg_002", "cfg_003"),
                   skill = c(0.615, 0.62, 0.55), n_params = c(2e5, 2e6, 5e4))
cmp <- merge(expand.grid(config_id = cfgs$config_id, fold = 1:3, seed = 1:2,
                         stringsAsFactors = FALSE), cfgs)
cmp$status  <- "success"
cmp$val_ccc <- cmp$skill + rnorm(nrow(cmp), 0, 0.02)
summarise_resamples(cmp, metrics = "val_ccc")

The forward and inverse functions for a named target transform.

Description

The forward and inverse functions for a named target transform.

Usage

target_transform_spec(name)

Arguments

Value

list(name, forward, inverse).

Examples

tr <- target_transform_spec("log1p")
tr$forward(10)
tr$inverse(tr$forward(10))