Find the best DLNM model structure using grouped validation

Evaluates combinations of exposure variables and cross-basis dimensions using grouped validation, then ranks candidate models with family-aware predictive performance metrics. Three grouped validation strategies are available:

Usage

find_bestfit(
  data,
  response = "y",
  group = "epi_id",
  time = "time",
  vars,
  max_lag,
  df_var_grid = c(3, 4, 5),
  df_lag_grid = c(3, 4, 5),
  min_vars = 1,
  max_vars = NULL,
  var_sets = NULL,
  fun_var = "ns",
  fun_lag = "ns",
  model_engine = "glmmTMB",
  family = "beta",
  random_effect = NULL,
  random_effect_prior = NULL,
  spatial_effect = NULL,
  spatial_structure = "matern",
  spatial_group = NULL,
  min_success = 2,
  rank_metric = NULL,
  threshold = 0.5,
  top_n = Inf,
  keep_fits = FALSE,
  verbose = TRUE,
  validation_method = c("LOOCV", "k-fold", "holdout"),
  k = 5,
  test_prop = 0.3,
  seed = NULL,
  ...
)

Arguments

data

Long-format data frame containing the response, grouping column, chronological time column, and candidate exposure variables.

response

Character scalar naming the response column in `data`. The response may be repeated over exposure-profile rows but must be constant within each complete validation group because one outcome is predicted per group.

group

Character scalar naming the independent grouped-validation unit. Every row belonging to one group remains together. Internally, `data_long$epi_id` is a canonical alias of `data_long[[group]]`; it does not redefine the user's grouping structure.

time

Character scalar naming the chronological time column. Time must be finite and unique within groups. Within every group observations must be complete and equally spaced, and the time step must be the same across groups. Every group must contain exactly `max_lag + 1` time points. Because all candidate exposure variables are finite columns of the same validated long-format rows, all variables necessarily use the same temporal support. Rows are ordered internally by `group` and `time`.

vars

Unique character vector of candidate exposure variables. Exposure variables must be distinct from `response`, `group`, and `time`.

max_lag

Non-negative integer maximum lag, expressed in observation steps. A two-element form such as `c(0, 85)` is accepted and normalized to its maximum. Because all groups must use the same time step, a lag index has the same temporal meaning across groups.

df_var_grid

Positive integer candidate dimensions for the exposure-response basis.

df_lag_grid

Positive integer candidate dimensions for the lag-response basis.

min_vars

Positive integer minimum number of exposure variables.

max_vars

Positive integer maximum number of exposure variables, or `NULL` to use `length(vars)`.

var_sets

Optional list of exact variable combinations. When supplied, `min_vars` and `max_vars` do not determine the combinations, although their basic argument validation is still applied.

fun_var

Character exposure-basis function passed to `define_exposures()`.

fun_lag

Character lag-basis function passed to `define_exposures()`.

model_engine

Character model engine supported by `fit_epidlnm()`: `"glm"`, `"glmmTMB"`, `"gam"`, `"gamm"`, `"gls"`, `"spamm"`, `"brms"`, `"inla"`, or `"bdlnm"`.

family

Character or family object passed to `prepare_response()` and `fit_epidlnm()`. Supported EpiExposure v1 canonical families are `"beta"`, `"binomial"`, `"poisson"`, `"gamma"`, `"gaussian"`, and `"negative_binomial"` (NB2). Ordinal outcomes and NB1 are not supported in EpiExposure v1. Family-object links are checked against the selected engine before validation and passed unchanged to `fit_epidlnm()`. To specify a non-default link, supply a supported family object, such as `family = stats::Gamma(link = "inverse")`; character family names use the EpiExposure default link.

random_effect

Optional character scalar naming one grouping column used as a conventional random intercept by `fit_epidlnm()`. It must be constant within each validation group. Random slopes, nested or crossed random-effect specifications, and arbitrary engine-specific random-effect expressions are not supported.

