A Complete Example

Pre-processing

Your data should be loaded into the R environment in the following format:

1. The data is in one or more data.frame class objects
2. The data is in wide form (one row per observation to cluster)
3. All data frames should have exactly one column that uniquely identifies each observation
4. All data should be complete (no missing values)

If you wish to use imputation to handle missingness in your data, you can take a look at the imputation vignette which outlines a basic workflow for meta clustering across multiple imputations of the same dataset.

The package comes with a few mock data frames based on real data from the Adolescent Brain Cognitive Development study:

• anxiety (anxiety scores from the CBCL)
• depress (depression scores from the CBCL)
• cort_t (cortical thicknesses)
• cort_sa (cortical surface areas in mm^2)
• subc_v (subcortical volumes in mm^3)
• income (household income on a 1-3 scale)
• pubertal (pubertal status on a 1-5 scale)

Here’s what the cortical thickness data looks like:

library(metasnf)

class(cort_t)
#> [1] "tbl_df"     "tbl"        "data.frame"

dim(cort_t)
#> [1] 188 152

str(cort_t[1:5, 1:5])
#> Classes 'tbl_df', 'tbl' and 'data.frame':    5 obs. of  5 variables:
#>  $ unique_id: chr  "NDAR_INV0567T2Y9" "NDAR_INV0GLZNC2W" "NDAR_INV0IZ157F8" "NDAR_INV0J4PYA5F" ...
#>  $ mrisdp_1 : num  2.6 2.62 2.62 2.6 2.53
#>  $ mrisdp_2 : num  2.49 2.85 2.29 2.67 2.76
#>  $ mrisdp_3 : num  2.8 2.78 2.53 2.68 2.83
#>  $ mrisdp_4 : num  2.95 2.85 2.96 2.94 2.99

cort_t[1:5, 1:5]
#>          unique_id mrisdp_1 mrisdp_2 mrisdp_3 mrisdp_4
#> 1 NDAR_INV0567T2Y9    2.601    2.487    2.801    2.954
#> 2 NDAR_INV0GLZNC2W    2.619    2.851    2.784    2.846
#> 3 NDAR_INV0IZ157F8    2.621    2.295    2.530    2.961
#> 4 NDAR_INV0J4PYA5F    2.599    2.670    2.676    2.938
#> 5 NDAR_INV0OYE291Q    2.526    2.761    2.829    2.986

The first column unique_id is the unique identifier (UID) for all observations in the data.

Here’s the household income data:

dim(income)
#> [1] 275   2

str(income[1:5, ])
#> Classes 'tbl_df', 'tbl' and 'data.frame':    5 obs. of  2 variables:
#>  $ unique_id       : chr  "NDAR_INV0567T2Y9" "NDAR_INV0GLZNC2W" "NDAR_INV0IZ157F8" "NDAR_INV0J4PYA5F" ...
#>  $ household_income: num  3 NA 1 2 1

income[1:5, ]
#>          unique_id household_income
#> 1 NDAR_INV0567T2Y9                3
#> 2 NDAR_INV0GLZNC2W               NA
#> 3 NDAR_INV0IZ157F8                1
#> 4 NDAR_INV0J4PYA5F                2
#> 5 NDAR_INV0OYE291Q                1

Putting everything in a list will help us get quicker summaries of all the data.

df_list <- list(
    anxiety,
    depress,
    cort_t,
    cort_sa,
    subc_v,
    income,
    pubertal
)

# The number of rows in each data frame:
lapply(df_list, dim)
#> [[1]]
#> [1] 275   2
#> 
#> [[2]]
#> [1] 275   2
#> 
#> [[3]]
#> [1] 188 152
#> 
#> [[4]]
#> [1] 188 152
#> 
#> [[5]]
#> [1] 174  31
#> 
#> [[6]]
#> [1] 275   2
#> 
#> [[7]]
#> [1] 275   2

