---
title: "Getting Started with NNS: Clustering and Regression"
author: "Fred Viole"
output:
  rmarkdown::html_vignette:
    check_title: false
vignette: >
  %\VignetteIndexEntry{07. Getting Started with NNS: Clustering and Regression}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(NNS)
options(mc.cores = 1)
RcppParallel::setThreadOptions(numThreads = 1)
Sys.setenv("OMP_THREAD_LIMIT" = 1)
```

```{r setup2, message=FALSE, warning=FALSE}
library(NNS)
require(knitr)
```


# Clustering and Regression
Below are some examples demonstrating unsupervised learning with NNS clustering and nonlinear regression using the resulting clusters.  As always, for a more thorough description and definition, please view the References.

## NNS Partitioning <span style="color:red">`NNS.part()`</span>
**`NNS.part`** is both a partitional and hierarchical clustering method.  `NNS` iteratively partitions the joint distribution into partial moment quadrants, and then assigns a quadrant identification (1:4) at each partition.

**`NNS.part`** returns a `data.frame` of observations along with their final quadrant identification.  It also returns regression points that summarize observations within occupied cells. Their location follows the chosen reduction rule; these points support the subsequent **`NNS.reg`** fit.
```{r linear}
x = seq(-5, 5, .05); y = x ^ 3

for(i in 1 : 4){NNS.part(x, y, order = i, Voronoi = TRUE, obs.req = 0)}
```


### X-only Partitioning
**`NNS.part`** offers a partitioning based on $x$ values only **`NNS.part(x, y, type = "XONLY", ...)`**, using the response observations within each predictor cell to form regression points. Consistency requires suitable shrinking cells and growing cell occupancy; a fixed order or occupancy cutoff alone does not establish those conditions.
```{r x part,results='hide'}
for(i in 1 : 4){NNS.part(x, y, order = i, type = "XONLY", Voronoi = TRUE)}
```

Note the partition identifications are limited to 1's and 2's (left and right of the partition respectively), not the 4 values per the $x$ and $y$ partitioning.
```{r res2, echo=FALSE}
NNS.part(x,y,order = 4, type = "XONLY")
```

## Clusters Used in Regression
The right column of plots shows the corresponding regression (plus endpoints and central point) for the order of `NNS` partitioning.
```{r depreg,results='hide'}
for(i in 1 : 3){NNS.part(x, y, order = i, obs.req = 0, Voronoi = TRUE, type = "XONLY") ; NNS.reg(x, y, order = i, ncores = 1)}
```


# NNS Regression `NNS.reg()`
**`NNS.reg`** estimates nonlinear relationships in univariate and multivariate settings using partition summaries and interpolation or distance weighting. Its fit depends on the observed support, selected resolution, and reduction rule. The examples below show its main outputs. Set `type = "REG"` when a numeric response should be treated as regression even if its values could trigger automatic classification.

## Univariate:
```{r nonlinear,fig.width=5,fig.height=3,fig.align = "center"}
NNS.reg(x, y, ncores = 1)
```

## Multivariate:
Multivariate regressions return a plot of $y$ and $\hat{y}$, as well as the regression points (`$RPM`) and partitions (`$rhs.partitions`) for each regressor.
```{r nonlinear multi,fig.width=5,fig.height=3,fig.align = "center"}
f = function(x, y) x ^ 3 + 3 * y - y ^ 3 - 3 * x
y = x ; z <- expand.grid(x, y)
g = f(z[ , 1], z[ , 2])
NNS.reg(z, g, order = "max", plot = FALSE, ncores = 1)
```

## Inter/Extrapolation
`NNS.reg` can inter- or extrapolate any point of interest.  The **`NNS.reg(x, y, point.est = ...)`** parameter permits any sized data of similar dimensions to $x$ and called specifically with **`NNS.reg(...)$Point.est`**.


## NNS Dimension Reduction Regression
`dim.red.method = "cor"` builds a signed synthetic predictor from the selected correlations. For a continuous response, `$equation` reports its coefficients and denominator. The response is fitted as a nonlinear function of that index; the index itself is not a linear prediction equation for the response.

```{r dimension_reduction,fig.width=5,fig.height=3,fig.align="center"}
regression_x <- iris[, 2:4]
regression_y <- iris$Sepal.Length
reduced <- NNS.reg(regression_x, regression_y, dim.red.method = "cor",
                   plot = FALSE, ncores = 1)
