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.

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.

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
#>         X1  X2  X3  X4  X5      y
#> ind_1  0.0 0.0 0.1 0.2 0.0   8.77
#> ind_2  0.5 0.1 0.2 0.0 0.1   9.09
#> ind_3  0.2 0.8 0.1 0.3 0.0   8.35
#> ind_4  0.9 0.9 0.2 0.2 0.1  25.00
#> ind_5  1.2 0.3 0.3 0.1 0.2 -10.00
#> ind_6  1.5 0.7 0.3 0.3 0.1  10.00
#> ind_7  4.8 5.0 5.0 4.7 4.9  12.47
#> ind_8  5.5 5.1 5.1 4.8 5.0  12.97
#> ind_9  5.2 5.8 4.9 5.3 5.1  12.31
#> ind_10 6.0 5.6 5.2 5.4 5.3  13.33

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

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.

model.lm <- lm(y ~ ., data = df)

summary(model.lm)
#> 
#> Call:
#> lm(formula = y ~ ., data = df)
#> 
#> Residuals:
#>    ind_1    ind_2    ind_3    ind_4    ind_5    ind_6    ind_7    ind_8 
#>   3.7388   5.8996  -8.0998  11.2753 -12.0956  -0.6438  -2.1699   0.9627 
#>    ind_9   ind_10 
#>  -2.9092   4.0419 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)   2.1986     9.5307   0.231    0.829
#> X1           -4.6526     8.0976  -0.575    0.596
#> X2           15.3786    14.9561   1.028    0.362
#> X3           27.4011    47.5770   0.576    0.596
#> X4            0.4626    32.3015   0.014    0.989
#> X5          -36.9997    63.5189  -0.582    0.591
#> 
#> Residual standard error: 10.23 on 4 degrees of freedom
#> Multiple R-squared:  0.3665, Adjusted R-squared:  -0.4254 
#> F-statistic: 0.4628 on 5 and 4 DF,  p-value: 0.7899

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.

y.hat <- model.lm$fitted.values

head(y.hat)
#>     ind_1     ind_2     ind_3     ind_4     ind_5     ind_6 
#>  5.031217  3.190372 16.449794 13.724683  2.095639 10.643766

summary(model.lm)$r.squared
#> [1] 0.3664888

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:

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

fit.umap.formula
#> Riemannian Linear Regression
#> -----------------------------
#> Similarity method: umap 
#> Riemannian center index lambda: 4 
#> Riemannian center y.lambda: 25 
#> Riemannian R-squared: 0.983437 
#> Solver: solve 
#> 
#> Riemannian coefficients:
#>         X1         X2         X3         X4         X5 
#>   3.707026  16.968374  43.124914   7.817295 -70.716339

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:

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.

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.

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:

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.

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

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

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:

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:

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:

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:

data.analysis <- df[, c("X1", "X2", "X3", "X4", "X5")]

y <- df$y

The Riemannian regression model can then be fitted as follows:

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

fit.umap
#> Riemannian Linear Regression
#> -----------------------------
#> Similarity method: umap 
#> Riemannian center index lambda: 4 
#> Riemannian center y.lambda: 25 
#> Riemannian R-squared: 0.983437 
#> Solver: solve 
#> 
#> Riemannian coefficients:
#>         X1         X2         X3         X4         X5 
#>   3.707026  16.968374  43.124914   7.817295 -70.716339

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

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

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
#>           y.R     yhat.R residual.R    weight
#> 1  -16.230000 -15.848717 -0.3812827 1.0000000
#> 2  -15.910000 -16.620968  0.7109684 1.0000000
#> 3    0.000000   0.000000  0.0000000 0.0000000
#> 4    0.000000   0.000000  0.0000000 1.0000000
#> 5   -6.028931  -2.172101 -3.8568297 0.1722552
#> 6    0.000000   0.000000  0.0000000 0.0000000
#> 7  -12.530000 -13.233280  0.7032798 1.0000000
#> 8  -12.030000 -10.918937 -1.1110628 1.0000000
#> 9  -12.690000 -11.941153 -0.7488472 1.0000000
#> 10 -11.670000 -12.793271  1.1232708 1.0000000

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

