LUCIDus
LUCID with Multiple Omics Data

Package index
Search the LUCIDus package
Vignettes
Functions
175
Source code
117
Man pages
25
Home
/
CRAN
/
LUCIDus
/
LUCID's Three Model Architectures: Early, Parallel, and Serial -- Continuous Outcome (HELIX Example)

LUCID's Three Model Architectures: Early, Parallel, and Serial -- Continuous 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.

What you will learn:

Important runtime note:

2) Data Objects and Statistical Roles

The HELIX example data (simulated_HELIX_data.rda) provides:

Model inputs used throughout:

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 | | R (in boot_lucid()) | number of bootstrap resamples | kept tiny (2-5) here for speed; use hundreds for a real analysis | | 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 17 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:

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 section 12.3 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. Section 7.2 shows what this value recovers;
# section 11 sweeps it, and sweeps 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

7.1 Penalized screening fit (Rho_* > 0)

Interpretation target:

set.seed(1101)

# Screening model: penalties help identify a parsimonious subset.
early_fit_pen <- estimate_lucid(
  lucid_model = "early",
  G = G,
  Z = Z_early_miss,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  family = "normal",
  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 = 1101,
  verbose = FALSE
)

# Summary shows model fit, missing profile, selection, and parameter tables.
summary(early_fit_pen)

7.1b Object structure: what's in the fitted object

summary() above is a printed report; early_fit_pen itself is a plain R list, and every field on it is documented on ?estimate_lucid's @return. This is what an early-integration fit actually holds:

names(early_fit_pen)

The fields fall into two groups:

cat("log-likelihood:", early_fit_pen$likelihood, "\n")
cat("exposures selected:", sum(get_selected_G(early_fit_pen)),
    "of", length(get_selected_G(early_fit_pen)), "\n")
str(early_fit_pen$em_control, max.level = 1)

7.2 Inspect selected features

early_selected_G <- names(which(get_selected_G(early_fit_pen)))

cat("exposures retained:", paste(early_selected_G, collapse = ", "), "\n\n")

# Score the selection against the planted truth from section 5: the three
# g_causal_* exposures drive the latent subgroup, the six g_noise_* do not.
data.frame(
  exposure = colnames(G),
  truth    = ifelse(colnames(G) %in% causal_exposures, "causal", "noise"),
  selected = colnames(G) %in% early_selected_G
)

cat("\ncausal exposures kept :", sum(causal_exposures %in% early_selected_G), "of 3\n")
cat("noise  exposures kept :",
    sum(!(early_selected_G %in% causal_exposures)), "of 6\n")

# The omics side is unpenalized in this screening fit, so every feature is
# retained. Section 11.3 penalizes the omics separately and shows why.
cat("omics features retained:", sum(get_selected_Z(early_fit_pen)),
    "of", length(get_selected_Z(early_fit_pen)), "(Rho_Z_Mu = 0 here)\n")

7.3 Zero-penalty selected-only refit

Why this step:

set.seed(1102)

# Build selected-only inputs.
early_inputs_refit <- prepare_early_selected_inputs(early_fit_pen, G, Z_early_miss)

# Refit with all penalties set to zero.
early_fit_refit <- refit_selected(
  "early",
  fit_pen = early_fit_pen,
  inputs = early_inputs_refit,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  seed = 1102,
  verbose = TRUE
)

summary(early_fit_refit)

Compare this summary against the penalized fit's in 7.1: the coefficients on the selected exposures/features should be of similar sign and rough magnitude (the model didn't change what it found, only stopped shrinking it), while any exposure/feature that was zeroed out in the screening step is now simply absent from the model rather than penalized toward zero. This refit -- not the penalized screening fit -- is the model bootstrapped next.

7.4 Bootstrap CI + summary

set.seed(1103)

# Bootstrap on zero-penalty refit for CI inference.
early_boot <- boot_lucid(
  G = early_inputs_refit$G,
  Z = early_inputs_refit$Z,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  model = early_fit_refit,
  R = 30,
  conf = 0.90
)

# Print summary with CI columns integrated into parameter tables.
summary(early_fit_refit, boot.se = early_boot)

Each coefficient table above now carries a normal-theory confidence interval, an odds-ratio view (OR, OR_lower, OR_upper -- the exp() of the coefficient-scale columns), and the sig column: "*" marks a row whose interval excludes 0, i.e. an effect the bootstrap distinguishes from no effect at all. R = 30 here keeps the vignette fast; a real analysis needs R in the hundreds, as noted in section 3's hyperparameter guide. A poorly identified coefficient -- the exposure-model intercept when the clusters separate cleanly, or a covariate nearly collinear with cluster membership -- can still show a very wide or off-centre interval; that is the bootstrap being honest, not a bug.

8) Parallel Model Tutorial

8.1 Penalized screening fit

