LUCID's Three Model Architectures: Early, Parallel, and Serial -- Binary Outcome (HELIX Example)

In LUCIDus: LUCID with Multiple Omics Data

1) Tutorial Goal and Scope

This tutorial is designed as a hands-on, end-to-end guide for fitting LUCID models on HELIX-style multi-omics data, with a binary outcome. Its companion vignette, lucid_3models_normal_outcome.Rmd, walks through the identical pipeline step by step for a continuous outcome, including the lucid() tuning wrapper, missing-data diagnostics, and prediction/g-computation -- material this vignette does not repeat. Visualization (Sankey diagram, cluster omics profiles) is shown here too, since neither depends on the outcome family.

What you will learn here:

• How to prepare common inputs (G, Z, Y_binary, CoG, CoY) from one dataset.
• How to fit and interpret, for a binary outcome:
• LUCID early integration
• LUCID in parallel
• LUCID in serial (all-early stages and mixed parallel+early stages)
• How to run the practical two-step strategy:
• Penalized screening fit (Rho_* > 0) for selection.
• Zero-penalty refit (Rho_* = 0) for stable inference + bootstrap CI.
• How predict_lucid()'s response argument controls whether a binary outcome comes back as class labels or as probabilities.

Important runtime note:

• Sections 7-10 fix K and do the screen-then-refit by hand, mirroring the continuous-outcome companion vignette's early/parallel/serial sections.
• We keep K fixed and use small bootstrap R so the vignette remains reproducible and fast.

2) Data Objects and Statistical Roles

The HELIX example data (simulated_HELIX_data.rda) provides:

• phenotype: exposure/covariate/outcome-style variables.
• three omics layers: methylome, transcriptome, miRNA.

Model inputs used throughout:

• G: exposure matrix (main predictors of latent cluster assignment). Built with a known signal structure -- three exposures correlated with the outcome and six pure-noise exposures -- so that feature selection can be checked against the truth rather than taken on trust. The outcome itself is deliberately kept out of G.
• Z: omics input (matrix for early, list for parallel, nested list for serial mixed).
• Y: continuous outcome vector; Y_binary is its median split, used to demonstrate the binary outcome model on the same subjects.
• CoG: covariates for G -> X model.
• CoY: covariates for X -> Y model.

3) Hyperparameter Guide (Practical)

| Hyperparameter | Meaning | Tutorial choice and rationale | |---|---|---| | K | number of latent clusters | fixed to small values for speed and interpretability | | Rho_G | penalty on G -> X coefficients | positive in screening fit, zero in inference refit | | Rho_Z_Mu | penalty on cluster-specific omics means | positive in screening fit, zero in inference refit | | Rho_Z_Cov | penalty on omics covariance matrices | positive in screening fit, zero in inference refit | | max_itr, max_tot.itr, tol | EM controls | modest values to balance speed/stability | | family | outcome model family | normal for this tutorial | | seed | reproducibility | fixed before each fit/bootstrap |

4) Setup and Source Package Code

# Keep knitting on error, so that one failing step reports itself and the rest of
# the tutorial still runs. The status table in section 13 records what happened.
knitr::opts_chunk$set(error = TRUE)

# Lightweight registry so the document can verify itself rather than relying on
# the reader to notice a missing output.
.reg <- new.env(parent = emptyenv()); .reg$rows <- list()
check_obj <- function(name, expected_class = NULL, section = "") {
  ok <- exists(name, envir = globalenv())
  cls <- if (ok) class(get(name, envir = globalenv()))[1] else NA_character_
  status <- if (!ok) "MISSING"
            else if (!is.null(expected_class) && !identical(cls, expected_class)) "unexpected class"
            else "ok"
  .reg$rows[[length(.reg$rows) + 1L]] <-
    data.frame(section = section, object = name, class = cls,
               status = status, stringsAsFactors = FALSE)
  invisible(NULL)
}

library(LUCIDus)

# The HELIX simulation object bundled with the package.
data(simulated_HELIX_data)

5) Build Modeling Inputs (With Missingness Injection)

This chunk creates a compact tutorial dataset and deliberately injects both:

• listwise missingness (whole row missing for a layer), and
• sporadic missingness (some cells missing),

so we can observe missing-data handling in summaries.

