Inference for State Occupation Probabilities and Time Summaries

Adds confidence intervals to SOP, average-time, and comparison objects using analytical delta, simulation-based, or bootstrap methods. The default is multivariate-normal coefficient simulation. Objects produced from rmsb::blrm() models already contain posterior uncertainty, so non-delta calls return them unchanged. method = "delta" is frequentist-only and errors for blrm objects.

Usage

inferences(
  x,
  method = "mvn",
  n_draws = 1000,
  vcov = NULL,
  cluster = NULL,
  workers = NULL,
  seed = NULL,
  conf_level = 0.95,
  conf_type = "auto",
  null = "auto",
  return_draws = TRUE,
  update_datadist = TRUE,
  use_coefstart = FALSE
)

Arguments

x

A markov_avg_sops object from avg_sops(), a markov_sops object from sops(), a markov_avg_comparisons object from avg_comparisons(), or a markov_avg_time object from avg_time().

method

Character. Inference method:

"mvn": Multivariate-normal coefficient draws.
"delta": Deterministic analytical delta-method inference.
"score_bootstrap": One-step score perturbation using fixed exponential cluster weights.
"bootstrap": Ordinary patient-level refit bootstrap.
"fwb": Fractional weighted refits using fixed exponential patient weights.

n_draws

Number of simulation draws (for simulation) or bootstrap iterations (for bootstrap). Default is 1000. For blrm SOP objects this argument is ignored; rerun sops()/avg_sops() with n_draws and seed to change posterior draws.

vcov

For method = "delta", "conditional" treats prediction or standardization covariates as fixed; "unconditional" includes their sampling variability and its cross-term with coefficient estimation. NULL (the default) selects unconditional inference for avg_sops() and avg_comparisons(), and conditional inference for sops(). Individual sops() results support only conditional inference. A complete, named coefficient covariance matrix selects conditional inference and overrides the model covariance. Character choices are available only for delta inference; other methods retain their existing matrix/NULL behavior. Averages whose variables set the starting state (the previous-state variable, usually yprev) support only conditional variance; see Details.

cluster

Optional patient-cluster specification. For analytical conditional and unconditional variance estimates, supply a vector aligned with the fitting rows or a one-sided formula selecting a stored fitting-data column. Otherwise the model's stored id_var is used; fitting rows are never treated as implicit clusters. For unconditional inference, the cluster labels must match the stored starting-profile patient IDs exactly. For score bootstrap with orm, this is the row-aligned cluster vector.

workers

Number of parallel workers. If NULL (default) or 1, uses sequential processing. If > 1, uses parallel processing with that many workers.

seed

Optional integer seed. The caller's complete RNG state is restored on exit.

conf_level

Confidence level for intervals (default 0.95). For blrm SOP objects this argument is ignored; rerun sops()/avg_sops() with conf_level to change posterior intervals.

conf_type

Type of frequentist confidence interval. "auto" uses percentile intervals for draw-based methods, componentwise logit-delta intervals for SOP probabilities, and identity-scale Wald intervals for analytical comparisons and average times:

"auto": Method-appropriate default.
"perc": Percentile-based intervals from the simulation distribution.
"wald": Uses simulation standard errors with normal quantiles.
"logit": Componentwise logit-delta limits for probabilities; available only with method = "delta" on SOP objects.

null

The value each estimate is tested against. A test adds three columns: the Wald statistic, (estimate - null) / std.error; its two-sided normal p.value; and the s.value, \(-\log_2(p)\), which expresses the evidence against the null in bits.

"auto" (the default): tests avg_comparisons() results against "no effect", which is 0 for comparison = "difference" and 1 for comparison = "ratio". SOPs and average times have no natural null value, so they are not tested.
A single finite number: tests every estimate against that value. Ratio comparisons cannot be tested against 0.
NULL: no test.

Ratios are tested on the ratio scale, using the ratio's own standard error. Bayesian blrm results are never tested; supplying a number for them gives a warning.

return_draws

Logical. If TRUE, stores individual simulation or bootstrap draws as an attribute. Extract with get_draws(). Delta inference never stores draws. Default is TRUE.

update_datadist

Logical. Whether to update datadist for rms models during bootstrap. Default is TRUE.

use_coefstart

Logical. If TRUE, vglm refits start from the original fit's coefficients, which can speed up convergence. Default is FALSE. See "Starting values in refits" in Details.

Value

The input object with added columns:

conf.low: Lower confidence bound
conf.high: Upper confidence bound
std.error: Standard error from analytical delta propagation, simulation, or bootstrap

When a null value is tested (by default, for avg_comparisons() results), the object also has statistic, p.value, and s.value columns; see null.

