Example Workflow

“Talk is cheap. Show me the code.”

-Linus Torvalds

The Question

This vignette is the main workflow argument for PreciseDist. The question is not “which clustering algorithm should I use?” The question comes earlier:

What numerical representation of object relatedness deserves to be carried downstream?

Every unsupervised analysis begins by turning observations into relationships. That step can look harmless because it often happens behind a default argument: Euclidean distance, correlation distance, cosine similarity, a nearest-neighbor graph, an embedding. But once that representation exists, every downstream operation inherits it. A clustering algorithm does not recover the hidden truth from raw data; it acts on the structure it is given.

PreciseDist exists to make that representation explicit, broad, typed, and inspectable. We will begin with simple baselines, then widen the search, compare metrics, choose a set worth fusing, project the result to a graph, and visualize the structure. The point is not to automate judgment. The point is to give judgment something concrete to inspect.

Workflow Map

The full spine used in this vignette is:

Make a direct raw-data UMAP sanity check before committing to any metric.
Build typed distance/similarity matrices with precise_dist().
Convert only when the conversion is explicit with precise_transform().
Triage candidate matrices with precise_diagnostics().
Inspect candidate views with precise_viz() and, when useful, precise_trellis().
Stress-test promising embeddings across a compact UMAP setting grid.
Measure metric agreement with precise_correlations().
Fuse selected matrices with precise_fusion().
Project relatedness into graph matrices with precise_graph().
Inspect the result with precise_viz() and ask whether the chosen metric set is fragile with precise_stability().

This is an intentionally iterative workflow. We should expect to find metrics that are unhelpful, metrics that are redundant, and visualizations that make us change our mind.

1. Load the Data

The bundled data_cell_cycle dataset is useful here because it has known labels we can hold back until interpretation. Those labels are not used to build the distance matrices. They are a visual check after the fact.

For the website build, we need a subset large enough to show real visual structure but small enough to keep examples fast and reproducible. This first block will establish that subset, preserve labels separately, and make row names stable so every downstream matrix carries meaningful observation names.

set.seed(1)

label_col <- data_cell_cycle$Cell_cycle
rows_by_label <- split(seq_len(nrow(data_cell_cycle)), label_col)
workflow_rows <- unlist(lapply(rows_by_label, head, 20), use.names = FALSE)
workflow_features <- seq_len(500)

cell_cycle <- as.matrix(data_cell_cycle[workflow_rows, -1, drop = FALSE][, workflow_features])
storage.mode(cell_cycle) <- "double"

cell_cycle_labels <- label_col[workflow_rows]
rownames(cell_cycle) <- sprintf("%s_%02d", cell_cycle_labels, seq_along(cell_cycle_labels))
cell_cycle_groups <- setNames(cell_cycle_labels, rownames(cell_cycle))
cell_cycle_palette <- c(G1 = "#0017FF", G2M = "#000000", S = "#FF0017")

dim(cell_cycle)
#> [1]  60 500
table(cell_cycle_labels)
#> cell_cycle_labels
#>  G1 G2M   S 
#>  20  20  20

The labels are now set aside as cell_cycle_labels. They will be useful for checking what we see, but the distance calculations below use only the numeric matrix cell_cycle.

2. Set a Baseline

Before searching widely, we need a baseline to exceed. Ordinary choices like Manhattan distance and correlation distance are the right opening move. A default metric may be perfectly defensible, but it is not self-justifying.

This section will build a small baseline collection, coerce it into a common distance space when the registry says that conversion is valid, and inspect the declared type rather than guessing from matrix values.

baseline_native <- precise_dist(
  cell_cycle,
  dists = c("manhattan", "correlation"),
  verbose = FALSE
)

baseline_native[, c("distance", "type")]
#> # A tibble: 2 × 2
#>   distance    type       
#>   <chr>       <chr>      
#> 1 manhattan   distance   
#> 2 correlation correlation

baseline_distances <- precise_transform(baseline_native, to = "distance")
baseline_distances[, c("distance", "type")]
#> # A tibble: 2 × 2
#>   distance    type    
#>   <chr>       <chr>   
#> 1 manhattan   distance
#> 2 correlation distance

This small table is already the core PreciseDist contract. Manhattan is born as a distance. Correlation is born as a correlation. The conversion to distance is an explicit operation, recorded by the typed object, not a silent assumption made later by a plotting or clustering function.

3. Look Before Believing

The first visual check should be modest. We are not trying to prove anything with a plot. We are asking whether the baseline representations make the available structure easy to see.

The most common first move is to run UMAP directly on the numeric data. That is useful as a sanity check, but it is not the same thing as asking which relatedness metric should drive the analysis. It hides the representation inside the embedding algorithm. By default, UMAP uses Euclidean geometry on the input features, so a direct raw-data UMAP and a UMAP from an explicitly computed Euclidean distance matrix are expected to agree when the remaining UMAP settings are held fixed. We will look at direct UMAP first, then compare it to typed PreciseDist baselines that ask different questions.

set.seed(42)
raw_umap <- uwot::umap(
  cell_cycle,
  n_neighbors = 15,
  min_dist = 0.01,
  spread = 1,
  metric = "euclidean",
  n_threads = 1,
  verbose = FALSE
)

raw_umap_data <- data.frame(
  d1 = raw_umap[, 1],
  d2 = raw_umap[, 2],
  observation = rownames(cell_cycle),
  label = cell_cycle_labels
)

raw_umap_plot <- plotly::plot_ly(
  raw_umap_data,
  x = ~d1,
  y = ~d2,
  type = "scatter",
  mode = "markers",
  color = ~label,
  colors = cell_cycle_palette,
  text = ~observation,
  hoverinfo = "text",
  marker = list(size = 8, opacity = 0.75)
)

raw_umap_plot <- plotly::layout(
  raw_umap_plot,
  title = "Direct UMAP on raw data",
  xaxis = list(title = "", showticklabels = FALSE, showgrid = FALSE, zeroline = FALSE),
  yaxis = list(title = "", showticklabels = FALSE, showgrid = FALSE, zeroline = FALSE),
  legend = list(title = list(text = ""))
)

plotly::config(
  raw_umap_plot,
  displaylogo = FALSE,
  modeBarButtonsToRemove = c("select2d", "lasso2d")
)

Now use precise_viz(views = "embedding") and precise_viz(views = "heatmap") on the baseline matrices. The visualization engine is explicit. For this article, standalone panels use Plotly when it is available; precise_viz() still returns the plot object and plotted data in a panel tibble. For embedding, UMAP is the default coordinate engine, MDS remains available as an explicit dependency-free fallback, and the held-back labels can be passed as color annotations.

baseline_embeddings <- precise_viz(
  baseline_distances,
  views = "embedding",
  params = list(
    embedding = list(
      method = "umap",
      n_neighbors = 15,
      min_dist = 0.01,
      spread = 1,
      engine = workflow_plot_engine,
      seed = 42,
      color = cell_cycle_groups,
      colors = cell_cycle_palette,
      show_labels = FALSE
    )
  ),
  verbose = FALSE
)

baseline_embeddings[, c("panel_id", "input", "metric", "view", "panel_class")]
#> # A tibble: 2 × 5
#>   panel_id               input       metric      view      panel_class
#>   <chr>                  <chr>       <chr>       <chr>     <chr>      
#> 1 manhattan__embedding   manhattan   manhattan   embedding plotly     
#> 2 correlation__embedding correlation correlation embedding plotly

manhattan_embedding <- baseline_embeddings$panel[[
  match("manhattan", baseline_embeddings$metric)
]]

correlation_embedding <- baseline_embeddings$panel[[
  match("correlation", baseline_embeddings$metric)
]]

Manhattan distance is a classic alternative to Euclidean distance. It measures absolute coordinate-wise differences rather than straight-line displacement, so it is a useful first reminder that even familiar distance metrics encode different assumptions.

manhattan_embedding

