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.
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.
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))
)
)
## Z Y g_computation result
## 1 supplied supplied FALSE works
## 2 supplied omitted FALSE works
## 3 omitted supplied FALSE error
## 4 omitted omitted FALSE error
## 5 omitted omitted TRUE works
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)
))
## Input data 'Z' is required for prediction. Omit it only with g_computation = TRUE, which predicts from the exposures alone.
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)")
)
## model block cluster n
## 1 early - 2 90
## 2 parallel (per layer) 1 1 45
## 3 parallel (per layer) 1 2 45
## 4 parallel (per layer) 2 1 47
## 5 parallel (per layer) 2 2 43
## 6 parallel (per layer) 3 1 48
## 7 parallel (per layer) 3 2 42
## 8 serial (per stage) 1 1 45
## 9 serial (per stage) 1 2 45
## 10 serial (per stage) 2 1 42
## 11 serial (per stage) 2 2 48
## 12 serial (per stage) 3 1 48
## 13 serial (per stage) 3 2 42
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")
)
## model n mean sd min median max
## 1 early 90 0.593 0.164 0.187 0.617 0.971
## 2 parallel 90 0.593 0.690 -0.334 0.205 1.546
## 3 serial 90 0.593 0.683 -0.321 0.178 1.567
## 4 observed outcome 90 0.593 1.182 -2.477 0.559 4.286
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.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]))
)
## scenario mean_pred_outcome prop_in_cluster_2
## 1 profile +1.5 SD 0.592603 1
## 2 profile -1.5 SD 0.592603 1
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")
## cluster outcome levels : 1.722 1.722
cat("gap between clusters :", round(diff(levels_by_cluster), 3), "\n")
## gap between clusters : 0
cat("shift in cluster-2 share :", round(shift_in_membership, 3), "\n")
## shift in cluster-2 share : 0
cat("implied outcome contrast :",
round(shift_in_membership * diff(levels_by_cluster), 4), "\n")
## implied outcome contrast : 0
cat("observed outcome contrast :",
round(mean(gc_hi$pred.y) - mean(gc_lo$pred.y), 4), "\n")
## observed outcome contrast : 0
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")
## predicted omics means differ between scenarios by (first 5 features):
print(round(head(colMeans(gc_hi$pred.z) - colMeans(gc_lo$pred.z), 5), 4))
## cg_GRHL3 cg_BTF3L4 cg_AL358472.7 cg_HDGF cg_TDRD5
## 0 0 0 0 0
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.
