Example Dataset

With the rCISSVAE package comes a sample dataset with predetermined clusters. We will use this for the sake of the tutorial.

The dataset, df_missing contains an index column as well as the following:

Age, Salary and ZipCode columns represent demographic data with no missingness and columns Y1t-Y5t represent biomarker data obtained at different timepoints t.

Explanation of run_cissvae() Parameters

The run_cissvae() function is a comprehensive wrapper for all basic steps in running the CISS-VAE model, including dataset preparation, optional clustering, and running the imputation model.

Running the imputation

To run the imputation model, first load your data. You can use the cluster_summary() function to visualize the missingness by cluster. The cluster_summary() function builds off of {gtsummary}.

  library(tidyverse)
  library(reticulate)
  library(rCISSVAE)
  library(gtsummary)


## Set correct virtualenv
reticulate::use_virtualenv("./cissvae_environment", required = TRUE)

## Load the data
data(df_missing)
data(clusters) ## actual cluster labels in clusters$clusters (other column is index)

cluster_summary(
  data = df_missing, 
  clusters = clusters$clusters, 
  include =setdiff(names(df_missing), "index"),
  statistic = list(
    all_continuous() ~ "{mean} ({sd})",
    all_categorical() ~ "{n} / {N}\n ({p}%)"), 
  missing = "always")

Then, plug your data and clusters into the run_cissvae() function.

## Run the imputation model. 
dat = run_cissvae(
  data = df_missing,
  index_col = "index",
  val_proportion = 0.1, ## pass a vector for different proportions by cluster
  columns_ignore = c("Age", "Salary", "ZipCode10001", "ZipCode20002", "ZipCode30003"), ## If there are columns in addition to the index you want to ignore when selecting validation set, list them here. In this case, we ignore the 'demographic' columns because we do not want to remove data from them for validation purposes. 
  clusters = clusters$clusters, ## we have precomputed cluster labels so we pass them here
  epochs = 500,
  return_silhouettes = FALSE,
  return_history = TRUE,  # Get detailed training history
  verbose = FALSE,
  return_model = TRUE, ## Allows for plotting model schematic
  device = "cpu",  # Explicit device selection
  layer_order_enc = c("unshared", "shared", "unshared"),
  layer_order_dec = c("shared", "unshared", "shared")
)

## Retrieve results
imputed_df <- dat$imputed
silhouette <- dat$silhouettes
training_history <- dat$history  # Detailed training progress

## Plot training progress
if (!is.null(training_history)) {
  plot(training_history$epoch, training_history$loss, 
       type = "l", main = "Training Loss Over Time", 
       xlab = "Epoch", ylab = "Loss")
}

plot_vae_architecture(model = dat$model, save_path = "test_plot_arch.png")

print(head(dat$imputed_dataset))

Cluster dataset:
 ClusterDataset(n_samples=8000, n_features=30, n_clusters=4)
  • Original missing: 61800 / 200000 (30.90%)
  • Validation held-out: 13783 (9.97% of non-missing)
  • .data shape:     (8000, 30)
  • .masks shape:    (8000, 30)
  • .val_data shape: (8000, 30)