Interpretation target:

set.seed(1201)

parallel_fit_pen <- estimate_lucid(
  lucid_model = "parallel",
  G = G,
  Z = Z_parallel_miss,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  family = "normal",
  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 = 1201,
  verbose = FALSE
)

summary(parallel_fit_pen)

8.1b Object structure: what's in the fitted object

names(parallel_fit_pen)

Compared to the early fit in section 7.1b, two fields are new and two common fields change shape:

cat("log-likelihood:", parallel_fit_pen$likelihood, "\n")
cat("sample size (N):", parallel_fit_pen$N, "\n")
cat("exposures selected per layer:\n")
print(sapply(seq_along(parallel_fit_pen$K), function(i) sum(get_selected_G(parallel_fit_pen, layer = i))))

8.2 Inspect selected features (union + per-layer)

# Exposure selection union across layers.
parallel_selected_G_union <- names(which(get_selected_G(parallel_fit_pen)))

# Exposure selection per layer.
parallel_selected_G_layer <- lapply(seq_along(parallel_fit_pen$K), function(i) {
  names(which(get_selected_G(parallel_fit_pen, layer = i)))
})

# Omics selection per layer (get_selected_Z already collapses any
# vector/matrix selection encoding to one logical value per feature).
parallel_selected_Z_layer <- lapply(get_selected_Z(parallel_fit_pen), function(s) names(which(s)))

parallel_selected_G_union
parallel_selected_G_layer
parallel_selected_Z_layer

8.3 Zero-penalty selected-only refit

set.seed(1202)

parallel_inputs_refit <- prepare_parallel_selected_inputs(parallel_fit_pen, G, Z_parallel_miss)
parallel_fit_refit <- refit_selected(
  "parallel",
  fit_pen = parallel_fit_pen,
  inputs = parallel_inputs_refit,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  seed = 1202,
  verbose = TRUE
)

summary(parallel_fit_refit)

Unlike the early refit in 7.3, this summary prints one G -> X and one X -> Z block per layer -- each layer's cluster variable is its own model, refit unpenalized on that layer's own selected features.

8.4 Bootstrap CI + summary

set.seed(1203)

parallel_boot <- boot_lucid(
  G = parallel_inputs_refit$G,
  Z = parallel_inputs_refit$Z,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  model = parallel_fit_refit,
  R = 30,
  conf = 0.90
)

summary(parallel_fit_refit, boot.se = parallel_boot)

As in 7.4, the sig column flags rows whose normal-theory interval excludes 0 in each layer's coefficient table; with layer-specific cluster variables, a feature can be significant in one layer's table and not another's, since each layer's bootstrap resamples and refits that layer's own model.

9) Serial Model Tutorial A: All-Early Stages

Serial model logic here:

Each stage uses previous latent information as designed by the serial pipeline.

9.1 Penalized screening fit

set.seed(1301)

# 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]]
)

serial_all_early_pen <- estimate_lucid(
  lucid_model = "serial",
  G = G,
  Z = Z_serial_all_early,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  family = "normal",
  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 = 1301,
  verbose = FALSE
)

summary(serial_all_early_pen)
serial_selection_report(serial_all_early_pen)

9.1b Object structure: what's in the fitted object

names(serial_all_early_pen)

Two fields are serial-only, beyond what early/parallel already showed:

likelihood and select at the top level are present for serial too, the same as early/parallel -- but unlike those two, a serial fit has no single joint EM loop, so both are aggregates over stages rather than values from one fit:

cat("top-level likelihood (sum over stages):", serial_all_early_pen$likelihood, "\n")
cat("per-stage likelihoods:",
    paste(round(sapply(serial_all_early_pen$submodel, `[[`, "likelihood"), 2),
          collapse = ", "), "\n")
cat("n_stages:", length(serial_all_early_pen$submodel), "\n")

The fit's top-level select is stage 1's own selection only, not one combined across all three stages. That is deliberate: stage 1 is the only stage whose "G" is the cohort's actual exposures. From stage 2 onward, "G" is the previous stage's latent cluster-membership probabilities -- there are no real exposures left to select among there, so Rho_G is forced to 0 for every stage after the first regardless of what penalty was requested, and that stage's select$selectG is not a meaningful exposure-selection result. The complete per-stage record (including every stage's real omics selection, selectZ) is always available via serial_selection_report() above, or directly at fit$submodel[[i]]$select.

# Top-level select == stage 1's select
identical(serial_all_early_pen$select, serial_all_early_pen$submodel[[1]]$select)

cat("Stage 1 exposures kept:",
    paste(names(which(get_selected_G(serial_all_early_pen))), collapse = ", "), "\n")

# Rho_G is 0 from stage 2 on -- there is nothing there to select
sapply(serial_all_early_pen$submodel, function(sm) sm$Rho$Rho_G)

