Fit a VGAM Markov Model

vglm_markov() is the recommended VGAM fitting entrypoint for mostr SOP workflows. It follows VGAM::vglm() closely, stores separate likelihood, refit, and designated-start profile data for downstream SOP inference, and marks the fit so the SOP prediction code can use package-native Markov prediction. For inline registered RMS spline terms such as rcs(time, 4), lsp(time, 5), or their namespace-qualified forms, the internal VGAM assignment metadata is split by generated spline column. This makes it possible to give individual spline basis columns separate constraint matrices for partial proportional odds models. Inline splines may be used for time or for numeric previous-state effects such as rcs(yprev, 6).

Usage

vglm_markov(
  formula,
  family = stop("argument 'family' needs to be assigned"),
  data = list(),
  id_var = NULL,
  type = c("HC0", "HC1"),
  cadjust = NULL,
  weights = NULL,
  subset = NULL,
  na.action,
  etastart = NULL,
  mustart = NULL,
  coefstart = NULL,
  control = VGAM::vglm.control(...),
  offset = NULL,
  model = FALSE,
  x.arg = TRUE,
  y.arg = TRUE,
  contrasts = NULL,
  constraints = NULL,
  extra = list(),
  form2 = NULL,
  qr.arg = TRUE,
  smart = TRUE,
  time_var = "time",
  first_followup_time = NULL,
  ...
)

Arguments

formula

a symbolic description of the model to be fit. The RHS of the formula is applied to each linear predictor. The effect of different variables in each linear predictor can be controlled by specifying constraint matrices—see constraints below.

family

A VGAM family object, e.g. VGAM::cumulative().

data

an optional data frame containing the variables in the model. By default the variables are taken from environment(formula), typically the environment from which vglm is called.

id_var

Optional character scalar naming the patient or cluster ID column in data. When supplied, robcov_vglm() is applied automatically.

type

HC correction passed to robcov_vglm() when id_var is supplied: "HC0" (the default) or "HC1".

cadjust

Optional logical cluster correction passed to robcov_vglm() when id_var is supplied. NULL uses the clustered default, which applies correction \(n / (n - 1)\), where \(n\) is the number of patient clusters.

weights

Optional positive prior weights. A vector supplies one weight per observation. For families with multiple responses, a matrix can supply one column per response; a vector is recycled across response columns. See VGAM::vglm() for family-specific weighting behavior. Weights are evaluated once when the model is fitted, so they can be a column of data or a vector defined outside it; they must be finite and non-negative. Each row keeps its weight in the data stored on the fit, so bootstrap refits resample rows together with their weights.

subset

an optional logical vector specifying a subset of observations to be used in the fitting process.

na.action

a function which indicates what should happen when the data contain NAs. The default is set by the na.action setting of options, and is na.fail if that is unset. The “factory-fresh” default is na.omit which is known as complete case analysis and applied to both sides of the formula.

etastart

Optional starting linear predictors: a matrix with one row per observation and \(L\) columns, where \(L\) is the number of linear predictors (\(K - 1\) for an ordinal model with \(K\) states). A vector is allowed when \(L = 1\).

mustart

Optional starting fitted values, with the same structure as fitted(fit). If a matrix, it must have one row per observation. Whether these values can initialize all linear predictors depends on the family; many families do not use this argument.

coefstart

optional starting values for the coefficient vector. The length and order must match that of coef(fit).

control

a list of parameters for controlling the fitting process. See vglm.control for details.

offset

Not supported. Supplying offsets raises an error.

model

a logical value indicating whether the model frame should be assigned in the model slot.

x.arg, y.arg

logical values indicating whether the LM matrix and response vector/matrix used in the fitting process should be assigned in the x and y slots. Note that the model matrix is the LM matrix; to get the VGLM matrix type model.matrix(vglmfit) where vglmfit is a vglm object.

contrasts

an optional list. See the contrasts.arg of model.matrix.default.

constraints

Optional VGAM constraints list. For inline registered RMS basis terms, names should match the column-level constraint names in a full proportional odds fit returned by vglm_markov().

extra