# Whether or not each data frame has missing values:
lapply(df_list,
    function(x) {
        any(is.na(x))
    }
)
#> [[1]]
#> [1] TRUE
#> 
#> [[2]]
#> [1] TRUE
#> 
#> [[3]]
#> [1] FALSE
#> 
#> [[4]]
#> [1] FALSE
#> 
#> [[5]]
#> [1] FALSE
#> 
#> [[6]]
#> [1] TRUE
#> 
#> [[7]]
#> [1] TRUE

Some of the data has missing values and not all of the data frames have the same number of participants. SNF can only be run with complete data, so you’ll need to either use complete case analysis (removal of observations with any missing values) or impute the missing values to proceed with the clustering. As mentioned above, metasnf can be used to visualize changes in clustering results across different imputations of the data.

For now, we’ll just examine the simpler complete-case analysis approach by reducing our data frames to only common and complete observations. This can be made easier using the get_complete_uids function.

complete_uids <- get_complete_uids(df_list, uid = "unique_id")

print(length(complete_uids))
#> [1] 87

# Reducing data frames to only common observations with no missing data
anxiety <- anxiety[anxiety$"unique_id" %in% complete_uids, ]
depress <- depress[depress$"unique_id" %in% complete_uids, ]
cort_t <- cort_t[cort_t$"unique_id" %in% complete_uids, ]
cort_sa <- cort_sa[cort_sa$"unique_id" %in% complete_uids, ]
subc_v <- subc_v[subc_v$"unique_id" %in% complete_uids, ]
income <- income[income$"unique_id" %in% complete_uids, ]
pubertal <- pubertal[pubertal$"unique_id" %in% complete_uids, ]

Generating the data list

The data_list class is a structured list of data frames (like the one already created), but with some additional metadata about each data frame. It should only contain the data frames we want to directly use as inputs for the clustering. Let’s say we are working in a context where the anxiety and depression data are especially important outcomes and we want to know if we can find subtypes using the other data which still do a good job of separating out observations by their anxiety and depression scores.

We’ll start by initializing a data list that stores our input features.

# Note that you do not need to explicitly name every single named element
# (data = ..., name = ..., etc.)
input_dl <- data_list(
    list(
        data = cort_t,
        name = "cortical_thickness",
        domain = "neuroimaging",
        type = "continuous"
    ),
    list(
        data = cort_sa,
        name = "cortical_surface_area",
        domain = "neuroimaging",
        type = "continuous"
    ),
    list(
        data = subc_v,
        name = "subcortical_volume",
        domain = "neuroimaging",
        type = "continuous"
    ),
    list(
        data = income,
        name = "household_income",
        domain = "demographics",
        type = "continuous"
    ),
    list(
        data = pubertal,
        name = "pubertal_status",
        domain = "demographics",
        type = "continuous"
    ),
    uid = "unique_id"
)

This process removes any observations that did not have complete data across all provided input data frames. The structure of the data list is a nested list tracking the data, the name of the data frame, what domain (broader source of information) the data frame belongs to, and the type of feature stored in that data frame. Options for feature type include “continuous”, “discrete”, “ordinal”, “categorical”, and “mixed”.

The uid parameter is the name of the column in the data frames that uniquely identifies each observation. Upon data list creation, the UID will be converted to "uid" and the UIDs themselves will be prefixed with "_uid" for ease of management across other functions in the package.

We can get a summary of our constructed data list with the summary function:

summary(input_dl)
#>                    name       type       domain length width
#> 1    cortical_thickness continuous neuroimaging     87   151
#> 2 cortical_surface_area continuous neuroimaging     87   151
#> 3    subcortical_volume continuous neuroimaging     87    30
#> 4      household_income continuous demographics     87     1
#> 5       pubertal_status continuous demographics     87     1

Each input data frame now has the same 87 observations with complete data. The width refers to the number of features in each data frame (not including the UID column).

data_list now stores all the features we intend on using for the clustering. We’re interested in knowing if any of the clustering solutions we generate can distinguish children apart based on their anxiety and depression scores. To do this, we’ll also create a data list storing only features that we’ll use for evaluating our cluster solutions and not for the clustering itself. We’ll refer to this as the target data list.

