# [![](reference/figures/clav1.png)](https://github.com/jbryer/clav) Cluster Analysis Validation (`clav`)

[Slides from New York Open Statistical Programming Meetup talk given on
March 10,
2026](https://github.com/jbryer/clav/blob/master/slides/clav_nyhackr_2026.pdf)

[Slides from Joint Statistical Meeting (JSM)
2025](https://github.com/jbryer/clav/blob/master/slides/clav_jsm_2025.pdf)

The `clav` package provides utilities for conducting cluster (profile)
analysis with an emphasis on the validating the stability of the
profiles both within a given data set as well as across data sets.
Unlike supervised models where the known class is measured, validation
of unsupervised models where there is no known class can be challenging.
The approach implemented here attempts to compare the cluster results
across many random samples.

## Installation

You can install the development version of clav like so:

``` r
remotes::install_github('jbryer/clav')
```

## Development

The following commands are useful for working with the package source
locally.

``` r
# Prep the PISA data set. This will take a while to run the first time.
source('data-raw/data-prep-pisa-2015.R')
# Generate the package documentation
usethis::use_tidy_description()
devtools::document()
devtools::check_man()
# Install the package
devtools::install()
devtools::install(build_vignettes = TRUE)
# Run CRAN check
devtools::check(cran = TRUE)

# Build the pkgdown site
pkgdown::build_site()
```

## Clustering Functions

This package is designed to work with nearly any clustering methods.
However, there is no consistency in how clustering functions are
implemented with regard to parameters and returned values. The functions
`clav` package expect a format similar to the `kmeans` function.
Specifically, there are two parameters passed:

1.  A data frame with the variables used to do the clustering.
2.  The number of clusters, k.

The parameters are passed unnamed so the order of the parameters
matters. Additionally, the returned object should have an element
`cluster` that provides the cluster membership. The length of this
vector should be equal to the number of rows in the data frame.

For example, to use hierarchical clustering with the `clav` package, we
need to implement a wrapper (note that the `hclust2` function below is
included in the package).

``` r
hclust2 <- function(x, k, ...) {
    result <- stats::hclust(dist(x), ...)
    result$k <- k
    result$data <- x
    result$cluster <- stats::cutree(result, k = k)
    class(result) <- c('hclust2', class(result))
    return(result)
}
```

Note that the first parameter, `x` is the data frame of variables used
for the clustering. The
[`stats::hclust()`](https://rdrr.io/r/stats/hclust.html) function expect
a dissimilarity structure as it’s first parameter (the result of
[`stats::dist()`](https://rdrr.io/r/stats/dist.html)). Since `hclust`
determines cluster membership for all *n* - 1 possibilities, we use the
[`stats::cutree()`](https://rdrr.io/r/stats/cutree.html) function for
the desired k. Other parameters for `hclust` are passed using the `...`
(dots) operator. Lastly, we add the `cluster` vector and modify the
class so that we could implement S3 methods if necessary.

Additionally, an [S3
method](http://adv-r.had.co.nz/OO-essentials.html#s3) for `predict`
needs to be implemented. The following example implements
[`predict()`](https://rdrr.io/r/stats/predict.html) for
[`hclust2()`](reference/hclust2.md).

## Example

``` r
library(clav)
data(pisa2015, package = 'clav')

cluster_vars <- c('interest', 'enjoyment', 'motivation', 'efficacy', "belonging")
outcome_vars <- c('science_score', 'principals')

pisa_usa <- pisa2015 |> dplyr::filter(country == 'UNITED STATES')
```

Finding the optimal number of clusters.

``` r
optimal <- optimal_clusters(pisa_usa[,cluster_vars], max_k = 5)
optimal
#>   k  wss silhoutte  gap calinski_harabasz davies_bouldin rand_index
#> 1 1 9232        NA 0.83               NaN            NaN         NA
#> 2 2 6746      0.24 0.87              1708            1.6       0.50
#> 3 3 5869      0.20 0.86              1332            1.8       0.63
#> 4 4 5259      0.20 0.86              1160            1.7       0.74
#> 5 5 4658      0.20 0.86              1125            1.6       0.68
optimmal_plots <- plot(optimal, ncol = 2)
names(optimmal_plots)
#> [1] "davies_bouldin" "calinski"       "wss"            "silhouette"    
#> [5] "gap"            "rand"
# optimmal_plots[['wss']]
optimmal_plots
```

![](reference/figures/README-optimal-clusters-1.png)

Validating cluster profiles using random samples of 50%. Out-of-bag uses
the remaining 50% to predict cluster membership.

``` r
pisa_cv_random <- pisa_usa |> 
    dplyr::select(dplyr::all_of(cluster_vars)) |>
    clav::cluster_validation(
        n_clusters = 3,
        sample_size = 0.5 * nrow(pisa_usa),
        replace = FALSE,
        n_samples = 100,
        seed = 42
)
plot(pisa_cv_random)
```

![](reference/figures/README-random-validation-1.png)

Re-estimate the clusters using the OOB sample instead of predicting
using the in sample model.

``` r
pisa_cv_random2 <- pisa_usa |> 
    dplyr::select(dplyr::all_of(cluster_vars)) |>
    clav::cluster_validation(
        n_clusters = 3,
        oob_predict_fun = function(fit, newdata) {
            newfit <- stats::kmeans(newdata, 3)
            newfit$cluster
        },
        sample_size = 0.5 * nrow(pisa_usa),
        replace = FALSE,
        n_samples = 100,
        seed = 42
)
plot(pisa_cv_random2)
```

![](reference/figures/README-oob-reestimate-validation-1.png)

Bootstrap approach to validation.

``` r
pisa_cv_bootstrap <- pisa_usa |> 
    dplyr::select(dplyr::all_of(cluster_vars)) |>
    clav::cluster_validation(
        n_clusters = 3,
        sample_size = nrow(pisa_usa),
        replace = TRUE,
        n_samples = 100,
        seed = 42
)
plot(pisa_cv_bootstrap)
```

![](reference/figures/README-bootstrap-validation-1.png)

Using latent profile analysis for estimating clusters.

``` r
library(tidyLPA)
lpa <- pisa_usa |> 
    dplyr::select(dplyr::all_of(cluster_vars)) |>
    tidyLPA::estimate_profiles(3)
# lpa_predict <- predict(lpa, pisaUSA15[sample(nrow(pisaUSA15), 100),])
# lpa_estimates <- get_estimates(lpa)
lpa_data <- get_data(lpa)

plot_profiles(lpa)
clav::profile_plot(pisa_usa[,cluster_vars],
                   clusters = lpa_data$Class,
                   df_dep = pisa_usa[,c('science_score')])

lpa_cv_random <- cluster_validation(
    pisaUSA15,
    n_clusters = 3,
    cluster_fun = estimate_profiles,
    oob_predict_fun = function(fit, newdata) {
        estimate_profiles(newdata, n_clusters)
    },
    sample_size = 0.5 * nrow(pisaUSA15),
    replace = FALSE,
    n_samples = 50,
    seed = 42
)
plot(lpa_cv_random)
```

## Profile Plot

``` r
fit <- pisa_usa |> 
    dplyr::select(dplyr::all_of(cluster_vars)) |>
    stats::kmeans(centers = 3)
clav::profile_plot(pisa_usa[,cluster_vars],
                   clusters = fit$cluster,
                   df_dep = pisa_usa[,outcome_vars],
                   center_band = 0.33,
                   cluster_order = cluster_vars)
```

![](reference/figures/README-profile-plot-1.png)

# Package index

## All functions

- [`n_cluster_message()`](clav-shiny.md)
  [`cluster_method_input()`](clav-shiny.md)
  [`n_clusters_input()`](clav-shiny.md)
  [`cluster_variable_input()`](clav-shiny.md)
  [`n_cluster_plot_output()`](clav-shiny.md)
  [`cluster_size_bar_plot_output()`](clav-shiny.md)
  [`profile_plot_output()`](clav-shiny.md)
  [`cluster_pairs_plot_output()`](clav-shiny.md)
  [`bivariate_cluster_plot_output()`](clav-shiny.md)
  [`discriminant_projection_plot_output()`](clav-shiny.md)
  [`dependent_variable_input()`](clav-shiny.md)
  [`dependent_variable_plot_output()`](clav-shiny.md)
  [`dependent_variable_table_output()`](clav-shiny.md)
  [`dependent_null_hypothesis_output()`](clav-shiny.md)
  [`optimal_clusters_plot_output()`](clav-shiny.md)
  [`cluster_valdiation_plot_output()`](clav-shiny.md)
  [`cluster_validation_distribution_plot_output()`](clav-shiny.md)
  [`cluster_module()`](clav-shiny.md) : Output for printing status
  messages from the Shiny module.
- [`cluster_explainer_shiny()`](clav_shiny.md)
  [`cluster_shiny()`](clav_shiny.md) [`clav_shiny_ui()`](clav_shiny.md)
  [`clav_shiny_server()`](clav_shiny.md) : Run the Cluster Explainer
  Shiny Application
- [`cluster_agreement_fit()`](cluster_agreement_fit.md)
  [`plot(`*`<cluster_agreement_fit>`*`)`](cluster_agreement_fit.md)
  [`hist(`*`<cluster_agreement_fit>`*`)`](cluster_agreement_fit.md)
  [`print(`*`<cluster_agreement_fit>`*`)`](cluster_agreement_fit.md)
  [`summary(`*`<cluster_agreement_fit>`*`)`](cluster_agreement_fit.md)
  [`cluster_agreement()`](cluster_agreement_fit.md) : Cluster Agreement
  Fit
- [`cluster_overlap_fit()`](cluster_overlap_fit.md)
  [`plot(`*`<cluster.overlap>`*`)`](cluster_overlap_fit.md) : Cluster
  analysis fit statistic using bootstrap distributions
- [`cluster_validation()`](cluster_validation.md)
  [`plot(`*`<clustervalidation>`*`)`](cluster_validation.md)
  [`print(`*`<clustervalidation>`*`)`](cluster_validation.md)
  [`summary(`*`<clustervalidation>`*`)`](cluster_validation.md)
  [`plot(`*`<cv_summary>`*`)`](cluster_validation.md)
  [`fix_cluster_labels()`](cluster_validation.md)
  [`plot_distributions()`](cluster_validation.md) : Cluster profile
  validation
- [`daacs`](daacs.md) : Diagnostic Assessment and Achievement of College
  Skills
- [`describe_by()`](describe_by.md) : Basic summary statistics by group
- [`optimal_clusters()`](fitmeasures.md)
  [`plot(`*`<optimalclusters>`*`)`](fitmeasures.md)
  [`print(`*`<optimalclustersplot>`*`)`](fitmeasures.md)
  [`euclidean_distance()`](fitmeasures.md)
  [`get_centers()`](fitmeasures.md) [`wss()`](fitmeasures.md)
  [`silhouette_score()`](fitmeasures.md)
  [`calinski_harabasz()`](fitmeasures.md)
  [`davies_bouldin()`](fitmeasures.md) [`rand_index()`](fitmeasures.md)
  : Utility function used to determine the optimal number of clusters.
- [`hclust2()`](hclust2.md) [`predict(`*`<hclust2>`*`)`](hclust2.md) :
  Hierarchical Clustering
- [`kmeans_iterative()`](kmeans_iterative.md) : Function to perform
  K-means step by step and save iterations
- [`pisa2015`](pisa2015.md) : Student questionnaire and science
  performance data from the 2015 PISA
- [`predict(`*`<kmeans>`*`)`](predict.kmeans.md) : Predict cluster class
  for k-means cluster
- [`print_p_value()`](print_p_value.md) : Utility function to print
  p-values
- [`profile_plot()`](profile_plot.md) : Profile plot for cluster
  analysis.
- [`scale_this()`](scale_this.md) : Convert vector to z-scores

# Articles

### All vignettes

- [Bootstrap Overlap Fit for Determining the Optimal Number of
  Clusters](bootstrap_overlap_fit.md):
- [Cluster Analysis Validation](clav.md):
- [Cluster Agreement Fit Statistic](cluster_agreement_fit.md):
- [College Readiness Profiles](daacs.md):