---
title: "Classical and Riemannian Linear Regression with UMAP, ISOMAP, and DBSCAN"
author: "Oldemar Rodríguez Rojas & Jennifer Lobo Vásquez"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Classical and Riemannian Linear Regression with UMAP, ISOMAP, and DBSCAN}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.align = "center"
)
```

# Overview

This example shows how to use the `riemannianStats` package to fit a
Riemannian linear regression model.

The objective is to compare a classical linear regression model with
Riemannian linear regression models. The Riemannian models use a similarity
matrix to incorporate the local geometric structure of the observations.

In this example, three Riemannian similarity methods are illustrated:

* UMAP similarities.
* ISOMAP similarities.
* DBSCAN similarities.

The `riem.lm()` function provides two interfaces for specifying the response
and explanatory variables:

* A standard R formula together with a data frame.
* A predictor data frame or matrix and a response vector passed separately.

# 1. Load the package

The first step is to load the `riemannianStats` package.

```{r load-package}
library(riemannianStats)
```

# 2. Load the data set

The data set contains 10 observations, identified as ind_1 through ind_10. Each observation has five numeric explanatory variables, named X1, X2, X3, X4, and X5, as well as a numeric response variable named y. The individual identifiers are stored as the row names of the data frame.

```{r create-data}
df <- data.frame(
  X1 = c(0.0, 0.5, 0.2, 0.9, 1.2, 1.5, 4.8, 5.5, 5.2, 6.0),
  X2 = c(0.0, 0.1, 0.8, 0.9, 0.3, 0.7, 5.0, 5.1, 5.8, 5.6),
  X3 = c(0.1, 0.2, 0.1, 0.2, 0.3, 0.3, 5.0, 5.1, 4.9, 5.2),
  X4 = c(0.2, 0.0, 0.3, 0.2, 0.1, 0.3, 4.7, 4.8, 5.3, 5.4),
  X5 = c(0.0, 0.1, 0.0, 0.1, 0.2, 0.1, 4.9, 5.0, 5.1, 5.3),
  y = c(
    8.77, 9.09, 8.35, 25.00, -10.00,
    10.00, 12.47, 12.97, 12.31, 13.33
  ),
  row.names = c(
    "ind_1", "ind_2", "ind_3", "ind_4", "ind_5",
    "ind_6", "ind_7", "ind_8", "ind_9", "ind_10"
  )
)

df
```

The number of neighbors used by the Riemannian similarity methods is defined
before fitting the models.

```{r model-parameters}
n.neighbors <- 3
```

# 3. Fit a classical linear regression model

As a reference, a classical linear regression model is fitted using the
`lm()` function. This model assumes the usual Euclidean structure of the
data.

```{r lm-classical}
model.lm <- lm(y ~ ., data = df)

summary(model.lm)
```

The expression `y ~ .` indicates that `y` is the response variable and that
all the remaining columns in `df` should be used as explanatory variables.

The fitted values and the classical coefficient of determination are
extracted from the model.

```{r lm-classical-results}
y.hat <- model.lm$fitted.values

head(y.hat)

summary(model.lm)$r.squared
```

# 4. Fit a Riemannian linear regression model using a formula

The `riem.lm()` function supports the standard R formula interface. In a
formula, the response variable is written on the left-hand side of `~`, while
the explanatory variables are written on the right-hand side.

For example, the following formula uses `y` as the response variable and all
the remaining columns in `df` as explanatory variables:

```{r lm-riem-umap-formula}
fit.umap.formula <- riem.lm(
  y ~ .,
  data = df,
  n.neighbors = n.neighbors,
  similarities = "umap",
  min.dist = 0.1,
  metric = "euclidean"
)