# Use a smaller subset for vignette speed while preserving model behavior.
idx <- 1:90
ph <- simulated_HELIX_data$phenotype[idx, ]
n <- nrow(ph)

set.seed(2026)

# ---------------------------------------------------------------------------
# A tutorial dataset with a KNOWN answer.
#
# The HELIX omics matrices are real simulated data with their own structure, and
# the exposures shipped with them have no relationship to it. That is fine for
# demonstrating that code runs, but it makes feature selection impossible to
# judge: there is no right answer to compare against. So we plant one.
#
# The generating story, which is the DAG LUCID assumes:
#
#     causal exposures  ->  latent subgroup  ->  omics profile
#                                            ->  outcome
#
# Three exposures carry the subgroup signal with graded strength; six are pure
# noise. Half the features of each omics layer are shifted by subgroup
# membership; the rest are left as they came. Selection therefore has an
# unambiguous target, and the tutorial can check its answer instead of asserting
# it.
# ---------------------------------------------------------------------------

# The true latent subgroup. Retained so every selection claim below can be
# checked against it.
x_true <- rbinom(n, 1, 0.5)

# Exposures. g_causal_* predict subgroup membership; g_noise_* do not.
# Effect sizes are deliberately moderate. Stronger exposures make selection
# look better but drive the G -> X model to saturation, where every subject sits
# at posterior probability 1 and no counterfactual shift can move anything --
# which would make the g-computation demonstration in the continuous-outcome
# companion vignette vacuous.
G <- cbind(
  g_causal_1 =  1.0 * (x_true - 0.5) + rnorm(n, sd = 0.8),   # strongest
  g_causal_2 = -0.8 * (x_true - 0.5) + rnorm(n, sd = 0.8),   # moderate, negative
  g_causal_3 =  0.6 * (x_true - 0.5) + rnorm(n, sd = 0.8),   # weakest
  g_noise_1 = rnorm(n), g_noise_2 = rnorm(n), g_noise_3 = rnorm(n),
  g_noise_4 = rnorm(n), g_noise_5 = rnorm(n), g_noise_6 = rnorm(n)
)
G <- as.matrix(scale(G))

causal_exposures <- c("g_causal_1", "g_causal_2", "g_causal_3")

# Exposure penalty used throughout. The continuous-outcome companion vignette
# shows what this value recovers and sweeps it, along with the omics penalty,
# separately.
RHO_G <- 0.05

# Covariates for G->X (CoG) and X->Y (CoY).
# Here we use age-related and sex covariates from phenotype.
CoG <- cbind(
  hs_child_age_yrs_None = as.numeric(ph$hs_child_age_yrs_None),
  sex_male = as.numeric(ph$e3_sex_None == "male")
)
CoY <- CoG

# Two outcomes on the SAME subjects, so the normal and binary results below are
# directly comparable: the only thing that changes between them is the outcome
# model, not the sample, the omics, or the injected missingness.
#
# Continuous outcome: the real CK-18 measurement, plus a subgroup effect so the
# cluster -> outcome arm of the model has something to estimate.
Y <- as.numeric(ph$ck18_scaled) + 1.2 * x_true

# Binary outcome: median split. The median is used rather than a higher
# threshold because it splits these 90 subjects 45/45, and a balanced outcome
# gives the K = 2 outcome model the most to work with at this sample size.
Y_binary <- as.integer(Y > median(Y))
cat("binary outcome balance:\n"); print(table(Y_binary))

# Construct three omics layers and standardize each, then plant the subgroup
# signal in the first three features of every layer. The remaining seven per
# layer are left as they came and act as omics noise.
meth <- scale(simulated_HELIX_data$methylome[idx, 1:10, drop = FALSE])
tran <- scale(simulated_HELIX_data$transcriptome[idx, 1:10, drop = FALSE])
mir  <- scale(simulated_HELIX_data$miRNA[idx, 1:10, drop = FALSE])

signal_features <- 1:3
omics_shift <- 3.0
meth[, signal_features] <- meth[, signal_features] + omics_shift * x_true
tran[, signal_features] <- tran[, signal_features] - omics_shift * x_true
mir[,  signal_features] <- mir[,  signal_features] + omics_shift * x_true

