Multivariate Sequential Ordinal Models

• Multivariate sequential ordinal models are statistical frameworks for ordered categorical data that integrate cross-response and sequential dependence using various constructions.
• They encompass methods such as hazard-style (continuation-ratio) formulations, dynamic multivariate time-series, and longitudinal latent-process models to address different types of ordinal dependence.
• These models support applications from invasive species abundance to unemployment and athlete wellness, often simplifying computation through binary likelihood recoding and advanced estimation strategies.

Multivariate sequential ordinal models are statistical models for ordered categorical data in which both cross-response dependence and an ordered sequence must be represented. Across the recent literature, the expression does not denote a single canonical specification. It is used for at least three distinct constructions: sequential category models based on continuation or hazard-style probabilities, dynamic multivariate ordinal time-series or panel models that specify a conditional law given the past, and longitudinal repeated-measures models with serial latent dependence or lagged covariate effects but without explicit observed-state transitions (Hosack et al., 2 Sep 2025, Nalpantidi et al., 13 Feb 2026, Nikoloulopoulos et al., 2016, Schliep et al., 2020, Vana-Gür, 2024).

1. Conceptual scope and terminology

The term ordinal refers to outcomes that are categorical with a natural ordering. The term multivariate can refer to several jointly modeled ordinal responses, as in multivariate ordinal time series or panel data, but in some sequential-ordinal formulations it refers more narrowly to the multicategory representation of a single ordinal response as a multinomial object recast into multiple binary contributions. The term sequential is equally overloaded. In one strand it refers to category progression through conditional stopping probabilities; in another it refers to temporal sequencing through conditioning on a filtration such as $\mathcal F_{t-1}$; in a third it refers only to repeated measurements with serially structured latent variables (Hosack et al., 2 Sep 2025, Nalpantidi et al., 13 Feb 2026, Vana-Gür, 2024).

This terminological heterogeneity is methodologically consequential. A dynamic spatio-temporal sequential ordinal model can be “sequential” because it parameterizes

$$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$

whereas a multivariate ordinal time-series model can be “sequential” because it specifies the conditional joint distribution of $\mathbf Z_t$ given the full past $\mathcal F_{t-1}$, and a repeated-measures ordinal regression can be called sequential even when the observed categories do not enter through a Markov transition law (Hosack et al., 2 Sep 2025, Nalpantidi et al., 13 Feb 2026, Schliep et al., 2020).

Several adjacent literatures are relevant but not identical. A Bayesian group sequential schema for ordinal endpoints concerns interim monitoring in a two-arm trial with a single ordinal endpoint, not multivariate ordinal dependence (Zhong et al., 2021). Bayesian nonparametric multivariate ordinal regression and hierarchical marginal models with latent uncertainty provide rich multivariate ordinal structure, but they do not develop sequential or temporal ordinal dependence in the sense of time-series or transition modeling (DeYoreo et al., 2014, Colombi et al., 2016). A Bayesian hierarchical multivariate probit latent class model for mixed dichotomous and ordinal tests is multivariate and ordered, but not sequential in the continuation-ratio or longitudinal-transition sense (Cerullo et al., 2021).

2. Sequential ordinal likelihoods and continuation-ratio constructions

In the continuation-ratio or hazard-style tradition, the central object is the conditional probability of stopping at category $c$ given that all lower categories have been bypassed: $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right) = \delta_{ic}(t)=g^{-1}(\eta_{ic}(t)),\qquad c=1,\dots,q.$$ The implied unconditional category probabilities are

$$\pi_{ic}(t)=\delta_{ic}(t)\prod_{\tilde c=1}^{c-1}\left(1-\delta_{i\tilde c}(t)\right),\qquad c=1,\dots,q,$$

with the top category probability

$$P(z_i(t)=C)=\prod_{\tilde c=1}^q (1-\delta_{i\tilde c}(t)).$$

This factorization yields the exact sequential or hazard-style interpretation: an observation “survives” lower categories and then stops at $c$ (Hosack et al., 2 Sep 2025).