fit.umap$beta.R
#>         X1         X2         X3         X4         X5 
#>   3.707026  16.968374  43.124914   7.817295 -70.716339

fit.umap$R2.R
#> [1] 0.9834373

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.

orig.umap <- riem.original.scale(fit.umap)

orig.umap$table
#>    y.original yhat.original residual.original y.c.lambda yhat.c.lambda
#> 1        8.77      9.151283        -0.3812827     -16.23   -15.8487173
#> 2        9.09      8.379032         0.7109684     -15.91   -16.6209684
#> 3        8.35     24.249116       -15.8991164     -16.65    -0.7508836
#> 4       25.00     25.000000         0.0000000       0.00     0.0000000
#> 5      -10.00     12.390212       -22.3902115     -35.00   -12.6097885
#> 6       10.00     28.924762       -18.9247617     -15.00     3.9247617
#> 7       12.47     11.766720         0.7032798     -12.53   -13.2332798
#> 8       12.97     14.081063        -1.1110628     -12.03   -10.9189372
#> 9       12.31     13.058847        -0.7488472     -12.69   -11.9411528
#> 10      13.33     12.206729         1.1232708     -11.67   -12.7932708
#>       weight
#> 1  1.0000000
#> 2  1.0000000
#> 3  0.0000000
#> 4  1.0000000
#> 5  0.1722552
#> 6  0.0000000
#> 7  1.0000000
#> 8  1.0000000
#> 9  1.0000000
#> 10 1.0000000

orig.umap$yhat.original
#>  [1]  9.151283  8.379032 24.249116 25.000000 12.390212 28.924762 11.766720
#>  [8] 14.081063 13.058847 12.206729

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:

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

fit.isomap
#> Riemannian Linear Regression
#> -----------------------------
#> Similarity method: isomap 
#> Riemannian center index lambda: 6 
#> Riemannian center y.lambda: 10 
#> Riemannian R-squared: 0.940558 
#> Solver: solve 
#> 
#> Riemannian coefficients:
#>         X1         X2         X3         X4         X5 
#>   0.608941  -0.112204   6.704858   3.996546 -10.029041

The separate interface is also valid:

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.

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
#>            y.R      yhat.R   residual.R     weight
#> 1  -0.17062230 -0.21814548  0.047523177 0.13871732
#> 2  -0.08597998 -0.22780609  0.141826114 0.09448349
#> 3  -0.17067994 -0.11801853 -0.052661412 0.10344239
#> 4   0.76014321 -0.07388316  0.834026367 0.05067621
#> 5  -0.85658435 -0.08308927 -0.773495082 0.04282922
#> 6   0.00000000  0.00000000  0.000000000 1.00000000
#> 7   1.87128474  1.88285078 -0.011566046 0.75760516
#> 8   2.35307946  2.35113116  0.001948301 0.79228265
#> 9   1.90186169  1.94361461 -0.041752916 0.82331675
#> 10  2.86807133  2.82122639  0.046844934 0.86128268

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

fit.isomap$beta.R
#>          X1          X2          X3          X4          X5 
#>   0.6089405  -0.1122042   6.7048582   3.9965459 -10.0290406

fit.isomap$R2.R
#> [1] 0.9405576

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.

orig.isomap <- riem.original.scale(fit.isomap)

orig.isomap$table
#>    y.original yhat.original residual.original y.c.lambda yhat.c.lambda
#> 1        8.77      8.427410       0.342590074      -1.23     -1.572590
#> 2        9.09      7.588932       1.501067636      -0.91     -2.411068
#> 3        8.35      8.859089      -0.509089295      -1.65     -1.140911
#> 4       25.00      8.542054      16.457945570      15.00     -1.457946
#> 5      -10.00      8.059986     -18.059986262     -20.00     -1.940014
#> 6       10.00     10.000000       0.000000000       0.00      0.000000
#> 7       12.47     12.485267      -0.015266588       2.47      2.485267
#> 8       12.97     12.967541       0.002459098       2.97      2.967541
#> 9       12.31     12.360713      -0.050713066       2.31      2.360713
#> 10      13.33     13.275610       0.054389732       3.33      3.275610
#>        weight
#> 1  0.13871732
#> 2  0.09448349
#> 3  0.10344239
#> 4  0.05067621
#> 5  0.04282922
#> 6  1.00000000
#> 7  0.75760516
#> 8  0.79228265
#> 9  0.82331675
#> 10 0.86128268