Correlation distance asks a different question: whether observations move together by pattern rather than raw magnitude. The object below is correlation after explicit conversion into distance form.

correlation_embedding

baseline_heatmaps <- precise_viz(
  baseline_distances,
  views = "heatmap",
  params = list(
    heatmap = list(
      engine = workflow_plot_engine
    )
  ),
  verbose = FALSE
)

baseline_heatmaps[, c("panel_id", "input", "metric", "view", "panel_class")]
#> # A tibble: 2 × 5
#>   panel_id             input       metric      view    panel_class
#>   <chr>                <chr>       <chr>       <chr>   <chr>      
#> 1 manhattan__heatmap   manhattan   manhattan   heatmap plotly     
#> 2 correlation__heatmap correlation correlation heatmap plotly

manhattan_heatmap <- baseline_heatmaps$panel[[
  match("manhattan", baseline_heatmaps$metric)
]]

correlation_heatmap <- baseline_heatmaps$panel[[
  match("correlation", baseline_heatmaps$metric)
]]

manhattan_heatmap

correlation_heatmap

The labels in the row names are only a visual aid. They were not used to compute the matrices. At this stage we should ask a simple question: do these ordinary baseline choices make the structure easy to see? If the answer is “not really,” that is not a failure of Manhattan distance or correlation distance. It is a reason to stop treating either one as an unquestioned default.

4. Broaden the Metric Search

Once the baseline is visible, we widen the search. The goal is not to run every possible metric for its own sake. The goal is to admit that the correct relatedness metric is rarely known beforehand, then make a controlled attempt to see which metrics produce coherent and interpretable structure.

The first thing to do is ask what the package can run. A user should not have to guess the available metric names, and the package should make the breadth of the search visible before any narrowing happens.

Because cell_cycle is an ordinary numeric matrix, the most natural broad family is static_dists: the registered metrics meant for a plain numeric sweep. Families such as probability, binary, nominal, transport, and time-series metrics are still available, but they have more specific input assumptions.

available_static_distances <- precise_dist_list("static_dists")

length(available_static_distances)
#> [1] 53
available_static_distances |>
  unlist() |>
  as.matrix()
#>       [,1]                
#>  [1,] "euclidean"         
#>  [2,] "supremum"          
#>  [3,] "manhattan"         
#>  [4,] "minkowski_0.5"     
#>  [5,] "minkowski_1.5"     
#>  [6,] "minkowski_2.5"     
#>  [7,] "sorensen"          
#>  [8,] "gower"             
#>  [9,] "gower_phil"        
#> [10,] "soergel"           
#> [11,] "soergel_phil"      
#> [12,] "kulczynski_d"      
#> [13,] "canberra"          
#> [14,] "canberra_phil"     
#> [15,] "lorentzian"        
#> [16,] "wave"              
#> [17,] "wavehedges"        
#> [18,] "intersection"      
#> [19,] "non_intersection"  
#> [20,] "czekanowski"       
#> [21,] "tanimoto"          
#> [22,] "ruzicka"           
#> [23,] "inner_product"     
#> [24,] "harmonic_mean"     
#> [25,] "cosine"            
#> [26,] "fjaccard"          
#> [27,] "ejaccard"          
#> [28,] "edice"             
#> [29,] "dice"              
#> [30,] "squared_euclidean" 
#> [31,] "squared_chi"       
#> [32,] "divergence"        
#> [33,] "clark"             
#> [34,] "additive_symm"     
#> [35,] "correlation"       
#> [36,] "mahalanobis"       
#> [37,] "bray"              
#> [38,] "chord"             
#> [39,] "geodesic"          
#> [40,] "whittaker"         
#> [41,] "avg"               
#> [42,] "random_forest_sqrt"
#> [43,] "random_forest_two" 
#> [44,] "random_forest_half"
#> [45,] "laplace_1"         
#> [46,] "laplace_2"         
#> [47,] "laplace_3"         
#> [48,] "rbf_1"             
#> [49,] "rbf_2"             
#> [50,] "rbf_3"             
#> [51,] "tsne_5"            
#> [52,] "tsne_30"           
#> [53,] "tsne_50"

identical(precise_dist_list("static"), available_static_distances)
#> [1] TRUE

The last line shows a convenience alias: static is equivalent to static_dists. The explicit _dists form is still used in prose because it makes clear that this is a family keyword, not an individual metric.

For a full local analysis, the broad sweep is the first place where the parallel backend matters. parallel = TRUE is only a flag telling PreciseDist to use the active future plan; the user still controls CPU and memory through future. This keeps the package portable across a laptop, a workstation, or a larger scheduler-backed environment.

library(future)

future::plan(
  future::multisession,
  workers = min(4, future::availableCores())
)

options(future.globals.maxSize = 8 * 1024^3)

cell_cycle_static <- precise_dist(
  cell_cycle,
  dists = "static",
  parallel = TRUE,
  max_time_per_metric = 60,
  verbose = TRUE
)

future::plan(future::sequential)

In that block, workers is the CPU limit, future.globals.maxSize is the memory limit for exported objects, and max_time_per_metric is the per-metric time limit inside precise_dist(). A full sweep is exactly where these controls matter: some metrics finish quickly, some are expensive, and some may be inappropriate for the data — and every metric that errors or times out is recorded in the run ledger (precise_run_report()) rather than silently dropped.

If you want the run checkpointed on disk, add file=. The envelope is written by the master process, so it combines with parallel = TRUE, and it doubles as exact persistence of the finished object.

cache_file <- "cell_cycle_static.rds"

cell_cycle_static <- precise_dist(
  cell_cycle,
  dists = "static",
  file = cache_file,
  parallel = TRUE,
  max_time_per_metric = 60,
  verbose = TRUE
)

cell_cycle_static_reread <- precise_read(cache_file)
precise_run_report(cell_cycle_static_reread)

For the rendered website, we will run a real broad sweep drawn from the static family, but we will leave out the rows that are intentionally expensive or have more specific input assumptions. This is still much broader than the set we will ultimately fuse, and it gives us the evidence needed to justify that later narrowing.

The chunk below uses the future backend for an ordinary installed-package or website render. It falls back to a single-process run when rendered directly from source, because future workers cannot attach an uninstalled package, and during R CMD check, because package-check hosts are not required to permit worker sockets. The analytic object is the same in every case.

candidate_metrics <- c(
  "euclidean",
  "supremum",
  "manhattan",
  "canberra",
  "canberra_phil",
  "clark",
  "divergence",
  "wavehedges",
  "minkowski_0.5",
  "minkowski_1.5",
  "minkowski_2.5",
  "sorensen",
  "gower",
  "gower_phil",
  "soergel",
  "soergel_phil",
  "kulczynski_d",
  "lorentzian",
  "wave",
  "czekanowski",
  "tanimoto",
  "dice",
  "squared_euclidean",
  "squared_chi",
  "additive_symm",
  "bray",
  "chord",
  "whittaker",
  "avg",
  "cosine",
  "correlation",
  "laplace_1",
  "laplace_2",
  "rbf_1",
  "rbf_2"
)
candidate_metrics <- candidate_metrics[candidate_metrics %in% unlist(available_static_distances)]

candidate_metrics
#>  [1] "euclidean"         "supremum"          "manhattan"        
#>  [4] "canberra"          "canberra_phil"     "clark"            
#>  [7] "divergence"        "wavehedges"        "minkowski_0.5"    
#> [10] "minkowski_1.5"     "minkowski_2.5"     "sorensen"         
#> [13] "gower"             "gower_phil"        "soergel"          
#> [16] "soergel_phil"      "kulczynski_d"      "lorentzian"       
#> [19] "wave"              "czekanowski"       "tanimoto"         
#> [22] "dice"              "squared_euclidean" "squared_chi"      
#> [25] "additive_symm"     "bray"              "chord"            
#> [28] "whittaker"         "avg"               "cosine"           
#> [31] "correlation"       "laplace_1"         "laplace_2"        
#> [34] "rbf_1"             "rbf_2"