target_dl <- data_list(
    list(anxiety, "anxiety", "behaviour", "ordinal"),
    list(depress, "depressed", "behaviour", "ordinal"),
    uid = "unique_id"
)

summary(target_dl)
#>        name    type    domain length width
#> 1   anxiety ordinal behaviour     87     1
#> 2 depressed ordinal behaviour     87     1

Note that it is not necessary to make use of a partition of input and out-of-model measures in this way. If you’d like to have no target data list and instead use every single feature of interest for the clustering, you can stick to just using one data list.

The settings data frame

The first part is the settings data frame:

my_sc$"settings_df"
#>                            1    2    3    4    5    6    7    8    9   10
#> SNF hyperparameters:
#> alpha                    0.5  0.4  0.3  0.3  0.5  0.4  0.7  0.8  0.3  0.6
#> k                         29   26   44   43   29   26   36   21   29   35  
#> t                         20   20   20   20   20   20   20   20   20   20  
#> SNF scheme:
#>                            2    1    2    1    2    2    2    3    1    3  
#> Clustering functions:
#>                            1    1    2    1    2    1    1    2    1    1  
#> Distance functions:
#> CNT                        1    1    1    1    1    1    1    1    1    1  
#> DSC                        1    1    1    1    1    1    1    1    1    1  
#> ORD                        1    1    1    1    1    1    1    1    1    1  
#> CAT                        1    1    1    1    1    1    1    1    1    1  
#> MIX                        1    1    1    1    1    1    1    1    1    1  
#> Component dropout:
#> cortical_thickness         ✔    ✔    ✔    ✔    ✔    ✔    ✔    ✔    ✔    ✔  
#> cortical_surface_area      ✖    ✔    ✖    ✔    ✔    ✔    ✔    ✔    ✔    ✔  
#> subcortical_volume         ✔    ✔    ✖    ✖    ✔    ✔    ✔    ✔    ✔    ✔  
#> household_income           ✖    ✔    ✔    ✔    ✔    ✔    ✔    ✖    ✔    ✖  
#> pubertal_status            ✔    ✔    ✔    ✔    ✔    ✔    ✔    ✔    ✔    ✖  
#> …and settings defined to create 10 more cluster solutions.

The settings data frame stores information about:

• SNF hyperparameters: alpha, k, and t hyperparameters that are directly used by the SNF algorithm.
• SNF scheme: the specific way in which input data gets collapsed into a final fused network (discussed further in the SNF schemes vignette).
• Clustering functions: Which clustering algorithm will be applied to the final fused network produced by SNF.
• Distance function: Which distance metric will be used when converting continuous, discrete, ordinal, categorical, or mixed type data into intermediate distance matrices (discussed further in the distance metrics vignette).
• Component dropout - Information on whether or not a particular input data frame in the data list will be in included or excluded from the corresponding SNF run (discussed further in the SNF config vignette).

The underlying structure of a settings data frame can be viewed by converting it from a settings_df class object into a regular data.frame class object:

head(as.data.frame(my_sc$"settings_df"))
#>   solution alpha  k  t snf_scheme clust_alg cnt_dist dsc_dist ord_dist cat_dist
#> 1        1   0.5 29 20          2         1        1        1        1        1
#> 2        2   0.4 26 20          1         1        1        1        1        1
#> 3        3   0.3 44 20          2         2        1        1        1        1
#> 4        4   0.3 43 20          1         1        1        1        1        1
#> 5        5   0.5 29 20          2         2        1        1        1        1
#> 6        6   0.4 26 20          2         1        1        1        1        1
#>   mix_dist inc_cortical_thickness inc_cortical_surface_area
#> 1        1                      1                         0
#> 2        1                      1                         1
#> 3        1                      1                         0
#> 4        1                      1                         1
#> 5        1                      1                         1
#> 6        1                      1                         1
#>   inc_subcortical_volume inc_household_income inc_pubertal_status
#> 1                      1                    0                   1
#> 2                      1                    1                   1
#> 3                      0                    1                   1
#> 4                      0                    1                   1
#> 5                      1                    1                   1
#> 6                      1                    1                   1