#>   index       Age   Salary ZipCode10001 ZipCode20002 ZipCode30003         Y11
#> 0     0 11.044449 6.366204            0            1            0  -4.0495372
#> 1     1  9.727260 5.912558            1            0            0   0.5461677
#> 2     2 11.383020 6.636472            0            1            0  -1.2134339
#> 3     3 13.560905 5.896255            0            0            1 -10.6082144
#> 4     4  9.542490 6.128326            1            0            0   0.3575883
#> 5     5  9.542521 6.393217            1            0            0   4.7617960
#>         Y12         Y13        Y14        Y15        Y21       Y22       Y23
#> 0 -14.39339  -0.5147629 -14.369148 -17.564449 -10.736261 -18.53887 -35.77263
#> 1 -19.02272 -12.1895180  -7.722473   3.319988  -7.470250 -20.88160 -25.92436
#> 2 -19.03144 -20.3589058 -15.126495 -17.251385 -18.448421 -21.01385 -34.40086
#> 3 -22.24773  -7.1759834 -14.207619 -21.339748 -21.971752 -24.32133 -40.18794
#> 4 -16.48769 -11.3127708   7.535458  10.184475  -7.576005 -27.63498 -15.74972
#> 5 -18.96558 -12.2694435   5.667511  -9.094982  -4.120708 -29.16836 -25.07210
#>          Y24       Y25       Y31        Y32        Y33         Y34        Y35
#> 0 -28.098907 -30.24259 -1.627203  -6.557133 -16.769653 -10.6946259 -13.900185
#> 1 -17.231422 -18.69529  4.355482  -9.500225 -10.927032  -5.8868866  -6.088921
#> 2 -27.250603 -28.83982 -2.169048 -10.735180 -17.193771 -10.4940796 -12.286621
#> 3 -26.304344 -33.35109 -7.478534 -13.641594 -25.310104  -0.5022125 -15.429726
#> 4  -1.223648 -18.55935  8.103012 -13.563183  -9.832379  11.2815590  -2.962196
#> 5  -4.624435 -19.35896  7.710739  -8.304111 -10.381554   7.6654510  -4.992691
#>          Y41       Y42       Y43       Y44       Y45      Y51         Y52
#> 0 -0.9049605 -2.512167  3.527925 -3.694853 -5.680294 2.587588  1.84293461
#> 1  2.6245856 -4.218277 -5.776196 -1.379498 -2.329605 6.080512  0.04265404
#> 2  0.1033239 -4.733656 -7.215717 -3.350798 -6.895340 2.531148  0.31131554
#> 3 -2.8251922 -6.273870 -8.293898 -2.398297  3.204130 0.138556 -0.79631424
#> 4  3.6173897 -4.619453 -3.856878  7.120106 -1.495659 2.406390 -3.34118080
#> 5  2.1276770 -4.269014 -4.685389  5.764248 -2.990555 2.464692 -3.37204909
#>         Y53        Y54        Y55
#> 0 -4.681194 -2.2484055 -2.6790791
#> 1 -2.290062 -0.8873978  0.5625324
#> 2 -5.427431 -1.3301620 -2.3243809
#> 3 -5.729860 -1.6395874 -4.4457073
#> 4 -1.915190  6.5168343  0.1034508
#> 5 -2.237415  5.5294094 -1.0737858

Using Pre-computed Missingness Proportion Matrix

To create clusters based on proportion of missingness across all timepoints for a given feature, you can provide a pre-computed missingness proportion matrix (by using create_missingness_prop_matrix() or manually) directly to run_cissvae():

## Standardize df_missing column names to feature_timepoint format
colnames(df_missing) = c('index', 'Age', 'Salary', 'ZipCode10001', 'ZipCode20002', 'ZipCode30003', 'Y1_1', 'Y1_2', 'Y1_3', 'Y1_4', 'Y1_5', 'Y2_1', 'Y2_2', 'Y2_3', 'Y2_4', 'Y2_5', 'Y3_1', 'Y3_2', 'Y3_3', 'Y3_4', 'Y3_5', 'Y4_1', 'Y4_2', 'Y4_3', 'Y4_4', 'Y4_5', 'Y5_1', 'Y5_2', 'Y5_3', 'Y5_4', 'Y5_5')

# Create and examine missingness proportion matrix
prop_matrix <- create_missingness_prop_matrix(df_missing, 
index_col = "index", 
cols_ignore = c('Age', 'Salary', 'ZipCode10001', 'ZipCode20002', 'ZipCode30003'),
repeat_feature_names = c("Y1", "Y2", "Y3", "Y4", "Y5"))


cat("Missingness proportion matrix dimensions:\n")
cat(dim(prop_matrix), "\n")
cat("Sample of proportion matrix:\n")
print(head(prop_matrix[, 1:5]))

# Use proportion matrix with scaling for better clustering
advanced_result <- run_cissvae(
  data = df_missing,
  index_col = "index",
  clusters = NULL,  # Let function cluster using prop_matrix
  columns_ignore = c('Age', 'Salary', 'ZipCode10001', 'ZipCode20002', 'ZipCode30003'), 
  missingness_proportion_matrix = prop_matrix,
  scale_features = TRUE,  # Standardize features before clustering
  n_clusters = 4,
  leiden_resolution = 0.1,  
  epochs = 5,
  return_history = TRUE,
  return_silhouettes = TRUE,
  device = "cpu",
  verbose = FALSE,
  return_clusters = TRUE
)

print("Clustering quality:")
print(paste("Silhouette score:", round(advanced_result$silhouette_width, 3)))

## Plotting imputation loss by epoch 

ggplot2::ggplot(data = advanced_result$training_history, aes(x = epoch, y = imputation_error)) + geom_point() + labs(y = "Imputation Loss", x = "Epoch") + 
  theme_classic()

