Bayesian Unit-level Modeling of Categorical Survey Data with a Longitudinal Design

In the context of survey data applications, modeling individual records (i.e., “unit-level modeling”) rather than area-level aggregates is often necessary. For ordinal-valued responses, modeling the individual survey responses allows us to capture dependence among the response categories, which is not possible in area-level models that model group averages. Likewise, for longitudinal studies, modeling at the unit-level is necessary in order to link individual responses and capture the serial correlation across a given respondent’s survey responses. In small area estimation applications, unit-level models have the additional benefit of allowing for arbitrary aggregation of data and may increase the precision of estimates (Hidiroglou and You,, 2016). Yet, despite the ubiquity of non-Gaussian survey data collected according to a longitudinal design, there is scarce unit-level methodology for such settings. Notably, there are few existing methods for the combination of ordinal and longitudinal data that are able to incorporate spatio-temporal dependence and account for complex survey designs.

In general, more attention has been paid to binary and nominal responses in the complex survey literature (Skinner,, 2018). In contrast, in the spatial literature for complex surveys, categorical data are typically modeled at the area level, which requires calculating weighted proportions first. When the data follow a natural ordering, it is common to ignore this structure and instead model the ordinal proportions as though they were nominal. For example, Bauder et al., (2021) model multinomial probabilities as multivariate normal after employing a logit transformation. Although it is possible to induce an ordinal structure by modifying the covariance matrix (Agresti,, 2010), this does not capture the dependence a unit-level approach would.

Meanwhile, most existing unit-level methodology for categorical survey data does not fully take into account the survey design or applies only to cross-sectional surveys. For example, Machini et al., (2022) present an ordinal model for spatial data, but disregard the survey design, Beltrán-Sánchez et al., (2024) model ordinal survey data with spatial dependence, but rely on the restrictive assumption that all design information is available in the form of covariates. Methods that do account for a longitudinal survey design generally assume Gaussian, binary, or occasionally nominal data. One exception is the longitudinal model of Kunihama et al., (2019) that is able to handle multivariate responses and incorporates the weights into a Bayesian nonparametric (BNP) sampling algorithm. However, the response types considered are continuous, binary, or nominal, but not ordinal. Sutradhar and Kovacevic, (2000) consider ordinal, longitudinal, complex design, but use GEE which does not easily extend to the spatial setting or allow for straightforward uncertainty quantification.

In this paper, we propose an ordinal model within a unit-level modeling framework that incorporates survey weights, spatial, and temporal dependence. We do so by extending the work of Parker et al., (2022) and Vedensky et al., (2023). The former models cross-sectional binary and nominal complex survey data at the unit level while the latter builds on this approach to cover longitudinal data with Gaussian or binary responses. By employing a particular stick-breaking representation of a multinomial likelihood (Linderman et al.,, 2015)—also known as a continuation ratio factorization (Fienberg,, 2007)—along with Pólya-Gamma data augmentation (Polson et al.,, 2013), we obtain conjugate sampling algorithms, with Bayesian, ordinal logistic regression as a special case. When data are especially high-dimensional, MCMC algorithms may prove insufficient, and we introduce Variational Bayes (VB) algorithms (Blei et al.,, 2017) as an alternative to Gibbs sampling for such applications.

Recent work has made use of this combination of the continuation ratio factorization and Pólya-Gamma data augmentation. Kang and Kottas, 2024a develop a BNP model for cross-sectional data and Kang and Kottas, 2024b extend this to functional, longitudinal data, modeling continuous curves. In the psychological literature, Jimenez et al., (2023) cover the special case of an item-response model with ordinal values. However, modifying these approaches to handle complex survey data with a spatial component is nontrivial, and to our knowledge, the link to Bayesian logistic regression has not been pointed out elsewhere.

To illustrate the efficacy of our proposed methodology, we apply our models to data from the Household Pulse Survey (HPS), a high-dimensional complex survey that follows a longitudinal design. This includes an empirical simulation study and a full data analysis.

This article proceeds as follows. In Section 2 we describe our general approach to modeling ordinal data. Section 3 describes how the ordinal and nominal models may be extended for use with survey data. In Section 4 variational Bayes algorithms for each of the models are introduced. We then present empirical results in Section 5 and conclude with a discussion in Section 6.

