# Step-by-Step Riemannian PCA with UMAP on Synthetic Data

#### Oldemar Rodríguez Rojas & Jennifer Lobo Vásquez

# Overview

This example explores a data set composed of 250 observations and 10 numeric variables. It also includes a column named `cluster`.

The numeric variables are used to perform the Riemannian analysis, while the `cluster` column is used only for visualization. In particular, this column allows the individuals to be colored according to the group they belong to.

In this data set, the first two variables, `x` and `y`, contain the most relevant structural information. The remaining variables correspond to added noise, included to increase the dimensionality of the problem.

The goal of this example is to explore the general structure of the data using Riemannian principal components and, using the known clusters, visualize possible patterns or characteristics among the individuals.

    library(riemannianStats)

# 1. Load the data set

    data.path
    #> [1] "/tmp/RtmpSeNnfL/Rinst1b08444218c7/riemannianStats/extdata/Data10D_250.csv"

    original.data<- read.csv(
      data.path, # It must be replaced with the path to the CSV file.
      sep = ",",
      dec = "."
    )

    original.data$cluster<- as.factor(original.data$cluster)

    str(original.data)
    #> 'data.frame':    250 obs. of  11 variables:
    #>  $ x      : num  2.06 1.57 4.8 -4.62 -3.68 ...
    #>  $ y      : num  -4.79 4.57 -1.37 0.56 -3 ...
    #>  $ var1   : num  -0.7546 -1.4854 1.7311 0.0109 1.6029 ...
    #>  $ var2   : num  -0.0093 -1.5914 0.9597 0.6607 1.3018 ...
    #>  $ var3   : num  -0.605 -0.158 -0.545 0.565 -1.083 ...
    #>  $ var4   : num  -1.1174 -0.9845 0.412 -1.7042 0.0538 ...
    #>  $ var5   : num  -0.196 -0.165 0.653 -1.71 1.863 ...
    #>  $ var6   : num  -0.563 1.128 0.225 0.307 1.202 ...
    #>  $ var7   : num  -0.5868 -0.0476 0.0693 0.2137 -0.6494 ...
    #>  $ var8   : num  -1.0984 1.3314 0.0857 0.3045 -0.7006 ...
    #>  $ cluster: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...

# 2. Prepare the data

Since the Riemannian analysis is applied only to numeric variables, the `cluster` column is first separated from the data set used in the analysis. This column is kept for later visualizations, for example, to color individuals according to the group they belong to.

Additionally, since this method is based on calculating distances between observations, it is recommended to scale the numerical variables before performing the analysis. This allows all variables to contribute in a comparable way and prevents those with larger magnitudes from dominating the results.

    clusters <- original.data$cluster

    data.analysis <- original.data[, setdiff(names(original.data), "cluster"), drop = FALSE]

    data.analysis.scaled <- scale(data.analysis)

    data.analysis.scaled <- as.data.frame(data.analysis.scaled)

    head(data.analysis.scaled)
    #>            x          y        var1        var2       var3        var4
    #> 1 -0.2074543 0.01870356 -0.67566177 -0.06441618 -0.6576387 -1.23307459
    #> 2 -0.2927010 1.51629287 -1.37317174 -1.71329539 -0.1935010 -1.09528575
    #> 3  0.2671849 0.56521880  1.69668460  0.94548345 -0.5950688  0.35293661
    #> 4 -1.3650547 0.87449934  0.05487369  0.63385730  0.5575895 -1.84154703
    #> 5 -1.2027221 0.30493810  1.57433004  1.30202055 -1.1536887 -0.01852609
    #> 6  0.2386114 0.40461734 -0.34518373 -0.17405326 -2.5224754 -1.44460491
    #>         var5       var6         var7       var8
    #> 1 -0.1809308 -0.4277941 -0.494942852 -1.0760166
    #> 2 -0.1505560  1.2854218 -0.004294505  1.2195707
    #> 3  0.6385524  0.3699252  0.102121416  0.0426612
    #> 4 -1.6416776  0.4533572  0.233465191  0.2494652
    #> 5  1.8060741  1.3601701 -0.551946857 -0.7001162
    #> 6 -2.1897655 -0.3961836  0.004922900  1.2618115

# 3. Choose the number of neighbors