For simulation and bootstrap methods, return_draws = TRUE also stores a "draws" attribute containing all individual draws. Delta results instead store a low-rank "analytical" attribute and never store draws. For ungrouped sops() objects evaluated on the stored prediction cohort, score-bootstrap and FWB draws include the draw-specific score_weight or fwb_weight column.

Details

Analytical Delta Method

method = "delta" differentiates the first-order full proportional-odds SOP recursion on the model's complete raw-coefficient scale. With vcov = "conditional", covariance is propagated as \(G V_\beta G^\top\), where \(\beta\) contains all model coefficients, \(G\) contains derivatives of the reported estimates with respect to \(\beta^\top\), evaluated at the fitted coefficients \(\hat\beta\), and \(V_\beta\) is their estimated covariance matrix. Rows of \(G\) correspond to estimates and columns to coefficients. For a simple patient average, \(G = \bar G = n^{-1}\sum_{i=1}^n G_i\), where \(G_i\) contains patient \(i\)'s prediction derivatives and \(n\) is the number of patients. An explicit covariance matrix passed as vcov overrides the model covariance for conditional inference.

vcov = "conditional" accounts for uncertainty in the estimated model coefficients while treating the patients' starting states and covariates used for prediction or averaging as given. vcov = "unconditional" is the default for averaged results and is available only for the patients stored with the fitted model. It also accounts for variation in which patients are sampled from the population, including its association with coefficient estimation from those same patients. Thus it estimates uncertainty in the population average. The calculation combines each patient's deviation from the average prediction with their contribution to coefficient estimation, then divides the sample covariance of these combined contributions by the number of patients. The sample covariance uses the usual divisor of one less than the number of patients; the model's type and cadjust corrections are not applied again. User-supplied prediction cohorts require vcov = "conditional"; unavailable unconditional inference errors rather than silently falling back. A custom covariance matrix always selects conditional inference because unconditional inference requires fitted-model scores and sensitivity. Every included patient must contribute at least one usable likelihood transition and have exactly one complete model-stored starting profile. The starting profile may have a missing first transition response; its ID, predictors, scheduled starting time, and previous state must be observed. Weighted and penalized ORM fits do not support unconditional variance because their score/sensitivity contracts have not been established. They remain eligible for conditional delta inference when their fitted covariance is valid.

Analytical covariance is patient-cluster robust only when valid patient IDs are supplied explicitly or retained as the fitted model's id_var. Observation rows are not an automatic clustering fallback, and multiway or otherwise non-patient cluster structures are not currently supported. Patient-cluster robustness protects variance estimation against within-patient score correlation; it does not remove bias from a misspecified transition model, informative observation process, or Markov assumption. blrm models are not supported by the analytical method.

Simulation Method

The simulation method works as follows:

Extract the fitted coefficient vector \(\hat\beta\) and its estimated covariance matrix \(V_\beta\).
Generate n_draws coefficient vectors using the selected method:
- method = "mvn": Multivariate-normal draws from \(\mathcal{N}(\hat\beta, V_\beta)\).
- method = "score_bootstrap": one-step score perturbation with cluster-level exponential multipliers.
For each draw, replace model coefficients and compute SOPs
Compute confidence intervals from the empirical distribution

Works for both individual-level (sops()) and averaged (avg_sops()) SOPs for method = "mvn"
method = "score_bootstrap" supports avg_sops() and sops() with robcov_vglm models and with orm models when cluster is supplied. When the prediction rows are the stored patients, the same cluster-level weights used for the score perturbation are used for every empirical averaging step. With by, weights are normalized within each subgroup; without by, they are normalized over the full set of stored patients.

Bootstrap Method

The bootstrap method refits the model on bootstrap-weighted data:

Sample patient IDs with replacement (engine = "standard") or draw exponential patient weights normalized to mean 1 (engine = "fwb")
Refit model on bootstrap sample (handles missing states through releveling)
Compute SOPs using G-computation
Compute percentile-based confidence intervals Bootstrap refits always resample the full longitudinal data (all time points) that orm_markov() or vglm_markov() stored on the fitted model. These are the rows used for fitting, after any subset, together with their fitting weights, so every refit uses the same data rules as the original fit. Models fitted without these wrappers cannot be used for refit bootstraps. The prediction rows are separate: supplying newdata to sops() or avg_sops() changes which patient profiles are predicted, not which data are refitted. Standard bootstrap is intentionally limited to marginal avg_sops() objects because ordinary resampling can drop outcome-state support needed by fixed individual prediction rows. Fractional weighted bootstrap keeps every row in the refit data with positive patient-level weights, so it is available for individual sops() objects. For FWB, the same patient weights used for refitting are used for marginal or subgroup averaging. Ungrouped stored-data sops() draws expose the patient weight as fwb_weight; grouped sops() and avg_sops() use the weights internally and return already averaged draws. Prediction data for such ungrouped draw output must not already contain a column named fwb_weight or score_weight, because those names are reserved for bootstrap draw weights.

