T1FF: Classification, Regression, Forecasting, Tuning, and Nested Validation

Nihat Tak

2026-09-05

1 Overview

T1FF fits Type-1 Fuzzy Function models for binary classification, numeric regression, and time-series forecasting with user-supplied temporal predictors. The method partitions the predictor space by fuzzy C-means, estimates a local statistical model for every fuzzy cluster, and combines local predictions using observation-specific membership degrees.

The package provides a complete modeling workflow:

  1. T1FF() fits a classification or regression model.
  2. predict() returns final predictions, local predictions, or memberships.
  3. tune.T1FF() selects the cluster count and fuzziness parameter.
  4. evaluate.T1FF() calculates task-appropriate performance measures.
  5. benchmark.T1FF() performs repeated nested cross-validation against a standard logistic or linear regression baseline.

The regression workflow can also be used for forecasting after temporal dependence has been encoded with lagged, seasonal, or trend predictors and the training, validation, and test partitions have been kept in chronological order.

The current release supports binary classification and numeric regression. Multiclass classification is outside the scope of this release.

2 Method

Let \(x_i\) denote the predictor vector for observation \(i\). Fuzzy C-means estimates \(c\) cluster centers and memberships \(\mu_{ik}\), where

\[ 0 \leq \mu_{ik} \leq 1, \qquad \sum_{k=1}^{c}\mu_{ik}=1. \]

For each cluster \(k\), the local model uses the original predictors together with three membership transformations:

\[ \mu_{ik}, \qquad \log\{\max(\mu_{ik},\epsilon)\}, \qquad \exp(\mu_{ik}). \]

By default, binary classification uses a local logistic model and regression uses a local ordinary least-squares model. Alternatively, local_model = "svm" uses probabilistic support-vector classifiers or epsilon-support-vector regression within the fuzzy clusters. If \(\widehat{p}_{ik}\) is the positive-class probability from local classifier \(k\), the final probability is

\[ \widehat{p}_i = \sum_{k=1}^{c}\mu_{ik}\widehat{p}_{ik}. \]

For regression, local responses are combined in the same way:

\[ \widehat{y}_i = \sum_{k=1}^{c}\mu_{ik}\widehat{y}_{ik}. \]

Predictors are standardized by default. All preprocessing parameters are estimated from training data and reused for new observations. During cross-validation, preprocessing is re-estimated inside each training fold to avoid information leakage.

3 Installation and package contents

After installing the source package, load it with:

library(T1FF)
packageVersion("T1FF")
#> [1] '0.1.0'

List the exported functions and open the package help index with:

ls("package:T1FF")
help(package = "T1FF")

4 Binary classification

The following example removes setosa to create a two-class outcome. virginica is explicitly defined as the positive class.

data(iris)
iris_binary <- droplevels(subset(iris, Species != "setosa"))
table(iris_binary$Species)
#> 
#> versicolor  virginica 
#>         50         50

Both a column interface and a formula interface are available.

classification_fit <- T1FF(
  Species ~ .,
  data = iris_binary,
  c = 2,
  m = 2,
  task = "classification",
  positive_class = "virginica",
  seed = 1
)
classification_fit
#> Type-1 Fuzzy Functions (T1FF)
#> ============================
#> Task            : classification 
#> Target          : Species 
#> Negative class  : versicolor 
#> Positive class  : virginica 
#> Local converged : 2 / 2 
#> Local engines  : ridge, ridge 
#> Separation     : 2 cluster(s) detected
#> GLM warnings    : 2 captured
#> Prob. clip     : 1e-06 
#> Local learner    : linear 
#> Clusters (c)    : 2 
#> Fuzziness (m)   : 2 
#> Scaled features : TRUE 
#> Features        : Sepal.Length, Sepal.Width, Petal.Length, Petal.Width 
#> Local models    : 2

The equivalent column interface is:

classification_fit <- T1FF(
  da = iris_binary,
  target_col = "Species",
  c = 2,
  m = 2,
  task = "classification",
  positive_class = "virginica",
  seed = 1
)

4.1 Prediction types

type = "prob" returns the final membership-weighted probability of the positive class. type = "class" applies a classification threshold.