candidate_parallel <- !isTRUE(precisedist_source_loaded) &&
  !nzchar(Sys.getenv("_R_CHECK_PACKAGE_NAME_")) &&
  requireNamespace("future", quietly = TRUE) &&
  requireNamespace("future.apply", quietly = TRUE)

if (candidate_parallel) {
  future::plan(
    future::multisession,
    workers = min(2, future::availableCores())
  )
  options(future.globals.maxSize = 2 * 1024^3)
}

candidate_run_log <- capture.output({
  candidate_native <- precise_dist(
    cell_cycle,
    dists = candidate_metrics,
    parallel = candidate_parallel,
    max_time_per_metric = 20,
    verbose = TRUE
  )
})

if (candidate_parallel) {
  future::plan(future::sequential)
}

candidate_missing <- setdiff(candidate_metrics, candidate_native$metric)

data.frame(
  requested = length(candidate_metrics),
  completed = nrow(candidate_native),
  not_completed = length(candidate_missing),
  used_future_backend = candidate_parallel
)
#>   requested completed not_completed used_future_backend
#> 1        35        35             0                TRUE

candidate_missing
#> character(0)

candidate_native[, c("distance", "metric", "type", "time_taken_seconds")]
#> # A tibble: 35 × 4
#>    distance      metric        type     time_taken_seconds
#>    <chr>         <chr>         <chr>                 <dbl>
#>  1 euclidean     euclidean     distance            0.0100 
#>  2 supremum      supremum      distance            0.00200
#>  3 manhattan     manhattan     distance            0.00300
#>  4 canberra      canberra      distance            0.00400
#>  5 canberra_phil canberra_phil distance            0.100  
#>  6 clark         clark         distance            0.00300
#>  7 divergence    divergence    distance            0.0310 
#>  8 wavehedges    wavehedges    distance            0.00300
#>  9 minkowski_0.5 minkowski_0.5 distance            0.0180 
#> 10 minkowski_1.5 minkowski_1.5 distance            0.0170 
#> # ℹ 25 more rows

candidate_distances <- candidate_native |>
  precise_transform(
    to = "distance",
    normalize = "range01",
    diagonal = 0
  )

candidate_distances[, c("distance", "type")]
#> # A tibble: 35 × 2
#>    distance      type    
#>    <chr>         <chr>   
#>  1 euclidean     distance
#>  2 supremum      distance
#>  3 manhattan     distance
#>  4 canberra      distance
#>  5 canberra_phil distance
#>  6 clark         distance
#>  7 divergence    distance
#>  8 wavehedges    distance
#>  9 minkowski_0.5 distance
#> 10 minkowski_1.5 distance
#> # ℹ 25 more rows

The completed set is deliberately mixed. Most rows are native distances, but cosine enters as a similarity, correlation enters as a correlation, and the kernel rows enter as affinities. They are allowed into the common distance space only because precise_transform() has registered conversion rules for them. If a metric cannot be converted honestly, the transform step should refuse it rather than fabricate a number.

The normalization is also explicit. normalize = "range01" puts each converted matrix on a common [0, 1] scale, and diagonal = 0 restores the conventional zero self-distance after the transformation. This is useful before comparing or fusing many matrices, but it is still a modeling choice. precise_transform() does not assume the user’s intent; it applies only the operations requested and records those operations in the output history.

Before looking at plots, we should ask for a quantitative readout of each matrix. precise_diagnostics() does not cluster the data and does not choose a winner. It gives us scalar summaries that can pressure-test what we think we see: spread, ties, nearest-neighbour contrast, hubness, graph fragmentation, and when labels are available, label separation.

In this tutorial we pass the known cell-cycle labels because they make the diagnostic columns easier to interpret. In a genuinely unsupervised analysis, leave labels = NULL and focus on the geometry-only columns.

candidate_diagnostics <- precise_diagnostics(
  candidate_distances,
  labels = cell_cycle_groups,
  k = c(5, 10, 15),
  verbose = FALSE
)

diagnostic_columns <- c(
  "distance",
  "metric",
  "diagnostic_score",
  "label_within_between_ratio",
  "label_silhouette",
  "label_pair_auc",
  "knn_k5_label_agreement",
  "knn_k5_components",
  "knn_k5_hubness_gini",
  "value_prop_tied",
  "nn_contrast_median"
)

candidate_diagnostics$rank_auc <- rank(
  -candidate_diagnostics$label_pair_auc,
  ties.method = "average",
  na.last = "keep"
)
candidate_diagnostics$rank_k5 <- rank(
  -candidate_diagnostics$knn_k5_label_agreement,
  ties.method = "average",
  na.last = "keep"
)
candidate_diagnostics$rank_ratio <- rank(
  candidate_diagnostics$label_within_between_ratio,
  ties.method = "average",
  na.last = "keep"
)
candidate_diagnostics$rank_silhouette <- rank(
  -candidate_diagnostics$label_silhouette,
  ties.method = "average",
  na.last = "keep"
)
candidate_diagnostics$rank_hubness <- rank(
  candidate_diagnostics$knn_k5_hubness_gini,
  ties.method = "average",
  na.last = "keep"
)

candidate_diagnostics$diagnostic_score <-
  candidate_diagnostics$rank_auc +
  candidate_diagnostics$rank_k5 +
  candidate_diagnostics$rank_ratio +
  candidate_diagnostics$rank_silhouette +
  0.5 * candidate_diagnostics$rank_hubness

candidate_diagnostics[
  order(candidate_diagnostics$diagnostic_score),
  diagnostic_columns
]
#> # A tibble: 35 × 11
#>    distance      metric diagnostic_score label_within_between…¹ label_silhouette
#>    <chr>         <chr>             <dbl>                  <dbl>            <dbl>
#>  1 minkowski_0.5 minko…              7                    0.921           0.0441
#>  2 lorentzian    loren…             24.5                  0.954           0.0261
#>  3 gower         gower              26                    0.948           0.0241
#>  4 manhattan     manha…             29.8                  0.950           0.0221
#>  5 gower_phil    gower…             29.8                  0.950           0.0221
#>  6 avg           avg                37                    0.951           0.0219
#>  7 wavehedges    waveh…             41                    0.964           0.0209
#>  8 soergel       soerg…             46                    0.954           0.0211
#>  9 squared_chi   squar…             52.5                  0.945           0.0135
#> 10 canberra_phil canbe…             57                    0.966           0.0166
#> # ℹ 25 more rows
#> # ℹ abbreviated name: ¹label_within_between_ratio
#> # ℹ 6 more variables: label_pair_auc <dbl>, knn_k5_label_agreement <dbl>,
#> #   knn_k5_components <int>, knn_k5_hubness_gini <dbl>, value_prop_tied <dbl>,
#> #   nn_contrast_median <dbl>

These numbers are not a substitute for seeing the matrices. They are a way to make the visual search more honest. A matrix with clean-looking structure but extreme hubness, weak nearest-neighbour contrast, or many tied values deserves skepticism. A matrix with good diagnostics but poor visual separation may still be useful downstream, but it should not be accepted blindly.

5. Browse Candidate Views

The diagnostic table gives us a quantitative first pass. It does not show what each metric does to the observations. The next step is to look at each candidate distance directly, first as an embedding and then as a heatmap.

For each candidate distance, we build one UMAP and one heatmap, then browse those panels with Trelliscope. precise_viz() builds the panel tibble, and precise_trellis() turns that tibble into a browsable display.

For the website build, the candidate set is small enough to keep the embedded Trelliscope browser responsive. In a full local analysis, use the full static sweep object in the same calls.

candidate_embedding_panels <- precise_viz(
  candidate_distances,
  views = "embedding",
  params = list(
    embedding = list(
      method = "umap",
      n_neighbors = 15,
      min_dist = 0.01,
      spread = 1,
      engine = workflow_plot_engine,
      seed = 42,
      color = cell_cycle_groups,
      colors = cell_cycle_palette,
      show_labels = FALSE
    )
  ),
  diagnostics = candidate_diagnostics,
  parallel = FALSE,
  verbose = FALSE
)

