10.1 Standard prediction
# Use lightweight no-covariate fits for robust prediction demo.
set.seed(205)
fit_early_pred <- estimate_lucid(
lucid_model = "early",
G = d$G,
Z = d$Z_early,
Y = d$Y_normal,
family = "normal",
K = 2,
max_itr = 6,
max_tot.itr = 20,
tol = 1e-2,
seed = 205
)
## Fitting LUCID early model (K = 2)...
## Finished LUCID early model.
set.seed(206)
fit_parallel_pred <- estimate_lucid(
lucid_model = "parallel",
G = d$G,
Z = d$Z_parallel,
Y = d$Y_normal,
family = "normal",
K = c(2, 2),
max_itr = 6,
max_tot.itr = 20,
tol = 1e-2,
seed = 206
)
## Fitting LUCID parallel model (2 layers)...
## Finished LUCID parallel model.
set.seed(207)
fit_serial_pred <- estimate_lucid(
lucid_model = "serial",
G = d$G,
Z = d$Z_serial_mixed,
Y = d$Y_normal,
family = "normal",
K = list(list(2, 2), 2),
max_itr = 6,
max_tot.itr = 24,
tol = 1e-2,
seed = 207
)
## Fitting LUCID serial model (2 stages)...
## Stage 1/2 (parallel) finished: log-likelihood = -611.077.
## Stage 2/2 (early) finished: log-likelihood = -383.375.
## Finished LUCID serial model.
pred_early <- predict_lucid(
model = fit_early_pred,
G = d$G,
Z = d$Z_early,
Y = d$Y_normal
)
pred_parallel <- predict_lucid(
model = fit_parallel_pred,
G = d$G,
Z = d$Z_parallel,
Y = d$Y_normal
)
pred_serial <- predict_lucid(
model = fit_serial_pred,
G = d$G,
Z = d$Z_serial_mixed,
Y = d$Y_normal
)
# Cluster assignment. pred.x is a vector for early and a list -- by layer for
# parallel, by stage for serial -- so the blocks are summarised in one place.
# 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)
}))
}
rbind(
cluster_sizes(pred_early$pred.x, "early"),
cluster_sizes(pred_parallel$pred.x, "parallel"),
cluster_sizes(pred_serial$pred.x, "serial")
)
## model block cluster n
## 1 early - 1 32
## 2 early - 2 48
## 3 parallel 1 1 44
## 4 parallel 1 2 36
## 5 parallel 2 1 43
## 6 parallel 2 2 37
## 7 serial 1.1 1 44
## 8 serial 1.1 2 36
## 9 serial 1.2 1 36
## 10 serial 1.2 2 44
## 11 serial 2 1 42
## 12 serial 2 2 38
The table above reports how many subjects fall in each cluster, per
model and (for parallel/serial) per layer/stage block. Roughly balanced
counts here are expected for this simulation’s roughly 50/50 latent
split; a wildly unbalanced split can be a sign that K is
larger than the data actually supports.
predict_lucid()’s pred.x above required
re-running the E-step on G/ Z. When the same
fitted-model hard assignment is wanted directly, without supplying data
again, get_cluster_assignment() reads it straight off the
fitted object’s own posterior (inclusion.p):
identical(as.numeric(get_cluster_assignment(fit_early_pred)),
as.numeric(pred_early$pred.x))
## [1] TRUE
The two agree here because pred_early was computed on
the same G/Z the model was fit on.
get_cluster_assignment() is the right tool whenever the
question is simply “what did this fitted model assign its subjects to,”
while predict_lucid() is for scoring genuinely new data or
getting a predicted outcome as well as an assignment.
Predicted outcomes, on the scale of the outcome that was
modelled:
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
)
}
rbind(
outcome_summary(pred_early$pred.y, "early"),
outcome_summary(pred_parallel$pred.y, "parallel"),
outcome_summary(pred_serial$pred.y, "serial")
)
## model n mean sd min median max
## 1 early 80 0.511 0.702 -0.356 1.114 1.122
## 2 parallel 80 0.506 0.691 -0.180 0.372 1.337
## 3 serial 80 0.512 0.673 -0.123 -0.090 1.252
The predicted outcome is a posterior-weighted average of the
cluster-specific outcome levels, so its spread is narrower than the
observed outcome’s: it carries no residual variation, only the
between-cluster differences. Comparing its range against
summary(d$Y_normal) shows how much of the outcome the
latent structure accounts for.
rbind(
outcome_summary(pred_early$pred.y, "predicted (early)"),
outcome_summary(d$Y_normal, "observed")
)
## model n mean sd min median max
## 1 predicted (early) 80 0.511 0.702 -0.356 1.114 1.122
## 2 observed 80 0.504 1.186 -2.024 0.306 3.470
10.2 G-computation mode (early, parallel, and serial)
g_computation = TRUE switches to a different question:
instead of scoring observed (G, Z) pairs, it asks what the
outcome distribution would be for a given exposure profile,
marginalizing over the fitted cluster and omics model – a simple causal
“what if this subject’s exposure had been X” calculation built on top of
the fitted parameters. This is why it needs only G (no
Z, no Y) and why it alone returns
pred.z, the implied omics profile under that hypothetical
exposure.
pred_early_g <- predict_lucid(
model = fit_early_pred,
G = d$G,
Z = NULL,
Y = NULL,
g_computation = TRUE
)
pred_parallel_g <- try(
predict_lucid(
model = fit_parallel_pred,
G = d$G,
Z = NULL,
Y = NULL,
g_computation = TRUE
),
silent = TRUE
)
names(pred_early_g)
## [1] "inclusion.p" "pred.x" "pred.z" "pred.y"
if (inherits(pred_parallel_g, "try-error")) {
"parallel g_computation returned a try-error on this demo object; code pattern is shown above."
} else {
names(pred_parallel_g)
}
## [1] "inclusion.p" "pred.x" "pred.z" "pred.y"
pred_serial_g <- predict_lucid(
model = fit_serial_pred,
G = d$G,
Z = NULL,
Y = NULL,
g_computation = TRUE
)
names(pred_serial_g)
## [1] "inclusion.p" "pred.x" "pred.z" "pred.y"
length(pred_serial_g$pred.z)
## [1] 2
10.3 Argument Modes: What to Pass and What You Get Back
Z is required for every model type. Only
g_computation = TRUE relaxes it. Y is always
optional; omitting it makes the prediction unsupervised.
| supplied |
supplied |
FALSE |
Supervised: outcome informs the posterior |
| supplied |
omitted |
FALSE |
Unsupervised: clusters from G and Z
only |
| omitted |
either |
FALSE |
Error naming Z and g_computation |
| omitted |
omitted |
TRUE |
Counterfactual prediction from G alone |
The posterior is formed from the exposure, omics and outcome
likelihood terms. Dropping Y removes one term and leaves a
well-defined posterior over the rest. Dropping Z removes
the term the clusters are defined by, leaving nothing to condition on.
g_computation = TRUE is a separate estimator that uses only
the exposure path, which is why it accepts Z = NULL and why
it alone returns pred.z; supplied Z and
Y are ignored in that mode.
Returned components, and their shape by model type:
inclusion.p |
Posterior cluster probabilities |
N x K matrix |
list by layer |
list by stage |
pred.x |
Cluster labels, 1..K since 3.1.0 |
vector |
list by layer |
list by stage |
pred.y |
Predicted outcome |
vector |
vector |
vector |
pred.z |
Implied omics profile, g_computation only |
matrix |
list by layer |
list by stage |
Reusing the fits from 10.1, with no refitting:
pred_unsup <- predict_lucid(
model = fit_early_pred,
G = d$G,
Z = d$Z_early
)
missing_z <- try(
predict_lucid(
model = fit_early_pred,
G = d$G,
Z = NULL,
Y = d$Y_normal
),
silent = TRUE
)
data.frame(
mode = c("Y omitted (unsupervised)", "Z omitted, no g-computation"),
result = c(
paste(names(pred_unsup), collapse = ", "),
if (inherits(missing_z, "try-error")) "error, as documented" else "unexpectedly succeeded"
)
)
## mode result
## 1 Y omitted (unsupervised) inclusion.p, pred.x, pred.y
## 2 Z omitted, no g-computation error, as documented
cat(if (inherits(missing_z, "try-error")) attr(missing_z, "condition")$message else "")
## Input data 'Z' is required for prediction. Omit it only with g_computation = TRUE, which predicts from the exposures alone.
A serial model must have at least two stages to be predicted. A
single-stage serial model is an equivalent early or parallel model and
should be fitted as one; prediction declines it with a message to that
effect.