predict(classification_fit, iris_binary[1:5, ], type = "prob")
#> [1] 5.312693e-05 1.419160e-04 4.080346e-03 1.902361e-04 1.227020e-03
predict(classification_fit, iris_binary[1:5, ], type = "class")
#> [1] versicolor versicolor versicolor versicolor versicolor
#> Levels: versicolor virginica

The local cluster probabilities and fuzzy memberships are available for model inspection.

predict(classification_fit, iris_binary[1:3, ], type = "cluster_prob")
#>         cluster_1    cluster_2
#> [1,] 4.753114e-05 6.154443e-05
#> [2,] 1.280531e-04 1.536250e-04
#> [3,] 3.908832e-03 4.422204e-03
predict(classification_fit, iris_binary[1:3, ], type = "membership")
#>      cluster_1 cluster_2
#> [1,] 0.6006802 0.3993198
#> [2,] 0.4578856 0.5421144
#> [3,] 0.6659061 0.3340939

The classification threshold can be changed without refitting:

predict(classification_fit, iris_binary[1:5, ],
        type = "class", threshold = 0.40)
#> [1] versicolor versicolor versicolor versicolor versicolor
#> Levels: versicolor virginica

4.2 Separation and stable probabilities

Cluster-specific logistic regression can encounter complete or quasi-complete separation. The default logistic_method = "auto" first fits an ordinary GLM, detects non-convergence, extreme coefficients, or extreme fitted probabilities, and refits affected local models using dependency-free ridge logistic regression.

stable_fit <- T1FF(
  Species ~ ., iris_binary,
  c = 2,
  logistic_method = "auto",
  ridge_lambda = 0.01,
  probability_clip = 1e-6,
  positive_class = "virginica",
  seed = 1
)
summary(stable_fit)
#> Type-1 Fuzzy Functions (T1FF)
#> ============================
#> Task            : classification 
#> Target          : Species 
#> Negative class  : versicolor 
#> Positive class  : virginica 
#> Local converged : 2 / 2 
#> Local engines  : ridge, ridge 
#> Separation     : 2 cluster(s) detected
#> GLM warnings    : 2 captured
#> Prob. clip     : 1e-06 
#> Local learner    : linear 
#> Clusters (c)    : 2 
#> Fuzziness (m)   : 2 
#> Scaled features : TRUE 
#> Features        : Sepal.Length, Sepal.Width, Petal.Length, Petal.Width 
#> Local models    : 2

Available strategies are:

probability_clip prevents exact zero and one predictions. This protects log-loss calculations without changing the predicted class at ordinary thresholds. The penalty and clipping value should be reported in reproducible analyses.

4.3 Support-vector local models

Set local_model = "svm" to replace the local logistic or linear models with support-vector models. Classification uses probability-enabled classifiers; regression uses epsilon-support-vector regression. The fuzzy memberships and their transformations remain local-model inputs, and final predictions are still combined by observation-specific memberships.

svm_fit <- T1FF(
  Species ~ ., iris_binary,
  c = 2,
  m = 2,
  local_model = "svm",
  svm_kernel = "rbfdot",
  svm_C = 1,
  positive_class = "virginica",
  seed = 8
)
predict(svm_fit, iris_binary[1:5, ], type = "prob")
#> [1] 0.02847482 0.01747250 0.14774417 0.03345157 0.03639519

When svm_sigma = NULL, the RBF scale is estimated separately in every local model. For regression, svm_epsilon controls the epsilon-insensitive loss.

5 Hyperparameter tuning

tune.T1FF() searches combinations of the cluster count c and fuzziness parameter m. Binary classification uses stratified K-fold cross-validation by default; regression uses ordinary K-fold cross-validation.