Let $\mathcal{U}=\{1,\ldots,N\}$ index a population, from which we observe a sample $\mathcal{S}\subseteq\mathcal{U}$ of size $|\mathcal{S}|=n$. Associated with each element of the sample are a set of covariates $\bm{x}_{i}=(x_{i1},\ldots,x_{iq})$ as well as an ordinal-valued response $y_{i}\in\{1,\ldots,K\}$, with $K>2$ possible categories. We also define the one-hot encoded vector representation of the response, $\widetilde{\bm{y}}_{i}=(\mathbbm{1}(y_{i}=1),\mathbbm{1}(y_{i}=2),\ldots,\mathbbm{1}(y_{i}=K))$, where $\mathbbm{1}(\cdot)$ denotes the indicator function. The response vectors for the entire sample may be collated into the $nK$-dimensional vector $\widetilde{\bm{y}}=(\widetilde{\bm{y}}_{1},\widetilde{\bm{y}}_{2},\ldots,\widetilde{\bm{y}}_{n}).$

It is common to model such data with a cumulative model (McCullagh,, 1980), where latent, ordered cutpoints $0\leq\gamma_{1}\leq\cdots\leq\gamma_{K-1}$ are introduced and it is assumed

with $F$ being any strictly monotonic distribution function and ${\bm{\beta}}$ a set of regression coefficients. This implies the probability of observing a given category to be

In the Bayesian setting, Albert and Chib, (1993) propose a popular data augmentation algorithm for cumulative, ordinal probit regression. They take $F$ to be the standard normal cdf $\Phi$ and add Gaussian latent variables to produce an easily implementable Gibbs sampler. This is the most widely used Bayesian formulation for ordinal models (Agresti,, 2010), and has been adapted to spatial models, as well (Higgs and Hoeting,, 2010; Schliep and Hoeting,, 2015; Carter et al.,, 2024).

Tutz, (1990, 1991) specify an alternative latent variable formulation for a sequential ordinal model, which has received relatively less attention in the literature. However, a number of authors have noted benefits, such as more flexible modeling of response curves and the possibility for category-specific covariates (Tutz et al.,, 2005; Bürkner and Vuorre,, 2019; Boes and Winkelmann,, 2006; Kang and Kottas, 2024a, ). In contrast to the cumulative approach, the sequential approach assumes no ordering for $\bm{\gamma}=(\gamma_{1},\ldots,\gamma_{K-1})$, and that the conditional probability for observing a given category $k<K$ is

and

for the last category. In this way, the probability of observing a particular category may be viewed as a sequential process, where category $k$ can only be attained if all previous categories, $j<k$ are not observed. The final category, $K$, is observed only if all other categories fail to be attained.

This formulation implies an independent binomial model for each response category, conditional on $\bm{\beta}$ and $\bm{\gamma}$, via the continuation-ratio parametrization of Fienberg, (1980), which is also referred to as a “stick-breaking” representation (Linderman et al.,, 2015). Specifically, a multinomial pmf parameterized by counts $n$ and a vector of probabilities $\bm{\pi}_{i}=(\pi_{i1},\ldots,\pi_{iK})$, can be expressed as the product

where $\pi_{ik}=\widetilde{\pi}_{ik}\prod_{j<k}(1-\widetilde{\pi}_{ij})$ and $n_{ik}=n_{i}-\sum_{j<k}y_{j}$ for $k=1,\ldots,K-1$. That is, if we let $\pi_{ik}=P(y_{i}=k|\bm{\beta},\bm{\gamma})$ and $\widetilde{\pi}_{ik}=F(\gamma_{k}-\bm{x}_{i}\bm{\beta})$, with $F$ the inverse logit function, we have that (2.1) and (2.2) define a multinomial likelihood over the ordinal categories.

Placing Gaussian priors on $\bm{\beta}$ and $\bm{\gamma}$, we may then specify the full ordinal logistic regression model hierarchically as a Bayesian GLM

For nominal responses, where there is no ordering of response categories, model (2) simplifies to estimating $\text{logit}(\widetilde{\pi}_{i})=\bm{x}_{i}\bm{\beta}_{k}$ with category-specific coefficients, $\bm{\beta}_{k}$ and with no need for cutpoints $\gamma_{k}$ (Parker et al.,, 2022)

Because both (2) and (2) have logistic likelihoods, it is possible to introduce latent variables that follow a Pólya-Gamma (PG) distribution (Polson et al.,, 2013) to obtain computationally efficient Gibbs samplers. In the case of (2) this yields a conjugate sampler for logistic ordinal regression.

A random variable $X$ is said to be distributed PG$(b,c)$ with shape parameter $b>0$ and scale parameter $c\in\mathbb{R}$ if it is equal in distribution to

where $g_{k}\overset{ind.}{\sim}\text{Gamma}(b,1)$. Polson et al., (2013) also show that