# Column positions of the signal features once the layers are stacked for the
# early model, so selection can be scored against them later.
signal_cols_early <- c(signal_features,
                       ncol(meth) + signal_features,
                       ncol(meth) + ncol(tran) + signal_features)

# Early model uses one combined Z matrix.
Z_early <- cbind(meth, tran, mir)

# Parallel model uses list-of-layers.
Z_parallel <- list(methylome = meth, transcriptome = tran, miRNA = mir)

# Inject listwise + sporadic missingness for demonstration.
Z_early_miss <- Z_early
Z_early_miss[1, ] <- NA      # listwise row
Z_early_miss[2:4, 1] <- NA   # sporadic block
Z_early_miss[5, 3] <- NA     # sporadic cell

Z_parallel_miss <- Z_parallel
Z_parallel_miss[[1]][1, ] <- NA  # listwise in layer 1
Z_parallel_miss[[2]][2, 2] <- NA # sporadic in layer 2
Z_parallel_miss[[3]][3, 1] <- NA # sporadic in layer 3

# Quick structural sanity check.
str(list(
  G = G,
  CoG = CoG,
  CoY = CoY,
  Y = Y,
  Z_early = Z_early_miss,
  Z_parallel = Z_parallel_miss
), max.level = 1)

6) Helper Functions for Selected-Feature Refit

These helpers implement a robust refit pipeline:

1. Read feature-selection indicators from penalized fit.
2. Build selected-only G/Z inputs.
3. Refit with all penalties set to zero for bootstrap inference.

# get_selected_G()/get_selected_Z() (from the package itself) already return a
# well-shaped, aligned logical mask straight from the fitted object -- no
# length mismatch is possible, since they derive it from the model's own
# recorded fields. The one thing left for a tutorial to decide is what to do
# if a penalty happened to deselect EVERY feature: refitting on zero columns
# would fail, so this keeps everything instead in that one edge case.
keep_or_all <- function(mask) if (any(mask, na.rm = TRUE)) mask else rep(TRUE, length(mask))

# Build selected-only inputs for early model.
prepare_early_selected_inputs <- function(fit_pen, G, Z) {
  list(
    G = as.matrix(G[, keep_or_all(get_selected_G(fit_pen)), drop = FALSE]),
    Z = as.matrix(Z[, keep_or_all(get_selected_Z(fit_pen)), drop = FALSE])
  )
}

# Build selected-only inputs for parallel model.
prepare_parallel_selected_inputs <- function(fit_pen, G, Z) {
  keep_g <- keep_or_all(get_selected_G(fit_pen))
  Z_sel <- lapply(seq_along(Z), function(i) {
    zi <- as.matrix(Z[[i]])
    zi[, keep_or_all(get_selected_Z(fit_pen, layer = i)), drop = FALSE]
  })
  names(Z_sel) <- names(Z)
  list(
    G = as.matrix(G[, keep_g, drop = FALSE]),
    Z = Z_sel
  )
}

# Serial stage>1 uses latent-cluster-derived "G" internally.
# We therefore subset stage-1 original G and each stage's Z where applicable.
prepare_serial_selected_inputs <- function(fit_pen, G, Z) {
  G_refit <- as.matrix(G)
  keep_g1 <- get_selected_G(fit_pen)
  if (length(keep_g1) == ncol(G_refit)) {
    G_refit <- G_refit[, keep_or_all(keep_g1), drop = FALSE]
  }

  selected_z <- get_selected_Z(fit_pen)
  Z_refit <- Z
  for (i in seq_along(fit_pen$submodel)) {
    sm <- fit_pen$submodel[[i]]
    if (inherits(sm, "early_lucid")) {
      zi <- as.matrix(Z_refit[[i]])
      Z_refit[[i]] <- zi[, keep_or_all(selected_z[[i]]), drop = FALSE]
    } else if (inherits(sm, "lucid_parallel")) {
      zi_list <- Z_refit[[i]]
      for (j in seq_along(zi_list)) {
        zij <- as.matrix(zi_list[[j]])
        zi_list[[j]] <- zij[, keep_or_all(selected_z[[i]][[j]]), drop = FALSE]
      }
      Z_refit[[i]] <- zi_list
    }
  }

  list(G = G_refit, Z = Z_refit)
}