#> Missingness proportion matrix dimensions:
#> 8000 5
#> Sample of proportion matrix:
#>    Y1  Y2  Y3  Y4  Y5
#> 1 0.4 0.4 0.2 0.4 0.2
#> 2 0.4 0.2 0.2 0.2 0.2
#> 3 0.4 0.2 0.2 0.4 0.2
#> 4 0.4 0.2 0.4 0.4 0.2
#> 5 0.4 0.4 0.2 0.2 0.4
#> 6 0.2 0.2 0.4 0.2 0.4

#> Clustering quality:
#> Silhouette score: 0.656

Understanding Parameter Types

For hyperparameter optimization in autotune_cissvae(), parameters can be specified as:

• Fixed value: beta = 0.01 → parameter remains constant across trials
• Categorical choice: c(64, 128, 256) → Optuna selects from the provided options
• Float range: reticulate::tuple(1e-4, 1e-3) → Optuna suggests floats in the specified range

Layer Placement Strategies

The layer arrangement strategies control how shared and unshared layers are positioned:

• "at_end": Places shared layers at the end of the encoder or start of the decoder
• "at_start": Places shared layers at the start of the encoder or end of the decoder
  
• "alternating": Distributes shared layers evenly throughout the architecture
• "random": Uses random placement of shared layers (with reproducible seed)

library(tidyverse)
library(reticulate)
library(rCISSVAE)
reticulate::use_virtualenv("./cissvae_environment", required = TRUE)

data(df_missing)
data(clusters)

aut <- autotune_cissvae(
  data = df_missing,
  index_col = "index",
  clusters = clusters$clusters,
  save_model_path = NULL,
  save_search_space_path = NULL,
  n_trials = 3, ## Using low number of trials for demo
  study_name = "comprehensive_vae_autotune",
  device_preference = "cpu",
  show_progress = FALSE,  # Set true for Rich progress bars with training visualization
  optuna_dashboard_db = "sqlite:///optuna_study.db",  # Save results to database
  load_if_exists = FALSE, ## Set true to load and continue study if it exists
  seed = 42, 
  verbose = FALSE,
  
  # Search strategy options
  constant_layer_size = FALSE,     # Allow different sizes per layer
  evaluate_all_orders = FALSE,     # Sample layer arrangements efficiently
  max_exhaustive_orders = 100,     # Limit for exhaustive search
  
  ## Hyperparameter search space
  num_hidden_layers = c(2, 5),     # Try 2-5 hidden layers
  hidden_dims = c(64, 512),        # Layer sizes from 64 to 512
  latent_dim = c(10, 100),         # Latent dimension range
  latent_shared = c(TRUE, FALSE),
  output_shared = c(TRUE, FALSE),
  lr = 0.01,  # Learning rate range
  decay_factor = 0.99,
  beta = 0.01,  # KL weight range
  num_epochs = 500,                # Fixed epochs for demo
  batch_size = c(1000, 4000),     # Batch size options
  num_shared_encode = c(0, 1, 2, 3),
  num_shared_decode = c(0, 1, 2, 3),
  
  # Layer placement strategies - try different arrangements
  encoder_shared_placement = c("at_end", "at_start", "alternating", "random"),
  decoder_shared_placement = c("at_start", "at_end", "alternating", "random"),
  
  refit_patience = 2,        # Early stopping patience
  refit_loops = 100,                # Fixed refit loops
  epochs_per_loop = 100,   # Epochs per refit loop
  reset_lr_refit = c(TRUE, FALSE)
)

# Analyze results
imputed <- aut$imputed
best_model <- aut$model
study <- aut$study
results <- aut$results

# View best hyperparameters
print("Trial results:")
results %>% kable() %>%
  kable_styling(font_size=12)

# Plot model architecture
plot_vae_architecture(best_model, title = "Optimized CISSVAE Architecture")

#> [1] "Trial results:"

trial_number	imputation_error	num_hidden_layers	hidden_dim_0	hidden_dim_1	hidden_dim_2	hidden_dim_3	hidden_dim_4	latent_dim	latent_shared	output_shared	batch_size	num_shared_encode	num_shared_decode	encoder_shared_placement	decoder_shared_placement	layer_order_enc_used	layer_order_dec_used
0	29.62110	5	512	64	512	64	64	100	FALSE	FALSE	1000	1	1	alternating	at_start	S,U,U,U,U	S,U,U,U,U
1	33.00819	5	64	64	64	512	64	10	TRUE	TRUE	1000	1	0	alternating	alternating	S,U,U,U,U	U,U,U,U,U
2	58.28946	2	512	64	NaN	NaN	NaN	10	TRUE	FALSE	4000	0	3	alternating	random	U,U	S,S

Autotuned VAE Architecture

For more information on using optuna dashboard, see the Optuna Dashboard tutorial