Each row of the settings data frame corresponds to a set of hyperparameters that can be used for generating a cluster solution. Without specifying any additional parameters, the snf_config function randomly populates these columns and ensures that no generated rows are identical. Further customization of the settings_df and the other parts of the SNF config will enable you to generate a broader space of cluster solutions from your data and ideally get you closer to finding a more useful solution for your context. More on settings_df customization can be found in the SNF config vignette.

Setting a seed prior to calling snf_config() ensures that the same random settings are generated each time we run the code.

Other parts of the SNF config

There are three other parts to an SNF config object: a distance functions list (dist_fns_list class object) that stores distance metric functions, a clustering functions list (clust_fns_list class object) that stores clustering algorithm functions, and a weights matrix (weights_matrix) class object that stores feature weights.

names(my_sc)
#> [1] "settings_df"    "dist_fns_list"  "clust_fns_list" "weights_matrix"

Customizing each of these parts is explained in more detail in separate vignettes:

• distance metrics
• clustering functions list
• weights matrix

For now, we’ll leave them as the default values generated during our snf_config call.

Calculating associations between cluster solutions and initial data

We start by running the extend_solutions function, which will calculate p-values representing the strength of the association between cluster membership (treated as a categorical feature) and any feature present in a provided data list and/or target_list.

extend_solutions also adds summary p-value measures (min, mean, and max) of any features present in the target list. This function can be quite slow depending on how large your data is; you can monitor its progress by setting verbose = TRUE.

ext_sol_df <- extend_solutions(
    mc_sol_df,
    dl = input_dl,
    target_dl = target_dl
)

ext_sol_df

#> 20 cluster solutions, 87 observations, and p-values for 336 features.
#> Cluster assignment columns:
#> solution nclust mc uid_NDAR_INV0567T2Y9 uid_NDAR_INV0J4PYA5F 
#>        1      5 E                     5                    2 
#>        2      3 D                     3                    3 
#>        3      9 C                     1                    7 
#>        4      2 D                     1                    2 
#>        5      8 C                     1                    6 
#>        6      2 B                     1                    2 
#>        7      4 C                     1                    4 
#>        8      5 E                     3                    4 
#>        9      3 D                     3                    3 
#>       10      4 A                     2                    4 
#> Association p-value columns:
#> solution mrisdp_1_pval mrisdp_2_pval mrisdp_3_pval mrisdp_4_pval
#>        1 4.0778e-01    8.7353e-01    1.8227e-01    6.4527e-01   
#>        2 2.9963e-01    3.7253e-01    6.0899e-01    2.9392e-02   
#>        3 2.4552e-01    2.2052e-01    1.1472e-01    5.2578e-01   
#>        4 2.0034e-01    8.0738e-01    3.5087e-01    4.7769e-01   
#>        5 1.8555e-01    4.5421e-01    1.0940e-01    4.5550e-01   
#>        6 7.9908e-01    3.7914e-02    4.8240e-01    7.0521e-01   
#>        7 4.7354e-01    3.7364e-01    1.0889e-01    8.1730e-01   
#>        8 4.0837e-01    6.6968e-01    1.6713e-01    6.7935e-01   
#>        9 2.5205e-01    3.4546e-01    3.6226e-01    1.4636e-02   
#>       10 5.2398e-01    4.1312e-01    2.9576e-01    1.9681e-01   
#> Summary p-value columns:
#> solution min_pval  mean_pval max_pval  
#>        1 7.263e-01 7.304e-01 7.344e-01 
#>        2 3.238e-01 4.941e-01 6.644e-01 
#>        3 3.961e-01 5.389e-01 6.817e-01 
#>        4 4.844e-01 6.204e-01 7.564e-01 
#>        5 2.257e-01 5.076e-01 7.895e-01 
#>        6 1.325e-01 1.946e-01 2.566e-01 
#>        7 1.230e-01 2.031e-01 2.832e-01 
#>        8 7.112e-01 7.380e-01 7.648e-01 
#>        9 5.110e-01 6.856e-01 8.602e-01 
#>       10 1.034e-02 1.833e-02 2.631e-02 
#> Summaries calculated from 2 features. Use `summary_features(x)` to see them.
#> 10 solutions and 332 features not shown.
#> Use `print(n = ...)` to change the number of rows printed.
#> Use `t()` to view compact cluster solution format.