candidate_embedding_panels[
  ,
  c(
    "panel_id",
    "input",
    "metric",
    "view",
    "panel_class",
    "label_within_between_ratio",
    "knn_k5_label_agreement",
    "knn_k5_components"
  )
]
#> # A tibble: 35 × 8
#>    panel_id                input metric view  panel_class label_within_between…¹
#>    <chr>                   <chr> <chr>  <chr> <chr>                        <dbl>
#>  1 euclidean__embedding    eucl… eucli… embe… plotly                       0.971
#>  2 supremum__embedding     supr… supre… embe… plotly                       0.981
#>  3 manhattan__embedding    manh… manha… embe… plotly                       0.950
#>  4 canberra__embedding     canb… canbe… embe… plotly                       0.965
#>  5 canberra_phil__embeddi… canb… canbe… embe… plotly                       0.966
#>  6 clark__embedding        clark clark  embe… plotly                       0.980
#>  7 divergence__embedding   dive… diver… embe… plotly                       0.961
#>  8 wavehedges__embedding   wave… waveh… embe… plotly                       0.964
#>  9 minkowski_0.5__embeddi… mink… minko… embe… plotly                       0.921
#> 10 minkowski_1.5__embeddi… mink… minko… embe… plotly                       0.961
#> # ℹ 25 more rows
#> # ℹ abbreviated name: ¹label_within_between_ratio
#> # ℹ 2 more variables: knn_k5_label_agreement <dbl>, knn_k5_components <int>

candidate_heatmap_panels <- precise_viz(
  candidate_distances,
  views = "heatmap",
  params = list(
    heatmap = list(
      cluster = TRUE,
      method = "complete",
      engine = workflow_plot_engine
    )
  ),
  diagnostics = candidate_diagnostics,
  parallel = FALSE,
  verbose = FALSE
)

candidate_heatmap_panels[
  ,
  c(
    "panel_id",
    "input",
    "metric",
    "view",
    "panel_class",
    "label_silhouette",
    "label_pair_auc",
    "value_prop_tied"
  )
]
#> # A tibble: 35 × 8
#>    panel_id       input metric view  panel_class label_silhouette label_pair_auc
#>    <chr>          <chr> <chr>  <chr> <chr>                  <dbl>          <dbl>
#>  1 euclidean__he… eucl… eucli… heat… plotly               0.00361          0.640
#>  2 supremum__hea… supr… supre… heat… plotly              -0.00457          0.598
#>  3 manhattan__he… manh… manha… heat… plotly               0.0221           0.691
#>  4 canberra__hea… canb… canbe… heat… plotly              -0.00162          0.591
#>  5 canberra_phil… canb… canbe… heat… plotly               0.0166           0.636
#>  6 clark__heatmap clark clark  heat… plotly               0.00440          0.613
#>  7 divergence__h… dive… diver… heat… plotly               0.00856          0.613
#>  8 wavehedges__h… wave… waveh… heat… plotly               0.0209           0.656
#>  9 minkowski_0.5… mink… minko… heat… plotly               0.0441           0.727
#> 10 minkowski_1.5… mink… minko… heat… plotly               0.0104           0.660
#> # ℹ 25 more rows
#> # ℹ 1 more variable: value_prop_tied <dbl>

The full browsing version is the same workflow inside a Trelliscope browser. When the optional trelliscopejs package is installed, this vignette embeds the browser directly. The supporting app files are written to a temporary directory during rendering, so the repository is not filled with generated browsing artifacts.

Because the panels were built with diagnostics = candidate_diagnostics, Trelliscope also receives the diagnostic columns as sortable and filterable cognostics. This lets the user look at an embedding or heatmap while keeping the quantitative matrix summaries visible in the same browser.

Trelliscope is an application-style viewer embedded inside a webpage. In this article layout, long label and sorting controls can be visually clipped in the embedded view. Use the Trelliscope fullscreen button when you want to inspect labels carefully.

precise_trellis(
  candidate_embedding_panels,
  key_col = "input",
  path = "trelliscope/cell_cycle_umap",
  name = "cell_cycle_umap",
  group = "common",
  desc = "UMAP embeddings from candidate distance matrices.",
  height = 500,
  width = 500
)

precise_trellis(
  candidate_heatmap_panels,
  key_col = "input",
  path = "trelliscope/cell_cycle_heatmap",
  name = "cell_cycle_heatmap",
  group = "common",
  desc = "Clustered heatmaps from candidate distance matrices.",
  height = 500,
  width = 500
)

These views are useful because they put pressure on the metrics in different ways. An embedding asks whether a distance creates a readable low-dimensional separation. A heatmap asks whether the full pairwise structure contains coherent blocks. Neither view proves that clusters are real, but both views can expose structure that summary statistics alone would flatten.

Here, the diagnostics and visual panels start to tell the same story. Metrics such as minkowski_0.5, additive_symm, lorentzian, and wavehedges show stronger local label agreement and clearer visual structure than the ordinary baselines in this example. We should still ask one more question before fusing them: are these metrics telling the same story, complementary stories, or unrelated stories?

6. Tune the Embedding Lens

The first UMAP pass used one ordinary setting. That is useful for browsing, but it is not enough for choosing. UMAP has its own modeling choices. A metric that only looks good at one narrow setting is less convincing than a metric whose structure persists across several reasonable settings.

Here we take the strongest distance candidates from the diagnostic table and render them across a compact UMAP grid. This is not a search for the single “best” UMAP. It is a way to ask whether the visual conclusion is robust to neighborhood size, local compression, and global spread.

diagnostic_order <- candidate_diagnostics$metric[
  order(candidate_diagnostics$diagnostic_score)
]

visual_anchor_metrics <- c(
  "minkowski_0.5",
  "additive_symm",
  "lorentzian",
  "wavehedges",
  "canberra_phil",
  "clark",
  "divergence"
)

umap_search_metrics <- unique(c(
  head(diagnostic_order, 3),
  visual_anchor_metrics
))

umap_search_metrics <- umap_search_metrics[
  umap_search_metrics %in% candidate_distances$metric
]

umap_search_metrics <- head(umap_search_metrics, 8)

umap_search_distances <- candidate_distances[
  candidate_distances$metric %in% umap_search_metrics,
]

umap_search_distances <- umap_search_distances[
  match(umap_search_metrics, umap_search_distances$metric),
]

umap_search_grid <- expand.grid(
  n_neighbors = c(5L, 15L, 30L),
  min_dist = c(0.001, 0.05),
  spread = c(1, 3),
  KEEP.OUT.ATTRS = FALSE
)

umap_search_grid$umap_setting <- sprintf(
  "nn%s_md%s_sp%s",
  umap_search_grid$n_neighbors,
  umap_search_grid$min_dist,
  umap_search_grid$spread
)

umap_search_metrics
#> [1] "minkowski_0.5" "lorentzian"    "gower"         "additive_symm"
#> [5] "wavehedges"    "canberra_phil" "clark"         "divergence"
umap_search_grid
#>    n_neighbors min_dist spread     umap_setting
#> 1            5    0.001      1  nn5_md0.001_sp1
#> 2           15    0.001      1 nn15_md0.001_sp1
#> 3           30    0.001      1 nn30_md0.001_sp1
#> 4            5    0.050      1   nn5_md0.05_sp1
#> 5           15    0.050      1  nn15_md0.05_sp1
#> 6           30    0.050      1  nn30_md0.05_sp1
#> 7            5    0.001      3  nn5_md0.001_sp3
#> 8           15    0.001      3 nn15_md0.001_sp3
#> 9           30    0.001      3 nn30_md0.001_sp3
#> 10           5    0.050      3   nn5_md0.05_sp3
#> 11          15    0.050      3  nn15_md0.05_sp3
#> 12          30    0.050      3  nn30_md0.05_sp3