For `brms`, the random intercept is fitted as `(1 | group)`. For the INLA-backed `inla` and `bdlnm` engines, it is fitted as `INLA::f(group, model = "iid")`. EpiExposure v1 does not expose structured INLA latent effects such as `"rw1"`, `"rw2"`, `"ar1"`, `"besag"`, `"bym"`, `"bym2"`, or SPDE effects.

Random effects may contribute to model fitting, but validation predictions are always population-level and therefore set fitted random effects to zero. `glm` and `gls` do not accept `random_effect`; `gamm` requires it under the current EpiExposure v1 fitting contract.

random_effect_prior

Optional named list defining the hyperprior for the precision of the IID random intercept when `model_engine = "inla"` or `"bdlnm"`. The list is passed unchanged to the `hyper` argument of the internally generated `INLA::f(..., model = "iid")` term. For example, `list(prec = list(prior = "pc.prec", param = c(1, 0.01)))`.

`NULL` retains the engine default. A non-`NULL` value requires `random_effect` and is currently supported only for INLA-backed engines. The argument customizes the hyperprior of the IID structure; it does not select another latent model. The same prior is used in every validation training fit and in retained full-data refits.

spatial_effect

`NULL` (default), or two distinct numeric coordinate column names for a Matérn term fitted by `model_engine = "spamm"` only. Coordinates are constant within each epidemic but may be identical across different epidemics; they never redefine the exposure-profile unit.

spatial_structure

Only `"matern"` is supported (case-insensitive) when spatial coordinates are supplied; otherwise this setting is inert.

spatial_group

`NULL` for one shared field or a factor, character, or integer grouping column such as `"year"` for independent spatial fields. Requires spatial coordinates and must be constant within each epidemic. Spatial effects and `random_effect` are distinct and may coexist.

min_success

Positive integer of at least 2 giving the minimum number of **evaluation groups with a finite out-of-sample validation prediction** required for a candidate to be eligible for metric calculation and, when requested, full-data refitting. Under LOOCV and k-fold, all groups are evaluation groups. Under holdout, only the fixed test groups count toward `n_success`, `n_failed`, and `min_success`; holdout training groups are never counted as failed predictions. Candidate diagnostics retain `n_success_folds` and `n_failed_folds` as complete-validation-iteration counters: for holdout there is one iteration, successful only when every designated test group receives a finite prediction.

rank_metric

Optional character scalar selecting the primary ranking metric. If `NULL`, the default is `"CCC"` for non-binary outcomes and `"ROC_AUC"` for binary outcomes. Non-binary metrics are `"CCC"`, `"Cb"`, `"rho"`, `"RMSE"`, and `"MAE"`. Binary metrics are `"ROC_AUC"`, `"Brier"`, `"LogLoss"`, `"Accuracy"`, `"Balanced_Accuracy"`, `"Sensitivity"`, `"Specificity"`, `"F1"`, `"MCC"`, and `"Precision"`. Metric direction is resolved automatically.

threshold

Numeric probability strictly between 0 and 1 used only for binomial classification metrics. The default is `0.5`. Predicted probabilities are retained; thresholding is used only to create the auxiliary 0/1 class required by classification metrics. ROC AUC, Brier score, and Log Loss do not depend on the threshold.

top_n

Positive integer or `Inf`; number of ranked candidates returned.

keep_fits

Logical. Validation metrics are always calculated first and exclusively from out-of-sample predictions produced by the selected validation strategy. If `TRUE`, every candidate that reaches `min_success` is **then** refitted on the complete data using a full-data basis. These full-data refits do not generate, replace, or alter validation predictions and do not participate in metric calculation or ranking. After ranking and `top_n` filtering, only retained candidate fits are stored in `attr(result, "fits")`. They are intended for downstream visualization, epidemiological interpretation, exposure-lag-response surfaces, effect summaries, lag decomposition, scenario simulation, and later prediction.