The `n.neighbors` parameter determines how many observations are considered in the local neighborhood of each point when constructing the UMAP similarity graph. As a practical rule, when a certain number of groups is known or expected in the data, it can be defined as the total number of observations divided by the expected number of groups.

    expected.groups <- 5

    n.neighbors <- as.integer(nrow(data.analysis) / expected.groups)
    n.neighbors
    #> [1] 50

# 4. Compute the UMAP similarities

This matrix represents the local similarity between observations.

    umap.similarities <- riem.similarities.umap(
      data = data.analysis,
      n.neighbors = n.neighbors,
      min.dist = 0.1,
      metric = "euclidean"
    )
    umap.similarities[1:5, 1:5]
    #>           [,1]        [,2]      [,3]        [,4]      [,5]
    #> [1,] 0.0000000 0.000000000 0.1027273 0.000000000 0.0390701
    #> [2,] 0.0000000 0.000000000 0.0000000 0.006585082 0.0000000
    #> [3,] 0.1027273 0.000000000 0.0000000 0.000000000 0.0000000
    #> [4,] 0.0000000 0.006585082 0.0000000 0.000000000 0.1910163
    #> [5,] 0.0390701 0.000000000 0.0000000 0.191016320 0.0000000

# 5. Compute the Rho matrix

The Rho matrix is computed as `1 - umap.similarities`. This matrix transforms local similarity into a local separation measure.

    rho <- riem.rho(umap.similarities)

    rho[1:5, 1:5]
    #>           [,1]      [,2]      [,3]      [,4]      [,5]
    #> [1,] 1.0000000 1.0000000 0.8972727 1.0000000 0.9609299
    #> [2,] 1.0000000 1.0000000 1.0000000 0.9934149 1.0000000
    #> [3,] 0.8972727 1.0000000 1.0000000 1.0000000 1.0000000
    #> [4,] 1.0000000 0.9934149 1.0000000 1.0000000 0.8089837
    #> [5,] 0.9609299 1.0000000 1.0000000 0.8089837 1.0000000

# 6. Compute the Riemannian vector differences

The Riemannian vector differences are stored in a three-dimensional array. The entry `riemannian.diff[i, j, ]` contains the weighted difference between observations `i` and `j`.

    riemannian.diff <- riem.diff(
      data = data.analysis,
      rho = rho
    )

    riemannian.diff[1, 2, ]
    #>           x           y        var1        var2        var3        var4
    #>  0.49167881 -9.35504083  0.73083244  1.58205801 -0.44714571 -0.13286946
    #>        var5        var6        var7        var8
    #> -0.03147282 -1.69081857 -0.53916487 -2.42982678

# 7. Compute the UMAP-based Riemannian distance matrix

The distance matrix summarizes each Riemannian difference vector into a single distance value.

    umap.distance.matrix <- riem.dist(riemannian.diff)

    umap.distance.matrix[1:5, 1:5]
    #>           [,1]      [,2]      [,3]      [,4]      [,5]
    #> [1,]  0.000000 10.003388  5.058490  9.031236  6.965314
    #> [2,] 10.003388  0.000000  8.229759  8.123170 10.642235
    #> [3,]  5.058490  8.229759  0.000000 10.333392  8.867062
    #> [4,]  9.031236  8.123170 10.333392  0.000000  4.961074
    #> [5,]  6.965314 10.642235  8.867062  4.961074  0.000000

# 8. Riemannian correlation matrix

The Riemannian correlation matrix is computed from the Riemannian covariance matrix.

    correlation.matrix <- riem.cor(
      data = data.analysis,
      rho = rho,
      umap.distance.matrix = umap.distance.matrix
    )

    correlation.matrix[1:5, 1:5]
    #>                x           y        var1        var2        var3
    #> x     1.00000000  0.45273855  0.13331328 -0.04184126  0.09649943
    #> y     0.45273855  1.00000000 -0.10806913 -0.02340368 -0.08123159
    #> var1  0.13331328 -0.10806913  1.00000000  0.07218516  0.31681632
    #> var2 -0.04184126 -0.02340368  0.07218516  1.00000000  0.02365621
    #> var3  0.09649943 -0.08123159  0.31681632  0.02365621  1.00000000

# 10. Riemannian covariance matrix