reduced$equation
```

With retained coefficients $c_j$, the index is
$X^\star = \sum_j c_j X_j / \sum_j \mathbf{1}\{|c_j| > 0\}$,
and the regression estimates a function $\hat y = g(X^\star)$.

### Threshold

`threshold` retains coefficients whose absolute association meets the cutoff. This controls which coordinates contribute to the synthetic index; it does not establish a significance level.

```{r dimension_reduction_threshold,fig.width=5,fig.height=3,fig.align="center"}
reduced_threshold <- NNS.reg(regression_x, regression_y,
                            dim.red.method = "cor", threshold = .75,
                            point.est = regression_x[1:10, ],
                            plot = FALSE, ncores = 1)
reduced_threshold$equation
reduced_threshold$Point.est
```

# Classification

Set `type = "CLASS"` for nominal classes. Factor, character, logical, and numeric `0/1` responses are classified automatically. Numeric class values need not start at one. Factor and character predictions retain numeric codes with the mapping returned in `$class.levels`.

```{r class,fig.width=5,fig.height=3,fig.align="center",message=FALSE}
classification <- NNS.reg(iris[, 1:4], iris$Species, type = "CLASS",
                          point.est = iris[1:10, 1:4], plot = FALSE, ncores = 1)
classification$class.levels[classification$Point.est]
classification$Point.probabilities
```

With dimension reduction, classification uses one signed index per class indicator. Its `$equation` is therefore a named list, with one equation per class. Class probabilities are combined before labels are chosen; the numeric labels are never interpolation coordinates.

```{r classification_indices}
class_reduced <- NNS.reg(iris[, 1:4], iris$Species, dim.red.method = "cor",
                         plot = FALSE, ncores = 1)
class_reduced$equation
```

These examples demonstrate fitted outputs at training observations. Independent classification accuracy requires heldout observations, as shown in the classification vignette.


# Cross-Validation `NNS.stack()`
The **`NNS.stack`** routine cross-validates for a given objective function the `n.best` parameter in the multivariate **`NNS.reg`** function as well as the `threshold` parameter in the dimension reduction **`NNS.reg`** version.  **`NNS.stack`** can be used for classification:

**`NNS.stack(..., type = "CLASS", ...)`** 

or continuous dependent variables:

**`NNS.stack(..., type = "REG", ...)`**.

Any objective function `obj.fn` can be called using `expression()` with the terms `predicted` and `actual`, even from external packages such as `Metrics`.

**`NNS.stack(..., obj.fn = expression(Metrics::mape(actual, predicted)), objective = "min")`**.


```{r stack,fig.width=5,fig.height=3,fig.align="center",message=FALSE}
set.seed(123)
stack_fit <- NNS.stack(IVs.train = iris[, 1:4],
                       DV.train = iris$Species,
                       IVs.test = iris[1:10, 1:4],
                       dim.red.method = "cor",
                       obj.fn = expression(mean(predicted == actual)),
                       objective = "max", type = "CLASS",
                       folds = 1, ncores = 1, status = FALSE)
stack_fit[c("OBJfn.reg", "OBJfn.dim.red", "NNS.reg.n.best",
            "NNS.dim.red.count", "weights", "class.levels")]