verbose

Logical. Print the validation strategy, candidate grid, execution mode, concise warning summaries, and completion information.

validation_method

Character grouped-validation strategy. `"LOOCV"` leaves one complete group out per iteration. `"k-fold"` assigns complete groups to `k` approximately balanced folds. `"holdout"` creates one fixed grouped proportional training/test split and uses the test groups as the validation set for every candidate. All strategies operate on complete groups and never split temporal rows from one group across partitions.

k

Positive integer number of folds used only when `validation_method = "k-fold"`. It must satisfy `2 <= k < number of groups`. Fold sizes differ by at most one group. For binomial outcomes, allocation is additionally stratified by the group-level 0/1 response and each response class must contain at least `k` groups.

test_prop

Numeric scalar strictly between 0 and 1 giving the requested approximate proportion of complete groups reserved for the test partition when `validation_method = "holdout"`. The default `0.30` requests approximately 70 percent of groups for training and 30 percent for testing. It is validated for every call but ignored by LOOCV and k-fold. For binary outcomes, the effective holdout size may be adjusted to the closest feasible value that keeps both outcomes 0 and 1 in both training and test partitions.

seed

Optional non-negative integer controlling grouped allocation for `validation_method = "k-fold"` and the grouped train/test split for `validation_method = "holdout"`. The same data, group order, validation method, `test_prop`, and seed reproduce the same allocation. The caller's existing global `.Random.seed` is restored after allocation. `seed` does not alter LOOCV membership and is not passed to `fit_epidlnm()`.

...

Named additional arguments passed to `fit_epidlnm()`. Core arguments managed by `find_bestfit()` (`data`, `model_engine`, `family`, `random_effect`, `random_effect_prior`, `spatial_effect`, `spatial_structure`, `spatial_group`, `epiexposure_spec`, and `basis_objects`) cannot be supplied again through `...`. For INLA-backed engines, likelihood availability, `control.family$control.link$model` consistency, and `control.compute$config = TRUE` compatibility are checked before validation. `fit_epidlnm()` supplies required INLA controls during each fit.

Value

A data frame ranked by `rank_metric`. Non-binary families report `CCC`, `Cb`, `rho`, `RMSE`, and `MAE`; binomial models report `ROC_AUC`, `Brier`, `LogLoss`, `Accuracy`, `Balanced_Accuracy`, `Sensitivity`, `Specificity`, `F1`, `MCC`, and `Precision`.

Candidate diagnostics include `n_groups_total`, `n_predictions`, `n_success`, `n_failed`, `n_folds`, `n_success_folds`, and `n_failed_folds`. `n_predictions` is the number of groups assigned to out-of-sample evaluation: all groups for LOOCV/k-fold, but only test groups for holdout. The `*_folds` fields are retained as validation-iteration counters; holdout therefore has one iteration.

Diagnostic attributes include `"predictions"`, `"predictions_by_model"`, `"failures"`, `"warnings"`, `"warning_summary"`, `"timing"`, and `"candidate_times"`. With `keep_fits = TRUE`, `"fits"` contains retained full-data refits and `"retained_fit_scope"` is `"full_data_refit_after_validation"`; otherwise `"retained_fit_scope"` is `NA_character_`.

Grouped-validation metadata include:

- `"validation_method"`: `"LOOCV"`, `"k-fold"`, or `"holdout"`; - `"validation_scheme"`: normalized descriptive scheme; - `"k"`: number of LOOCV iterations for LOOCV, supplied `k` for k-fold, and `NA_integer_` for holdout; - `"validation_seed"`: supplied k-fold/holdout seed or `NA_integer_`; - `"validation_stratified"`: whether binomial stratification was used; - `"test_prop"`: requested holdout proportion, otherwise `NA_real_`; - `"n_groups_total"` and `"n_evaluation_groups"`; - `"training_groups"`, `"test_groups"`, and `"evaluation_groups"`; - `"validation_assignments"`: data frame with `group`, `partition`, `fold`, and `observed`; holdout uses fixed `train`/`test` partitions, whereas LOOCV/k-fold use `partition = "evaluation"` and the evaluation fold; - `"validation_balance"`: fold-level counts for LOOCV/k-fold or train/test counts for holdout, including 0/1 class counts for binomial outcomes.