A major computational consequence is that the sequential ordinal likelihood can be rewritten as a product of Bernoulli terms with shared latent effects. For each ordinal observation, one defines binary pseudo-observations corresponding to the run of rejected lower categories and the accepted stopping category. Stacking across observations produces a likelihood that is formally binary conditional on regression coefficients and latent random effects. This is the mechanism by which dynamic generalized linear model machinery and INLA can be used for high-dimensional spatio-temporal ordinal data (Hosack et al., 2 Sep 2025).

The sequential formulation differs sharply from cumulative ordinal models. In cumulative models, threshold parameters must satisfy ordering constraints across categories. In the sequential formulation developed for spatio-temporal invasive-weed abundance, the model does not require strictly ordered cutpoints or thresholds across categories, which simplifies Bayesian estimation in high-dimensional latent-variable settings. With complementary log-log link and global effects, the model also admits a grouped-Cox or proportional-hazards interpretation and preserves collapsibility under aggregation of ordinal categories (Hosack et al., 2 Sep 2025).

An important limitation is that this use of “multivariate” is narrower than a full joint multiresponse ordinal model. In the invasive-weed application, four species are modeled separately; there is no full joint cross-species ordinal distribution, no multivariate cross-species covariance matrix, and no shared latent process across species in the fitted application. The multivariate contribution lies primarily in the multicategory ordinal likelihood, its binary recoding, and shared latent spatio-temporal random effects within an ordinal response (Hosack et al., 2 Sep 2025).

3. Dynamic multivariate ordinal time series and panel Markov models

The most direct realization of a multivariate sequential ordinal model in the temporal sense specifies the conditional distribution of an ordinal vector given the past. For a $K$-variate ordinal time series,

$$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$0

with filtration

$$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$1

the dependence structure can be decomposed into serial dependence, contemporaneous dependence, and cross-series lag dependence. One approach models each conditional marginal $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$2 by an ordinal autoregressive logit model of order $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$3,

$$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$4

while coupling the contemporaneous pair $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$5 by a copula conditional on $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$6 (Nalpantidi et al., 13 Feb 2026).

In this copula-based construction, temporal effects are pushed into the conditional marginals and lag-0 dependence is handled by the copula. The discrete joint pmf for a pair is obtained by a rectangle probability, that is, a finite difference of the copula cdf. Higher-dimensional inference is then approximated through conditional pairwise likelihood rather than a full $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$7-variate discrete copula likelihood. The proposed estimation strategy maximizes each bivariate model independently and synthesizes duplicated parameter estimates by a weighted mean derived from the Hessian matrices of the pairwise models (Nalpantidi et al., 13 Feb 2026).

A second, more fully specified panel formulation models each ordinal component with a copula-based first-order Markov structure over time and then joins the conditional distributions at each time point through a multivariate $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$8 copula. For subject $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right),$$9, outcome component $\mathbf Z_t$0, and time $\mathbf Z_t$1, the observed response is

$\mathbf Z_t$2

with ordinal margins generated by a cumulative-link latent-threshold construction. Serial dependence is introduced through a bivariate copula linking adjacent observations within each component, and contemporaneous dependence across components at time $\mathbf Z_t$3 is introduced by a multivariate $\mathbf Z_t$4 copula applied to the conditional distributions $\mathbf Z_t$5 (Nikoloulopoulos et al., 2016).

Conceptually, this yields

$\mathbf Z_t$6

where $\mathbf Z_t$7, but the conditional law is not a fully general vector Markov chain. It is assembled from $\mathbf Z_t$8 componentwise first-order Markov transitions and a time-$\mathbf Z_t$9 multivariate $\mathcal F_{t-1}$0-copula coupling. This separation of within-outcome persistence and contemporaneous cross-outcome dependence is one of the defining features of the model (Nikoloulopoulos et al., 2016).

These time-series and panel copula models are the closest direct matches to the phrase multivariate sequential ordinal model. They are multivariate because several ordinal outcomes are modeled jointly, sequential because the conditional law at time $\mathcal F_{t-1}$1 depends on previous observations, and ordinal because the observed states are ordered categories rather than continuous responses (Nalpantidi et al., 13 Feb 2026, Nikoloulopoulos et al., 2016).

4. Longitudinal multivariate ordinal models with serial latent structure

A separate but closely related literature models repeated multivariate ordinal outcomes through latent serial processes rather than explicit observed-state transitions. In a three-dimensional panel setting with subjects $\mathcal F_{t-1}$2, outcomes $\mathcal F_{t-1}$3, and times $\mathcal F_{t-1}$4, one formulation uses latent variables