The rendered grid is intentionally small: at most eight metrics, three neighborhood sizes, two local compression levels, and two spread values. That is enough to reveal fragile views without turning the vignette into a benchmark.

umap_search_panels <- lapply(seq_len(nrow(umap_search_grid)), function(i) {
  setting <- umap_search_grid[i, ]

  panels <- precise_viz(
    umap_search_distances,
    views = "embedding",
    params = list(
      embedding = list(
        method = "umap",
        n_neighbors = setting$n_neighbors,
        min_dist = setting$min_dist,
        spread = setting$spread,
        engine = workflow_plot_engine,
        seed = 42,
        color = cell_cycle_groups,
        colors = cell_cycle_palette,
        show_labels = FALSE
      )
    ),
    diagnostics = candidate_diagnostics,
    parallel = FALSE,
    verbose = FALSE
  )

  panels$n_neighbors <- setting$n_neighbors
  panels$min_dist <- setting$min_dist
  panels$spread <- setting$spread
  panels$umap_setting <- setting$umap_setting
  panels$panel_id <- paste(panels$metric, panels$umap_setting, sep = "__")
  panels
})

umap_search_panels <- dplyr::bind_rows(umap_search_panels)

umap_search_panels[
  ,
  c(
    "panel_id",
    "metric",
    "n_neighbors",
    "min_dist",
    "spread",
    "label_pair_auc",
    "knn_k5_label_agreement",
    "label_silhouette"
  )
]
#> # A tibble: 96 × 8
#>    panel_id                    metric n_neighbors min_dist spread label_pair_auc
#>    <chr>                       <chr>        <int>    <dbl>  <dbl>          <dbl>
#>  1 minkowski_0.5__nn5_md0.001… minko…           5    0.001      1          0.727
#>  2 lorentzian__nn5_md0.001_sp1 loren…           5    0.001      1          0.720
#>  3 gower__nn5_md0.001_sp1      gower            5    0.001      1          0.696
#>  4 additive_symm__nn5_md0.001… addit…           5    0.001      1          0.570
#>  5 wavehedges__nn5_md0.001_sp1 waveh…           5    0.001      1          0.656
#>  6 canberra_phil__nn5_md0.001… canbe…           5    0.001      1          0.636
#>  7 clark__nn5_md0.001_sp1      clark            5    0.001      1          0.613
#>  8 divergence__nn5_md0.001_sp1 diver…           5    0.001      1          0.613
#>  9 minkowski_0.5__nn15_md0.00… minko…          15    0.001      1          0.727
#> 10 lorentzian__nn15_md0.001_s… loren…          15    0.001      1          0.720
#> # ℹ 86 more rows
#> # ℹ 2 more variables: knn_k5_label_agreement <dbl>, label_silhouette <dbl>

Browse this object the same way as the first candidate screen. In Trelliscope, filter by metric, then move across n_neighbors, min_dist, and spread. The goal is to find metrics whose separation is not an artifact of one UMAP setting.

precise_trellis(
  umap_search_panels,
  key_col = "panel_id",
  path = "trelliscope/cell_cycle_umap_settings",
  name = "cell_cycle_umap_settings",
  group = "common",
  desc = "UMAP parameter sweep over the strongest candidate distance matrices.",
  height = 500,
  width = 500
)

This is the point where the visual workflow becomes disciplined. We are not only asking “which distance looks good?” We are asking “which distance keeps showing useful structure when the embedding lens changes?”

7. Treat Metrics as Objects

Before fusing anything, ask whether the candidate matrices agree. Strongly agreeing metrics may be redundant. Weakly agreeing metrics may add information, or they may be noise. The agreement matrix does not make the decision for us, but it changes the decision from a hunch into an object we can inspect.

This section will run precise_correlations() and render the same information with the agreement view in precise_viz().

metric_agreement <- precise_correlations(
  candidate_distances,
  method = "pearson",
  permutations = 99,
  parallel = FALSE,
  verbose = FALSE
)

agreement_matrix <- metric_agreement$statistic
agreement_pair_index <- upper.tri(agreement_matrix)

agreement_pairs <- data.frame(
  metric_1 = rownames(agreement_matrix)[row(agreement_matrix)[agreement_pair_index]],
  metric_2 = colnames(agreement_matrix)[col(agreement_matrix)[agreement_pair_index]],
  pearson = agreement_matrix[agreement_pair_index],
  row.names = NULL
)

agreement_pairs <- agreement_pairs[!is.na(agreement_pairs$pearson), ]

strongest_agreements <- agreement_pairs[
  order(-agreement_pairs$pearson),
]

weakest_agreements <- agreement_pairs[
  order(agreement_pairs$pearson),
]

head(strongest_agreements, 10)
#>         metric_1    metric_2   pearson
#> 81     manhattan  gower_phil 1.0000000
#> 183     sorensen czekanowski 1.0000000
#> 206 soergel_phil    tanimoto 1.0000000
#> 312     sorensen        bray 1.0000000
#> 320  czekanowski        bray 1.0000000
#> 433        chord      cosine 1.0000000
#> 528    laplace_1   laplace_2 1.0000000
#> 381    manhattan         avg 0.9999900
#> 392   gower_phil         avg 0.9999900
#> 595        rbf_1       rbf_2 0.9998546

The full statistic matrix is usually too large to print usefully in a vignette. The pair table is a more readable first pass: it shows which metrics are nearly telling the same story. The least-agreeing pairs are also useful, because a metric that disagrees with everything else may be complementary, noisy, or capturing a genuinely different geometry.

head(weakest_agreements, 10)
#>     metric_1      metric_2    pearson
#> 278 supremum additive_symm 0.06280504
#> 17  supremum    divergence 0.06425962
#> 12  supremum         clark 0.06515364
#> 353 supremum     whittaker 0.06935708
#> 8   supremum canberra_phil 0.08752872
#> 5   supremum      canberra 0.08955204
#> 408 supremum        cosine 0.09623093
#> 327 supremum         chord 0.09623093
#> 23  supremum    wavehedges 0.10582265
#> 437 supremum   correlation 0.15016751

For a larger run, the same agreement calculation can use the registered parallel backend. This is usually more relevant when permutations is high or the candidate set contains many matrices.

library(future)

future::plan(
  future::multisession,
  workers = min(4, future::availableCores())
)

metric_agreement <- precise_correlations(
  candidate_distances,
  method = "pearson",
  permutations = 99,
  parallel = TRUE,
  seed = 1,
  verbose = TRUE
)

future::plan(future::sequential)

The statistic matrix is now a relatedness object about the metrics themselves. Here, some metrics move together tightly while others tell a different story. That does not make one block “right” and another block “wrong.” It tells us where fusion would be redundant, where fusion would be diverse, and where a metric may need closer visual inspection before it is trusted.

agreement_panels <- precise_viz(
  candidate_distances,
  views = "agreement",
  params = list(
    agreement = list(
      engine = workflow_plot_engine
    )
  ),
  agreement = metric_agreement,
  verbose = FALSE
)

agreement_panels$panel[[1]]

The vignette uses 99 permutations to keep the website fast. For a real analysis, increase permutations so the p-value floor is appropriate for the decision being made.

8. Choose a Metric Set

Now we have three kinds of evidence in front of us:

scalar diagnostics from precise_diagnostics();
embeddings, heatmaps, and UMAP-setting sweeps browsed through precise_viz() and Trelliscope;
pairwise agreement among the candidate matrices.

The selected set below is not a universal recommendation. It is the choice this particular run supports. minkowski_0.5, additive_symm, and lorentzian score well in the diagnostics, while wavehedges preserves useful visual structure and adds a different but still related view of the observations.

selected_metrics <- c("minkowski_0.5", "additive_symm", "lorentzian", "wavehedges")

selected_metric_evidence <- candidate_diagnostics[
  candidate_diagnostics$metric %in% selected_metrics,
  diagnostic_columns
]