Standard downstream metadata are retained in `"family"`, `"outcome_type"`, `"rank_metric"`, `"threshold"`, `"prediction_level"`, `"prediction_estimand"`, `"prediction_contract"`, `"random_effect"`, `"random_effect_model"`, `"random_effect_prior"`, `"spatial_effect"`, `"spatial_structure"`, `"spatial_group"`, `"spatial_term"`, `"has_spatial_effect"`, `"basis_training_only"`, `"max_lag"`, `"history_length"`, `"history_contract"`, and `"time_step"`.

Details

- `"LOOCV"`: grouped leave-one-out cross-validation; - `"k-fold"`: grouped k-fold cross-validation; - `"holdout"`: one grouped proportional holdout validation split.

The validation unit is always the complete column named by `group`. The function creates the internal canonical alias `epi_id` from that column, but never splits a group's exposure-profile rows across training and evaluation. With the default `group = "epi_id"`, complete epidemics are therefore kept intact in every strategy.

The statistical target used for every validation prediction is the same target used by the EpiExposure v1 prediction layer: the population/fixed-component expected response, with fitted random effects set to zero. Model ranking uses deterministic predictions from the harmonized central parameter estimate and does not propagate coefficient, posterior, residual, or future-observation uncertainty.

## Grouped validation and cross-validation

Validation membership is constructed once before candidate evaluation and is reused unchanged for every candidate. This guarantees fair comparison on the same out-of-sample groups. LOOCV produces one evaluation iteration per group. Grouped k-fold randomizes complete groups into approximately balanced folds; for binomial outcomes, allocation is stratified so both classes occur in every fold whenever the documented class-count requirement is satisfied.

## Grouped proportional holdout

With `validation_method = "holdout"`, the function creates exactly one grouped training/test split before candidate evaluation. `test_prop` determines the requested approximate fraction of complete groups assigned to testing. For non-binary outcomes, test groups are sampled directly. For binary outcomes, sampling is stratified at the **group-level response**, not at longitudinal row level, and both outcomes 0 and 1 are required in both partitions. At least two complete groups in each class are therefore required.

Every candidate uses the same fixed holdout partition. `define_exposures()` is called only on the holdout training groups, so exposure-basis knots, boundary knots, effective dimensions, `argvar`, `arglag`, templates, and related cross-basis attributes are learned exclusively from training data. The stored training template is then transported unchanged to the holdout test profiles by `.build_design_from_templates()`. Test values never redefine the basis.

Holdout metrics are calculated exclusively from predictions for the fixed test groups. Because that same test subset is used to compare and rank candidate structures, it is a **validation** holdout, not an untouched final external test set. A final independent external evaluation requires another data set or partition that was never used for candidate selection. If `keep_fits = TRUE`, eligible candidates are refitted on all available data only after their validation metrics have already been calculated; these refits are downstream models and never replace holdout predictions.

## Training-only basis construction

For every validation iteration, `define_exposures()` receives only the corresponding training groups. The effective `argvar`, `arglag`, lag range, spline knots, boundary knots, and other returned basis attributes are reused to transform the evaluation profiles. For highly skewed predictors or many repeated values, `splines::ns()` may warn when an interior knot coincides with a boundary knot; the effective training basis returned after that adjustment is the template transported to evaluation data.