9.2 Zero-penalty selected-input refit

set.seed(1302)

serial_all_early_inputs_refit <- prepare_serial_selected_inputs(
  serial_all_early_pen, G, Z_serial_all_early
)
serial_all_early_refit <- refit_selected(
  "serial",
  fit_pen = serial_all_early_pen,
  inputs = serial_all_early_inputs_refit,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  seed = 1302,
  verbose = TRUE
)

summary(serial_all_early_refit)

As with the other architectures, this refit's coefficients are unshrunk versions of the screening fit's; a serial model additionally reports one X -> Z/G -> X block per stage, plus res_Delta, the transition coefficients linking each stage's cluster to the next stage's.

9.3 Bootstrap CI + summary

set.seed(1303)

serial_all_early_boot <- boot_lucid(
  G = serial_all_early_inputs_refit$G,
  Z = serial_all_early_inputs_refit$Z,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  model = serial_all_early_refit,
  R = 30,
  conf = 0.90
)

summary(serial_all_early_refit, boot.se = serial_all_early_boot)

Each stage's coefficient tables get their own sig markers, independently of the other stages -- a stage-2 effect being significant says nothing about stage 3, since each stage's bootstrap resamples and refits that stage's own conditional model.

10) Serial Model Tutorial B: Mixed Parallel + Early

Serial mixed architecture used here:

10.1 Penalized screening fit

set.seed(1401)

# 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]]
)

serial_mixed_pen <- estimate_lucid(
  lucid_model = "serial",
  G = G,
  Z = Z_serial_mixed,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  family = "normal",
  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 = 1401,
  verbose = FALSE
)

summary(serial_mixed_pen)
serial_selection_report(serial_mixed_pen)

10.2 Zero-penalty selected-input refit

set.seed(1402)

serial_mixed_inputs_refit <- prepare_serial_selected_inputs(serial_mixed_pen, G, Z_serial_mixed)
serial_mixed_refit <- refit_selected(
  "serial",
  fit_pen = serial_mixed_pen,
  inputs = serial_mixed_inputs_refit,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  seed = 1402,
  verbose = FALSE
)

summary(serial_mixed_refit)

Stage 1's block here is itself a parallel-model summary (one G -> X/X -> Z per layer, as in section 8), nested inside the serial report; stage 2's block looks like the single-layer early summary from section 7. This is the same mix-and-match structure estimate_lucid() fit -- the summary just mirrors it.

10.3 Bootstrap CI + summary

set.seed(1403)

serial_mixed_boot <- boot_lucid(
  G = serial_mixed_inputs_refit$G,
  Z = serial_mixed_inputs_refit$Z,
  Y = Y,
  CoG = CoG,
  CoY = CoY,
  model = serial_mixed_refit,
  R = 30,
  conf = 0.90
)

summary(serial_mixed_refit, boot.se = serial_mixed_boot)

With this, every architecture the package supports (early, parallel, a serial chain of early stages, and a serial chain mixing a parallel stage with an early stage) has been fit, refit unpenalized, bootstrapped, and summarized with significance markers -- the full pipeline this tutorial set out to demonstrate.

11) Missing-Data Handling: Diagnostics and Verification

Sections 7-10 injected missingness and let summary() display the resulting profile. That shows the model noticed the missing values; it does not show it handled them correctly. This section checks the handling.

LUCID classifies every subject, per omics layer, into one of three states:

| code | state | meaning | |---|---|---| | 1 | complete | every feature observed | | 2 | sporadic | some features observed, some missing | | 3 | listwise | no feature observed for that layer |

The distinction matters: imputation is triggered by sporadic subjects only. A subject missing an entire layer has no within-layer information to condition on, so the omics term drops out of that subject's likelihood and they contribute through the exposure and outcome arms instead. Inventing values for them would be fabricating data.

na_early <- check_na(Z_early_miss, lucid_model = "early")

cat("subjects by missingness code (1 complete, 2 sporadic, 3 listwise):\n")
print(table(na_early$indicator_na))

cat("\nimpute_flag (TRUE only when sporadic subjects exist):",
    na_early$impute_flag, "\n")

analyze_missing_pattern() gives the per-feature and per-subject view. The count of distinct patterns is the number worth watching: the observed-data likelihood is evaluated once per distinct missingness pattern, so pattern diversity, not sparsity alone, is what drives fitting cost.

pat <- analyze_missing_pattern(Z_early_miss)

cat("distinct missingness patterns :", pat$n_patterns, "\n")
cat("complete subjects             :", pat$n_complete, "\n")
cat("overall proportion missing    :", round(pat$total_missing, 4), "\n")

11.1 Verify the handling, do not just display it

The two claims above are checkable against the fitted object, so we check them rather than assert them in prose. Recall the injection at section 5: row 1 of Z_early was blanked entirely (listwise), and a few cells in rows 2-5 were blanked individually (sporadic).