selected_metric_evidence[
  match(selected_metrics, selected_metric_evidence$metric),
]
#> # A tibble: 4 × 11
#>   distance      metric  diagnostic_score label_within_between…¹ label_silhouette
#>   <chr>         <chr>              <dbl>                  <dbl>            <dbl>
#> 1 minkowski_0.5 minkow…              7                    0.921           0.0441
#> 2 additive_symm additi…             90                    0.943           0.0160
#> 3 lorentzian    lorent…             24.5                  0.954           0.0261
#> 4 wavehedges    wavehe…             41                    0.964           0.0209
#> # ℹ abbreviated name: ¹label_within_between_ratio
#> # ℹ 6 more variables: label_pair_auc <dbl>, knn_k5_label_agreement <dbl>,
#> #   knn_k5_components <int>, knn_k5_hubness_gini <dbl>, value_prop_tied <dbl>,
#> #   nn_contrast_median <dbl>

9. Fuse, Graph, and Inspect the Selected Metrics

We learn a lot from the embedding, heatmap, and agreement views above. Several distances show structure by themselves, and the diagnostics help us avoid choosing by eye alone. None of that proves these are the “right” distances. It does give us a concrete, inspectable reason to isolate them and ask what their fused structure looks like.

This is the point of the package. PreciseDist does not pretend there is a universal rule for choosing a distance. It gives the user enough parallel views of the same objects to make a defensible choice.

First, extract those distances from the typed object precise_dist() gave us. The output is a tibble, so ordinary row filtering is clear and explicit. We filter by metric, not distance, because metric is the canonical registry name.

selected_distances <- candidate_distances[
  candidate_distances$metric %in% selected_metrics,
]

selected_distances <- selected_distances[
  match(selected_metrics, selected_distances$metric),
]

selected_distances[, c("distance", "metric", "type")]
#> # A tibble: 4 × 3
#>   distance      metric        type    
#>   <chr>         <chr>         <chr>   
#> 1 minkowski_0.5 minkowski_0.5 distance
#> 2 additive_symm additive_symm distance
#> 3 lorentzian    lorentzian    distance
#> 4 wavehedges    wavehedges    distance

This step is intentionally plain R subsetting. The important object is still the typed tibble returned by the package; the choice of rows is the user’s analytic judgment.

Now we will fuse the distances. precise_fusion() can run several fusion rules and return one typed row per method. That matters for the same reason metric breadth matters. There is rarely a single obviously correct way to combine relatedness matrices, so we should look before choosing.

selected_fusions <- precise_fusion(
  selected_distances,
  methods = "all",
  verbose = FALSE
)

selected_fusions[, c("method", "type", "time_taken_seconds")]
#> # A tibble: 8 × 3
#>   method        type       time_taken_seconds
#>   <chr>         <chr>                   <dbl>
#> 1 mean          distance                0    
#> 2 weighted_mean distance                0    
#> 3 median        distance                0.105
#> 4 trimmed_mean  distance                0.087
#> 5 rank_mean     distance                0.002
#> 6 analogue_fuse distance                0.001
#> 7 snf           similarity              0.032
#> 8 distatis      distance                0.013
selected_fusions$inputs[[1]]
#> [1] "minkowski_0.5" "additive_symm" "lorentzian"    "wavehedges"

On this machine, the optional fusion backends are installed, so "all" includes the base reducers, analogue’s max-scaled fuse, Similarity Network Fusion, and DiSTATIS. On a lean installation, "all" still means all compatible methods whose backends are available. If you request methods = "snf" or methods = "distatis" explicitly and the backend is missing, precise_fusion() will error before doing partial work.

The fused matrices are still just relatedness matrices, but each row carries a different claim: not “this metric is correct,” but “this rule turns these selected views into this consensus representation.” Those claims can now be inspected like any other.

First, prepare heatmaps for the fused matrices, then compare them in Trelliscope. This view includes both distance and similarity outputs, because a heatmap can display either type honestly when the type is declared.

selected_fusion_heatmaps <- precise_viz(
  selected_fusions,
  views = "heatmap",
  params = list(
    heatmap = list(
      cluster = TRUE,
      method = "complete",
      engine = workflow_plot_engine
    )
  ),
  verbose = FALSE
)

selected_fusion_heatmaps[, c("panel_id", "input", "input_type", "view")]
#> # A tibble: 8 × 4
#>   panel_id               input         input_type view   
#>   <chr>                  <chr>         <chr>      <chr>  
#> 1 mean__heatmap          mean          distance   heatmap
#> 2 weighted_mean__heatmap weighted_mean distance   heatmap
#> 3 median__heatmap        median        distance   heatmap
#> 4 trimmed_mean__heatmap  trimmed_mean  distance   heatmap
#> 5 rank_mean__heatmap     rank_mean     distance   heatmap
#> 6 analogue_fuse__heatmap analogue_fuse distance   heatmap
#> 7 snf__heatmap           snf           similarity heatmap
#> 8 distatis__heatmap      distatis      distance   heatmap

precise_trellis(
  selected_fusion_heatmaps,
  key_col = "input",
  path = "trelliscope/cell_cycle_fusion_heatmap",
  name = "cell_cycle_fusion_heatmap",
  group = "common",
  desc = "Heatmaps from selected-metric fusion methods.",
  height = 500,
  width = 500
)

Next, prepare embeddings for the fusion outputs that are distances, then browse them together. SNF returns a similarity, so it belongs in the heatmap comparison above and in graph-style similarity workflows below; it is not coerced into a distance just to satisfy an embedding view.

selected_distance_fusions <- selected_fusions[
  selected_fusions$type == "distance",
]

selected_fusion_embeddings <- precise_viz(
  selected_distance_fusions,
  views = "embedding",
  params = list(
    embedding = list(
      method = "umap",
      n_neighbors = 15,
      min_dist = 0.01,
      spread = 1,
      engine = workflow_plot_engine,
      seed = 42,
      color = cell_cycle_groups,
      colors = cell_cycle_palette,
      show_labels = FALSE
    )
  ),
  verbose = FALSE
)

selected_fusion_embeddings[, c("panel_id", "input", "input_type", "view")]
#> # A tibble: 7 × 4
#>   panel_id                 input         input_type view     
#>   <chr>                    <chr>         <chr>      <chr>    
#> 1 mean__embedding          mean          distance   embedding
#> 2 weighted_mean__embedding weighted_mean distance   embedding
#> 3 median__embedding        median        distance   embedding
#> 4 trimmed_mean__embedding  trimmed_mean  distance   embedding
#> 5 rank_mean__embedding     rank_mean     distance   embedding
#> 6 analogue_fuse__embedding analogue_fuse distance   embedding
#> 7 distatis__embedding      distatis      distance   embedding

precise_trellis(
  selected_fusion_embeddings,
  key_col = "input",
  path = "trelliscope/cell_cycle_fusion_umap",
  name = "cell_cycle_fusion_umap",
  group = "common",
  desc = "UMAP embeddings from selected-metric distance fusion methods.",
  height = 500,
  width = 500
)

After looking across the fusion views, we will carry forward analogue’s max-scaled fuse, represented as method == "analogue_fuse". If that optional backend is unavailable, the code falls back to the simple mean so the workflow remains executable, but the intended tutorial path uses analogue_fuse.

selected_graph_fusion_method <- if ("analogue_fuse" %in% selected_fusions$method) {
  "analogue_fuse"
} else {
  "mean"
}

selected_fusion <- selected_fusions[
  selected_fusions$method == selected_graph_fusion_method,
]

selected_fusion[, c("method", "type", "time_taken_seconds")]
#> # A tibble: 1 × 3
#>   method        type     time_taken_seconds
#>   <chr>         <chr>                 <dbl>
#> 1 analogue_fuse distance              0.001