The difference between the data list passed into the dl parameter and the one passed into the target_dl parameter is that the target_dl features are the only ones used for generating the p-value summary columns.

# No features are used to calculate summary p-values
ext_sol_df_no_summaries <- extend_solutions(
    mc_sol_df,
    dl = c(input_dl, target_dl)
)

# Every features are used to calculate summary p-values
ext_sol_df_all_summaries <- extend_solutions(
    mc_sol_df,
    target_dl = c(input_dl, target_dl)
)

Visualizing feature associations with meta clustering results

This ext_sol_df we created can be converted into a standard data frame and passed into the meta_cluster_heatmap function to easily visualize the level of separation on each of our features for each of our cluster solutions.

annotated_ari_hm <- meta_cluster_heatmap(
    sol_aris,
    order = meta_cluster_order,
    split_vector = split_vec,
    data = as.data.frame(ext_sol_df),
    top_hm = list(
        "Depression p-value" = "cbcl_depress_r_pval",
        "Anxiety p-value" = "cbcl_anxiety_r_pval",
        "Overall outcomes p-value" = "mean_pval"
    ),
    bottom_bar = list(
        "Number of Clusters" = "nclust"
    ),
    annotation_colours = list(
        "Depression p-value" = colour_scale(
            ext_sol_df$"cbcl_depress_r_pval",
            min_colour = "purple",
            max_colour = "black"
        ),
        "Anxiety p-value" = colour_scale(
            ext_sol_df$"cbcl_anxiety_r_pval",
            min_colour = "green",
            max_colour = "black"
        ),
        "Overall outcomes p-value" = colour_scale(
            ext_sol_df$"mean_pval",
            min_colour = "lightblue",
            max_colour = "black"
        )
    )
)

The meta_cluster_heatmap function wraps around a more generic similarity_matrix_heatmap function in this package, which itself is a wrapper around ComplexHeatmap::Heatmap(). Consequently, any parameter that can be used with ComplexHeatmap::Heatmap() package can be used here. This also makes the documentation for the ComplexHeatmap package one of the best places to learn more about what can be done for further customizing these heatmaps.

The data for the annotations don’t necessarily need to come from functions in metasnf. For example, to highlight the solutions where two very specific observations happened to cluster together, we can easily add that information in as another annotation. By converting an ext_solutions_df class object into a data frame with the parameter keep_attributes = TRUE, the resulting data frame will retain the p-value information from the ext_solutions_df side of things, cluster assignment columns from the solutions_df side of things, and settings_df and weights_matrix columns from the embedded original SNF config.

Then, we can use dplyr::mutate() to create a new column flagging when two observations co-cluster.

annotation_data <- ext_sol_df |>
    as.data.frame(keep_attributes = TRUE) |>
    dplyr::mutate(
        key_subjects_cluster_together = dplyr::case_when(
            uid_NDAR_INVLF3TNDUZ == uid_NDAR_INVLDQH8ATK ~ TRUE,
            TRUE ~ FALSE
        )
    )

annotated_ari_hm2 <- meta_cluster_heatmap(
    sol_aris,
    order = meta_cluster_order,
    split_vector = split_vec,
    data = annotation_data,
    top_hm = list(
        "Depression p-value" = "cbcl_depress_r_pval",
        "Anxiety p-value" = "cbcl_anxiety_r_pval",
        "Key Subjects Clustered Together" = "key_subjects_cluster_together"
    ),
    bottom_bar = list(
        "Number of Clusters" = "nclust"
    ),
    annotation_colours = list(
        "Depression p-value" = colour_scale(
            ext_sol_df$"cbcl_depress_r_pval",
            min_colour = "purple",
            max_colour = "black"
        ),
        "Anxiety p-value" = colour_scale(
            ext_sol_df$"cbcl_anxiety_r_pval",
            min_colour = "purple",
            max_colour = "black"
        ),
        "Key Subjects Clustered Together" = c(
            "TRUE" = "blue",
            "FALSE" = "black"
        )
    )
)