# Zero-penalty refit, for any model type.
#
# The three model types previously had three byte-identical wrappers differing
# only in `lucid_model` and whether `useY` was forwarded; they are one function
# here. Everything about the model -- family, K, initialization, EM controls --
# is carried over from the screening fit, so the ONLY difference between the
# screening fit and this one is that the penalties are zero. That is what makes
# the refit estimates unshrunk and therefore suitable for bootstrap inference.
refit_selected <- function(model_type, fit_pen, inputs, Y,
                           CoG = NULL, CoY = NULL, seed = 1, verbose = FALSE) {
  args <- list(
    lucid_model = model_type,
    G = inputs$G,
    Z = inputs$Z,
    Y = Y,
    CoG = CoG,
    CoY = CoY,
    family = fit_pen$family,
    K = fit_pen$K,
    init_omic.data.model = fit_pen$init_omic.data.model,
    init_impute = fit_pen$init_impute,
    init_par = fit_pen$init_par,
    Rho_G = 0,
    Rho_Z_Mu = 0,
    Rho_Z_Cov = 0,
    max_itr = fit_pen$em_control$max_itr,
    max_tot.itr = fit_pen$em_control$max_tot.itr,
    tol = fit_pen$em_control$tol,
    seed = seed,
    verbose = verbose
  )
  # Every fitted class records useY, so it is carried over for all three model
  # types. The original three wrappers omitted it on the early path, which meant
  # an unsupervised screening fit would have been silently refitted supervised.
  args$useY <- fit_pen$useY
  do.call(estimate_lucid, args)
}

# Dispatcher for the three input-preparation helpers above.
prepare_selected_inputs <- function(model_type, fit_pen, G, Z) {
  switch(model_type,
    early    = prepare_early_selected_inputs(fit_pen, G, Z),
    parallel = prepare_parallel_selected_inputs(fit_pen, G, Z),
    serial   = prepare_serial_selected_inputs(fit_pen, G, Z),
    stop("unknown model_type: ", model_type)
  )
}


# Compact stage-wise feature-selection report for serial fits, built entirely
# from get_selected_G()/get_selected_Z() -- no per-stage dispatch of its own.
serial_selection_report <- function(fit_serial_pen) {
  selected_z <- get_selected_Z(fit_serial_pen)
  out <- vector("list", length(fit_serial_pen$submodel))
  for (i in seq_along(fit_serial_pen$submodel)) {
    sm <- fit_serial_pen$submodel[[i]]
    if (inherits(sm, "early_lucid")) {
      out[[i]] <- list(
        stage = i,
        model = "early",
        selected_G = if (i == 1) sum(get_selected_G(fit_serial_pen)) else NA,
        total_G = if (i == 1) length(get_selected_G(fit_serial_pen)) else NA,
        selected_Z = sum(selected_z[[i]]),
        total_Z = length(selected_z[[i]])
      )
    } else {
      out[[i]] <- list(
        stage = i,
        model = "parallel",
        selected_G = if (i == 1) sum(get_selected_G(fit_serial_pen)) else NA,
        total_G = if (i == 1) length(get_selected_G(fit_serial_pen)) else NA,
        selected_Z_by_layer = sapply(selected_z[[i]], sum),
        total_Z_by_layer = sapply(selected_z[[i]], length)
      )
    }
  }
  out
}

7) Early Model Tutorial: Binary Outcome

Y_binary is the median split of the continuous outcome built in section 5, on the same subjects, with the same injected missingness. Fitting proceeds in the same three explicit steps as the continuous-outcome companion vignette: penalized screening fit, zero-penalty refit on the survivors, then bootstrap. Only the outcome model changes.

7.1 Penalized screening fit

set.seed(1105)

early_pen_bin <- estimate_lucid(
  lucid_model = "early",
  G = G,
  Z = Z_early_miss,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  family = "binary",
  K = 2,
  Rho_G = RHO_G,
  Rho_Z_Mu = 0,
  Rho_Z_Cov = 0,
  max_itr = 15,
  max_tot.itr = 40,
  tol = 1e-2,
  seed = 1105,
  verbose = FALSE
)

summary(early_pen_bin)

7.2 Zero-penalty selected-only refit

set.seed(1106)

early_inputs_bin <- prepare_early_selected_inputs(early_pen_bin, G, Z_early_miss)

