Calculate Averaged State Occupation Probabilities (Marginal Effects)

Averages state occupation probabilities over patients' starting states and covariates. With variables = NULL, uses each patient's observed covariates. Supply variables to set the specified covariates to common values before averaging, for example to estimate probabilities under each treatment.

Usage

avg_sops(
  model,
  newdata = NULL,
  variables = NULL,
  by = NULL,
  times,
  y_levels = NULL,
  absorb = 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 data frame of standardization profiles. When supplied, every row supplies a separate starting state and set of covariates. The internal rowid column is regenerated. If NULL, uses data stored by orm_markov(), blrm_markov(), or vglm_markov() and requires exactly one complete designated starting profile per fitted patient.
variables: Optional named list specifying covariates to set before averaging. E.g., list(tx = c(0, 1)) estimates probabilities if every patient received control or treatment. The default, NULL, averages with each patient's observed covariates, including their observed treatment.
by: Optional character vector of additional variables to group by after standardization. Rows missing any grouping value are omitted; an error is raised if none remain.
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.
id_var: Name of the patient ID variable. Required for bootstrap inference and for blrm random-effect prediction. If NULL, defaults to "id"; for wrapper-fitted models, it is inferred from stored metadata. For user-supplied newdata, id_var is required only for blrm random effects. Refit bootstrap inference resamples patients by this column in the data stored on the wrapper-fitted model. For blrm models with include_re = TRUE, model$clusterInfo$name is used before the final "id" fallback when wrapper metadata is absent.
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 precomputed time-basis columns. Registered inline terms such as rms::rcs(time, 4) and rms::lsp(time, c(3, 7)) can be used without it.
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().
...: Additional arguments reserved for future extensions.

Value

A data frame of class markov_avg_sops with columns:

time: Time point
state: State level
(variables): Values requested in variables, when supplied
estimate: Average probability across individuals
conf.low, conf.high, std.error: For blrm fits, posterior uncertainty summaries computed directly from marginalized SOP draws

Details

With variables = NULL, this function averages patients' SOPs using their observed covariates, returning one estimate per time and state. Supply by to average separately within groups.

Supplying variables implements G-computation (standardization):

Counterfactual Creation: For each value in variables, creates a copy of newdata with that variable set to the specified value.
Individual Prediction: Calls sops() to compute individual-level SOPs for each patient under each counterfactual scenario.
Marginalization: Averages individual SOPs across patients within each time-state-treatment combination, yielding population-average (marginal) probabilities.

With variables = list(tx = c(0, 1)), the result represents the expected SOP if the entire population received treatment or control, averaged over the observed covariate distribution. For automatic wrapper-stored standardization, that distribution contains exactly the fitted patients with a complete designated starting profile and at least one usable likelihood transition. A missing response at the designated row alone does not exclude a patient who contributes a later transition. Patients with no usable fitted transition are excluded, and a fitted patient without a complete designated profile is an error. For rmsb::blrm() models, patient-level posterior arrays are chunked and marginalized by draw before summarizing to reduce memory use. 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.

For supported frequentist models, inferences(method = "delta") adds analytical confidence intervals. Its default vcov = "unconditional" accounts for coefficient estimation and sampling the patients used in the average. Use vcov = "conditional" to treat those patients' starting states and covariates as given and account only for coefficient estimation. User-supplied newdata requires the conditional choice. See inferences() for model restrictions.

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{
fit <- vglm_markov(
  ordered(y) ~ rms::rcs(time, 4) + tx + yprev,
  family = VGAM::cumulative(reverse = TRUE, parallel = TRUE),
  data = data,
  id_var = "id"
)

# Wrapper-fitted models can use their stored designated starting profiles.
observed <- avg_sops(fit, times = 1:60, y_levels = 1:6, absorb = 6)

# Set treatment to each requested value before averaging.
result <- avg_sops(
  model = fit,
  variables = list(tx = c(0, 1)),
  times = 1:60,
  y_levels = 1:6,
  absorb = 6
) |> inferences(method = "mvn", n_draws = 500)

# Refit bootstrap resamples the patients stored on the fitted model.
result_boot <- avg_sops(
  model = fit,
  variables = list(tx = c(0, 1)),
  times = 1:60,
  y_levels = 1:6,
  absorb = 6
) |> inferences(method = "fwb", n_draws = 500)
} # }