The function validates basis transport explicitly. Requested `max_lag`, the lag stored in training cross-bases, and the lag stored in the exposure specification must agree. Reconstructed training and evaluation cross-bases must have matching dimensions and native column names/order. Canonical `cb_<variable>_<index>` names are assigned deterministically. For `model_engine = "bdlnm"`, an epidemic-level matrix-form cross-basis is built from the same training parameters, numerically checked against the canonical design, and exposure-prefixed so internal INLA names remain globally unique.

## Temporal requirements

Every complete group must contain exactly `max_lag + 1` ordered, equally spaced observations. Profiles are never truncated, padded, or silently realigned. All candidate exposures are finite columns of the same validated long-format rows and therefore share identical temporal support.

## Random-effect scope

Candidate models may include one conventional random intercept through `random_effect`. The same random-intercept specification is used in every validation training fit and, when `keep_fits = TRUE`, in the full-data refits retained after validation.

For Bayesian engines, EpiExposure v1 supports `(1 | group)` through `brms` and `f(group, model = "iid")` through `inla` and `bdlnm`. The candidate grid does not search over alternative latent-effect structures, and `find_bestfit()` does not support INLA models such as `"rw1"`, `"rw2"`, `"ar1"`, `"besag"`, `"bym"`, `"bym2"`, or SPDE effects.

This restriction keeps the validation target and population-level prediction contract harmonized across candidates and engines. All fitted random-effect contributions are set to zero when out-of-sample validation predictions are generated.

## Optional spatial covariance

Only `spaMM` accepts `spatial_effect = c("x_coord", "y_coord")` through the EpiExposure v1 interface. This adds a Matérn term during fitting; `spatial_group` may define independent fields. This spaMM-specific interface is distinct from INLA/SPDE or areal latent-effect models, which are not currently exposed by `find_bestfit()`.

Conventional random effects remain distinct. All fitted random and spatial effects are set to zero for population-level validation predictions.

## Validation prediction target, warnings, and ranking

Every successful evaluation group receives one deterministic out-of-sample expected-response prediction from the harmonized central fixed/population parameter estimate. For LOOCV and k-fold these are out-of-fold predictions; for holdout they are out-of-sample holdout predictions. Coefficient/posterior uncertainty and residual/future-observation noise are excluded from ranking.

Warnings raised during fitting/prediction are captured and classified as `"knot"`, `"convergence"`, `"hessian"`, or `"other"`; warnings do not by themselves make a validation iteration fail. Candidate ranking uses `rank_metric` first, with the remaining family-appropriate metrics as deterministic tie-breakers in canonical order.

Candidate-level parallelism is preserved. When the caller configures a `future` plan with more than one worker, one future owns each candidate and that candidate's validation iterations remain sequential. The function does not change the caller's global `future::plan()`.

Examples

# \donttest{
best_models <- find_bestfit(
  data = epi_data,
  response = "y",
  group = "epi_id",
  time = "time",
  vars = "wetness",
  max_lag = 85,
  df_var_grid = 2,
  df_lag_grid = 2,
  var_sets = list("wetness"),
  model_engine = "glmmTMB",
  family = "beta",
  validation_method = "holdout",
  test_prop = 0.20,
  seed = 123,
  min_success = 2,
  top_n = 1,
  keep_fits = FALSE,
  verbose = FALSE
)

head(best_models)
#>   rank model_id df_var df_lag    vars n_vars       CCC        Cb       rho
#> 1    1        1      2      2 wetness      1 0.6431441 0.8898498 0.7227558
#>        RMSE       MAE n_groups_total n_folds n_success_folds n_failed_folds
#> 1 0.1976386 0.1707208            520       1               1              0
#>   n_predictions n_success n_failed n_warning_folds warning_rate
#> 1           104       104        0               0            0
#>   n_warning_events n_knot_warning_folds n_convergence_warning_folds
#> 1                0                    0                           0
#>   n_hessian_warning_folds n_other_warning_folds n_knot_warning_events
#> 1                       0                     0                     0
#>   n_convergence_warning_events n_hessian_warning_events n_other_warning_events
#> 1                            0                        0                      0
# }