---
title: "T1FF: Classification, Regression, Forecasting, Tuning, and Nested Validation"
author: "Nihat Tak"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    number_sections: true
vignette: >
  %\VignetteIndexEntry{T1FF: Classification, Regression, Forecasting, Tuning, and Nested Validation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)
library(T1FF)
```

# 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.

# 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.

# Installation and package contents

After installing the source package, load it with:

```{r load}
library(T1FF)
packageVersion("T1FF")
```

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

```{r help-index, eval=FALSE}
ls("package:T1FF")
help(package = "T1FF")
```

# Binary classification

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

```{r classification-data}
data(iris)
iris_binary <- droplevels(subset(iris, Species != "setosa"))
table(iris_binary$Species)
```

Both a column interface and a formula interface are available.

```{r classification-fit}
classification_fit <- T1FF(
  Species ~ .,
  data = iris_binary,
  c = 2,
  m = 2,
  task = "classification",
  positive_class = "virginica",
  seed = 1
)
classification_fit
```

The equivalent column interface is:

```{r column-interface, eval=FALSE}
classification_fit <- T1FF(
  da = iris_binary,
  target_col = "Species",
  c = 2,
  m = 2,
  task = "classification",
  positive_class = "virginica",
  seed = 1
)
```

## Prediction types

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

```{r classification-prediction}
predict(classification_fit, iris_binary[1:5, ], type = "prob")
predict(classification_fit, iris_binary[1:5, ], type = "class")
```

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

```{r local-prediction}
predict(classification_fit, iris_binary[1:3, ], type = "cluster_prob")
predict(classification_fit, iris_binary[1:3, ], type = "membership")
```

The classification threshold can be changed without refitting:

```{r threshold}
predict(classification_fit, iris_binary[1:5, ],
        type = "class", threshold = 0.40)
```

## 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.

```{r stable-logistic}
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)
```

Available strategies are:

* `"auto"`: detect unstable GLMs and use ridge only when needed;
* `"glm"`: always retain ordinary local logistic models;
* `"ridge"`: always use local ridge-logistic models.

`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.

## 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.

```{r svm-fit}
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")
```

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

# 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.

```{r tuning}
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)
```

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

```{r tuning-components}
tuned_fit$results
c(c = tuned_fit$best_c, m = tuned_fit$best_m)
predict(tuned_fit, iris_binary[1:5, ], type = "prob")
```

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

```{r threshold-tuning}
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
predict(threshold_tuned_fit, iris_binary[1:5, ], type = "class")
```

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.

# 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.

```{r evaluation}
classification_evaluation <- evaluate.T1FF(
  tuned_fit,
  iris_binary,
  truth = "Species",
  threshold = 0.5
)
classification_evaluation
```

Components can be extracted for custom tables:

```{r evaluation-components}
classification_evaluation$metrics
classification_evaluation$confusion_matrix
```

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

# 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.

```{r benchmark}
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
```

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

```{r benchmark-components}
classification_benchmark$fold_results
classification_benchmark$selected_parameters
classification_benchmark$timing_summary
```

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.

# Regression

Regression requires a finite numeric outcome.

```{r regression-fit}
regression_fit <- T1FF(
  mpg ~ wt + hp + disp,
  data = mtcars,
  c = 2,
  m = 2,
  task = "regression",
  seed = 4
)
regression_fit
```

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

```{r regression-prediction}
predict(regression_fit, mtcars[1:5, ], type = "response")
predict(regression_fit, mtcars[1:3, ], type = "cluster_response")
predict(regression_fit, mtcars[1:3, ], type = "membership")
```

```{r regression-evaluation}
regression_evaluation <- evaluate.T1FF(
  regression_fit, mtcars, truth = "mpg"
)
regression_evaluation
```

A nested regression benchmark uses linear regression as the baseline:

```{r regression-benchmark, eval=FALSE}
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
)
```

# 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.

```{r forecasting-data}
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)
range(test_time$period)
```

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

```{r forecasting-fit}
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")])
```

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

```{r forecasting-evaluation}
forecast_evaluation <- evaluate.T1FF(
  forecast_fit,
  test_time,
  truth = "y"
)
forecast_evaluation$metrics
```

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.

# 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.

```{r missing-values}
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")
```

# Reproducibility and reporting

For reproducible analyses:

* set `seed` explicitly;
* report the candidate `c` and `m` grids;
* report inner and outer resampling designs;
* identify the positive class and classification threshold;
* report `local_model` and, for SVMs, the kernel, cost, sigma, and epsilon;
* report `logistic_method`, `ridge_lambda`, and `probability_clip`;
* retain fold-level benchmark results rather than only means;
* report package and R versions.

```{r reproducibility}
packageVersion("T1FF")
sessionInfo()
```

# 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.