an optional list with any extra information that might be needed by the VGAM family function.

form2

the second (optional) formula. If argument xij is used (see vglm.control) then form2 needs to have all terms in the model. Also, some VGAM family functions such as micmen use this argument to input the regressor variable. If given, the slots @Xm2 and @Ym2 may be assigned. Note that smart prediction applies to terms in form2 too.

qr.arg

logical value indicating whether the slot qr, which returns the QR decomposition of the VLM model matrix, is returned on the object.

smart

logical value indicating whether smart prediction (smartpred) will be used.

time_var

Character scalar naming the modeled time column used to identify the designated starting-profile row.

first_followup_time

First scheduled post-baseline outcome time used to select the starting-profile row. For numeric time, NULL uses 1. Numeric schedules may not contain values below 1, and an explicit value must be the earliest observed time. Factor and character time require an explicit matching value.

...

further arguments passed into vglm.control.

Value

A fitted S4 vglm object with internal Markov marker attributes, or a robcov_vglm object when id_var is supplied.

Details

vglm_markov() also makes the registered rms formula helpers rcs(), lsp(), and %ia% available while building model frames and prediction design matrices. This lets formulas use time %ia% yprev for a linear-time by previous-state interaction without requiring library(rms) in the user's session.

Non-spline terms, factor terms, and workflows that use explicit precomputed basis columns keep ordinary VGAM assignment behavior. When id_var is supplied, the returned object is wrapped with robcov_vglm() using that patient/cluster ID. When id_var is omitted, a plain S4 vglm object is returned; if data contains an id column, a warning reminds the user that automatic cluster-robust covariance and FWB refits need id_var.

Offsets are not supported in mostr SOP workflows, so vglm_markov() rejects both offset() terms and the offset argument.

Starting profiles are retained before response-driven model-frame omission. Every fitted patient must have exactly one complete profile at the cohort-wide first_followup_time, including ID, model predictors, time, and the previous state; the transition response on that row may be missing. A patient is included only when at least one usable likelihood transition is fitted somewhere. Patients with no fitted transition are excluded, and a fitted patient without a complete designated profile is an error. Later likelihood rows are not substituted for that profile.

Examples

if (FALSE) { # \dontrun{
library(VGAM)
library(rms)

fit_po <- vglm_markov(
  ordered(y) ~ rcs(time, 4) + tx + age + yprev,
  family = cumulative(reverse = TRUE, parallel = TRUE),
  data = data,
  id_var = "id"
)

fit_numeric_prev <- vglm_markov(
  ordered(y) ~ rcs(time, 4) + tx + age + rcs(yprev, 6),
  family = cumulative(reverse = TRUE, parallel = TRUE),
  data = prepare_markov_data(data, factor_previous = FALSE)
)

fit_linear_time_by_prev <- vglm_markov(
  ordered(y) ~ rcs(time, 6) + tx + age + yprev + time %ia% yprev,
  family = cumulative(reverse = TRUE, parallel = TRUE),
  data = data
)

# In contrast, rcs(time, 6) %ia% yprev follows rms semantics and interacts
# the spline basis columns with yprev when yprev is categorical.
fit_spline_basis_by_prev <- vglm_markov(
  ordered(y) ~ rcs(time, 6) + tx + age + yprev + rcs(time, 6) %ia% yprev,
  family = cumulative(reverse = TRUE, parallel = TRUE),
  data = data
)

constraints <- fit_po$vglm_fit@constraints
linear_time_term <- grep("^rcs\\(time, 4\\).*time$", names(constraints),
  value = TRUE
)[[1]]
n_thresholds <- nrow(constraints[[linear_time_term]])
constraints[[linear_time_term]] <- cbind(
  PO_effect = rep(1, n_thresholds),
  linear_deviation = seq_len(n_thresholds) - 1
)

fit_ppo <- vglm_markov(
  ordered(y) ~ rcs(time, 4) + tx + age + yprev,
  family = cumulative(reverse = TRUE, parallel = FALSE),
  data = data,
  constraints = constraints
)
} # }

Usage

Arguments

Value

Details

See also

Examples