The Riemannian covariance matrix is computed using the data centered with respect to the Riemannian mean observation.

    covariance.matrix <- riem.cov(
      data = data.analysis,
      rho = rho,
      umap.distance.matrix = umap.distance.matrix
    )

    covariance.matrix[1:5, 1:5]
    #>               x          y        var1        var2        var3
    #> x    34.2801633 16.4339501  0.96098014 -0.23074306  0.63621468
    #> y    16.4339501 38.4367006 -0.82488661 -0.13666573 -0.56709460
    #> var1  0.9609801 -0.8248866  1.51579026  0.08370858  0.43922250
    #> var2 -0.2307431 -0.1366657  0.08370858  0.88716620  0.02509025
    #> var3  0.6362147 -0.5670946  0.43922250  0.02509025  1.26798750

# 10. Riemannian principal components

The Riemannian principal components are computed from the original data, the Riemannian correlation matrix, the Rho matrix, and the UMAP-based Riemannian distance matrix.

    components <- riem.ind.coord(
      data = data.analysis,
      correlation.matrix = correlation.matrix,
      rho = rho,
      umap.distance.matrix = umap.distance.matrix
    )

    components[1:5, 1:5]
    #>      Component.1 Component.2 Component.3 Component.4 Component.5
    #> [1,]   0.0000000   0.0000000   0.0000000  0.00000000   0.0000000
    #> [2,]  -1.1200072  -1.1488369   1.3831043 -0.57534001  -0.9438186
    #> [3,]  -2.0943296  -0.3782564  -1.0144378 -0.06380339  -0.5234230
    #> [4,]  -0.9573561   0.5963693   0.2614249  1.66878216  -0.7702065
    #> [5,]  -1.7154839   0.5926663  -2.1134278 -1.19228986  -1.2825122
    dim(components)
    #> [1] 250  10

Each row corresponds to an observation, and each column corresponds to a Riemannian principal component.

# 11. Explained inertia

Explained inertia measures the proportion of the total Riemannian variance represented by the selected components.

In R, indices start at 1. Therefore, the first two components are selected using `component1 = 1` and `component2 = 2`.

    inertia <- riem.inertia(
      correlation.matrix = correlation.matrix,
      component1 = 1,
      component2 = 2
    ) * 100

    inertia
    #> [1] 43.61026

The result is multiplied by 100 to express it as a percentage.

# 12. Correlations between variables and components

The `riem.loadings()` function computes the Riemannian correlation between each original variable and the first two Riemannian principal components.

# Estan mal los resultado******************* Dan mal los negativos

    correlations <- riem.var.coord(
      data = data.analysis,
      components = components,
      rho = rho,
      umap.distance.matrix = umap.distance.matrix
    )

    correlations
    #>      Component.1  Component.2
    #> x    -0.26606396 -0.814317121
    #> y     0.01744553 -0.854836997
    #> var1 -0.67413322  0.091342678
    #> var2 -0.12043951  0.153521772
    #> var3 -0.67769815  0.115575190
    #> var4 -0.82429868 -0.023679400
    #> var5 -0.28379687 -0.088618277
    #> var6 -0.43762930  0.031409815
    #> var7 -0.60903855  0.116612496
    #> var8 -0.75561732  0.002370848

These correlations help interpret the axes of the Riemannian principal plane. Variables with large absolute values are strongly associated with the corresponding component.

# 13. Visualizations

## 13.1 Principal plane

The principal plane is constructed from the Riemannian component matrix and represents the observations in the space generated by the first two principal components.

In this case, since the data set includes a column with cluster labels, each individual can be visualized on the principal plane and colored according to the cluster it belongs to. This makes it easier to visually explore the structure of the data and identify possible patterns, groupings, or separations between the different groups.

    riem.plot(
      data = data.analysis,
      choix = "ind",
      components = components,
      clusters = clusters,
      explained.inertia = inertia,
      show.labels = TRUE
    )

    riem.plot(
      data = data.analysis,
      choix = "ind",
      components = components,
      clusters = clusters,
      explained.inertia = inertia,
      show.labels = TRUE,
      title = "Data10D_250",
      interactive = TRUE
    )

## 13.2 Correlation circle