fitted_Z <- early_fit_pen$Z

# The listwise subject must still be NA after fitting.
stopifnot(all(is.na(fitted_Z[1, ])))

# The sporadic cells must have been filled.
stopifnot(!anyNA(fitted_Z[2:5, ]))

# Observed values must never have been altered.
obs <- !is.na(Z_early_miss)
stopifnot(isTRUE(all.equal(fitted_Z[obs], Z_early_miss[obs])))

cat("verified: listwise row still NA, sporadic cells filled,",
    "observed values unchanged\n")

One asymmetry to know about when you go looking for imputed values: fit$Z holds the fitted omics for early and parallel models, but for a serial model it holds the omics as supplied, because imputation happens inside each stage. The imputed data for serial stage i is at fit$submodel[[i]]$Z.

cat("early  fit$Z still has NAs? ", anyNA(early_fit_pen$Z), "\n")
cat("serial fit$Z still has NAs? ", anyNA(serial_all_early_pen$Z[[1]]),
    "  (imputed data lives on submodel[[i]]$Z)\n")

11.2 Single-value imputation, and scoring it

safe_impute() is the cruder alternative: fill each feature with one summary of its observed values. It is what LUCID uses to initialize before the EM algorithm takes over with model-based imputation, and it is exported so you can inspect that starting point. check_imputation_quality() scores any imputation by asking whether the filled values have shifted the centre or spread of the data relative to what was observed.

Z_mean_filled <- safe_impute(Z_early_miss, method = "mean")

qual <- check_imputation_quality(Z_early_miss, Z_mean_filled)
cat("mean-filled  : valid =", qual$is_valid,
    " mean shift =", round(qual$mean_diff, 3),
    " sd ratio =", round(qual$sd_ratio, 3), "\n")

An sd_ratio well below 1 is the signature of mean-filling: every imputed value sits exactly at the centre, so the imputed cells carry no variation. The model-based imputation inside the EM loop conditions on each subject's observed coordinates and on their cluster, which is why it recovers the withheld values more accurately than this baseline.

11.3 Penalizing the omics: what the two penalties actually do

Rho_G and Rho_Z_Mu are both called LASSO penalties, but they act on different parts of the model and they live on completely different scales. This is the single easiest thing to get wrong.

Rho_G penalizes the exposure coefficients in the G -> X model. Values in the region of 0.01-0.1 are meaningful, which is why this tutorial uses RHO_G = 0.05.

Rho_Z_Mu penalizes the cluster-specific omics means, and it is compared against those means directly. Our signal features are separated by roughly three standard deviations, so a penalty of 0.02 -- the value this tutorial used before -- is far too small to threshold anything at all.

sweep_rho_z <- function(rz) {
  invisible(capture.output(
    f <- estimate_lucid(
      lucid_model = "early", G = G, Z = Z_early_miss, Y = Y, CoG = CoG, CoY = CoY,
      family = "normal", K = 2, init_omic.data.model = "EEV",
      Rho_G = 0, Rho_Z_Mu = rz, Rho_Z_Cov = if (rz > 0) 0.02 else 0,
      max_itr = 30, max_tot.itr = 90, tol = 1e-2, seed = 1101
    )
  ))
  kept <- which(get_selected_Z(f))
  acc <- mean(get_cluster_assignment(f) == x_true + 1); acc <- max(acc, 1 - acc)
  data.frame(
    Rho_Z_Mu       = rz,
    features_kept  = length(kept),
    signal_kept    = sum(kept %in% signal_cols_early),
    converged      = isTRUE(f$em_control$converged),
    subgroup_recovery = round(acc, 2)
  )
}

do.call(rbind, lapply(c(0, 2, 10), sweep_rho_z))

Two things to take from that table, and the second is the reason this tutorial leaves the omics unpenalized.

The penalty does eventually select. At 0.02 -- the value this tutorial used before -- nothing is dropped at all. It takes a penalty two to three orders of magnitude larger before the feature count moves.

But look at signal_kept as it does. The penalty is not preferentially keeping the features that carry the subgroup signal; it drops those alongside the rest. On this data it is shrinking indiscriminately rather than selecting.

And it costs cluster recovery. Watch the subgroup_recovery column against the known x_true. Unpenalized, the model recovers the planted subgroup well above chance. As the omics penalty rises, the cluster means are shrunk toward each other, the subgroups become harder to tell apart, and recovery falls back toward 0.5 while the fit stops converging.

That is not a bug; it is the bias-variance trade the penalty exists to make, and here the variance saving is not worth the bias. Two practical rules follow:

12) Prediction and G-Computation

Fitting is only half of it. predict_lucid() assigns clusters and predicts outcomes, either on new data or -- which is how you extract the fitted cluster memberships -- on the training data.