tuned_fit <- tune.T1FF(
  Species ~ ., iris_binary,
  task = "classification",
  c_values = 2:3,
  m_values = c(1.5, 2),
  metric = "logloss",
  resampling = "stratified_kfold",
  folds = 3,
  positive_class = "virginica",
  seed = 2,
  verbose = FALSE
)
summary(tuned_fit)
#> T1FF tuning results
#> ===================
#> Task       : classification 
#> Metric     : logloss (minimize) 
#> Resampling : stratified_kfold 
#> Folds      : 3 
#> Repeats    : 1 
#> Best c     : 2 
#> Best m     : 2 
#> Best threshold: 0.5 
#> Best score : 0.1122412 
#> 
#>  c   m threshold mean_score   sd_score valid_splits total_splits
#>  2 2.0       0.5  0.1122412 0.06848711            3            3
#>  3 1.5       0.5  0.1149799 0.08385234            3            3
#>  2 1.5       0.5  0.1239898 0.05662326            3            3
#>  3 2.0       0.5  0.1253452 0.06523663            3            3

The tuning result stores the full grid, the selected parameters, and a final model refitted on all supplied observations.

tuned_fit$results
#>   c   m threshold mean_score   sd_score valid_splits total_splits
#> 1 2 1.5       0.5  0.1239898 0.05662326            3            3
#> 2 3 1.5       0.5  0.1149799 0.08385234            3            3
#> 3 2 2.0       0.5  0.1122412 0.06848711            3            3
#> 4 3 2.0       0.5  0.1253452 0.06523663            3            3
c(c = tuned_fit$best_c, m = tuned_fit$best_m)
#> c m 
#> 2 2
predict(tuned_fit, iris_binary[1:5, ], type = "prob")
#> [1] 5.308239e-05 1.418200e-04 4.078445e-03 1.901028e-04 1.226182e-03

For threshold-dependent objectives, the decision threshold can be tuned jointly with c and m without refitting the same fold models for every threshold.

threshold_tuned_fit <- tune.T1FF(
  Species ~ ., iris_binary,
  c_values = 2,
  m_values = 2,
  threshold_values = seq(0.3, 0.7, by = 0.1),
  metric = "balanced_accuracy",
  folds = 3,
  positive_class = "virginica",
  seed = 7,
  verbose = FALSE
)
threshold_tuned_fit$best_threshold
#> [1] 0.3
predict(threshold_tuned_fit, iris_binary[1:5, ], type = "class")
#> [1] versicolor versicolor versicolor versicolor versicolor
#> Levels: versicolor virginica

The tuned threshold is automatically used by predict() for class prediction. An explicit threshold supplied to predict() overrides it. Multiple candidate thresholds apply only to accuracy, balanced accuracy, F1, sensitivity, and specificity; probability and ranking metrics do not depend on a decision threshold.

Classification tuning metrics are logloss, brier, roc_auc, pr_auc, accuracy, balanced_accuracy, f1, sensitivity, and specificity. Regression metrics are rmse, mse, mae, mape, smape, and r2. MAPE excludes zero actual values; both MAPE and SMAPE are reported as percentages. Use probability-based metrics such as log-loss when calibrated probabilities are important and threshold metrics such as F1 when a particular decision rule is central.

6 Model evaluation

evaluate.T1FF() accepts either a fitted model or a tuned model. For classification it reports a confusion matrix and discrimination, calibration, and threshold-dependent metrics. It first obtains predictions internally with predict() (type = "prob" for classification and type = "response" for regression), then compares them with the supplied truth. Therefore, call predict() directly when predictions alone are needed; call evaluate.T1FF() when known outcome values should be used to assess those predictions.

classification_evaluation <- evaluate.T1FF(
  tuned_fit,
  iris_binary,
  truth = "Species",
  threshold = 0.5
)
classification_evaluation
#> T1FF model evaluation
#> =====================
#> Task        : classification 
#> Truth       : Species 
#> Rows used   : 100 / 100 
#> Positive    : virginica 
#> Threshold   : 0.5 
#> 
#> Confusion matrix:
#>             Predicted
#> Actual       versicolor virginica
#>   versicolor         49         1
#>   virginica           1        49
#> 
#> Metrics:
#>             metric      value
#>           accuracy 0.98000000
#>  balanced_accuracy 0.98000000
#>        sensitivity 0.98000000
#>        specificity 0.98000000
#>          precision 0.98000000
#>             recall 0.98000000
#>                 f1 0.98000000
#>            roc_auc 0.99840000
#>             pr_auc 0.99846792
#>            logloss 0.05311984
#>              brier 0.01681536

Components can be extracted for custom tables:

classification_evaluation$metrics
#>               metric      value
#> 1           accuracy 0.98000000
#> 2  balanced_accuracy 0.98000000
#> 3        sensitivity 0.98000000
#> 4        specificity 0.98000000
#> 5          precision 0.98000000
#> 6             recall 0.98000000
#> 7                 f1 0.98000000
#> 8            roc_auc 0.99840000
#> 9             pr_auc 0.99846792
#> 10           logloss 0.05311984
#> 11             brier 0.01681536
classification_evaluation$confusion_matrix
#>             Predicted
#> Actual       versicolor virginica
#>   versicolor         49         1
#>   virginica           1        49

Evaluating on training data describes apparent fit and can be optimistic. Use an untouched test set or benchmark.T1FF() when estimating generalization performance.

7 Nested cross-validation benchmark

Nested cross-validation separates model selection from model assessment. For each outer split, benchmark.T1FF() tunes c and m using only the outer training observations. It then assesses the selected model on the untouched outer test observations. A standard logistic regression baseline is evaluated on the same splits.

The small example below is intended to keep vignette build time short. Applied studies should generally use at least 5 outer folds, 5 inner folds, and multiple outer repeats when sample size permits.

classification_benchmark <- benchmark.T1FF(
  Species ~ ., iris_binary,
  task = "classification",
  c_values = 2,
  m_values = 2,
  tune_metric = "logloss",
  metrics = c("roc_auc", "pr_auc", "logloss", "brier", "f1"),
  outer_folds = 2,
  inner_folds = 2,
  repeats = 1,
  positive_class = "virginica",
  seed = 3,
  verbose = FALSE
)
classification_benchmark
#> T1FF nested cross-validation benchmark
#> ======================================
#> Task        : classification 
#> Outer folds : 2 
#> Inner folds : 2 
#> Repeats     : 1 
#> Tune metric : logloss 
#> Confidence  : 95% 
#> 
#>     model  metric n       mean          sd           se   conf_low conf_high
#>  Baseline   brier 2 0.05999988 0.028284215 0.0199999600 0.00000000 0.3141235
#>      T1FF   brier 2 0.04505955 0.010389587 0.0073465473 0.00000000 0.1384063
#>  Baseline      f1 2 0.93910256 0.031729150 0.0224358974 0.65402746 1.0000000
#>      T1FF      f1 2 0.95015699 0.016095837 0.0113814757 0.80554163 1.0000000
#>  Baseline logloss 2 0.82893157 0.390761620 0.2763101912 0.00000000 4.3397854
#>      T1FF logloss 2 0.15666365 0.015990773 0.0113071840 0.01299225 0.3003350
#>  Baseline  pr_auc 2 0.81903515 0.014639503 0.0103516916 0.68750444 0.9505659
#>      T1FF  pr_auc 2 0.99092531 0.001007275 0.0007122507 0.98187531 0.9999753
#>  Baseline roc_auc 2 0.94000000 0.028284271 0.0200000000 0.68587591 1.0000000
#>      T1FF roc_auc 2 0.99120000 0.001131371 0.0008000000 0.98103504 1.0000000

Detailed results include one row per outer split, model, and metric; parameter selections from every inner search; and elapsed-time summaries.

classification_benchmark$fold_results
#>    repeat_id fold    model  metric      value
#> 1          1    1     T1FF roc_auc 0.99200000
#> 2          1    1     T1FF  pr_auc 0.99163756
#> 3          1    1     T1FF logloss 0.16797083
#> 4          1    1     T1FF   brier 0.05240610
#> 5          1    1     T1FF      f1 0.93877551
#> 6          1    1 Baseline roc_auc 0.92000000
#> 7          1    1 Baseline  pr_auc 0.82938684
#> 8          1    1 Baseline logloss 1.10524176
#> 9          1    1 Baseline   brier 0.07999984
#> 10         1    1 Baseline      f1 0.91666667
#> 11         1    2     T1FF roc_auc 0.99040000
#> 12         1    2     T1FF  pr_auc 0.99021306
#> 13         1    2     T1FF logloss 0.14535646
#> 14         1    2     T1FF   brier 0.03771301
#> 15         1    2     T1FF      f1 0.96153846
#> 16         1    2 Baseline roc_auc 0.96000000
#> 17         1    2 Baseline  pr_auc 0.80868346
#> 18         1    2 Baseline logloss 0.55262138
#> 19         1    2 Baseline   brier 0.03999992
#> 20         1    2 Baseline      f1 0.96153846
classification_benchmark$selected_parameters
#>   repeat_id fold c m inner_score
#> 1         1    1 2 2   0.1067597
#> 2         1    2 2 2   0.2019673
classification_benchmark$timing_summary
#>      model mean_seconds   sd_seconds
#> 1 Baseline       0.0015 0.0007071068
#> 2     T1FF       0.0255 0.0021213203