Characterizing individual solutions representative of each meta cluster

Now that we’ve visually delineated our meta clusters, we can get a quick summary of what sort of separation exists across our input and held out features by taking a closer look at one representative cluster solution from each meta cluster.

This can be achieved by the get_representative_solutions function, which extracts one cluster solution per meta cluster based on having the highest average ARI with the other solutions in that meta cluster.

rep_solutions <- get_representative_solutions(sol_aris, ext_sol_df)

mc_manhattan <- mc_manhattan_plot(
    rep_solutions,
    dl = input_dl,
    target_dl = target_dl,
    hide_x_labels = TRUE,
    point_size = 2,
    text_size = 12,
    threshold = 0.05,
    neg_log_pval_thresh = 5
)

neg_log_pval_thresh sets a threshold for the negative log of the p-values to be displayed. At a value of 5, any p-value that is smaller than 1e-5 will be truncated to 1e-5.

The plot as it is is a bit unwieldy plot given how many neuroimaging ROIs are present. Let’s take out the cortical thickness and surface area measures to make the plot a little clearer.

We’ll also be able to read some feature measures more clearly when we dial the number of features plotted back a bit, as well.

rep_solutions_no_cort <- dplyr::select(rep_solutions, -dplyr::contains("mrisdp"))

mc_manhattan2 <- mc_manhattan_plot(
    ext_sol_df = rep_solutions_no_cort,
    dl = input_dl,
    target_dl = target_dl,
    point_size = 4,
    threshold = 0.01,
    text_size = 12,
    domain_colours = c(
        "neuroimaging" = "cadetblue",
        "demographics" = "purple",
        "behaviour" = "firebrick"
    )
)

Note that the Manhattan plot automatically uses a vertical line to separate features in the data_list argument from those in the target_list. Vertical boundaries can be controlled with the xints parameter.

Relating results to metasnf hyperparameters

If you see something interesting in your heatmap, you may be curious to know how that corresponds to the settings that were in your settings data frame.

You could certainly stack on each setting to the ARI heatmap as annotations, but that may be a bit cumbersome given how many settings there are. Another option is to start by taking a first look at entire settings_df, sorted by the meta cluster results, through the settings_df_heatmap function. The settings can be passed in through either an snf_config or settings_df class object.

config_hm <- config_heatmap(
    my_sc,
    order = meta_cluster_order,
    hide_fixed = TRUE
)

This heatmap rescales all the columns in the settings_df to have a maximum value of 1. The purpose of re-ordering the settings data frame in this way is to see if any associations exist between certain settings values and pairwise cluster solution similarities. If there are any particular important settings, you can simply add them into your adjusted rand index heatmap annotations. Recall that the solutions_df class (and, by extension, the ext_solutions_df class that contains a solutions_df object as an attribute) contain the SNF config as an attribute, so no further data manipulation is needed to add a setting as a heatmap annotation.

Give it a try with the code below:

annotation_data <- ext_sol_df |>
    as.data.frame(keep_attributes = TRUE) |>
    dplyr::mutate(
        key_subjects_cluster_together = dplyr::case_when(
            uid_NDAR_INVLF3TNDUZ == uid_NDAR_INVLDQH8ATK ~ TRUE,
            TRUE ~ FALSE
        )
    )

annotation_data$"clust_alg" <- as.factor(annotation_data$"clust_alg")