The correlation circle shows how the original variables are related to the first two Riemannian components.

    riem.plot(
      data = data.analysis,
      choix = "var",
      correlations = correlations,
      explained.inertia = inertia,
      title = "Data10D_250"
    )

Variables pointing in similar directions are positively related in the component space. Variables pointing in opposite directions have a negative association. Longer arrows indicate a stronger correlation with the plotted components.

## 13.3 Biplot

The biplot combines, in a single visualization, the representation of individuals on the principal plane and the information associated with the variables through arrows.

This type of plot makes it possible to jointly analyze the position of the observations and the contribution of the variables to the principal components. In other words, it helps interpret which variables are most related to the separation, grouping, or dispersion of individuals in the plane.

    riem.biplot(
      data = data.analysis,
      components = components,
      correlations = correlations,
      clusters = clusters,
      explained.inertia = inertia,
      title = "Data10D_250"
    )

    riem.biplot(
      data = data.analysis,
      components = components,
      correlations = correlations,
      clusters = clusters,
      show.ind.labels = FALSE,
      show.var.labels = TRUE,
      var.color = "red",
      interactive = TRUE
    )

## 13.4 3D plot

Before generating the three-dimensional plot, a new `data.frame` is created by combining the original variables, the column with cluster labels, and the previously computed Riemannian components.

In particular, the first three Riemannian components are added. These components allow each observation to be represented in a three-dimensional space.

The goal of this visualization is to explore the dispersion of individuals in the space defined by the Riemannian components. In addition, by coloring the points according to the `cluster` variable, it is possible to identify whether the known groups show differentiated patterns, separations, or structures in three dimensions.

    data.viz <- original.data

    data.viz$Riemannian.Component1 <- components[, 1]
    data.viz$Riemannian.Component2 <- components[, 2]
    data.viz$Riemannian.Component3 <- components[, 3]

    head(data.viz)
    #>           x          y        var1         var2       var3        var4
    #> 1  2.058716 -4.7859969 -0.75455676 -0.009300336 -0.6053267 -1.11741505
    #> 2  1.567038  4.5690440 -1.48538920 -1.591358346 -0.1581809 -0.98454559
    #> 3  4.796300 -1.3720619  1.73112481  0.959672891 -0.5450474  0.41197187
    #> 4 -4.617992  0.5599311  0.01087891  0.660675462  0.5654121 -1.70416361
    #> 5 -3.681704 -2.9979661  1.60292495  1.301761249 -1.0832163  0.05377125
    #> 6  4.631496 -2.3752964 -0.40829064 -0.114494351 -2.4018919 -1.32139325
    #>         var5       var6        var7        var8 cluster Riemannian.Component1
    #> 1 -0.1961988 -0.5627749 -0.58675363 -1.09844412       1            0.00000000
    #> 2 -0.1647260  1.1280436 -0.04758876  1.33138266       1           -1.12000718
    #> 3  0.6529101  0.2245156  0.06934983  0.08565077       1           -2.09432960
    #> 4 -1.7097540  0.3068569  0.21368121  0.30454811       1           -0.95735615
    #> 5  1.8626395  1.2018148 -0.64939433 -0.70056215       1           -1.71548394
    #> 6 -2.2776562 -0.5315776 -0.03745991  1.37609350       1           -0.02195544
    #>   Riemannian.Component2 Riemannian.Component3
    #> 1             0.0000000             0.0000000
    #> 2            -1.1488369             1.3831043
    #> 3            -0.3782564            -1.0144378
    #> 4             0.5963693             0.2614249
    #> 5             0.5926663            -2.1134278
    #> 6            -0.3783834             0.7487315

The first three Riemannian components are visualized using a three-dimensional plot. This representation makes it possible to analyze the spatial distribution of the observations and facilitates the interpretation of possible groupings when a cluster variable is used to color the points.

    riem.plot3d(
      data = data.viz,
      x.col = "Riemannian.Component1",
      y.col = "Riemannian.Component2",
      z.col = "Riemannian.Component3",
      cluster.col = "cluster",
      title = "Data10D_250 - Riemannian Components",
      explained.inertia = inertia
    )

This example shows how to use `riemannianStats` with a higher-dimensional synthetic data set through a functional workflow. Starting from numeric data, it builds a UMAP-based Riemannian structure, extracts components, interprets variable contributions, and visualizes the resulting geometry.