Higher values are better for AUC, accuracy, F1, sensitivity, specificity, and \(R^2\). Lower values are better for log-loss, Brier score, RMSE, MSE, MAE, MAPE, and SMAPE. Repeated cross-validation folds are dependent; confidence intervals summarize variation across resamples but should not be interpreted as if all fold values were independent experimental observations.

8 Regression

Regression requires a finite numeric outcome.

regression_fit <- T1FF(
  mpg ~ wt + hp + disp,
  data = mtcars,
  c = 2,
  m = 2,
  task = "regression",
  seed = 4
)
regression_fit
#> Type-1 Fuzzy Functions (T1FF)
#> ============================
#> Task            : regression 
#> Target          : mpg 
#> Local learner    : linear 
#> Clusters (c)    : 2 
#> Fuzziness (m)   : 2 
#> Scaled features : TRUE 
#> Features        : wt, hp, disp 
#> Local models    : 2

Regression predictions include the combined response, cluster-specific responses, and memberships.

predict(regression_fit, mtcars[1:5, ], type = "response")
#> [1] 23.64968 22.69738 25.49542 20.36035 18.50277
predict(regression_fit, mtcars[1:3, ], type = "cluster_response")
#>      cluster_1 cluster_2
#> [1,]  23.10288  23.65734
#> [2,]  23.50167  22.67083
#> [3,]  24.15973  25.51036
predict(regression_fit, mtcars[1:3, ], type = "membership")
#>       cluster_1 cluster_2
#> [1,] 0.01382260 0.9861774
#> [2,] 0.03195720 0.9680428
#> [3,] 0.01106263 0.9889374
regression_evaluation <- evaluate.T1FF(
  regression_fit, mtcars, truth = "mpg"
)
regression_evaluation
#> T1FF model evaluation
#> =====================
#> Task        : regression 
#> Truth       : mpg 
#> Rows used   : 32 / 32 
#> 
#> Metrics:
#>  metric    value
#>    rmse 2.386298
#>     mse 5.694420
#>     mae 1.866707
#>    mape 9.601167
#>   smape 9.738728
#>      r2 0.838176

A nested regression benchmark uses linear regression as the baseline:

regression_benchmark <- benchmark.T1FF(
  mpg ~ ., mtcars,
  task = "regression",
  c_values = 2:4,
  m_values = c(1.5, 2, 2.5),
  tune_metric = "rmse",
  metrics = c("rmse", "mse", "mae", "mape", "smape", "r2"),
  outer_folds = 5,
  inner_folds = 5,
  repeats = 5,
  seed = 5
)

9 Time-series forecasting

T1FF treats the supplied rows as a predictive data set; it does not infer a time index or create temporal features automatically. To use the regression workflow for forecasting, temporal dependence should first be represented by response lags, seasonal variables, and/or trend terms. Training must use earlier observations, while validation and testing must use later observations.

The built-in tune.T1FF() resampling options assign observations randomly. They are suitable for ordinary cross-validation but do not constitute time-aware validation. Forecasting hyperparameters should therefore be chosen with an external chronological validation set or a rolling-origin loop.

The following example uses the base R AirPassengers series. It creates one-month and twelve-month lags, adds seasonal and trend terms, and reserves the latest 20% of observations as an untouched test period. Because each test row’s lag1 uses the most recently observed passenger count, this is a rolling one-step-ahead evaluation rather than a recursive multi-step forecast.

data(AirPassengers)
passengers <- as.numeric(AirPassengers)
period <- as.numeric(time(AirPassengers))
month <- as.numeric(cycle(AirPassengers))
n_time <- length(passengers)