early_bin <- list(fit_pen = early_pen_bin, inputs = early_inputs_bin)
early_bin$fit_refit <- refit_selected(
  "early",
  fit_pen = early_pen_bin,
  inputs = early_inputs_bin,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  seed = 1106
)

summary(early_bin$fit_refit)

7.3 Bootstrap CI + summary

set.seed(1107)

early_bin$boot <- boot_lucid(
  G = early_inputs_bin$G,
  Z = early_inputs_bin$Z,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  model = early_bin$fit_refit,
  R = 30,
  conf = 0.90
)

summary(early_bin$fit_refit, boot.se = early_bin$boot)

Every bootstrap CI table now also carries an odds-ratio view (OR, OR_lower, OR_upper) alongside the coefficient-scale columns and sig, and the exposure/assignment (3) E table includes its intercept and covariate rows. R = 30 keeps the vignette fast; a real analysis needs R in the hundreds. A poorly identified coefficient -- the exposure-model intercept, or a covariate nearly collinear with cluster membership -- can show a very wide or off-centre interval; that is the bootstrap being honest, not a bug.

Three differences are worth noting in that output whenever an outcome switches from continuous to binary -- they are the whole payload of the outcome-family axis, and apply identically to the parallel and serial models fitted below:

• Outcome effects are on the log-odds scale. For a continuous outcome the cluster terms are means; here they are log-odds, and the printed exposure effects carry odds ratios alongside the coefficients.
• There is no residual standard deviation. A Bernoulli outcome has no dispersion parameter to estimate, so the sigma row that appears for a normal outcome is simply absent.
• The omics and exposure parts are unchanged in form. res_Mu and res_Beta mean the same thing regardless of family -- only the X -> Y arm differs.

8) Parallel Model Tutorial: Binary Outcome

8.1 Penalized screening fit

set.seed(1205)

parallel_pen_bin <- estimate_lucid(
  lucid_model = "parallel",
  G = G,
  Z = Z_parallel_miss,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  family = "binary",
  K = c(2, 2, 2),
  Rho_G = RHO_G,
  Rho_Z_Mu = 0,
  Rho_Z_Cov = 0,
  max_itr = 15,
  max_tot.itr = 40,
  tol = 1e-2,
  seed = 1205,
  verbose = FALSE
)

summary(parallel_pen_bin)

8.2 Zero-penalty selected-only refit

set.seed(1206)

parallel_inputs_bin <- prepare_parallel_selected_inputs(parallel_pen_bin, G, Z_parallel_miss)

parallel_bin <- list(fit_pen = parallel_pen_bin, inputs = parallel_inputs_bin)
parallel_bin$fit_refit <- refit_selected(
  "parallel",
  fit_pen = parallel_pen_bin,
  inputs = parallel_inputs_bin,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  seed = 1206
)

summary(parallel_bin$fit_refit)

8.3 Bootstrap CI + summary

set.seed(1207)

parallel_bin$boot <- boot_lucid(
  G = parallel_inputs_bin$G,
  Z = parallel_inputs_bin$Z,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  model = parallel_bin$fit_refit,
  R = 30,
  conf = 0.90
)

summary(parallel_bin$fit_refit, boot.se = parallel_bin$boot)

One structural difference between the families is specific to the parallel model and easy to miss. For a normal outcome, the early model estimates a per-cluster residual standard deviation -- one value per latent cluster -- whereas the parallel model estimates a single pooled standard deviation across the joint cluster configuration. For a binary outcome neither exists. So if you are comparing dispersion across model types, compare like with like.

9) Serial Model Tutorial A (All-Early Stages): Binary Outcome

9.1 Penalized screening fit

# Serial structure: list of early-stage matrices.
Z_serial_all_early <- list(
  methylome = Z_parallel_miss[[1]],
  transcriptome = Z_parallel_miss[[2]],
  miRNA = Z_parallel_miss[[3]]
)

set.seed(1305)

serial_ae_pen_bin <- estimate_lucid(
  lucid_model = "serial",
  G = G,
  Z = Z_serial_all_early,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  family = "binary",
  K = list(2, 2, 2),
  Rho_G = RHO_G,
  Rho_Z_Mu = 0,
  Rho_Z_Cov = 0,
  max_itr = 15,
  max_tot.itr = 40,
  tol = 1e-2,
  seed = 1305,
  verbose = FALSE
)