$\mathcal F_{t-1}$5

with temporal persistence introduced by an autoregressive model on the latent errors,

$\mathcal F_{t-1}$6

In the implemented model, $\mathcal F_{t-1}$7, so each outcome has its own persistence parameter $\mathcal F_{t-1}$8, while contemporaneous residual dependence at a given time is encoded by correlations $\mathcal F_{t-1}$9 in $c$0 (Vana-Gür, 2024).

This is sequential only in a latent-process sense. The current observed ordinal response is not modeled by conditioning directly on the previous observed ordinal response. Instead, temporal dependence is induced by AR(1) structure in the latent errors. The paper therefore distinguishes between correlations at a single point in time and persistence over time, and treats the model as a multivariate longitudinal or panel ordinal regression with serial dependence rather than a classical response-driven ordinal Markov model (Vana-Gür, 2024).

An even more explicit boundary case is a multivariate ordinal longitudinal latent factor distributed lag model for athlete wellness. There, six daily ordinal wellness outcomes are jointly modeled through one or two latent Gaussian factors, and the sequential component enters through lagged training and recovery covariates: $c$1 or, in the two-factor version,

$c$2

The ordinal observation model is cumulative probit,

$c$3

with thresholds mapping latent Gaussian responses to five ordered categories (Schliep et al., 2020).

This formulation is multivariate, ordinal, and longitudinal, but it is not a Markov or transition model. There is no autoregressive term in $c$4 and no direct dependence of the form $c$5. The sequential aspect comes from the distributed history of workload and recovery over prior days. The paper explicitly frames the model as a distributed lag longitudinal latent factor model rather than a hidden Markov or ordinal transition model (Schliep et al., 2020).

5. Dependence structures, identifiability, and computation

Across these model classes, dependence can be introduced through several technically distinct mechanisms. Copula-based time-series models distinguish serial dependence within each response component from contemporaneous dependence across components and, in the pairwise-copula construction, from cross-series lag dependence encoded in the marginal regressions (Nalpantidi et al., 13 Feb 2026, Nikoloulopoulos et al., 2016). Latent-factor longitudinal models induce cross-metric dependence through shared latent wellness factors, while hierarchical subject-specific lag coefficients represent between-subject heterogeneity with partial pooling (Schliep et al., 2020). Latent-AR panel models distinguish persistence parameters $c$6 from same-time cross-outcome correlations $c$7 (Vana-Gür, 2024).

Identifiability constraints differ sharply by formulation. In cumulative ordinal models, ordered thresholds are central. In the athlete-monitoring model, the first threshold is fixed and factor loadings are anchored by fixing the first loading to $c$8; threshold spacings are reparameterized on a log scale to enforce ordering, and Dirichlet priors are placed on variance partitions for scale identification (Schliep et al., 2020). In the panel AR(1) model, the diagonal entries of the latent-error covariance are fixed to $c$9 for identifiability of the marginal ordinal models (Vana-Gür, 2024). In Bayesian nonparametric multivariate ordinal regression, fixed cut-offs are used rather than estimated thresholds, and kernel-level identifiability is established provided each ordinal response has more than two categories (DeYoreo et al., 2014). By contrast, the sequential ordinal model based on continuation probabilities explicitly avoids threshold-ordering constraints across categories (Hosack et al., 2 Sep 2025).