For marginal avg_sops() objects built from user-supplied newdata, the supplied patients' starting states and covariates are treated as given. This also applies to sops(newdata = ...). Score bootstrap and FWB use the data stored on the fitted model for coefficient or refit uncertainty, but do not attach or apply draw weights to the supplied prediction profiles because those rows cannot be assumed to align with the bootstrap clusters.

Starting values in refits

Starting values are the coefficient values from which the fitting algorithm begins its search. They can come from use_coefstart = TRUE (the original fit's coefficients) or from a coefstart argument you passed to vglm_markov(). A refit uses them only when its bootstrap sample contains every outcome state. When a state is missing, the refit has fewer intercepts than the original model, so the starting values no longer fit and the model is refitted without any. If a refit with starting values fails, for example because a covariate level is absent from the sample, it is retried once without them.

Choosing conditional or unconditional analytical variance changes the standard errors and confidence intervals, not the point estimates.

Setting the Starting State

When variables in avg_sops(), avg_time(), or avg_comparisons() sets the previous-state variable, for example variables = list(yprev = c("1", "2")), every patient starts in the set state in each scenario. Real-time summaries then start from that state rather than from the patients' observed starting states. Such results support only conditional variance: use method = "delta" with vcov = "conditional" (or a coefficient covariance matrix), or method = "mvn". Unconditional delta inference and the score-bootstrap, bootstrap, and FWB methods error, because they also account for which patients were sampled, and that calculation is not supported once the starting state is set.

Examples

if (FALSE) { # \dontrun{
# Step 1: Fit model with stored data and cluster-robust vcov
fit_robust <- vglm_markov(
  ordered(y) ~ time + tx + yprev,
  family = cumulative(reverse = TRUE, parallel = TRUE),
  data = data,
  id_var = "id"
)

# Step 2a: Simulation inference
result_sim <- avg_sops(
  model = fit_robust,
  variables = list(tx = c(0, 1)),
  times = 1:60,
  y_levels = factor(1:6),
  absorb = "6"
) |>
  inferences(method = "mvn", n_draws = 1000)

# Step 2a-bis: Score bootstrap simulation (requires robcov_vglm)
result_score <- avg_sops(
  model = fit_robust,
  variables = list(tx = c(0, 1)),
  times = 1:60,
  y_levels = factor(1:6),
  absorb = "6"
) |>
  inferences(method = "score_bootstrap", n_draws = 1000)

# orm_markov() stores the package-owned full robust covariance matrix.
dd <- rms::datadist(data)
options(datadist = "dd")
fit_orm <- orm_markov(y ~ time + tx + yprev, data = data, id_var = "id")

avg_orm <- avg_sops(
  model = fit_orm,
  variables = list(tx = c(0, 1)),
  times = 1:60,
  y_levels = factor(1:6),
  absorb = "6"
)
result_orm <- avg_orm |>
  inferences(method = "mvn", n_draws = 1000)

result_orm_score <- avg_orm |>
  inferences(
    method = "score_bootstrap",
    cluster = data$id,
    n_draws = 1000
  )

# Step 2b: Bootstrap inference reuses stored full data
result_boot <- avg_sops(
  model = fit_robust,
  variables = list(tx = c(0, 1)),
  times = 1:60,
  y_levels = factor(1:6),
  absorb = "6"
) |>
  inferences(method = "bootstrap", n_draws = 500)

# Individual-level SOPs with simulation inference
ind_result <- sops(
  model = fit_robust,
  times = 1:60,
  y_levels = factor(1:6),
  absorb = "6") |>
  inferences(n_draws = 500)

# Individual-level refit bootstrap uses FWB.
ind_boot <- sops(
  model = fit_robust,
  times = 1:60,
  y_levels = factor(1:6),
  absorb = "6") |>
  inferences(method = "fwb", n_draws = 500)

# Extract draws for custom analyses
get_draws(ind_result)

# Compute treatment effect on time in state
library(dplyr)
library(tidyr)
treatment_effect <- draws |>
  filter(state == 1) |>
  group_by(draw_id, tx) |>
  summarise(total_time = sum(draw), .groups = "drop") |>
  pivot_wider(names_from = tx, values_from = total_time) |>
  mutate(effect = `1` - `0`)
quantile(treatment_effect$effect, c(0.025, 0.5, 0.975))
} # }