stack_fit$class.levels[stack_fit$stack]
```

The displayed objectives describe internal tuning folds. The final predictions above are at supplied training rows and are not an independent test score. For nominal models, the selected dimension count is reported directly; a single scalar coefficient threshold cannot summarize the class-specific equations.

# Increasing Dimensions
A duplicated predictor changes its weight in the multivariate geometry and can be used to compare that prediction path with univariate interpolation. It adds no information. Choose this representation and the neighbor count `n.best` using heldout performance:

**`NNS.stack(IVs.train = cbind(x, x), DV.train = y, method = 1, ...)`**.

```{r stack2, message = FALSE,fig.width=5,fig.height=3,fig.align = "center",results='hide'}
set.seed(123)
x = rnorm(100); y = rnorm(100)

nns.params = NNS.stack(IVs.train = cbind(x, x),
                        DV.train = y,
                        method = 1, ncores = 1)
```

```{r stack2res, fig.width=5,fig.height=3,fig.align = "center",results='hide'}
duplicated_fit <- NNS.reg(cbind(x, x), y,
                          n.best = nns.params$NNS.reg.n.best,
                          point.est = cbind(x, x), plot = FALSE,
                          residual.plot = FALSE, ncores = 1,
                          confidence.interval = .95)
plot(x, y, pch = 19, col = "steelblue", main = "Tuned multivariate prediction path")
o <- order(x)
lines(x[o], duplicated_fit$Point.est[o], col = "red", lwd = 2)
```



# Smoothing Option
Smoothness is not required for curve fitting, but the `NNS.reg` function offers an optional smoothed fit. This feature applies a smoothing spline to regression points generated internally using the partitioning method described earlier.

```{r smooth, fig.width=5,fig.height=3,fig.align = "center",results='hide'}
NNS.reg(x, y, smooth = TRUE)
```


# Imputation
Imputation in `NNS` is a direct application of nearest neighbor regression. When values of $y$ are missing, we use the observed $(X,y)$ pairs as the training set and the predictors of the missing rows as `point.est`.

A key insight is that even in univariate regressions, `NNS.reg` benefits from the increasing dimensions trick: by duplicating the predictor into a multivariate form, e.g. `cbind(x, x)`, the distance function underlying `NNS.reg` operates in a 2-D space. This sharpened distance metric allows a more robust donor selection, effectively turning univariate imputation into a special case of multivariate nearest neighbor regression.

For multivariate predictors, the same form applies directly — supply the full set of observed predictors in $x$, the observed responses in $y$, and the incomplete rows in `point.est`. With `order = "max", n.best = 1`, the imputation is always 1-NN donor-based: each missing $y$ is filled in by the response of its closest donor under the `NNS` hybrid distance. This ensures imputations remain strictly within the support of the observed data.

**Categorical data** is handled analogously, only requiring `NNS.reg(..., type = "CLASS")` in the procedure.

## Univariate Imputation

```{r uniimpute, fig.width=7, fig.height=5}
set.seed(123)

# Univariate predictor with nonlinear signal
n <- 400
x <- sort(runif(n, -3, 3))
y <- sin(x) + 0.2 * x^2 + rnorm(n, 0, 0.25)

# Induce ~25% MCAR missingness in y
miss <- rbinom(n, 1, 0.25) == 1
y_mis <- y
y_mis[miss] <- NA

# ---- Increasing dimensions trick ----
# Duplicate x so the distance operates in a 2D space: cbind(x, x).
# This chooses the multivariate neighbor path; duplication adds no information.
x2_train <- cbind(x[!miss], x[!miss])
x2_miss  <- cbind(x[miss],  x[miss])

# 1-NN donor imputation with NNS.reg
y_hat_uni <- NNS::NNS.reg(
  x         = x2_train,             # predictors (duplicated x)
  y         = y[!miss],             # observed responses
  point.est = x2_miss,              # rows to impute
  order     = "max",                # maximum requested partition order
  n.best    = 1,                    # 1-NN donor
  plot      = FALSE
)$Point.est

# Fill back
y_completed_uni <- y_mis
y_completed_uni[miss] <- y_hat_uni

# Plot observed vs imputed (NNS 1-NN)
plot(x, y, pch = 1, col = "steelblue", cex = 1.5, lwd = 2,
     xlab = "x", ylab = "y", main = "NNS 1-NN Imputation")