# CoG and CoY must be supplied here exactly as they were at fitting time: the
# G -> X design includes the CoG columns, so omitting them makes the design
# matrix too narrow for the fitted coefficients.
# lucid_model is auto-detected from early_fit_refit's own class.
pred_train <- predict_lucid(
  model = early_fit_refit,
  G = early_inputs_refit$G, Z = early_inputs_refit$Z, Y = Y,
  CoG = CoG, CoY = CoY
)

# Predicting on the training data reproduces the fit's own assignment.
stopifnot(identical(
  as.numeric(pred_train$pred.x),
  as.numeric(get_cluster_assignment(early_fit_refit))
))

cat("cluster sizes:\n"); print(table(pred_train$pred.x))

Cluster labels run 1..K. Versions before 3.1.0 returned 0..K-1, so code written against the old convention needs its + 1 removed.

12.1 What you may omit, and why

predict_lucid() takes both Z and Y, and both look optional from the signature. Only one of them is.

| | Y supplied | Y omitted | |---|---|---| | Z supplied | works -- supervised | works -- unsupervised | | Z omitted | error | error | | Z omitted, g_computation = TRUE | works | works |

The asymmetry follows from what the E-step actually does. Cluster membership is a posterior formed from three likelihood terms: the exposures, the omics, and (optionally) the outcome. Dropping Y removes one term and leaves a perfectly well-defined posterior over the other two -- that is unsupervised prediction. Dropping Z removes the term the clusters are defined by, and there is nothing left to condition on.

g_computation = TRUE is not a way around that. It is a different estimator: it discards the omics and outcome terms entirely and derives the posterior from the exposure path alone, which is exactly why it can run with no omics data and why it is the only mode that returns pred.z. Passing Z or Y to it is accepted but they are ignored, with a printed notice to that effect.

Rather than take the table on trust, run it:

legal <- function(expr) {
  r <- try(suppressWarnings(suppressMessages(invisible(capture.output(force(expr))))),
           silent = TRUE)
  if (inherits(r, "try-error")) "error" else "works"
}

Gp <- early_inputs_refit$G
Zp <- early_inputs_refit$Z
call_p <- function(...) predict_lucid(model = early_fit_refit,
                                      G = Gp, CoG = CoG, CoY = CoY, ...)

data.frame(
  Z              = c("supplied", "supplied", "omitted", "omitted", "omitted"),
  Y              = c("supplied", "omitted",  "supplied", "omitted", "omitted"),
  g_computation  = c(FALSE, FALSE, FALSE, FALSE, TRUE),
  result         = c(
    legal(call_p(Z = Zp, Y = Y)),
    legal(call_p(Z = Zp)),
    legal(call_p(Z = NULL, Y = Y)),
    legal(call_p(Z = NULL)),
    legal(call_p(Z = NULL, g_computation = TRUE))
  )
)

And this is the message you get if you forget it, which names the one mode that relaxes the rule:

cat(tryCatch(
  call_p(Z = NULL, Y = Y),
  error = function(e) conditionMessage(e)
))

Two traps worth knowing before the sections below.

Pass CoG and CoY exactly as you fitted them. The G -> X design includes the CoG columns, so omitting them at prediction leaves the design matrix narrower than the fitted coefficients and you get non-conformable arguments -- an error that says nothing about the actual mistake.

A single-stage serial model cannot be predicted. It fits without complaint, but a one-stage chain is a fully equivalent early or parallel model, and prediction declines it with a message saying so. Fit it as the model it is.

12.2 Prediction across all three architectures

The rule above applies identically to every model type. Here it is exercised on the three models this tutorial has already fitted, reporting what was actually predicted rather than only that the call returned.

# pred.x nests: a vector for early, a list by layer for parallel, and for serial
# a list by stage whose elements are themselves lists when that stage is a
# parallel submodel. Flatten to the leaves so each cluster variable is counted on
# its own -- pooling a parallel stage's layers would report twice as many
# assignments as there are subjects.
flatten_blocks <- function(x, path = "") {
  if (!is.list(x)) return(stats::setNames(list(as.numeric(x)), path))
  out <- list()
  for (i in seq_along(x)) {
    nm <- if (nzchar(path)) paste0(path, ".", i) else as.character(i)
    out <- c(out, flatten_blocks(x[[i]], nm))
  }
  out
}

cluster_sizes <- function(pred_x, label) {
  blocks <- flatten_blocks(pred_x)
  nms <- names(blocks)
  # index by position: the single block of an early model is named "", and
  # blocks[[""]] does not select anything.
  do.call(rbind, lapply(seq_along(blocks), function(i) {
    tb <- table(factor(blocks[[i]]))
    data.frame(model = label,
               block = if (!nzchar(nms[i])) "-" else nms[i],
               cluster = names(tb), n = as.integer(tb), row.names = NULL)
  }))
}