fit.umap.formula
```

## Formula specification options

There are several ways to specify the explanatory variables in a formula.

### Use one explanatory variable

The following formula fits a simple Riemannian linear regression model using
only `X1` as an explanatory variable:

```r
riem.lm(
  y ~ X1,
  data = df,
  similarities = "umap"
)
```

### List the explanatory variables explicitly

The variables can be listed individually by separating them with `+`. This
option makes it clear which variables are included in the model.

```r
riem.lm(
  y ~ X1 + X2 + X3 + X4 + X5,
  data = df,
  similarities = "umap"
)
```

### Use all the remaining variables

The dot `.` represents all the columns in `data` except the response
variable. This is useful when all the remaining variables should be included
as predictors.

```r
riem.lm(
  y ~ .,
  data = df,
  similarities = "umap"
)
```

### Exclude a variable

A variable can be removed from the formula using the minus sign `-`. In the
following example, all the explanatory variables except `X5` are used:

```r
riem.lm(
  y ~ . - X5,
  data = df,
  similarities = "umap"
)
```

### Include an interaction

The operator `*` includes the main effects of both variables and their
interaction. Therefore, `X1 * X2` is equivalent to
`X1 + X2 + X1:X2`.

```r
riem.lm(
  y ~ X1 * X2,
  data = df,
  similarities = "umap"
)
```

If only the interaction is required, the `:` operator can be used:

```r
riem.lm(
  y ~ X1:X2,
  data = df,
  similarities = "umap"
)
```

### Include a mathematical transformation

Mathematical transformations can be included directly in the formula. For
example, `I()` can be used to include a quadratic term:

```r
riem.lm(
  y ~ X1 + I(X1^2),
  data = df,
  similarities = "umap"
)
```

Other functions, such as `log()` or `sqrt()`, can also be used when they are
appropriate for the data:

```r
riem.lm(
  y ~ X1 + log(X2 + 1),
  data = df,
  similarities = "umap"
)
```

### Include categorical variables

If an explanatory variable is a factor, the formula interface processes it
as a categorical variable. Internally, `model.matrix()` converts its levels
into indicator variables.

For example, if the data contained a factor named `group`, it could be
included as follows:

```r
riem.lm(
  y ~ X1 + group,
  data = df,
  similarities = "umap"
)
```

Consequently, the formula interface supports common R formula notation,
including explicit variable selection, `.`, exclusions, interactions,
transformations, and categorical variables.

# 5. Fit the model by passing the variables separately

Alternatively, the explanatory variables and the response variable can be
passed separately. In this interface, `data` contains the explanatory
variables and `y` contains the numeric response vector.

First, the explanatory variables and the response are extracted from the
data frame:

```{r prepare-riemannian-data}
data.analysis <- df[, c("X1", "X2", "X3", "X4", "X5")]

y <- df$y
```

The Riemannian regression model can then be fitted as follows:

```{r lm-riem-umap-separate}
fit.umap <- riem.lm(
  data = data.analysis,
  y = y,
  n.neighbors = n.neighbors,
  similarities = "umap",
  min.dist = 0.1,
  metric = "euclidean"
)

fit.umap
```

Therefore, the following two calls represent equivalent ways of specifying
the same response and explanatory variables:

```r
# Formula interface
riem.lm(
  y ~ X1 + X2 + X3 + X4 + X5,
  data = df,
  similarities = "umap"
)

# Separate data and response interface
riem.lm(
  data = df[, c("X1", "X2", "X3", "X4", "X5")],
  y = df$y,
  similarities = "umap"
)
```

The formula interface is convenient when the variables are stored in a
single data frame. The separate interface is useful when the predictor
matrix and response vector have already been prepared independently.

# 6. Inspect the Riemannian regression results with UMAP

The fitted model stores the transformed response, fitted values, residuals,
and Riemannian weights. These values are expressed on the Riemannian
transformed scale.

```{r lm-riem-umap-results}
riem.results.umap <- data.frame(
  y.R = fit.umap$y.R,
  yhat.R = fit.umap$yhat.R,
  residual.R = fit.umap$residuals.R,
  weight = fit.umap$weights
)

riem.results.umap
```

The Riemannian regression coefficients and the Riemannian coefficient of
determination can also be extracted directly from the fitted model.

```{r lm-riem-umap-coefficients}
fit.umap$beta.R

fit.umap$R2.R
```

# 7. Convert UMAP Riemannian results to the original scale

The `riem.original.scale()` function converts the fitted values from the
Riemannian transformed scale back to the original response scale.

```{r lm-riem-umap-original-scale}
orig.umap <- riem.original.scale(fit.umap)

orig.umap$table

orig.umap$yhat.original
```

# 8. Fit a Riemannian linear regression model using ISOMAP similarities

A second Riemannian regression model is fitted using ISOMAP similarities.
ISOMAP constructs a neighborhood graph and uses shortest-path distances to
approximate geodesic distances along the underlying data manifold.

The formula interface can be used with ISOMAP similarities as follows:

```{r lm-riem-isomap-formula}
fit.isomap <- riem.lm(
  y ~ .,
  data = df,
  n.neighbors = n.neighbors,
  similarities = "isomap",
  metric = "euclidean"
)

fit.isomap
```

The separate interface is also valid:

```r
fit.isomap <- riem.lm(
  data = data.analysis,
  y = y,
  n.neighbors = n.neighbors,
  similarities = "isomap",
  metric = "euclidean"
)
```

The `n.neighbors` argument controls the number of nearest neighbors used to
construct the neighborhood graph. The optional `max.neighbors` argument can be
used when the implementation needs to increase the neighborhood size to obtain
a connected graph.

# 9. Inspect the Riemannian regression results with ISOMAP

The transformed response, fitted values, residuals, and weights are extracted
from the ISOMAP-based Riemannian model.

```{r lm-riem-isomap-results}
riem.results.isomap <- data.frame(
  y.R = fit.isomap$y.R,
  yhat.R = fit.isomap$yhat.R,
  residual.R = fit.isomap$residuals.R,
  weight = fit.isomap$weights
)