points(x[miss], y_hat_uni, col = "red", pch = 15, cex = 1.3)

legend("topleft",
       legend = c("Observed", "Imputed (NNS 1-NN)"),
       col    = c("steelblue", "red"),
       pch    = c(1, 15),
       pt.lwd = c(2, NA),
       bty    = "n")
```


## Multivariate Imputation
```{r multiimpute, fig.width=7, fig.height=5}
set.seed(123)

# Multivariate predictors with nonlinear & interaction structure
n <- 800
X <- cbind(
  x1 = rnorm(n),
  x2 = runif(n, -2, 2),
  x3 = rnorm(n, 0, 1)
)

f <- function(x1, x2, x3) 1.1*x1 - 0.8*x2 + 0.5*x3 + 0.6*x1*x2 - 0.4*x2*x3 + 0.3*sin(1.3*x1)
y <- f(X[,1], X[,2], X[,3]) + rnorm(n, 0, 0.4)

# Induce ~30% MCAR missingness in y
miss <- rbinom(n, 1, 0.30) == 1
y_mis <- y
y_mis[miss] <- NA

# Training (observed) vs rows to impute
X_obs <- X[!miss, , drop = FALSE]
y_obs <- y[!miss]
X_mis <- X[ miss, , drop = FALSE]

# 1-NN donor imputation with NNS.reg
y_hat_mv <- NNS::NNS.reg(
  x         = X_obs,     # all observed predictors
  y         = y_obs,     # observed responses
  point.est = X_mis,     # rows to impute
  order     = "max",     # maximum requested partition order
  n.best    = 1,         # 1-NN donor
  plot      = FALSE
)$Point.est

# Completed vector
y_completed_mv <- y_mis
y_completed_mv[miss] <- y_hat_mv

# Plot observed vs imputed (multivariate, NNS 1-NN)
plot(seq_along(y), y, 
     pch = 1, col = "steelblue", cex = 1.5, lwd = 2,
     xlab = "Observation index", ylab = "y",
     main = "NNS 1-NN Multivariate Imputation")

# Overlay imputed values
points(which(miss), y_hat_mv, pch = 15, col = "red", cex = 1.2)

# Legend
legend("topleft",
       legend = c("Observed", "Imputed (NNS 1-NN)"),
       col    = c("steelblue", "red"),
       pch    = c(1, 15),
       pt.lwd = c(2, NA),
       bty    = "n")
```


## A Note on Uncertainty Propagation

A completed data set does not by itself represent imputation uncertainty. Bootstrap model refits can produce multiple completed data sets, but valid pooled inference also depends on the missingness assumptions, estimand, imputation design, and sampling calibration. Applying Rubin's pooling formula alone does not establish those conditions. A smaller pooled standard error is not, by itself, evidence of more accurate uncertainty.

[NNS Multiple Imputation vs MICE](https://github.com/OVVO-Financial/NNS/blob/NNS-Beta-Version/examples/NNS_MI_vs_MICE.md) supplies an illustrative simulation. A general comparison requires repeated samples and assessment of bias, interval coverage, and the missingness mechanism.


# References
If the user is so motivated, detailed arguments further examples are provided within the following:

-  [Nonlinear Nonparametric Statistics: Using Partial Moments](https://ovvo-financial.github.io/NNS/book/)

-  [Deriving Nonlinear Correlation Coefficients from Partial Moments](https://doi.org/10.2139/ssrn.2148522)

-  [Nonparametric Regression Using Clusters](https://doi.org/10.1007/s10614-017-9713-5)

-  [Clustering and Curve Fitting by Line Segments](https://doi.org/10.2139/ssrn.2861339)

-  [Classification Using NNS Clustering Analysis](https://doi.org/10.2139/ssrn.2864711)

-  [Partitional Estimation Using Partial Moments](https://doi.org/10.2139/ssrn.3592491)