outcome_summary <- function(pred_y, label) {
  v <- as.numeric(unlist(pred_y))
  data.frame(model = label, n = length(v),
             mean = round(mean(v), 3), sd = round(stats::sd(v), 3),
             min = round(min(v), 3), median = round(stats::median(v), 3),
             max = round(max(v), 3), row.names = NULL)
}

pred_early_m <- predict_lucid(
  model = early_fit_refit,
  G = early_inputs_refit$G, Z = early_inputs_refit$Z, Y = Y,
  CoG = CoG, CoY = CoY
)
pred_parallel_m <- predict_lucid(
  model = parallel_fit_refit,
  G = parallel_inputs_refit$G, Z = parallel_inputs_refit$Z, Y = Y,
  CoG = CoG, CoY = CoY
)
pred_serial_m <- predict_lucid(
  model = serial_all_early_refit,
  G = serial_all_early_inputs_refit$G, Z = serial_all_early_inputs_refit$Z, Y = Y,
  CoG = CoG, CoY = CoY
)

rbind(
  cluster_sizes(pred_early_m$pred.x,    "early"),
  cluster_sizes(pred_parallel_m$pred.x, "parallel (per layer)"),
  cluster_sizes(pred_serial_m$pred.x,   "serial (per stage)")
)

Predicted outcomes for the same three models:

rbind(
  outcome_summary(pred_early_m$pred.y,    "early"),
  outcome_summary(pred_parallel_m$pred.y, "parallel"),
  outcome_summary(pred_serial_m$pred.y,   "serial"),
  outcome_summary(Y,                      "observed outcome")
)

Two things are worth reading off that comparison.

The predicted outcome is much less variable than the observed one. It is a posterior-weighted average of the cluster-specific outcome levels, so it carries no residual variation at all -- only the between-cluster differences. The ratio of its spread to the observed spread is a rough indication of how much of the outcome the latent structure explains.

Cluster sizes are worth checking before anything else. A block where almost every subject lands in one cluster is a degenerate solution, and any downstream interpretation of that layer or stage is unsafe regardless of what the summary tables report.

12.3 Supervised versus unsupervised prediction

Supplying Y lets the outcome inform the posterior, exactly as it does during fitting. Omitting it predicts clusters from exposures and omics alone -- which is what you want when the outcome is unavailable, or when it must not influence the assignment.

pred_unsup <- predict_lucid(
  model = early_fit_refit,
  G = early_inputs_refit$G, Z = early_inputs_refit$Z,
  CoG = CoG, CoY = CoY
)

cat("agreement between supervised and unsupervised assignment:",
    round(mean(pred_train$pred.x == pred_unsup$pred.x), 3), "\n")

12.5 G-computation: what if the exposure were different?

G-computation predicts from the exposures alone, dropping the omics and outcome terms from the E-step. That makes it a counterfactual tool: hold the fitted model fixed, change G, and read off what the model implies. A single call demonstrates nothing on its own -- the contrast between two exposure scenarios is the point.

Here we shift the strongest signal exposure by one standard deviation in each direction.

G_ref <- early_inputs_refit$G

# Shift every retained causal exposure by 1.5 SD, each in the direction of its
# own effect, i.e. "what if this whole exposure profile were less favourable?"
kept_causal <- intersect(colnames(G_ref), causal_exposures)
direction <- c(g_causal_1 = 1, g_causal_2 = -1, g_causal_3 = 1)

shift_profile <- function(Gx, delta) {
  for (nm in kept_causal) Gx[, nm] <- Gx[, nm] + delta * direction[[nm]]
  Gx
}

gc_hi <- predict_lucid(model = early_fit_refit,
                       G = shift_profile(G_ref,  1.5), Z = NULL,
                       CoG = CoG, CoY = CoY, g_computation = TRUE)
gc_lo <- predict_lucid(model = early_fit_refit,
                       G = shift_profile(G_ref, -1.5), Z = NULL,
                       CoG = CoG, CoY = CoY, g_computation = TRUE)

data.frame(
  scenario          = c("profile +1.5 SD", "profile -1.5 SD"),
  mean_pred_outcome = c(mean(gc_hi$pred.y), mean(gc_lo$pred.y)),
  prop_in_cluster_2 = c(mean(gc_hi$inclusion.p[, 2]), mean(gc_lo$inclusion.p[, 2]))
)

Read that table as a decomposition. G-computation changes the outcome only by moving subjects between clusters, so the effect it can produce is governed by