annotated_ari_hm3 <- meta_cluster_heatmap(
    sol_aris,
    order = meta_cluster_order,
    split_vector = split_vec,
    data = annotation_data,
    top_hm = list(
        "Depression p-value" = "cbcl_depress_r_pval",
        "Anxiety p-value" = "cbcl_anxiety_r_pval",
        "Key Subjects Clustered Together" = "key_subjects_cluster_together"
    ),
    left_hm = list(
        "Clustering Algorithm" = "clust_alg" # from the original settings
    ),
    bottom_bar = list(
        "Number of Clusters" = "nclust" # also from the original settings
    ),
    annotation_colours = list(
        "Depression p-value" = colour_scale(
            ext_sol_df$"cbcl_depress_r_pval",
            min_colour = "purple",
            max_colour = "black"
        ),
        "Anxiety p-value" = colour_scale(
            ext_sol_df$"cbcl_anxiety_r_pval",
            min_colour = "purple",
            max_colour = "black"
        ),
        "Key Subjects Clustered Together" = c(
            "TRUE" = "blue",
            "FALSE" = "black"
        )
    )
)

annotated_ari_hm3

Quality measures

Quality metrics are another useful heuristic for the goodness of a cluster that don’t require any contextualization of results in the domain they may be used in. metasnf enables measures of silhouette scores, Dunn indices, and Davies-Bouldin indices. To calculate these values, we’ll need not only the cluster results but also the final fused network (the similarity matrices produced by SNF) that the clusters came from. These similarity matrices can be collected from the batch_snf using the return_sim_mats parameter:

sol_df <- batch_snf(dl = input_dl, sc = my_sc, return_sim_mats = TRUE)

Now, the solutions data frame contains a non-empty list of similarity matrices as one of its attributes. These similarity matrices are the final SNF-derived networks from each SNF run. Using that attribute, we can calculate the above mentioned quality metrics:

silhouette_scores <- calculate_silhouettes(sol_df)
dunn_indices <- calculate_dunn_indices(sol_df)
db_indices <- calculate_db_indices(sol_df)

The first function is a wrapper around cluster::silhouette while the second and third come from the clv package. clv isn’t set as a mandatory part of the installation, so you’ll ned to install it yourself to calculate these two metrics.

The original documentation on these functions can be helpful for interpreting and working with them:

1. cluster::silhouette documentation
2. clv::clv.Dunn documentation
3. clv::clv.Davies.Bouldin documentation

Stability measures

metasnf offers tools to evaluate two different measures of stability:

1. Pairwise adjusted Rand indices (across resamplings of the clustering, on average, how similar was every pair of solutions according to the adjusted Rand index?)
2. Fraction clustered together (what is the average fraction of times that observations who clustered together in the full results clustered together in resampled results?)

You can learn more about running stability calculations in the stability and coclustering vignette.

Evaluating separation across “target features” of importance

If you can specify a metric or objective function that may tell you how useful a clustering solution will be for your purposes in advance, that makes the cluster selection process much less arbitrary.

There are many ways to go about doing this, but this package offers one way through a target_list. The target_list contains data frames what we can examine our clustering results over through linear regression (continuous data), ordinal regression (ordinal data), or the Chi-squared test (categorical data).

Just like when generating the initial data_list, we need to specify the name of the column in the provided data frames that is originally being used to uniquely identify the different observations from each other with the uid parameter.

We will next extend our sol_df with p-values from regressing the target_list features onto our generated clusters.

ext_sol_df <- extend_solutions(sol_df, target_dl)

There is a heatmap for visualizing this too:

pval_hm <- pval_heatmap(ext_sol_df, order = meta_cluster_order)

pval_hm

These p-values hold no real meaning for the traditional hypothesis-testing context, but they are reasonable proxies of the magnitude of the effect size / separation of the clusters across the features in question. Here, they are just a tool to find clustering solutions that are well-separated according to the outcome measures you’ve specified. Finding a cluster solution like this is similar to a supervised learning approach, but where the optimization method is just random sampling. The risk for overfitting your data with this approach is considerable, so make sure you have some rigorous external validation before reporting your findings.

We recommend using label propagation (provided by the SNFtool package in the groupPredict function) for validation: take the top clustering solutions found in some training data, assign predicted clusters to some held out test observations, and then characterize those test observations to see how well the clustering solution seemed to have worked.