summary(serial_ae_pen_bin)

9.2 Zero-penalty selected-input refit

set.seed(1306)

serial_ae_inputs_bin <- prepare_serial_selected_inputs(serial_ae_pen_bin, G, Z_serial_all_early)

serial_ae_bin <- list(fit_pen = serial_ae_pen_bin, inputs = serial_ae_inputs_bin)
serial_ae_bin$fit_refit <- refit_selected(
  "serial",
  fit_pen = serial_ae_pen_bin,
  inputs = serial_ae_inputs_bin,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  seed = 1306
)

summary(serial_ae_bin$fit_refit)

9.3 Bootstrap CI + summary

set.seed(1307)

serial_ae_bin$boot <- boot_lucid(
  G = serial_ae_inputs_bin$G,
  Z = serial_ae_inputs_bin$Z,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  model = serial_ae_bin$fit_refit,
  R = 30,
  conf = 0.90
)

summary(serial_ae_bin$fit_refit, boot.se = serial_ae_bin$boot)

For a serial model the outcome family applies to the final stage only. Upstream stages are fitted unsupervised regardless of what you pass, because the outcome enters the chain once, at the end -- so the binary/normal distinction shows up in the last stage's report and nowhere else.

10) Serial Model Tutorial B (Mixed Parallel + Early): Binary Outcome

10.1 Penalized screening fit

# Nested list signals a parallel submodel at stage 1, followed by early stage 2.
Z_serial_mixed <- list(
  list(
    methylome = Z_parallel_miss[[1]],
    transcriptome = Z_parallel_miss[[2]]
  ),
  miRNA = Z_parallel_miss[[3]]
)

set.seed(1405)

serial_mixed_pen_bin <- estimate_lucid(
  lucid_model = "serial",
  G = G,
  Z = Z_serial_mixed,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  family = "binary",
  K = list(list(2, 2), 2),
  Rho_G = RHO_G,
  Rho_Z_Mu = 0,
  Rho_Z_Cov = 0,
  max_itr = 15,
  max_tot.itr = 40,
  tol = 1e-2,
  seed = 1405,
  verbose = FALSE
)

summary(serial_mixed_pen_bin)

10.2 Zero-penalty selected-input refit

set.seed(1406)

serial_mixed_inputs_bin <- prepare_serial_selected_inputs(serial_mixed_pen_bin, G, Z_serial_mixed)

serial_mixed_bin <- list(fit_pen = serial_mixed_pen_bin, inputs = serial_mixed_inputs_bin)
serial_mixed_bin$fit_refit <- refit_selected(
  "serial",
  fit_pen = serial_mixed_pen_bin,
  inputs = serial_mixed_inputs_bin,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  seed = 1406
)

summary(serial_mixed_bin$fit_refit)

10.3 Bootstrap CI + summary

set.seed(1407)

serial_mixed_bin$boot <- boot_lucid(
  G = serial_mixed_inputs_bin$G,
  Z = serial_mixed_inputs_bin$Z,
  Y = Y_binary,
  CoG = CoG,
  CoY = CoY,
  model = serial_mixed_bin$fit_refit,
  R = 30,
  conf = 0.90
)

summary(serial_mixed_bin$fit_refit, boot.se = serial_mixed_bin$boot)

This is the most general configuration the package supports: a serial chain whose first stage is itself a parallel model over two omics layers, fitted to a binary outcome, with both listwise and sporadic missingness present. If this runs and reports sane estimates, the combination space is covered.

11) Visualization: Sankey Diagram and Cluster Omics Profiles