(change in cluster membership) x (gap between the clusters' outcome levels)

The two agree exactly when the outcome model has no covariates; with CoY in the model each subject carries their own covariate offset, so the identity holds up to that variation.

# early_gamma_levels() is the internal accessor predict_lucid() itself uses to
# turn the fitted outcome parameters into one level per cluster -- not part of
# the public API, called here (via :::) only to show the arithmetic behind the
# g-computation contrast above. Reading res_Gamma$beta directly would give the
# stored parameterization, which is not the same thing.
levels_by_cluster <- LUCIDus:::early_gamma_levels(early_fit_refit$res_Gamma, early_fit_refit$K)
shift_in_membership <- mean(gc_hi$inclusion.p[, 2]) - mean(gc_lo$inclusion.p[, 2])

cat("cluster outcome levels     :", round(levels_by_cluster, 3), "\n")
cat("gap between clusters       :", round(diff(levels_by_cluster), 3), "\n")
cat("shift in cluster-2 share   :", round(shift_in_membership, 3), "\n")
cat("implied outcome contrast   :",
    round(shift_in_membership * diff(levels_by_cluster), 4), "\n")
cat("observed outcome contrast  :",
    round(mean(gc_hi$pred.y) - mean(gc_lo$pred.y), 4), "\n")

This is worth internalizing before applying g-computation to real data. If the contrast comes out near zero there are exactly two possible reasons, and they call for different responses: either the exposures barely move cluster membership, or the clusters barely differ in outcome. Printing both terms tells you which. A near-zero contrast with a large membership shift means the clusters are not outcome-relevant; a large gap with no membership shift means the exposures are not cluster-relevant.

A third possibility worth ruling out first: if prop_in_cluster_2 is 0 or 1 in both scenarios, the exposure model has saturated and no shift of any size will move it. That is a property of the fit, not a finding about the exposure.

pred.z -- the omics profile the model implies under each scenario -- is returned in this mode as well, so the counterfactual can be read on the omics layer and not only on the outcome.

cat("predicted omics means differ between scenarios by (first 5 features):\n")
print(round(head(colMeans(gc_hi$pred.z) - colMeans(gc_lo$pred.z), 5), 4))

These differences reflect the same underlying cluster-membership shift as the outcome contrast above, just read off the omics side instead of Y -- a feature with a larger implied shift here is one whose cluster-specific means are further apart, consistent with plot_cluster_omic_profile()'s ranking in section 15.

13) The lucid() Wrapper and Tuning K

Sections 7-10 deliberately did the screen-then-refit by hand, because seeing the steps is the point of a tutorial. In practice lucid() does it for you: it tunes over a grid, selects by BIC, and refits the winner on the selected features without a penalty.

set.seed(1501)
fit_tuned <- lucid(
  G = G, Z = Z_early_miss, Y = Y, CoG = CoG, CoY = CoY,
  lucid_model = "early", family = "normal",
  K = 2:3,
  Rho_G = c(0, 0.02),
  init_omic.data.model = NULL,
  max_itr = 15, max_tot.itr = 40, tol = 1e-2
)

cat("selected K :", fit_tuned$K, "\n")
cat("BIC        :", round(summary(fit_tuned, auto_print = FALSE)$BIC, 2), "\n")

When a penalty selects a strict subset, the fit carries a selection component recording what was dropped from the original inputs, while select describes the refit's own dimensions.

if (!is.null(fit_tuned$selection)) {
  cat("exposures retained:",
      paste(fit_tuned$selection$Gnames[fit_tuned$selection$selectG], collapse = ", "), "\n")
  cat("exposures dropped :",
      paste(fit_tuned$selection$Gnames[!fit_tuned$selection$selectG], collapse = ", "), "\n")
} else {
  cat("no variable was deselected at this penalty grid\n")
}

One trap worth stating plainly: fit$Rho records the tuned penalty as metadata -- the value that produced the selection -- not the penalty the final refit ran with. The refit is unpenalized, which is precisely why its estimates are not shrunk.

14) Visualizing the Model: Sankey Diagram

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. Node colour encodes variable type, link width is the magnitude of the estimated effect, and link colour is its sign.

plot(early_fit_refit)

Only features the model retained are drawn, so a penalized fit produces a correspondingly sparser diagram -- which makes this a quick visual check on whether selection did what you expected.

plot() on a parallel or serial fit currently raises an error by design; those diagrams are not implemented yet. That is a known limitation, not a bug.

15) Cluster Omics Profiles

The Sankey shows the path structure. It does not answer the question a reader of the results actually asks first: what are these clusters? LUCID's clusters are defined by their omics profiles, so that question is answered by looking at the fitted cluster means -- and plot_cluster_omic_profile() does exactly that, for whichever architecture you fitted.

15.1 The default view

prof_early <- plot_cluster_omic_profile(early_fit_refit, top_n = 10)
prof_early[[1]]

Read it the way you would a single-cell marker heatmap. Each row is an omics feature, each column a latent cluster, and the fill says whether that cluster sits high or low for that feature relative to the others.