and derive the integral identity,

where $\kappa=a-b/2$ and $\omega\sim PG(b,0)$. In the right-hand side of (2.7), after conditioning on $\omega$, we have a Gaussian density in terms of $\lambda$. If we condition on $\lambda$, then we have a PG density in $\omega$. Therefore, sampling the model only requires alternating between Gaussian and PG random draws.

Notation for the full sampling algorithm is greatly simplified if, following Albert and Chib, (2001)’s algorithm for probit regression, we augment the covariate vector for each response such that if $y_{i}=k$, then $\bm{\widetilde{x}}_{ij}=(0,\ldots,1,\ldots,0,-\bm{x}_{i}^{\prime})$, for $j=1,\ldots,\min(k,K-1)$, where the $1$ is in the $j$th position. Further, let $\bm{\widetilde{X}}_{i}=(\widetilde{\bm{x}}_{i1},\ldots,\widetilde{\bm{x}}_{ij})^{\prime}$ be the new covariate matrix for respondent $i$ and

the augmented design matrix for the entire sample. Letting $\bm{\delta}=(\bm{\gamma}^{\prime},\bm{\beta}^{\prime})^{\prime},$ we then have

so that the integral identity (2.7) allows us to write the multinomial data model in (2) as

with $\kappa_{ik}=y_{ik}-1/2.$

Assuming the prior distribution $\bm{\delta}\sim N(\bm{d},\bm{D})$, a Gibbs sampler for the model requires iterating over only two steps

1. 1.

   $\omega_{ik}|\bm{\delta}\sim\text{PG}(1,\bm{\widetilde{x}_{ik}}^{\prime}\bm{\delta})$ (for $i=1,\ldots,n$ and $k=1,\ldots,k_{i}$)

2. 2.

   $\bm{\delta}|\bm{\Omega}\sim\text{N}(\bm{\mu}_{\delta},\bm{\Sigma}_{\delta})$,

where $\bm{\Sigma}_{\delta}=(\bm{D}^{-1}+\bm{\widetilde{X}^{\prime}\Omega\widetilde{X}})^{-1}$ and $\bm{\mu}_{\delta}=\bm{\Sigma}_{\delta}(\bm{\widetilde{X}}^{\prime}\widetilde{\bm{\kappa}}+\bm{D}^{-1}\bm{d})$ with $\bm{\omega}=(\omega_{11},\omega_{12},\ldots,\omega_{nk_{n}})$ and $\bm{\Omega}=\text{diag}(\bm{\omega}).$ The vector $\widetilde{\bm{\kappa}}$ is constructed by the same process as (2.8), where $\widetilde{\bm{\kappa}}_{i}=(\kappa_{i1},\ldots,\kappa_{ij})^{\prime}$ and $\widetilde{\bm{\kappa}}=(\widetilde{\kappa}_{1}^{\prime},\ldots,\widetilde{\kappa}_{n}^{\prime})^{\prime},$

In order to apply the above models to unit-level survey data, it is important to consider a number of additional factors. First, complex surveys often collect data according to a design that may not be known to the data modeler (Gelman,, 2007). Even if the design is known, factors such as non-response, attrition (in the case of longitudinal surveys) (Thompson,, 2015), and informative sampling (Parker et al.,, 2023) may cause bias. These complications must be accounted for through proper use of survey weights. Large-scale surveys are also likely to exhibit spatial dependence and it is common to produce predictions for spatial (or demographic) domains with low sample sizes (i.e., small area estimation). Additionally, many surveys are conducted repeatedly over time, or employ a longitudinal structure, where the same people or households are re-interviewed, leading to temporal correlation. The Bayesian hierarchical modeling framework allows us to address each of these types of dependence by extending the data and process model components of (2).

Despite the conjugacy in the above models, the Gibbs samplers may be inadequate for extremely high-dimensional problems, which require sampling a large number of latent variables and random effects (e.g., large sample sizes and fine-scaled geographies). In this case, we also propose a set of variational Bayes algorithms (Bishop,, 2006; Blei et al.,, 2017) as an alternative to MCMC. These algorithms are much faster albeit at the cost of approximating the posterior distribution rather than sampling from the exact posterior.

Instead of sampling sequentially, VB methods minimize the Kullback-Leibler divergence (Kullback and Leibler,, 1951)

between the true posterior $p(\theta|x)$ and an approximation $q(\theta)$ selected from a given class of distributions $\mathcal{Q}$. That is, an optimal approximation is selected as