At this point, we will use precise_graph() to build a graph we can then visualize using precise_viz(). We could also input the fused distance directly into an embedding view, but graph projection often returns a more pleasing and more interpretable visual object. The graph step asks a clean question: which relationships are strong enough to keep as edges?

precise_graph() is deliberately simple: it is a typed matrix-to-graph-matrix transformer. It returns an adjacency/similarity matrix, not a materialized graph object, which keeps the core workflow composable.

selected_graph <- precise_graph(
  selected_fusion,
  methods = "knn",
  params = list(
    knn = list(k = 15)
  ),
  parallel = FALSE,
  verbose = FALSE
)

selected_graph[, c("distance", "type", "graph_method", "input", "input_type")]
#> # A tibble: 1 × 5
#>   distance           type       graph_method input         input_type
#>   <chr>              <chr>      <chr>        <chr>         <chr>     
#> 1 analogue_fuse__knn similarity knn          analogue_fuse distance
data.frame(
  edges = selected_graph$meta[[1]]$edges,
  components = selected_graph$meta[[1]]$components
)
#>   edges components
#> 1   646          1

The graph output is typed as a similarity because an adjacency matrix encodes connection strength: one means an edge is present, zero means it is absent. The meta fields give a quick sanity check before visualization.

Before choosing one representation, it is worth looking across the graph construction and layout choices. The next block rebuilds the k-nearest-neighbour graph with k = 5, k = 10, and k = 15, then renders each graph with the igraph layouts exposed by precise_viz(). All panels use the same threejs renderer so the comparison is about graph density and layout algorithm, not about the plotting engine. Some igraph layouts are natively 2-D; when requested in this 3-D browser, PreciseDist renders them as flat z = 0 layouts and records that fact in the flat_3d column.

graph_k_values <- c(5, 10, 15)
graph_layout_values <- c(
  "mds", "kk", "fr", "drl", "nicely", "sphere",
  "gem", "graphopt", "lgl", "dh"
)

graph_layout_grid <- expand.grid(
  graph_k = graph_k_values,
  igraph_layout = graph_layout_values,
  stringsAsFactors = FALSE
)

graph_layout_panels <- lapply(seq_len(nrow(graph_layout_grid)), function(i) {
  graph_i <- precise_graph(
    selected_fusion,
    methods = "knn",
    params = list(
      knn = list(k = graph_layout_grid$graph_k[[i]])
    ),
    parallel = FALSE,
    verbose = FALSE
  )

  panel_i <- precise_viz(
    graph_i,
    views = "graph_layout",
    params = list(
      graph_layout = list(
        layout = graph_layout_grid$igraph_layout[[i]],
        dimensions = 3,
        engine = "threejs",
        render = "graph",
        seed = 1,
        color = cell_cycle_groups,
        colors = c(
          G1 = "#0017FF",
          G2M = "#000000",
          S = "#FF0017"
        ),
        show_labels = FALSE,
        size = 1
      )
    ),
    verbose = FALSE
  )

  panel_i$graph_k <- graph_layout_grid$graph_k[[i]]
  panel_i$igraph_layout <- graph_layout_grid$igraph_layout[[i]]
  panel_i$edges <- graph_i$meta[[1]]$edges
  panel_i$components <- graph_i$meta[[1]]$components
  panel_i$native_dimensions <- panel_i$meta[[1]]$native_layout_dimensions
  panel_i$flat_3d <- panel_i$meta[[1]]$flat_3d
  panel_i$panel_id <- sprintf(
    "k_%02d__%s",
    graph_layout_grid$graph_k[[i]],
    graph_layout_grid$igraph_layout[[i]]
  )
  panel_i
})

graph_layout_panels <- dplyr::bind_rows(graph_layout_panels)

graph_layout_panels[
  ,
  c(
    "panel_id",
    "graph_k",
    "igraph_layout",
    "edges",
    "components",
    "native_dimensions",
    "flat_3d"
  )
]
#> # A tibble: 30 × 7
#>    panel_id  graph_k igraph_layout edges components native_dimensions flat_3d
#>    <chr>       <dbl> <chr>         <int>      <int>             <int> <lgl>  
#>  1 k_05__mds       5 mds             236          1                 3 FALSE  
#>  2 k_10__mds      10 mds             449          1                 3 FALSE  
#>  3 k_15__mds      15 mds             646          1                 3 FALSE  
#>  4 k_05__kk        5 kk              236          1                 3 FALSE  
#>  5 k_10__kk       10 kk              449          1                 3 FALSE  
#>  6 k_15__kk       15 kk              646          1                 3 FALSE  
#>  7 k_05__fr        5 fr              236          1                 3 FALSE  
#>  8 k_10__fr       10 fr              449          1                 3 FALSE  
#>  9 k_15__fr       15 fr              646          1                 3 FALSE  
#> 10 k_05__drl       5 drl             236          1                 3 FALSE  
#> # ℹ 20 more rows

precise_trellis(
  graph_layout_panels,
  key_col = "panel_id",
  path = "trelliscope/cell_cycle_graph_layouts",
  name = "cell_cycle_graph_layouts",
  group = "common",
  desc = "3-D graph-layout comparison across k-nearest-neighbour graph densities and igraph layouts.",
  height = 550,
  width = 550
)

Layout is only one part of the problem. The graph matrix itself controls which relationships are available for the layout algorithm to arrange. The next Trelliscope view fixes the layout search to two useful 3-D force layouts and varies the graph construction instead: binary kNN, weighted kNN, mutual kNN, and quantile thresholding. It also tries two graph normalizations, chua and laplacian, after the initial graph has been constructed.

igraph’s DrL backend can reject a particular graph when its density grid is exceeded. Those specific DrL combinations are omitted; any other error remains fatal so the comparison cannot silently hide a workflow defect.

graph_knn_specs <- expand.grid(
  graph_constructor = c("knn", "weighted_knn", "mutual_knn"),
  graph_parameter = c(5, 10, 15),
  stringsAsFactors = FALSE
)
graph_knn_specs$graph_parameter_name <- "k"

graph_threshold_specs <- data.frame(
  graph_constructor = "quantile_threshold",
  graph_parameter = c(0.03, 0.06, 0.10),
  graph_parameter_name = "prop",
  stringsAsFactors = FALSE
)

graph_constructor_specs <- rbind(graph_knn_specs, graph_threshold_specs)

graph_construction_grid <- merge(
  graph_constructor_specs,
  expand.grid(
    graph_normalization = c("none", "chua", "laplacian"),
    igraph_layout = c("fr", "drl"),
    stringsAsFactors = FALSE
  ),
  all = TRUE
)

graph_construction_panels <- lapply(seq_len(nrow(graph_construction_grid)), function(i) {
  constructor <- graph_construction_grid$graph_constructor[[i]]
  parameter_name <- graph_construction_grid$graph_parameter_name[[i]]
  parameter_value <- graph_construction_grid$graph_parameter[[i]]
  normalization <- graph_construction_grid$graph_normalization[[i]]
  layout <- graph_construction_grid$igraph_layout[[i]]

  constructor_params <- if (identical(parameter_name, "k")) {
    list(k = parameter_value)
  } else {
    list(prop = parameter_value)
  }
  constructor_params <- stats::setNames(list(constructor_params), constructor)

  graph_i <- precise_graph(
    selected_fusion,
    methods = constructor,
    params = constructor_params,
    parallel = FALSE,
    verbose = FALSE
  )

  graph_for_viz <- if (identical(normalization, "none")) {
    graph_i
  } else {
    precise_graph(
      graph_i,
      methods = normalization,
      parallel = FALSE,
      verbose = FALSE
    )
  }

  panel_i <- tryCatch(
    precise_viz(
      graph_for_viz,
      views = "graph_layout",
      params = list(
        graph_layout = list(
          layout = layout,
          dimensions = 3,
          engine = "threejs",
          render = "graph",
          seed = 1,
          color = cell_cycle_groups,
          colors = c(
            G1 = "#0017FF",
            G2M = "#000000",
            S = "#FF0017"
          ),
          show_labels = FALSE,
          size = 1
        )
      ),
      verbose = FALSE
    ),
    error = function(e) {
      known_drl_limit <- identical(layout, "drl") &&
        grepl("Exceeded density grid in DrL", conditionMessage(e), fixed = TRUE)
      if (!known_drl_limit) stop(e)
      message(
        "Skipping ", constructor, "/", normalization, "/", layout,
        " because the layout backend failed: ", conditionMessage(e)
      )
      NULL
    }
  )

  if (is.null(panel_i)) return(NULL)

  panel_i$graph_constructor <- constructor
  panel_i$graph_parameter <- paste0(parameter_name, "=", parameter_value)
  panel_i$graph_normalization <- normalization
  panel_i$igraph_layout <- layout
  panel_i$edges <- graph_for_viz$meta[[1]]$edges
  panel_i$components <- graph_for_viz$meta[[1]]$components
  panel_i$weighted <- isTRUE(graph_i$meta[[1]]$weighted)
  panel_i$panel_id <- paste(
    constructor,
    paste0(parameter_name, "_", parameter_value),
    normalization,
    layout,
    sep = "__"
  )
  panel_i
})