forecast_data <- data.frame(
  period = period[13:n_time],
  y = passengers[13:n_time],
  lag1 = passengers[12:(n_time - 1)],
  lag12 = passengers[1:(n_time - 12)],
  trend = seq_len(n_time - 12),
  season_sin = sin(2 * pi * month[13:n_time] / 12),
  season_cos = cos(2 * pi * month[13:n_time] / 12)
)

time_split <- floor(0.80 * nrow(forecast_data))
train_time <- forecast_data[seq_len(time_split), ]
test_time <- forecast_data[(time_split + 1):nrow(forecast_data), ]

range(train_time$period)
#> [1] 1950.000 1958.667
range(test_time$period)
#> [1] 1958.750 1960.917

Fit only on the earlier period and predict the later period:

forecast_fit <- T1FF(
  y ~ lag1 + lag12 + trend + season_sin + season_cos,
  data = train_time,
  c = 2,
  m = 2,
  task = "regression",
  seed = 9
)

test_time$forecast <- predict(forecast_fit, test_time, type = "response")
head(test_time[, c("period", "y", "forecast")])
#>       period   y forecast
#> 106 1958.750 359 386.1115
#> 107 1958.833 310 343.1084
#> 108 1958.917 337 358.0561
#> 109 1959.000 360 365.4618
#> 110 1959.083 342 353.3427
#> 111 1959.167 406 388.7863

The observed test outcomes are used only after prediction to calculate out-of-sample performance metrics:

forecast_evaluation <- evaluate.T1FF(
  forecast_fit,
  test_time,
  truth = "y"
)
forecast_evaluation$metrics
#>   metric      value
#> 1   rmse  21.768312
#> 2    mse 473.859421
#> 3    mae  17.856863
#> 4   mape   4.073175
#> 5  smape   4.106509
#> 6     r2   0.925397

For a genuine multi-step forecast beyond the last observed value, future lags must instead be updated recursively from earlier forecasts. That data-preparation step is not automated by the current package.

10 Categorical predictors and missing values

Formula and column interfaces encode factor and character predictors using a design matrix learned from the training data. New observations must have compatible factor levels.

The default na_action = "fail" stops when training data contain missing values. With na_action = "omit", incomplete training rows are removed. At prediction time, incomplete rows receive NA predictions while original row order and length are preserved.

iris_missing <- iris_binary
iris_missing$Sepal.Length[c(2, 7)] <- NA_real_

missing_fit <- T1FF(
  Species ~ ., iris_missing,
  c = 2,
  positive_class = "virginica",
  na_action = "omit",
  seed = 6
)
predict(missing_fit, iris_missing[1:8, ], type = "prob")
#> [1] 5.417852e-05           NA 4.109074e-03 2.033966e-04 1.283804e-03
#> [6] 5.489761e-05           NA 1.000000e-06

11 Reproducibility and reporting

For reproducible analyses:

packageVersion("T1FF")
#> [1] '0.1.0'
sessionInfo()
#> R version 4.5.0 (2025-04-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS 26.3.1
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
#> 
#> time zone: Europe/Istanbul
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] T1FF_0.1.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.56      
#>  [5] e1071_1.7-17    cachem_1.1.0    knitr_1.51      htmltools_0.5.9
#>  [9] rmarkdown_2.30  lifecycle_1.0.5 cli_3.6.5       sass_0.4.10    
#> [13] proxy_0.4-29    jquerylib_0.1.4 class_7.3-23    compiler_4.5.0 
#> [17] kernlab_0.9-33  tools_4.5.0     evaluate_1.0.5  bslib_0.10.0   
#> [21] yaml_2.3.12     rlang_1.1.7     jsonlite_2.0.0

12 Current limitations

The current implementation is intentionally focused. It supports binary rather than multiclass classification, fuzzy C-means as the clustering method, and cluster-specific linear/logistic or support-vector response models. Nested benchmarking uses a standard GLM/LM baseline; comparisons with additional algorithms should be performed externally or added in a future release. Forecast lags, recursive multi-step updates, and rolling-origin resampling must currently be constructed externally. These limitations should be stated when reporting empirical results.