Computation is correspondingly heterogeneous. The sequential spatio-temporal model exploits the binary recoding of the ordinal likelihood so that estimation and prediction can proceed using INLA for latent Gaussian binary models; the shared latent effects are assigned through the mapping matrix and, in SPDE implementations, the projection matrix (Hosack et al., 2 Sep 2025). The copula-based multivariate ordinal panel Markov model uses maximum simulated likelihood with multivariate $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right) = \delta_{ic}(t)=g^{-1}(\eta_{ic}(t)),\qquad c=1,\dots,q.$$0 rectangle probabilities evaluated by randomized quasi Monte Carlo methods of Genz and Bretz, and practical estimation proceeds through a staged strategy or IFM when full optimization is burdensome (Nikoloulopoulos et al., 2016). The multivariate ordinal time-series copula model uses conditional pairwise likelihood, separate pairwise maximization, Hessian-based weighted synthesis, and Godambe or sandwich covariance estimation (Nalpantidi et al., 13 Feb 2026). The multivariate repeated-measurements model uses pairwise composite likelihood with optional lag truncation, Godambe information, and the R package mvordflex; the paper reports that pairwise likelihood is competitive in terms of computation time and, in one simulation, much faster than a Bayesian rstan implementation (Vana-Gür, 2024). The athlete distributed-lag latent-factor model uses fully Bayesian MCMC with a hybrid Metropolis-within-Gibbs algorithm, Albert–Chib latent-variable augmentation, and Metropolis updates for ordered thresholds (Schliep et al., 2020).

6. Applications, misconceptions, and boundaries of the concept

The empirical range of these models is broad. Sequential likelihood models have been used for foliage cover abundance of four actively managed invasive alien species observed over space and time (Hosack et al., 2 Sep 2025). Copula-based multivariate ordinal time-series models have been illustrated on quarterly unemployment states for six EU countries (Nalpantidi et al., 13 Feb 2026). Copula-based panel Markov models have been used for seven well-being dimensions across survey waves (Nikoloulopoulos et al., 2016). Longitudinal latent-process models have been applied to athlete daily wellness, training, and recovery data across two Major League Soccer seasons and to panel data containing firm failure and credit ratings from S&P and Moody’s for US listed companies (Schliep et al., 2020, Vana-Gür, 2024).

Several recurring misconceptions follow directly from the diversity of the literature. First, a repeated multivariate ordinal model is not automatically an ordinal transition model. The distributed-lag athlete model and the AR(1) latent-error panel model both contain temporal structure, yet neither specifies explicit transition probabilities of the form $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right) = \delta_{ic}(t)=g^{-1}(\eta_{ic}(t)),\qquad c=1,\dots,q.$$1 (Schliep et al., 2020, Vana-Gür, 2024). Second, a model called “multivariate sequential ordinal” may not be multivariate in the sense of a fully joint covariance structure across several ordinal responses. The spatio-temporal sequential ordinal model for invasive weeds is multivariate primarily in the multinomial or multicategory sense and models species separately in the application (Hosack et al., 2 Sep 2025). Third, “sequential” can refer to trial monitoring rather than data-generating dependence; Bayesian group sequential designs for ordinal endpoints are sequential in the interim-analysis sense and are explicitly univariate with respect to the modeled endpoint (Zhong et al., 2021).

A further boundary concerns adjacent multivariate ordinal models that are not sequential. Hierarchical marginal models with latent uncertainty components separate aware from uncertain responses in multivariate ordinal ratings and allow local, global, continuation, or reverse continuation logits as parameterizations, but the core contribution is a latent uncertainty mixture rather than a sequential category-generation mechanism (Colombi et al., 2016). Bayesian nonparametric multivariate ordinal regression models the joint distribution of latent responses and covariates through a Dirichlet process mixture and induces flexible multivariate ordinal regression surfaces, but does not explicitly model sequential, longitudinal, time-series, or Markov dependence (DeYoreo et al., 2014). A latent class multivariate ordered probit model for diagnostic-test meta-analysis supplies multivariate thresholds, latent Gaussian dependence, and hierarchical pooling, but it is not sequential in either the continuation-ratio or time-series sense (Cerullo et al., 2021).

Taken together, these papers suggest that precision in technical labeling is essential. A plausible implication is that the phrase multivariate sequential ordinal model should be unpacked into one of three more specific meanings whenever possible: a continuation-ratio or hazard-style ordinal likelihood, a dynamic multivariate ordinal conditional law for $$P\!\left(z_i(t)=c \mid z_i(t)\ge c,\eta_{ic}(t)\right) = \delta_{ic}(t)=g^{-1}(\eta_{ic}(t)),\qquad c=1,\dots,q.$$2, or a longitudinal multivariate ordinal model with serial latent dependence. Without that distinction, closely related but substantively different constructions are easily conflated (Hosack et al., 2 Sep 2025, Nalpantidi et al., 13 Feb 2026, Nikoloulopoulos et al., 2016, Vana-Gür, 2024).