orig.isomap$yhat.original
#>  [1]  8.427410  7.588932  8.859089  8.542054  8.059986 10.000000 12.485267
#>  [8] 12.967541 12.360713 13.275610

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:

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

fit.dbscan
#> Riemannian Linear Regression
#> -----------------------------
#> Similarity method: dbscan 
#> Riemannian center index lambda: 6 
#> Riemannian center y.lambda: 10 
#> Riemannian R-squared: 0.732843 
#> Solver: solve 
#> 
#> Riemannian coefficients:
#>         X1         X2         X3         X4         X5 
#>   0.020087   1.088430  13.134867   7.220052 -19.915526

The separate interface is also valid:

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.

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
#>           y.R     yhat.R residual.R    weight
#> 1  -0.3641445 -0.6363519  0.2722075 0.2960524
#> 2  -0.1745784 -0.7966619  0.6220835 0.1918443
#> 3  -0.3318597 -0.1111610 -0.2206987 0.2011271
#> 4   2.1731498 -0.2651037  2.4382535 0.1448767
#> 5  -2.7235950 -0.5279636 -2.1956314 0.1361797
#> 6   0.0000000  0.0000000  0.0000000 1.0000000
#> 7   2.4353468  2.6168785 -0.1815317 0.9859704
#> 8   2.9301543  2.7831115  0.1470427 0.9865839
#> 9   2.2802246  2.5352275 -0.2550029 0.9871102
#> 10  3.2891279  3.0086889  0.2804390 0.9877261

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

fit.dbscan$beta.R
#>           X1           X2           X3           X4           X5 
#>   0.02008663   1.08843009  13.13486727   7.22005203 -19.91552618

fit.dbscan$R2.R
#> [1] 0.7328433

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.

orig.dbscan <- riem.original.scale(fit.dbscan)

orig.dbscan$table
#>    y.original yhat.original residual.original y.c.lambda yhat.c.lambda
#> 1        8.77      7.850543         0.9194570      -1.23    -2.1494570
#> 2        9.09      5.847353         3.2426470      -0.91    -4.1526470
#> 3        8.35      9.447310        -1.0973096      -1.65    -0.5526904
#> 4       25.00      8.170142        16.8298579      15.00    -1.8298579
#> 5      -10.00      6.123039       -16.1230390     -20.00    -3.8769610
#> 6       10.00     10.000000         0.0000000       0.00     0.0000000
#> 7       12.47     12.654115        -0.1841147       2.47     2.6541147
#> 8       12.97     12.820958         0.1490423       2.97     2.8209577
#> 9       12.31     12.568333        -0.2583327       2.31     2.5683327
#> 10      13.33     13.046076         0.2839238       3.33     3.0460762
#>       weight
#> 1  0.2960524
#> 2  0.1918443
#> 3  0.2011271
#> 4  0.1448767
#> 5  0.1361797
#> 6  1.0000000
#> 7  0.9859704
#> 8  0.9865839
#> 9  0.9871102
#> 10 0.9877261

orig.dbscan$yhat.original
#>  [1]  7.850543  5.847353  9.447310  8.170142  6.123039 10.000000 12.654115
#>  [8] 12.820958 12.568333 13.046076

14. Compare the models

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

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
#>                    Model        R2
#> 1           Classical LM 0.3664888
#> 2   Riemannian LM - UMAP 0.9834373
#> 3 Riemannian LM - ISOMAP 0.9405576
#> 4 Riemannian LM - DBSCAN 0.7328433

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:

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

The second passes the explanatory variables and response variable separately:

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.

Classical and Riemannian Linear Regression with UMAP, ISOMAP, and DBSCAN