Calculate Individual State Occupation Probabilities

Computes individual-level state occupation probabilities (SOPs) for each row in a dataset. Optionally aggregates results within strata defined by grouping variables.

Usage

sops(
  model,
  newdata = NULL,
  times,
  y_levels = NULL,
  absorb = NULL,
  by = NULL,
  id_var = NULL,
  time_var = "time",
  p_var = "yprev",
  p2_var = NULL,
  gap_var = NULL,
  time_covariates = NULL,
  include_re = FALSE,
  n_draws = 100L,
  seed = NULL,
  posterior_summary = c("mean", "median"),
  conf_level = 0.95,
  return_draws = FALSE,
  ...
)

Arguments

model: A fitted model object (e.g., vglm, orm, or blrm). For vglm models, the family must be cumulative-logit with reverse = TRUE. For orm models, the family must be logistic.
newdata: Optional. A data frame of prediction profiles. When supplied, every row supplies a separate starting state and set of covariates. The internal rowid column is regenerated. If NULL, uses the data stored by orm_markov(), blrm_markov(), or vglm_markov() and requires exactly one complete designated starting profile per fitted patient.
times: Required. The consecutive visits at which to calculate SOPs, in increasing order, for example times = 1:30 for daily visits up to day 30. Each value is one transition, so leaving out visits does not skip time; see the section "Choosing visits". Numeric time variables use numeric values; factor-valued visit indices use fitted visit levels.
y_levels: A vector of state levels in fitted outcome order. Its length must agree with the fitted model's threshold count. If NULL, attempts to infer the levels from model.
absorb: The absorbing state.
by: Optional character vector of variable names to stratify by. When provided, the function aggregates (averages) SOPs within each stratum defined by combinations of these variables. E.g., by = "ecog" aggregates within ECOG levels; by = c("ecog", "age_group") aggregates within each combination of ECOG and age group. Rows missing any grouping value are omitted; an error is raised if none remain. This is simple aggregation within observed strata, not G-computation standardization (use avg_sops() for that).
id_var: Character ID column used when include_re = TRUE. For blrm, NULL is inferred from wrapper metadata, then model$clusterInfo$name when available, otherwise "id". For frequentist models, id_var is used for stored starting-profile validation and refit bootstrap metadata, not for ordinary user-supplied prediction rows.
time_var: Name of the time variable in the model.
p_var: Name of the previous state variable in the model.
p2_var: Optional second previous-state variable. NULL uses a first-order Markov recursion; a non-NULL column name uses a second-order recursion.
gap_var: Optional name of the variable holding the time elapsed since the previous visit, if the model includes one. Each transition then uses the difference between consecutive values of times (the first value for the first transition), unless time_covariates contains this column, so times may skip visits. See the section "Choosing visits".
time_covariates: Optional time-varying covariate lookup table for explicit basis columns. Registered inline terms such as rms::rcs(time, 4) and rms::lsp(time, c(3, 7)) can be used without supplying time_covariates.
include_re: Logical. For rmsb::blrm() fits with cluster(), include fitted random-effect draws in posterior predictions for known IDs.
n_draws: Integer number of posterior draws to sample for blrm, or NULL to use all stored draws. Defaults to 100.
seed: Optional random seed for reproducible random draw sampling.
posterior_summary: Summary used for blrm posterior SOP draws: "mean" or "median".
conf_level: Confidence level for blrm posterior intervals.
return_draws: Logical. For blrm, store posterior SOP draws for extraction with get_draws(). Large draw objects are guarded to protect memory.
...: Additional arguments (currently unused).

Value

A data frame of class markov_sops containing:

rowid: Row identifier from newdata (omitted if by is specified)
time: Time point
state: State name
estimate: Probability of being in the state (individual-level or stratum average)
conf.low, conf.high, std.error: For blrm fits, posterior uncertainty summaries computed directly from SOP draws
(by variables): Stratification variables if by is specified

Plus all columns from newdata (for individual-level results).

Details

For supported frequentist models, pass the result to inferences(method = "delta", vcov = "conditional") for analytical confidence intervals. These account for coefficient estimation while treating the starting states and covariates used for prediction as given. sops() does not support unconditional analytical variance. See inferences() for model restrictions and other inference methods.

This function wraps soprob_markov() and converts its array output to a tidy data frame. The output contains one row per patient-time-state combination (or stratum-time-state if by is used).

Automatic prediction from a wrapper-fitted model uses the cohort-wide designated first_followup_time retained before response-driven model-frame omission. The transition response is not required on that row, but ID, modeled predictors, time information, and previous state must be complete. Each included patient must contribute at least one usable fitted transition somewhere; later rows are not substituted for a missing starting profile.

For rmsb::blrm() models, SOPs are computed on sampled posterior draws and then summarized. The point estimate is the requested posterior summary of the SOP draws, not a plug-in calculation from summarized model parameters. State-wise medians and interval bounds are not constrained to sum to one across states; use posterior_summary = "mean" when the displayed estimates themselves need to preserve total probability. Draw-level probabilities stored with return_draws = TRUE remain normalized within each draw.

Model Requirements:

For vglm models, only the cumulative-logit family with reverse = TRUE is supported. For orm models, only the logistic family is supported. This is because the package's Markov simulation logic expects higher-numbered states to represent worse outcomes, and uses reverse cumulative probabilities to model the probability of being in state k or worse.

Stratification:

When by is specified, the function:

Computes individual-level SOPs for all patients
Groups results by time, state, and the specified variables
Averages the SOPs within each group

This aggregation preserves heterogeneity across observed strata without creating counterfactual scenarios (unlike G-computation in avg_sops()).

For computing marginal/standardized SOPs (G-computation), use avg_sops() instead, which creates counterfactual datasets and averages over individuals.

Choosing visits

SOPs are calculated one visit at a time. The first value of times is the first transition from the starting state, and each later value is one more transition. The values set the time variable, and any time-varying covariates, for each transition; they do not say how much time passes between transitions. Leaving out visits therefore does not skip time. With daily visits, times = c(1, 7, 14) gives the probabilities after 1, 2, and 3 transitions, labelled 1, 7, and 14, instead of the probabilities on days 7 and 14. To get SOPs at selected visits, request every visit, such as times = 1:14, and keep the rows you need.

When newdata is omitted, each patient starts from their state before the first follow-up visit (first_followup_time of the fitting wrapper), so times must begin at that visit. Each row of a supplied newdata starts from its previous state at the first requested visit, which may be a later visit.

A warning is given when times leaves out visits that occur in the data the model was fitted to, is not in increasing order, or, when newdata is omitted, does not begin at the first follow-up visit. Times after the last fitted visit are allowed. There is no warning when gap_var is supplied, because each transition then uses its own elapsed time, so the visits may be spaced unevenly.

Examples

if (FALSE) { # \dontrun{
# Individual-level SOPs
sops_ind <- sops(fit, newdata = baseline_data, times = 1:30, y_levels = 1:6)

# Stratified by ECOG performance status
sops_ecog <- sops(fit, newdata = baseline_data, times = 1:30, y_levels = 1:6,
                  by = "ecog")

# Stratified by multiple variables
sops_strat <- sops(fit, newdata = baseline_data, times = 1:30, y_levels = 1:6,
                   by = c("ecog", "age_group"))

# Compatible with inferences()
sops_strat |> inferences(method = "mvn", n_draws = 1000)
} # }