A common choice for $\mathcal{Q}$ is the mean-field variational family (Blei et al.,, 2017) which consists of all densities $q$ that can be factored into mutually independent “variational factors” $q_{j}(\theta_{j})$ for each of $J$ latent variables so that

Since the KL divergence is not necessarily computable, the optimization problem may be reformulated as a maximization with respect to the evidence lower bound (ELBO) (Blei et al.,, 2017)

Durante and Rigon, (2019), building on earlier work of Jaakkola and Jordan, (2000), show that logistic models augmented by Pólya-Gamma latent variables have a particularly simple VB algorithm for performing this maximization, and we adapt this to the ordinal case.

Returning to the notation of the two-stage Gibbs sampler at the end of Section 2, note that the Gaussian and PG densities therein have the exponential family representation

and

respectively, with natural parameters $\bm{\eta}_{1}(\bm{y})=\bm{\widetilde{X}}^{\prime}\widetilde{\bm{\kappa}}+\bm{D}^{-1}\bm{d}$ and $\bm{\eta}_{2}(\bm{\omega})=-\frac{1}{2}\bm{D}^{-1}+\widetilde{\bm{X}}^{\prime}\bm{\Omega}\widetilde{\bm{X}}$ and $\eta_{ik}(\bm{\omega})=-\frac{1}{2}(\widetilde{\bm{x}}_{ik}^{\prime}\bm{\beta})^{2}.$ Durante and Rigon, (2019) show that the optimal CAVI parameter value updates at iteration $j$ are

1. 1.

   $\bm{\lambda}_{1}^{(j)}=\mathbb{E}_{q^{(j-1)(\bm{\omega})}}[\bm{\eta}_{1}(\bm{y})]$

2. 2.

   $\bm{\lambda}_{2}^{(j)}=\mathbb{E}_{q^{(j-1)}(\bm{\omega})}[\bm{\eta}_{2}(\bm{\omega}]$

3. 3.

   $\widecheck{\xi}_{ik}^{(j)}=\mathbb{E}_{q^{(j)}(\bm{\delta}_{ik})}[\eta_{ik}(\bm{\omega})]$ for $i=1,\ldots,n$ and $k=1,\ldots,k_{i}$,

where the $q^{(r)}(\bm{\delta}_{ik})$ densities are Gaussian and the $q^{(r)}(\omega_{ik})$ densities are PG.

Calculation of the expectation in Step 2 is straightforward due to the closed-form expression provided by (2.6). Convergence is measured by comparing the change in successive iterations of the ELBO and once this value drops below a predetermined threshold, the algorithm terminates. Independent samples are then drawn from the variational posterior N$(\widecheck{\bm{\mu}},\widecheck{\bm{\Sigma}})$, as an approximation to the true posterior, with $\widecheck{\bm{\mu}}=(-2\bm{\lambda}_{2})^{-1}\bm{\lambda}_{1}$ and $\widecheck{\bm{\Sigma}}=(-2\bm{\lambda}_{2})^{-1}$.

Parker et al., (2022) adapt this approach to cross-sectional binary and nominal pseudo-likelihood models with random effects. We further develop a cross-sectional, ordinal pseudo-likelihood sampling procedure in Algorithm 1 and a longitudinal variant in Algorithm 2. The matrices $\widetilde{G}$, $\widetilde{X}$ and $\widetilde{\Psi}$ are constructed in the same manner as (2.8), where if $y_{it}=k$, then we define $j_{it}=\min(k,K-1)$ and $\widetilde{G}_{it}=(\bm{u}_{it1},\ldots,\bm{u}_{itj_{it}})$ and let

and similarly for $\widetilde{X}=(\bm{x}_{it1},\ldots,\bm{x}_{itj_{it}})^{\prime}$ and $\widetilde{\Psi}=(\bm{\psi}_{it1},\ldots,\bm{\psi}_{itj_{it}})^{\prime}$.

In Algorithm 2, calculations for the variational parameters $\widetilde{\mu}_{\phi}$ depend on moments of the truncated normal distribution (Johnson et al.,, 1994). In that algorithm, $\Phi$ denotes the standard normal cdf and $\varphi$ the pdf. Note, too, the distinction between $(\widetilde{\mu}_{\phi})^{2}$ and $\widetilde{\mu}_{\phi^{2}}$. The former is the squared mean, while the latter is the variational variance parameter for $\phi$. Subscripts on the matrices, such as $X_{t}$, refer to the subset of rows corresponding to respondents at time $t$. A subscript such as, $\widecheck{\bm{\Sigma}}_{\eta_{-T}}$ refers to all entries of the matrix except those corresponding to the last time point $t=T$.