plot() renders an early-integration fit as a Sankey diagram: exposures flow into the latent clusters, and the clusters flow on into the omics features and the outcome. Everything about how to read it -- node colour, link width and sign -- is exactly as in the continuous-outcome case (see the companion vignette's section 14); a binary outcome only changes what the final cluster -> outcome link represents (a log-odds effect rather than a mean difference).

sankey_early_bin <- plot(early_bin$fit_refit)
sankey_early_bin

As in the continuous case, plot() on a parallel or serial fit currently raises an error by design -- not implemented yet.

plot_cluster_omic_profile() -- which shows what the clusters are, via their fitted omics means -- is unaffected by outcome family entirely: res_Mu and the separation/range/sd ranking behind it describe the X -> Z arm of the model, which a binary Y never touches. Every panel, per architecture:

prof_early_bin <- plot_cluster_omic_profile(early_bin$fit_refit, top_n = 10)
prof_early_bin[[1]]

prof_par_bin <- plot_cluster_omic_profile(
  parallel_bin$fit_refit,
  layer_names = c("methylome", "transcriptome", "miRNA"),
  top_n = 8
)
for (nm in names(prof_par_bin)) print(prof_par_bin[[nm]])

prof_ser_bin <- plot_cluster_omic_profile(serial_ae_bin$fit_refit, top_n = 8)
for (nm in names(prof_ser_bin)) print(prof_ser_bin[[nm]])

The parallel model's methylome panel and the serial model's stage-1 panel above look almost identical, and that isn't a rendering glitch: stage 1 of an all-early serial chain and the corresponding layer of a parallel fit are both fit as an early-integration model on the same methylome matrix, with the same K and the same (default, mclust-based) initialization. The exposure/outcome coupling that distinguishes them is comparatively weak, so both converge to nearly the same cluster solution -- which is itself a useful sanity check that the methylation clustering is robust to which architecture surfaces it.

For the full discussion of what the importance argument measures and why (separation vs. range vs. sd), and how to pull the ranking behind a plot out as a plain table, see lucid_3models_normal_outcome.Rmd's section 15 -- none of that changes here, so it isn't repeated.

12) Prediction: Labels or Probabilities

predict_lucid()'s response argument controls whether a binary outcome comes back as class labels or as probabilities. Both draw on the early-model fit from section 7.

pred_lab <- predict_lucid(model = early_bin$fit_refit,
                          G = early_bin$inputs$G, Z = early_bin$inputs$Z,
                          CoG = CoG, CoY = CoY, response = TRUE)
pred_prob <- predict_lucid(model = early_bin$fit_refit,
                           G = early_bin$inputs$G, Z = early_bin$inputs$Z,
                           CoG = CoG, CoY = CoY, response = FALSE)

cat("response = TRUE  ->", paste(head(pred_lab$pred.y, 8), collapse = " "), "(class labels)\n")
cat("response = FALSE ->", paste(round(head(pred_prob$pred.y, 8), 3), collapse = " "), "(probabilities)\n")

The two rows describe the same underlying prediction at different granularities: each label in the first row is simply the second row's probability rounded to whichever side of 0.5 it falls on. Use response = FALSE when the downstream use needs the actual predicted risk (e.g. computing a mean predicted probability, or a classification threshold other than 0.5), and response = TRUE when a hard label is what's needed.

13) Closing Notes and Session Info

• This vignette mirrors its continuous-outcome companion (lucid_3models_normal_outcome.Rmd) one-for-one: same data, same penalized-screen-then-refit pipeline, same three model architectures -- only Y_binary in place of Y.
• For a serial model the outcome family applies to the final stage only; upstream stages are always fitted unsupervised, regardless of outcome family.
• Seeds are fixed before each fit/bootstrap step for deterministic tutorial output.

Every model this vignette fits is registered below. This is the document checking itself: if a fit failed, its object would be missing or of the wrong class, and it would be listed here rather than passing unnoticed.

check_obj("early_bin",             "list", "7 early binary")
check_obj("parallel_bin",          "list", "8 parallel binary")
check_obj("serial_ae_bin",         "list", "9 serial all-early binary")
check_obj("serial_mixed_bin",      "list", "10 serial mixed binary")
check_obj("prof_early_bin",        "list", "11 omics profile (early)")
check_obj("prof_par_bin",          "list", "11 omics profile (parallel)")
check_obj("prof_ser_bin",          "list", "11 omics profile (serial)")
check_obj("pred_lab",              "list", "12 predict (labels)")
check_obj("pred_prob",             "list", "12 predict (probabilities)")

status <- do.call(rbind, .reg$rows)
print(status, row.names = FALSE)

cat(sprintf("\n%d of %d registered steps ok; %d not ok\n",
            sum(status$status == "ok"), nrow(status),
            sum(status$status != "ok")))

Session Info

sessionInfo()

LUCIDus documentation built on Sept. 3, 2026, 1:06 a.m.