The fill is a z-score across clusters, not the raw mean. Each feature is standardized over the clusters before colouring, which is what stops a feature with a large baseline from washing out every other row. The raw means are still available -- see 15.4. Set scale = FALSE to colour by the mean instead.

15.2 Why the ranking is not the spread of the means

Only the top_n most discriminating features are drawn, and how "most discriminating" is defined matters more than it first appears.

The obvious choice is the spread of a feature's cluster means: the further apart they are, the more that feature separates the clusters. That is available as importance = "range", and it is what you get by eyeballing res_Mu. But it ignores noise. A feature whose cluster means differ by two units tells you nothing if its within-cluster standard deviation is also two -- the clusters overlap almost completely on it.

The default, "separation", divides the spread of the means by the typical within-cluster spread. It is a between-to-within ratio, the same idea as an effect size, and it is the only one of the three that can tell a genuinely separating feature from a merely variable one. It is also scale-free, so features on different measurement scales can be ranked in the same panel.

rank_by <- function(measure) {
  d <- attr(plot_cluster_omic_profile(early_fit_refit, top_n = 5,
                                      importance = measure)[[1]], "profile_data")
  as.character(unique(d$feature)[order(-unique(d[, c("feature","score")])$score)])
}

data.frame(
  rank        = 1:5,
  separation  = rank_by("separation"),
  range       = rank_by("range")
)

Where those two columns disagree, the difference is entirely within-cluster variance.

15.3 Bars, and more than one layer

The same information as bars, with clusters as shades of the layer's colour. The shading is a sequential ramp rather than a categorical palette, so it keeps working as the number of clusters grows instead of running out and recycling.

plot_cluster_omic_profile(early_fit_refit, type = "bar", top_n = 10)[[1]]

A parallel or serial fit returns one plot per layer, named, so nothing is squeezed into a single figure. All three layers, in turn:

prof_par <- plot_cluster_omic_profile(
  parallel_fit_refit,
  layer_names = c("methylome", "transcriptome", "miRNA"),
  top_n = 8
)
names(prof_par)
for (nm in names(prof_par)) print(prof_par[[nm]])

plot_cluster_omic_profile(parallel_fit_refit, type = "bar",
                          layer_names = c("methylome", "transcriptome", "miRNA"),
                          top_n = 8)[["transcriptome"]]

For a serial fit there is one plot per stage, and one per layer within any stage that is itself a parallel sub-model. All three stages, in turn:

prof_ser <- plot_cluster_omic_profile(serial_all_early_refit, top_n = 8)
names(prof_ser)
for (nm in names(prof_ser)) print(prof_ser[[nm]])

plot_cluster_omic_profile(serial_all_early_refit, type = "bar", top_n = 8)[[2]]

15.4 The numbers behind the figure

Each plot carries the data it drew, so the ranking can be reported in a table without recomputing it:

pd <- attr(prof_early[[1]], "profile_data")
head(pd[order(-pd$score), c("feature", "cluster", "mean", "sd", "score")], 6)

mean is the fitted cluster mean on the data's own scale and sd the within-cluster standard deviation, so a claim made from the figure can be quoted with its actual magnitude rather than a colour.

16) How to Read the Main Outputs

When interpreting summary(...) output, focus on:

  1. Model specification block Confirms family, sample size, cluster structure.

  2. Missing-data profile Verifies listwise vs sporadic missingness handling in each model/layer.

  3. Feature-selection overview Shows selected counts and rates for G and Z.

  4. Detailed parameter estimates

  5. (1) Y model terms (including intercept and cluster effects)
  6. (2) Z cluster means

  7. (3) E/G effects on latent-cluster assignment (with OR)

  8. Bootstrap CI columns (when boot.se is supplied) Compare estimate with its normal-theory confidence interval, and use the sig column ("*" where that interval excludes 0) as a quick scan across a table before reading every interval individually.

17) Practical Notes and Reproducibility

18) Verification Summary and Session Info

Every model this tutorial 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_fit_refit",       "early_lucid",    "7 early normal")
check_obj("parallel_fit_refit",    "lucid_parallel", "8 parallel normal")
check_obj("serial_all_early_refit","lucid_serial",   "9 serial all-early normal")
check_obj("serial_mixed_refit",    "lucid_serial",   "10 serial mixed normal")
check_obj("na_early",              "list",           "11 missing diagnostics")
check_obj("pred_train",            "list",           "12 prediction")
check_obj("gc_hi",                 "list",           "12.3 g-computation")
check_obj("fit_tuned",             "early_lucid",    "13 lucid() wrapper")
check_obj("prof_early",            "list",           "15 omics profile (early)")
check_obj("prof_par",              "list",           "15.3 omics profile (parallel)")
check_obj("prof_ser",              "list",           "15.3 omics profile (serial)")

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()


Try the LUCIDus package in your browser

Any scripts or data that you put into this service are public.

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

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close