riem.results.isomap
```

The coefficients and coefficient of determination can also be extracted
directly from the fitted model.

```{r lm-riem-isomap-coefficients}
fit.isomap$beta.R

fit.isomap$R2.R
```

# 10. Convert ISOMAP Riemannian results to the original scale

The fitted values from the ISOMAP-based model are converted back to the
original response scale.

```{r lm-riem-isomap-original-scale}
orig.isomap <- riem.original.scale(fit.isomap)

orig.isomap$table

orig.isomap$yhat.original
```

# 11. Fit a Riemannian linear regression model using DBSCAN similarities

A third Riemannian regression model is fitted using DBSCAN similarities.
This allows the model to incorporate a density-based local structure.

The formula interface can also be used with DBSCAN similarities:

```{r lm-riem-dbscan-formula}
fit.dbscan <- riem.lm(
  y ~ .,
  data = df,
  n.neighbors = n.neighbors,
  similarities = "dbscan",
  metric = "euclidean"
)

fit.dbscan
```

The separate interface is also valid:

```r
fit.dbscan <- riem.lm(
  data = data.analysis,
  y = y,
  n.neighbors = n.neighbors,
  similarities = "dbscan",
  metric = "euclidean"
)
```

The choice between the formula and separate interfaces changes only how the
response and predictors are supplied. It does not change the selected
similarity method.

# 12. Inspect the Riemannian regression results with DBSCAN

The transformed response, fitted values, residuals, and weights are extracted
from the DBSCAN-based Riemannian model.

```{r lm-riem-dbscan-results}
riem.results.dbscan <- data.frame(
  y.R = fit.dbscan$y.R,
  yhat.R = fit.dbscan$yhat.R,
  residual.R = fit.dbscan$residuals.R,
  weight = fit.dbscan$weights
)

riem.results.dbscan
```

The Riemannian regression coefficients and the Riemannian coefficient of
determination are shown separately.

```{r lm-riem-dbscan-coefficients}
fit.dbscan$beta.R

fit.dbscan$R2.R
```

# 13. Convert DBSCAN Riemannian results to the original scale

The fitted values from the DBSCAN-based Riemannian model are converted back
to the original response scale.

```{r lm-riem-dbscan-original-scale}
orig.dbscan <- riem.original.scale(fit.dbscan)

orig.dbscan$table

orig.dbscan$yhat.original
```

# 14. Compare the models

The classical model and the Riemannian models can be compared using their
respective coefficients of determination.

```{r compare-models}
comparison <- data.frame(
  Model = c(
    "Classical LM",
    "Riemannian LM - UMAP",
    "Riemannian LM - ISOMAP",
    "Riemannian LM - DBSCAN"
  ),
  R2 = c(
    summary(model.lm)$r.squared,
    fit.umap$R2.R,
    fit.isomap$R2.R,
    fit.dbscan$R2.R
  )
)

comparison
```

The classical and Riemannian coefficients of determination are calculated on
different representations of the response. The classical value is calculated
on the original scale, whereas `R2.R` is calculated on the Riemannian
transformed scale. Therefore, the comparison should be interpreted with this
distinction in mind.

# 15. Summary

This example illustrates how to fit a classical linear regression model and
Riemannian linear regression models using the `riemannianStats` package.

The `riem.lm()` function supports two interfaces. The first uses a standard R
formula together with a data frame:

```r
riem.lm(
  y ~ X1 + X2 + X3,
  data = df
)
```

The second passes the explanatory variables and response variable separately:

```r
riem.lm(
  data = X,
  y = y
)
```

Several formula specifications are available:

* `y ~ X1` fits a model with one explanatory variable.
* `y ~ X1 + X2` lists the explanatory variables explicitly.
* `y ~ .` uses all the remaining variables in the data frame.
* `y ~ . - X5` uses all the remaining variables except `X5`.
* `y ~ X1 * X2` includes main effects and an interaction.
* `y ~ X1:X2` includes only the interaction.
* `y ~ X1 + I(X1^2)` includes a quadratic term.
* `y ~ X1 + group` can include a categorical variable.

The formula interface is useful when the response and explanatory variables
are stored in the same data frame. The separate interface is useful when the
predictor matrix and response vector have already been prepared.

The Riemannian models differ in how the local geometric structure is
constructed. The UMAP model uses a fuzzy neighborhood-based similarity
structure, the ISOMAP model approximates geodesic distances through shortest
paths in a neighborhood graph, and the DBSCAN model uses a density-based
similarity structure.

The fitted Riemannian models store their results on the transformed
Riemannian scale. The `riem.original.scale()` function can then be used to
convert the fitted values back to the original response scale.