graph_construction_panels <- dplyr::bind_rows(graph_construction_panels)

graph_construction_panels[
  ,
  c(
    "panel_id",
    "graph_constructor",
    "graph_parameter",
    "graph_normalization",
    "igraph_layout",
    "edges",
    "components",
    "weighted"
  )
]
#> # A tibble: 72 × 8
#>    panel_id  graph_constructor graph_parameter graph_normalization igraph_layout
#>    <chr>     <chr>             <chr>           <chr>               <chr>        
#>  1 knn__k_5… knn               k=5             none                fr           
#>  2 weighted… weighted_knn      k=5             none                fr           
#>  3 mutual_k… mutual_knn        k=5             none                fr           
#>  4 knn__k_1… knn               k=10            none                fr           
#>  5 weighted… weighted_knn      k=10            none                fr           
#>  6 mutual_k… mutual_knn        k=10            none                fr           
#>  7 knn__k_1… knn               k=15            none                fr           
#>  8 weighted… weighted_knn      k=15            none                fr           
#>  9 mutual_k… mutual_knn        k=15            none                fr           
#> 10 quantile… quantile_thresho… prop=0.03       none                fr           
#> # ℹ 62 more rows
#> # ℹ 3 more variables: edges <int>, components <int>, weighted <lgl>

precise_trellis(
  graph_construction_panels,
  key_col = "panel_id",
  path = "trelliscope/cell_cycle_graph_construction",
  name = "cell_cycle_graph_construction",
  group = "common",
  desc = "3-D graph comparison across graph constructors, graph normalizations, and force layouts.",
  height = 550,
  width = 550
)

Now feed it into precise_viz(). A 3-D force-directed graph is useful here because the extra dimension makes some neighbourhood structures easier to rotate and inspect. The choices are explicit: layout = "fr" requests the Fruchterman-Reingold layout, dimensions = 3 asks igraph for a 3-D layout, and engine = "threejs" renders the result as a browser widget.

selected_graph_layout_3d <- precise_viz(
  selected_graph,
  views = "graph_layout",
  params = list(
    graph_layout = list(
      layout = "fr",
      dimensions = 3,
      engine = "threejs",
      render = "graph",
      seed = 1,
      color = cell_cycle_groups,
      colors = c(
        G1 = "#0017FF",
        G2M = "#000000",
        S = "#FF0017"
      ),
      show_labels = FALSE,
      size = 1
    )
  ),
  verbose = FALSE
)

selected_graph_layout_3d$panel[[1]]

Move the 3-D graph with the mouse and use hover text to inspect individual observations. You can also zoom in and out with the scroll wheel. The edge matrix is the same selected_graph object as above; only the visual layout and renderer have changed.

This is the same discipline as before: the visualization is a view of a declared object. The matrix was built by selected metrics, fused by a named reducer, projected by a named graph constructor, and then rendered by a view that accepts similarity/adjacency input. Each step is inspectable.

As we can see, fusing four highly correlated distances that originally showed some structure does a good job of separating the held-back classes. What if we transform the input though?

For a stricter second graph representation, we can estimate a sparse graph from a similarity-like object with precise_graph(methods = "graphical_lasso") and the optional glasso backend. The question is whether the same broad organization survives regularized conditional-dependence estimation.

The k-nearest-neighbour graph asks which objects are close enough to retain as neighbours. A sparse precision graph asks which relationships survive after regularized conditional-dependence estimation. That question requires similarity input that is covariance/correlation-like, not a distance matrix. Here we use Similarity Network Fusion to fuse the selected distance matrices into one similarity, then pass that similarity into the graphical lasso graph constructor.

selected_similarity <- precise_fusion(
  selected_distances,
  methods = "snf",
  params = list(
    snf = list(
      K = 15,
      sigma = 0.5,
      t = 20
    )
  ),
  verbose = FALSE
)

selected_sparse_graph <- precise_graph(
  selected_similarity,
  methods = "graphical_lasso",
  params = list(
    graphical_lasso = list(
      rho = 0.022,
      edge_tol = 1e-8
    )
  ),
  verbose = FALSE
)

selected_sparse_graph[, c("distance", "type", "graph_method", "input", "input_type")]
#> # A tibble: 1 × 5
#>   distance             type       graph_method    input input_type
#>   <chr>                <chr>      <chr>           <chr> <chr>     
#> 1 snf__graphical_lasso similarity graphical_lasso snf   similarity
data.frame(
  edges = selected_sparse_graph$meta[[1]]$edges,
  components = selected_sparse_graph$meta[[1]]$components,
  precision_nonzero = selected_sparse_graph$meta[[1]]$precision_nonzero
)
#>   edges components precision_nonzero
#> 1   104         12               104

The penalty value is explicit because sparsity is an analytic choice, not a default to hide. Increasing rho removes more conditional-dependence edges; decreasing it keeps more. The goal is not to tune a clustering result. The goal is to inspect whether a stricter graph representation preserves a coherent structure.

selected_sparse_layout_3d <- precise_viz(
  selected_sparse_graph,
  views = "graph_layout",
  params = list(
    graph_layout = list(
      layout = "fr",
      dimensions = 3,
      engine = "threejs",
      render = "graph",
      seed = 1,
      color = cell_cycle_groups,
      colors = c(
        G1 = "#0017FF",
        G2M = "#000000",
        S = "#FF0017"
      ),
      show_labels = FALSE,
      size = 1
    )
  ),
  verbose = FALSE
)

selected_sparse_layout_3d$panel[[1]]

The sparse graph is not a proof of truth, but it is a useful second pressure test. If the same broad organization remains visible after moving from neighbourhood edges to regularized conditional edges, the structure has survived a more demanding representation.

What Comes Next

This main workflow stops after the first complete path through metric search, fusion, graph construction, and visualization. The follow-on vignettes take the same objects in three more focused directions:

Function Factory and Parameter Sweeps shows how to use precise_func_fact() when the question is not one metric, but a family of related metrics.
Metric Set Stability uses precise_stability() to ask whether the selected metric set is fragile.
Minkowski Parameter Neighborhood explores a local parameter sweep around Minkowski p = 0.5 and then uses graph-before-fusion.

What This Workflow Does Not Do

This vignette stops at typed relatedness objects, fusion, graph projection, and visual inspection. It does not choose a clustering algorithm. That is a deliberate boundary. Many tools can cluster a distance matrix, a similarity matrix, an adjacency matrix, or an embedding. PreciseDist’s responsibility is to help the user build and interrogate the representation those tools consume.

The result of this workflow is not an answer stamped “true.” It is a more honest object of analysis: a set of relatedness matrices, an agreement structure, a declared fusion, a graph projection, and visual evidence that can be argued with.

Brian Muchmore

2026-07-14