Bayesian Two-level Model for Repeated Partially Ordered Responses: Application to Adolescent Smoking Behavior Analysis

Abstract

Given the questionnaire design and the nature of the problem, partially ordered data that are neither completely ordered nor completely unordered are frequently encountered in social, behavioral, and medical studies. However, early developments in partially ordered data analysis are very limited and restricted only to cross-sectional data. In this study, we propose a Bayesian two-level regression model for analyzing repeated partially ordered responses in longitudinal data. The first-level model is defined for partially ordered observations of interest that are taken at each time point nested within individuals, while the second-level model is defined for individuals to assess the effects of their characteristics on the first-level model. A full Bayesian approach with the Markov chain Monte Carlo algorithm is developed for statistical inference. Simulation studies demonstrate the satisfactory performance of the developed methodology. The methodology is then applied to a longitudinal study on adolescent smoking behavior.

Keywords

two-level model partially ordered data adolescent smoking behavior longitudinal data MCMC methods

Introduction

Partially ordered data are categorical data that are frequently encountered in behavioral, social, psychological, and medical sciences. Unlike the well-studied nominal and ordinal data, partially ordered data that are neither completely ordered nor completely unordered can have highly complex data structures. A typical example of partially ordered data is the disability states of older adults in the Health, Aging, and Body Composition dataset described by Ip et al. (2013). As shown in Figure 1, different categories represent different disability states. Category 1 represents the nondisabled state; categories 2–5 indicate having one form of difficulty in performing daily activities, such as vertical mobility difficulty, horizontal mobility difficulty, and four other forms of difficulties; categories 6–8 indicate having two forms of difficulties; category 9 indicates having all six forms of difficulties captured by the study; and category D represents death. An arrow from one category to another indicates the dominance relationship. For example, the arrow between categories 1 and 4 shows that category 1 (no disability) dominates category 4 (horizontal mobility difficulty), whereas categories 4 and 5 (vertical mobility difficulty) are incomparable and are not linked by an arrow. In this typical example of partially ordered data, some categories follow rank order while others do not.

Figure 1.

The partially ordered structure of the 10 disability states of the older adults in the Health ABC data in Ip et al. (2013).

To cope with this special data structure, one may combine the incomparable categories at nearly the same level into one equivalent class and then treat the combined classes as completely ordered categories (Wilson 1992) to use the well-known ordinal models. In the aforementioned disability states example, applying this method leads to five ordinal categories, {1}, {2,3,4,5},{6,7,8},{9}, and {D}, which can be interpreted as categories of no disability, mild disability, severe disability, worst disability, and death, respectively. However, one consequence of this combination procedure is that the interpretation ability of the model is jeopardized because the combined categories are not differentiable after combination. For example, one cannot find the unique risk factors corresponding to category 4 after categories 2–5 have been combined. Another simple approach for analyzing partially ordered data is to regard all categories as nominal. However, such method ignores the valuable ordinal information and leads to information loss. Recently, researchers have proposed improved methods to model the dominance relationship for every category in the partially ordered structure without combining the incomparable categories. Meulders, Ip, and Boeck (2005) used a dominance matrix to describe the dominance relationship among the partially ordered categories and proposed latent variable models to analyze partially ordered responses to anger-related feelings. Ip et al. (2013) developed a partially ordered mixed hidden Markov model to analyze the hidden disability states of older adults with the structure of states shown in Figure 1. However, this approach does not model the relationships among incomparable categories, such as categories 2–5 in Figure 1, and is unable to manage highly complicated partially ordered data structures, such as the disjoint ones.

The ultimate goal of partially ordered data analysis is to model the relationships among the incomparable categories while retaining all ordinal information of the data structure. To this end, several studies have proposed alternative ways of modeling partially ordered data precisely. Zhang and Ip (2012) regarded the partially ordered structure as a hierarchy of nominal and ordinal substructures and developed partitioned conditional models through several classes of generalized linear models. Peyhardi, Trottier, and Guédon (2016) defined a class of partitioned conditional generalized linear models through a partition tree of categories to handle complex data structures with any level of hierarchy and applied this method on partially ordered responses. Ip, Chen, and Quandt (2016) developed an item response model for analyzing multiple disjoint partially ordered responses by using a two-stage estimation algorithm. However, the aforementioned developments have focused only on cross-sectional data, and none of them consider repeated partially ordered responses in the longitudinal setting. Moreover, a Bayesian method for analyzing partially ordered data is yet to be proposed.

Substantive research often encounters longitudinal data where the observed variables are measured repeatedly over time. Given that the classical assumption of independence between observations does not hold, those statistical methods that ignore the dependency structure can underestimate the variance of the estimated parameters and are therefore inapplicable. Statistical models, such as the random effects model, that consider the interdependence among repeated measurements of longitudinal data are popular in the statistical literature and have attracted wide attention in applications (see, e.g., Diggle 2002; Laird and Ware 1982; Verbeke and Molenberghs 2009). Longitudinal data can also be viewed as two-level data with repeated measurements (first level) nested within individuals (second level; Hedeker and Gibbons 2006). The two-level type model not only introduces random coefficients in the first-level regression to address the interdependence among the repeated measurements but also includes time-invariant covariates as explanatory variables in the second-level model to examine the influence of individual characteristics on the first-level model (Hox, Moerbeek, and van de Schoot 2010; Longford 1995). One of the primary advantages of two-level models is that they provide insights into the effects of time-invariant factors, such as gender and race, on the random coefficients in order to better reflect the relationships between the longitudinal covariates and the outcomes for each individual (Muñoz and Chang 2007; Raudenbush and Bryk 2002; Willms and Raudenbush 1989).

This research is motivated by a data set on adolescent smoking behavior taken from the National Longitudinal Surveys of Youth (NLSY). Adolescent smoking is one of the driving forces of many public health problems and has been linked to negative physical, psychological, and interpersonal consequences (Petraitis, Flay, and Miller 1995; Repetti, Taylor, and Seeman 2002a; Weschler et al. 1994). Adolescent smokers have different smoking behaviors, and each smoking behavior has unique risk factors (Schane, Ling, and Glantz 2010). Obtaining a comprehensive understanding of behavior-specific risk factors can inform treatment and prevention strategies for teenagers at various conditions and stages of smoking. Adolescent smokers can be classified into several ordered or not strictly ordered categories, namely, nonsmoker, former smoker, light-nonfrequent smoker, light-frequent smoker, heavy-nonfrequent smoker, and heavy-frequent smoker, and a comprehensive approach for dealing with such partially ordered data structure needs to be developed (Zhang and Ip 2012). Moreover, the NLSY is a longitudinal study that contains repeated partially ordered responses. This longitudinal feature introduces additional difficulties in modeling the partially ordered data structure. To our knowledge, this is the first study that precisely models partially ordered data in the longitudinal setting.

In this study, a Bayesian two-level model for analyzing repeated partially ordered adolescent smoking behavior data is proposed. In the first-level model, both the response and covariates are time varying, and the regression coefficients are regarded as random so that both the intercepts and slopes differ across individuals. We propose partitioned conditional models to formulate the accurate structure of the smoking behaviors without combining incomparable categories. A conditional coding algorithm is simultaneously developed to facilitate the formulation of the complex data structure. The model enables the identification of different risk factors for smoking behaviors from incomparable categories. In the second-level model, we model the variation of random coefficients through time-invariant covariates and random components at the individual level. This model also explains how the individual-level covariates affect the relationship between the time-varying covariates and the response in the first-level model. In the proposed model setting, the observed data likelihood involves intractable high-dimensional integrals, thereby making the existing maximum likelihood methods indirectly applicable for partially ordered data. This study develops a Bayesian approach, along with efficient Markov chain Monte Carlo (MCMC) techniques, for statistical inference. The Bayesian approach is proposed not only for its basic advantages, such as utilizing useful prior information and producing reliable results in small samples (Ansari and Jedidi 2000; Lee and Song 2004), but also for its power and efficiency in managing complicated models and data structures for which other methods are either unfeasible or work poorly (Congdon 2007; Skrondal and Rabe-Hesketh 2004). In this study, due to the complex dependence of the likelihood of the proposed model, many numerical optimization procedures are likely to fail in parameter estimation. By contrast, the sampling-based Bayesian approach with MCMC methods has been proven to be a powerful tool to overcome the various challenges confronting the conventional optimization procedures (e.g., Congdon 2007; Dunson 2000; Feng, Wu, and Song 2017; Song et al. 2013). Another appealing feature of MCMC is that various characteristics of the posterior distribution, such as the posterior mean, mode, and percentiles, can be investigated based on stationary simulated observations.

The rest of this article is organized as follows. Section 2 describes the partially ordered set theory and the model formulation. Section 3 briefly introduces the Bayesian analysis procedure. Section 4 evaluates the empirical performance of the proposed model through a simulation study. Section 5 presents a real data example concerning the risk factors of adolescent smoking. Section 6 concludes the article with a discussion. Technical details are provided in the Online Appendix (which can be found at http://smr.sagepub.com/supplemental/).

Model Description

Partially Ordered Set (Poset) Theory

The mathematical model for partially ordered data is a poset, which has been intensively studied as a mathematical object (Dushnik and Miller 1941; Birkhoff 1940). A poset is a set (P) equipped with a binary relation (≤). Any two distinct elements in a poset can be comparable or incomparable with respect to the binary relation. As a typical poset example, assume four elements e₁, e₂, e₃, and e₄, where $e_{1} \leq e_{2}$ , $e_{1} \leq e_{3}$ , $e_{2} \leq e_{4}$ , $e_{3} \leq e_{4}$ , and the other pairs are incomparable. The poset structure can be visualized by the Hasse diagram in Figure 2, where the nodes represent the elements and the arrows represent the dominance relationship under the binary relation. A chain in a poset is a totally ordered subset $C \subset P$ , while an antichain is a subset A of pairwise incomparable elements. For two subsets $S_{1}, S_{2} \subset P$ , either S₁ weakly dominates S₂ or S₁ and S₂ are weakly ordered if at least one element in S₂ is dominated by some element in S₁ and if no element in S₂ dominates any element in S₁. A set of subsets is called totally weekly ordered if pairwise subsets are weekly ordered. For instance, the poset in Figure 2, {{1},{2,3},{4}}, is totally weekly ordered.

Figure 2.

A simple partially ordered set structure with four categories.

The poset structure can be regarded as a hierarchy of nominal and ordinal substructures. The following ordered partition theorem, which is proven in Zhang and Ip (2012), provides a basis for specifying the hierarchical substructures of a poset starting from a set of given pairwise comparable and incomparable relationships among elements.

Theorem 1 (ordered partition): A finite poset can always be partitioned into totally weakly ordered antichains.

We take the poset mentioned in Zhang and Ip (2012) as an illustrative example. As shown in Figure 3, the poset ${a_{1}, \dots, a_{9}}$ has the following relations: $a_{2} > a_{4}, a_{3} > a_{4}, a_{3} > a_{5}, a_{3} > a_{6}, a_{5} > a_{7}, a_{6} > a_{7}$ , and $a_{8} > a_{9}$ . As a general procedure, we first partition the poset into disjoint connected components, and then for every connected component, we apply the ordered partition theorem to partition the poset into totally weakly ordered antichains. Afterward, we model the incomparable categories in the antichains. The detailed procedures are described as follows.

Figure 3.

The Hasse diagram of the partially ordered set of nine categories in the illustrative example. Note that there are three connected components.

Let Y be a random variable that takes the values in the poset categories from 1 to 9. The conditional modeling technique uses three random variables, namely, $Y^{(1)}$ , $Y^{(2)}$ , and $Y^{(3)}$ , to describe Y, where $Y^{(1)}$ , $Y^{(2)}$ , and $Y^{(3)}$ represent the disjoint connected components, totally weekly ordered antichains, and incomparable element in the antichains, respectively. The procedure for creating $Y^{(1)}$ , $Y^{(2)}$ , and $Y^{(3)}$ demonstrates how to partition the poset into a hierarchy of nominal and ordinal substructures.

First, we partition the poset into the disjoint connected components {1}, {2,3,4,5,6,7}, and {8,9} by using the breadth-first search algorithm in graph theory (Cormen 2009). The random variable $Y^{(1)}$ takes the values of 0, 1, and 2 corresponding to the three connected components, respectively. Given that the disjoint connected components are incomparable, $Y^{(1)}$ is nominally distributed with categories 0, 1, and 2 as follows:

Y^{(1)} = {\begin{array}{l} 0 & for Y = {1}, \\ 1 & for Y = {2,3,4,5,6,7}, \\ 2 & for Y = {8, 9} . \end{array}

Second, each connected component is partitioned into totally weekly ordered antichains. The first connected component {1} has only one category. We do not consider a model for it because the random variable is degenerate. The second connected component {2,3,4,5,6,7} corresponds to $Y^{(1)} = 1$ . Based on ordered partition theory, this component is partitioned into three totally weakly ordered antichains {2,3},{4,5,6}, and {7}, where the former weakly dominates the latter successively. Zhang and Ip (2012) describe the procedure of constructing antichains in detail. Given $Y^{(1)} = 1$ , let $Y^{(2)}$ take the values of 1, 2, and 3, which correspond to the three antichains, respectively. We have the following conditional distribution with ordinal categories:

[Y^{(2)} | Y^{(1)} = 1] = {\begin{array}{l} 1 & for Y = {2, 3}, \\ 2 & for Y = {4, 5, 6}, \\ 3 & for Y = {7} . \end{array}

Similarly, for the third connected component {8,9}, the ordinal random variable $[Y^{(2)} | Y^{(1)} = 2]$ can be defined as follows:

[Y^{(2)} | Y^{(1)} = 2] = {\begin{matrix} 1 & for Y = {8}, \\ 2 & for Y = {9} . \end{matrix}

Third, the elements in each antichain are pairwise incomparable and are thus modeled through nominal models. Taking the antichain {4,5,6} corresponding to ${Y^{(1)} = 1, Y^{(2)} = 2}$ as an example, we let $[Y^{(3)} | Y^{(1)} = 1, Y^{(2)} = 2]$ take the values of 0, 1, and 2 corresponding to the three elements, respectively. Then, $[Y^{(3)} | Y^{(1)} = 1, Y^{(2)} = 2]$ is nominally distributed as follows:

[Y^{(3)} | Y^{(1)} = 1, Y^{(2)} = 2] = {\begin{matrix} 0 & for Y = {4}, \\ 1 & for Y = {5}, \\ 2 & for Y = {6} . \end{matrix}

The conditional models for other antichains are similarly defined.

Through these steps, we can clearly demonstrate the partially ordered structure of the categories in a Hasse diagram as shown in Figure 3. Based on the definition of $Y^{(1)}$ , $Y^{(2)}$ , and $Y^{(3)}$ , the values of the triple correspond to each poset category (see Table 1). For example, { $Y^{(1)} = 1, Y^{(2)} = 1, Y^{(3)} = 0$ } corresponds to the poset category 2 with $Pr (Y = 2) = Pr (Y^{(1)} = 1, Y^{(2)} = 1, Y^{(3)} = 0)$ . Therefore, the partially ordered random variable Y can be fully described through the set of random variables { $Y^{(1)}, Y^{(2)}, Y^{(3)}$ } and can be conditionally modeled through the well-known ordinal and nominal models.

Table 1.

The One-to-one Mapping Between Y and $Y^{(1)}$ $Y^{(2)}$ $Y^{(3)}$ for the Partially Ordered Set (Poset) in the Illustrated Example.

Poset Category of Y	$Y^{(1)}$	$Y^{(2)}$	$Y^{(3)}$
1	0	NA	NA
2	1	1	0
3	1	1	1
4	1	2	0
5	1	2	1
6	1	2	2
7	1	3	NA
8	2	1	NA
9	2	2	NA

Note. NA = not applicable.

Two-level Model for Repeated Partially Ordered Responses

Consider a set of partially ordered observations ${y_{i t} : i = 1, \dots, N, t = 1, \dots, T}$ for individuals $i = 1, \dots, N$ that are measured at multiple time points $t = 1, \dots, T$ with the structure shown in Figure 3. Based on the poset partition technique described in Partially Ordered Set (Poset) Theory subsection, the modeling of the poset random variable Y corresponding to response $y_{i t}$ is converted to the conditional modeling of the triple { $Y^{(1)}, Y^{(2)}, Y^{(3)}$ }. Afterward, the complex partially ordered structure can be modeled through the conditional ordinal and nominal hierarchies.

The first-level models for the repeated partially ordered responses are defined through the following conditional models:

1. For $Y_{i t}^{(1)} = 0, 1, 2$ , which corresponds to the three disjoint components {1}, {2,3,4,5,6,7} and {8,9}, a random coefficient nominal model is considered. Using the first category $Y_{i t}^{(1)} = 0$ as a reference, the model for $Y_{i t}^{(1)}$ is defined as follows:

log (\frac{π_{i t k}^{(1)}}{π_{i t 0}^{(1)}}) = α_{i 1 k} + {β^{'}}_{i 1 k} x_{i t}, k = 1, 2,

where $x_{i t}$ is an $s \times 1$ vector of time-varying first-level covariates, $α_{i 1 k}$ and $β_{i 1 k}$ are the random intercepts and random slopes corresponding to individual i, $π_{i t k}^{(1)} = P r (Y_{i t}^{(1)} = k)$ , and $π_{i t 0}^{(1)} + π_{i t 1}^{(1)} + π_{i t 2}^{(1)} = 1$ .

2. For $[Y_{i t}^{(2)} = 1, 2, 3 | Y_{i t}^{(1)} = 1]$ that corresponds to the weekly ordered antichains ${2, 3}$ , ${4, 5, 6}$ , and ${7}$ , the conditional proportional odds model with random coefficients is defined for $Y_{i t}^{(2)}$ as follows:

log (\frac{\sum_{j \leq k} π_{i t j}^{(2.1)}}{\sum_{j > k} π_{i t j}^{(2.1)}}) = α_{i 2 k} + β'_{i 21} x_{i t}, k = 1, 2,

where $π_{i t k}^{(2.1)} = P r (Y_{i t}^{(2)} = k | Y_{i t}^{(1)} = 1)$ . The ordinal feature imposes the constraint that $α_{i 22}$ ≥ $α_{i 21}$ and the proportional odds assumption that $β_{i 21}$ is invariant to k.

For $[Y_{i t}^{(2)} = 1, 2 | Y_{i t}^{(1)} = 2]$ that corresponds to the weekly ordered antichains ${8}$ and ${9}$ , the conditional ordinal model with random coefficients is defined as:

log (\frac{π_{i t 1}^{(2.2)}}{π_{i t 2}^{(2.2)}}) = α_{i 31} + {β^{'}}_{i 31} x_{i t},

where $π_{i t k}^{(2.2)} = P r (Y_{i t}^{(2)} = k | Y_{i t}^{(1)} = 2)$ .

3. We regard the elements within antichains as nominal categories and propose conditional nominal models as follows:

For $[Y_{i t}^{(3)} = 0, 1 | Y_{i t}^{(1)} = 1, Y_{i t}^{(2)} = 1]$ that corresponds to the elements in the antichain ${2, 3}$ , the model is defined as follows:

log (\frac{π_{i t 1}^{(3.11)}}{π_{i t 0}^{(3.11)}}) = α_{i 41} + {β^{'}}_{i 41} x_{i t},

where $π_{i t k}^{(3.11)} = P r (Y_{i t}^{(3)} = k | Y_{i t}^{(1)} = 1, Y_{i t}^{(2)} = 1)$ .

For $[Y_{i t}^{(3)} = 0, 1, 2 | Y_{i t}^{(1)} = 1, Y_{i t}^{(2)} = 2]$ that corresponds to the antichain ${4, 5, 6}$ , the model is defined as follows:

log (\frac{π_{i t k}^{(3.12)}}{π_{i t 0}^{(3.12)}}) = α_{i 5 k} + {β^{'}}_{i 5 k} x_{i t}, k = 1, 2,

where $π_{i t k}^{(3.12)} = P r (Y_{i t}^{(3)} = k | Y_{i t}^{(1)} = 1, Y_{i t}^{(2)} = 2)$ .

The second-level model that examines the interdependence among the repeated measurements and the variation of random coefficients $α_{i r k}$ and $β_{i r k}$ is defined as follows:

\begin{matrix} α_{i r k} = μ_{α r k} + {γ^{'}}_{α} z_{i} + ε_{α i}, \\ β_{i r k} = μ_{β r k} + {γ^{'}}_{β} z_{i} + ε_{β i}, \end{matrix}

where $α_{i r k}$ and $β_{i r k}$ are the intercepts and regression coefficients involved in equations (1) –(5), $z_{i}$ is a $q \times 1$ vector of time-invariant individual-level covariates, $γ_{α}$ is a $q \times 1$ vector of coefficients, $γ_{β}$ is a $q \times s$ matrix of coefficients, $μ_{α r k}$ and $μ_{β r k}$ are the intercepts, and $ε_{α i}$ and $ε_{β i}$ are the scalar and $s \times 1$ vector of residual errors, respectively. $ε_{i} = (ε_{α i}, {ε^{'}}_{β i})^{'} \sim N (0, Ψ)$ , where $Ψ = diag (ψ_{α}, ψ_{β 1}, \dots, ψ_{β 5})$ is an $(s + 1) \times (s + 1)$ diagonal covariance matrix. The second-level model investigates the effects of the time-invariant covariates on the first-level regression parameters while addressing the dependence of the repeated measurements within individuals. Notably, we assume that $γ_{α}$ , $γ_{β}$ , and $Ψ$ in equation (6) are invariant to r and k and that $Ψ$ is a diagonal matrix. An extension to the case with $(r, k)$ -specific $γ_{α r k}$ , $γ_{β r k}$ , and nondiagonal $Ψ_{r k}$ is straightforward. However, given that the proposed model is already parameter-rich, such extension increases the number of parameters considerably and thus may cause difficulties in interpretation.

The two-level model defined above is particularly appropriate for research designs where data for participants are collected at a nested structure (e.g., repeated measurements nested within individuals). In the longitudinal study of adolescent smoking behavior discussed in the Introduction section, while potential predictors, such as age and educational enrollment status, can influence adolescent smoking behavior over time, these effects may exhibit heterogeneity due to their dependence on certain time-invariant individual characteristics, such as gender and race. A simple way to incorporate this into the regression model is to add additional independent binary variables to account for gender (male = 0) and race (nonwhite = 0). This would have the effect of shifting the log odds in equations (1 )–(5) up or down, but it still assumes that the effects of age and educational enrollment status on adolescents’ smoking behavior are the same regardless of their gender and race. This may not be the case in reality. Gender or racial gap is likely to cause all of the predictors to have different sorts of effects in males and females or whites and nonwhites. Moreover, compared with the fixed effect model that treats cluster-specific characteristics as nuisances and does not estimate them, the proposed two-level random effects model views the clustering of individuals’ repeated measurements as a feature of interest in its own right, and not just a nuisance to be adjusted for, thereby providing additional insights into the rationale behind the clustering.

Bayesian Inference

Let $μ$ be a vector that includes all the elements of $μ_{α r k}$ and $μ_{β r k}$ , $γ = (γ_{α}, γ_{β})^{'}$ , and $θ = (μ, γ, Ψ)$ . We assign prior distributions for the unknown parameters in $θ$ as follows: for $j = 1, \dots, s + 1$ ,

p (μ) \sim N (μ_{0}, Σ_{0}) p (ψ_{j}^{- 1}) \sim Gamma (a_{0 j}, b_{0 j}), p (γ_{j} | ψ_{j}) \sim N (γ_{0 j}, ψ_{j}^{- 1} H_{0 j}),

where $ψ_{j}$ and $γ_{j}$ are the jth diagonal element of $Ψ$ and the jth row of $γ$ , respectively, while $μ_{0}$ , $a_{0 j}$ , $b_{0 j}$ , $γ_{0 j}$ , and the positive definite matrices $Σ_{0}$ and $H_{0 j}$ are hyperparameters which values can be determined based on prior information.

Let $Y = {y_{i t} : i = 1, \dots, N, t = 1, \dots, T}$ be the set of partially ordered responses, $X = {x_{i t} : i = 1, \dots, N, t = 1, \dots, T}$ be the set of first-level covariates, $Z = {z_{i} : i = 1, \dots, N}$ be the set of second-level covariates, $β_{i}$ be a vector that includes all the random intercepts $α_{i r k}$ and the random coefficients $β_{i r k}$ in models (1)–(5) for the ith individual and $β = {β_{i} : i = 1, \dots, N}$ . The Bayesian inference is based on the posterior distribution of $θ$ given the observed data, that is, $p (θ | Y, X, Z)$ . However, given the presence of the random coefficients in $β$ , this posterior distribution involves intractable integrals. Following Tanner and Wong (1987), we augment the unobservable random coefficients in the posterior distribution and work on the joint distribution $p (θ, β | Y, X, Z)$ . Although this joint posterior distribution is relatively easy to manage, sampling directly from such distribution is complex and requires the use of MCMC methods, such as the Gibbs sampler (Geman and Geman 1984) and the Metropolis–Hastings algorithm (Hastings 1970; Metropolis et al. 1953), to sample $θ$ from $p (θ | Y, X, Z, β)$ and $β$ from $p (β | Y, X, Z, θ)$ iteratively. The Bayesian estimate of $θ$ is taken as the sample mean of the random observations generated. The Gibbs sampling procedure and the full conditional distributions involved are provided in Online Appendix A (which can be found at http://smr.sagepub.com/supplemental/).

Simulation Study

In this section, we conduct two simulations to evaluate the performance of the proposed method. Simulation 1 (Simulation 1 subsection) examines the performance of parameter estimation. Simulation 2 (Simulation 2 subsection) checks the robustness of the proposed model and MCMC method to prior inputs, the burn-in phase and length of MCMC, the initial values of the unknown parameters, and model misspecification.

Simulation 1

We consider the model defined by equations (1 )–(6) with $s = 5$ , $q = 2$ , and $T = 5$ . The response variable $y_{i t}$ follows the poset structure shown in Figure 3. For the first-level covariates $x_{i t j}, j = 1, \dots,5$ , $x_{i t 1}$ and $x_{i t 2}$ are independently generated from $N (0, 1)$ , while $x_{i t 3}, x_{i t 4}$ , and $x_{i t 5}$ are independently generated from $U (0, 1)$ , $t (3)$ , and the Bernoulli distribution with a success probability of .6, where $U (0, 1)$ denotes the uniform distribution in [0,1] and $t (3)$ denotes the student t-distribution with three degrees of freedom. The second-level covariates $z_{i 1}$ and $z_{i 2}$ are independently generated from $N (0, 1)$ and the Bernoulli distribution with a success probability of .6. In the above setting, models (1)–(5) include parameters ${(α_{i r k}, β_{i r k}), r = 1, k = 1, 2}$ , ${(α_{i r k}, β_{i r k}), r = 2, k = 1, 2}$ , ${(α_{i r k}, β_{i r k}), r = 3, k = 1}$ , ${(α_{i r k}, β_{i r k}), r = 4, k = 1}$ , and ${(α_{i r k}, β_{i r k}), r = 5, k = 1, 2}$ , respectively. Let $μ_{r k} = (μ_{α r k}, {μ^{'}}_{β r k})^{'}$ . Then, for each combination of $(r, k)$ , $μ_{r k}$ is a six-dimensional vector, in which the first element is the mean of the random intercept $α_{i r k}$ and the other five elements are the means of the random coefficients in $β_{i r k}$ . The population values of the model parameters are given as follows: $ψ_{α} = ψ_{β 1} = \dots = ψ_{β 5} = 0.36$ ,

(\begin{array}{l} μ_{11} & μ_{12} & μ_{21} & μ_{22} & μ_{31} & μ_{41} & μ_{51} & μ_{52} \\ 2.0 & 0.5 & - 1.0 & 1.5 & 0.5 & 0.0 & 0.0 & 0.0 \\ 0.8 & - 0.8 & 1.0 & 1.0 & 0.8 & 0.8 & 1.2 & - 1.2 \\ - 0.8 & 0.8 & 0.8 & 0.8 & - 1.0 & - 0.8 & 1.0 & 0.8 \\ 1.0 & 0.8 & 0.8 & 0.8 & 1.2 & 0.8 & 0.8 & - 0.8 \\ 1.2 & 1.0 & - 0.8 & - 0.8 & 0.8 & 0.8 & 0.8 & 0.8 \\ 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 & 0.8 \end{array})

γ' = [\begin{matrix} γ_{α 1} & γ_{β 11} & γ_{β 21} & γ_{β 31} & γ_{β 41} & γ_{β 51} \\ γ_{α 2} & γ_{β 12} & γ_{β 22} & γ_{β 32} & γ_{β 42} & γ_{β 52} \end{matrix}] = [\begin{matrix} 0.8 & 0.8 & - 1.0 & 1.0 & - 0.8 & 0.8 \\ - 0.8 & 1.2 & 0.8 & 0.8 & 0.8 & 0.8 \end{matrix}] ≜ [\begin{matrix} γ_{1} \\ γ_{2} \end{matrix}] .

To generate $y_{i t}$ , we generate $ε_{i}, i = 1, \dots, N$ from $N (0, Ψ)$ and compute $α_{i r k}$ and $β_{i r k}$ according to model (6). We then obtain the probability of each conditional nominal or ordinal categories of $Y^{(1)}$ , $[Y^{(2)} | Y^{(1)}]$ , and $[Y^{(3)} | Y^{(1)}, Y^{(2)}]$ according to models (1)–(5). The probability of each poset category can be calculated based on the conditional probabilities. For example, the probability of the poset category {2} can be calculated as follows:

\begin{array}{l} P r (y_{i t} = 2) = P r (Y_{i t}^{(1)} = 1, Y_{i t}^{(2)} = 1, Y_{i t}^{(3)} = 0), \\ = P r (Y_{i t}^{(1)} = 1) P r (Y_{i t}^{(2)} = 1 | Y_{i t}^{(1)} = 1) P r (Y_{i t}^{(3)} = 0 | Y_{i t}^{(1)} = 1, Y_{i t}^{(2)} = 1), \\ = π_{i t 1}^{(1)} π_{i t 1}^{(2.1)} π_{i t 0}^{(3.12)} . \end{array}

Afterward, we generate the partially ordered $y_{i t}$ shown in Figure 3 via the probabilities of the nine categories.

In the simulation study, the prior inputs for the hyperparameters in the prior distributions are assigned as follows: (prior 1) $a_{0 j} = 11, b_{0 j} = 3, j = 1, \dots, 6$ , the elements of $γ_{0 j}, j = 1, \dots,6$ are all set to 0.8, the elements of $μ_{0}$ are set twice their population values, and $Σ_{0}$ and $H_{0 j}$ are identity matrices with appropriate dimensions. Two sample sizes, N = 1,000 and N = 2,000, are considered in this simulation study. A total of 100 data sets are generated under each sample size. We first conduct a pilot study to determine the burn-in phase of MCMC. We randomly select a data set generated above and run three parallel sequences that start from different initial values. The MCMC algorithm converges if the three parallel sequences mix well. The burn-in phase is then determined as the number of iterations at convergence. In this simulation, the algorithm converges within 8,000 iterations. After a burn-in phase of 8,000 iterations, we collect additional 10,000 observations to obtain the Bayesian estimates of the unknown parameters and their standard error (SE) estimates. Let $θ$ denote an arbitrary element of $θ$ , $θ_{0}$ be the true value of $θ$ , and ${\hat{θ}}^{(j)}$ be the estimate of $θ$ at the j th replication. On the basis of 100 replications, the performance of the parameter estimation is evaluated through the bias (bias), root mean square error (RMS), and the ratio of SE/standard deviation (SD), where the bias of $\hat{θ} {= 100}^{- 1} \sum_{j = 1}^{100} {\hat{θ}}^{(j)} - θ_{0}$ , the RMS error of $\hat{θ} {= {100}^{- 1} \sum_{j = 1}^{100} {({\hat{θ}}^{(j)} - θ_{0})}^{2}}^{1 / 2}$ , and the SE/SD $\hat{θ} = M e a n o f {S E ({\hat{θ}}^{(j)}), j = 1, \dots,100} / S D o f {{\hat{θ}}^{(j)}, j = 1, \dots,100}$ , in which $S E ({\hat{θ}}^{(j)})$ stands for the standard error estimate of ${\hat{θ}}^{(j)}$ . Table 2 presents the estimation results obtained under $N = 1, 000$ and $N = 2, 000$ . All the Bias and RMS values are reasonably small and most of SD/SE values are close to 1, thereby indicating that the proposed method performs satisfactorily in terms of both parameter estimates and their SE estimates. As expected, the performance is improved when the sample size increases.

Table 2.

Estimation Results in Simulation 1.

Para	$N = 1, 000$			$N = 2, 000$			Para	$N = 1, 000$			$N = 2, 000$
Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD	Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD
$μ_{α 11}$	.028	.147	0.979	.045	.068	1.053	$μ_{α 12}$	$- .008$	.151	1.088	.026	.140	0.977
$μ_{β 111}$	.021	.065	1.052	.012	.058	0.942	$μ_{β 121}$	.021	.077	0.961	.016	.067	1.069
$μ_{β 112}$	.024	.090	1.120	−. $008$	.065	0.961	$μ_{β 122}$	.022	.098	1.069	$- .009$	.070	1.169
$μ_{β 113}$	.009	.088	1.004	.001	.065	0.910	$μ_{β 123}$	.041	.085	1.002	.004	.065	0.976
$μ_{β 114}$	.037	.230	1.102	.019	.176	1.030	$μ_{β 124}$	$- .025$	.246	1.050	.024	.194	0.970
$μ_{β 115}$	.021	.144	0.966	−. $032$	.066	0.932	$μ_{β 125}$	.002	.153	0.902	$- .030$	.072	0.986
$μ_{α 21}$	.012	.118	0.985	.006	.111	1.047	$μ_{α 22}$	.060	.124	1.073	.024	.108	1.008
$μ_{β 211}$	.019	.073	0.939	.006	.050	0.987	—	—	—	—	—	—	—
$μ_{β 212}$	.028	.077	0.972	.000	.056	1.046	—	—	—	—	—	—	—
$μ_{β 213}$	.039	.097	1.002	−. $001$	.071	0.994	—	—	—	—	—	—	—
$μ_{β 214}$	.056	.155	1.143	.025	.158	1.051	—	—	—	—	—	—	—
$μ_{β 215}$	.022	.104	1.032	$- .049$	.095	1.097	—	—	—	—	—	—	—
$μ_{α 31}$	$- .046$	.244	1.094	.026	.183	1.055	$μ_{α 41}$	.029	.180	0.977	.062	.129	1.076
$μ_{β 311}$	.005	.135	1.102	.005	.092	1.086	$μ_{β 411}$	.021	.081	1.127	.002	.062	0.935
$μ_{β 312}$	.030	.118	1.130	$- .001$	.093	1.032	$μ_{β 412}$	.034	.091	1.042	.007	.072	1.070
$μ_{β 313}$	$- .016$	.139	1.063	.005	.096	0.978	$μ_{β 413}$	$- .061$	.109	1.082	.003	.075	1.045
$μ_{β 314}$	.014	.332	1.052	$- .013$	.221	0.974	$μ_{β 414}$	.098	.259	1.054	.012	.195	1.043
$μ_{β 315}$	$- .023$	.178	1.106	$- .015$	.143	1.113	$μ_{β 415}$	.047	.157	1.050	$- .027$	.110	1.098
$μ_{α 51}$	.035	.159	1.103	.019	.138	1.096	$μ_{α 52}$	.055	.166	1.103	.084	.150	1.080
$μ_{β 511}$	.008	.092	1.062	$- .004$	.073	0.998	$μ_{β 521}$	.014	.099	1.062	.001	.069	1.065
$μ_{β 512}$	.032	.106	1.130	$- .005$	.081	1.071	$μ_{β 522}$	.026	.101	1.130	$- .004$	.071	1.030
$μ_{β 513}$	.043	.115	1.067	$- .004$	.086	1.067	$μ_{β 523}$	.047	.122	1.067	$- .009$	.085	1.027
$μ_{β 514}$	$- .015$	.255	0.942	$- .009$	.230	1.046	$μ_{β 524}$	.068	.258	0.942	.014	.224	0.922
$μ_{β 515}$	$- .020$	.187	1.068	$- .072$	.131	1.025	$μ_{β 525}$	.017	.154	1.068	$- .022$	.147	1.026
$γ_{α 1}$	.065	.079	0.955	.039	.066	0.916	$γ_{α 2}$	$- .057$	.112	0.897	$- .056$	.110	0.942
$γ_{β 11}$	.039	.056	0.995	.008	.039	0.963	$γ_{β 12}$	.033	.074	.0933	.003	.060	0.998
$γ_{β 21}$	.051	.049	1.064	.014	.035	0.967	$γ_{β 22}$	.018	.087	0.902	.012	.060	0.946
$γ_{β 31}$	.050	.057	0.903	.012	.032	1.067	$γ_{β 32}$	.008	.105	1.013	.012	.066	0.999
$γ_{β 41}$	.029	.130	0.914	.002	.103	0.982	$γ_{β 42}$	$- .035$	.169	0.888	.019	.139	0.962
$γ_{β 51}$	.037	.086	1.069	$- .015$	.059	0.920	$γ_{β 52}$	.044	.127	0.961	.030	.093	0.975
$ψ_{α}$	.060	.068	1.080	.018	.054	0.947
$ψ_{β 1}$	$- .011$	.053	1.226	.003	.044	0.909
$ψ_{β 2}$	$- .007$	.049	0.992	$- .001$	.045	1.062
$ψ_{β 3}$	.045	.057	1.110	.012	.040	0.931
$ψ_{β 4}$	.016	.045	1.515	.016	.056	1.167
$ψ_{β 5}$	.033	.071	1.488	.009	.063	0.964

Note. $N = 1, 000$ and $N = 2, 000$ . RMS = root mean square; SE/SD = standard error/standard deviation.

The computer code is written in R. The R code for this simulation is provided in a freely accessable website (http://www.sta.cuhk.edu.hk/xysong/BTLM).

Simulation 2

To check the sensitivity of Bayesian results to prior inputs, we disturb prior 1 by modifying the hyperparameters as follows: (prior 2) $a_{0 j} = 5$ , $b_{0 j} = 4$ , $j = 1, \dots,6$ , the elements of $γ_{0 j}$ , $j = 1, \dots,6$ , and $μ_{0}$ are all set to 0, $Σ_{0}$ and $H_{0 j}$ are $4 I$ , where $I$ represents the identity matrix with an appropriate dimension. The 100 data sets that are simulated with $N = 1, 000$ in simulation 1 are reanalyzed under prior 2. Table 3 reports the estimation results obtained under prior 2, along with those obtained under prior 1 for ease comparison. No substantial difference exists between the parameter estimates under priors 1 and 2.

Table 3.

Estimation Results Under Different Priors in Simulation 2.

Para	Prior 1			Prior 2			Para	Prior 1			Prior 2
Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD	Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD
$μ_{α 11}$	.028	.147	0.979	.058	.171	0.948	$μ_{α 12}$	$- .008$	.151	1.088	.021	.182	0.976
$μ_{β 111}$	.021	.065	1.052	.024	.091	0.954	$μ_{β 121}$	.021	.077	0.961	.041	.094	1.010
$μ_{β 112}$	.024	.090	1.120	.023	.094	1.013	$μ_{β 122}$	.022	.098	1.069	.038	.099	1.083
$μ_{β 113}$	.009	.088	1.004	.009	.102	0.982	$μ_{β 123}$	.041	.085	1.002	.027	.115	1.008
$μ_{β 114}$	.037	.230	1.102	.041	.228	0.903	$μ_{β 124}$	$- .025$	.246	1.050	.016	.191	0.909
$μ_{β 115}$	.021	.144	0.966	$- .038$	.164	0.958	$μ_{β 125}$	.002	.153	0.902	$- .024$	.168	0.966
$μ_{α 21}$	.012	.118	0.985	.034	.131	1.043	$μ_{α 22}$	.060	.124	1.073	.071	.150	0.939
$μ_{β 211}$	.019	.073	0.939	.022	.081	0.959	—	—	—	—	—	—	—
$μ_{β 212}$	.028	.077	0.972	.026	.097	0.904	—	—	—	—	—	—	—
$μ_{β 213}$	.039	.097	1.002	.016	.117	1.104	—	—	—	—	—	—	—
$μ_{β 214}$	.056	.155	1.143	$- .025$	.107	0.987	—	—	—	—	—	—	—
$μ_{β 215}$	.022	.104	1.032	.013	.126	0.952	—	—	—	—	—	—	—
$μ_{α 31}$	$- .046$	.244	1.094	.029	.279	1.092	$μ_{α 41}$	.029	.180	0.977	.058	.167	1.089
$μ_{β 311}$	.005	.135	1.102	.018	.122	1.056	$μ_{β 411}$	.021	.081	1.127	.043	.089	1.116
$μ_{β 312}$	.030	.118	1.130	$- .013$	.124	1.084	$μ_{β 412}$	.034	.091	1.042	.039	.098	1.051
$μ_{β 313}$	$- .016$	.139	1.063	.023	.161	1.130	$μ_{β 413}$	$- .061$	.109	1.082	.023	.121	1.088
$μ_{β 314}$	.014	.332	1.052	$- .012$	.208	0.998	$μ_{β 414}$	.098	.259	1.054	.042	.245	1.005
$μ_{β 315}$	$- .023$	.178	1.106	$- .022$	.159	1.129	$μ_{β 415}$	.047	.157	1.050	$- .015$	.168	1.107
$μ_{α 51}$	.035	.159	1.103	.053	.200	1.086	$μ_{α 52}$	.055	.166	1.103	.075	.185	1.071
$μ_{β 511}$	.008	.092	1.062	.030	.094	1.081	$μ_{β 521}$	.014	.099	1.062	.040	.095	1.066
$μ_{β 512}$	.032	.106	1.130	.017	.108	1.095	$μ_{β 522}$	.026	.101	1.130	.031	.104	1.132
$μ_{β 513}$	.043	.115	1.067	.037	.127	1.112	$μ_{β 523}$	.047	.122	1.067	.025	.123	1.088
$μ_{β 514}$	$- .015$	.255	0.942	.034	.219	1.046	$μ_{β 524}$	.068	.258	0.942	$- .028$	.278	0.994
$μ_{β 515}$	$- .020$	.187	1.068	.027	.172	0.987	$μ_{β 525}$	.017	.154	1.068	$- .018$	.190	1.077
$γ_{α 1}$	.065	.079	0.955	.096	.083	1.043	$γ_{α 2}$	$- .057$	.112	0.897	$- .076$	.132	0.960
$γ_{β 11}$	.039	.056	1.005	.061	.051	1.119	$γ_{β 12}$	.033	.074	0.933	.015	.096	0.976
$γ_{β 21}$	.051	.049	1.064	.063	.056	1.025	$γ_{β 22}$	.018	.087	0.902	.013	.102	0.916
$γ_{β 31}$	.050	.057	0.903	.055	.052	1.019	$γ_{β 32}$	.008	.105	1.013	.022	.104	1.078
$γ_{β 41}$	.029	.130	0.914	$- .049$	.137	0.941	$γ_{β 42}$	$- .035$	.169	0.888	$- .027$	.213	0.919
$γ_{β 51}$	.037	.086	1.069	.038	.109	0.913	$γ_{β 52}$	.044	.127	0.961	.061	.139	0.979
$ψ_{α}$	.060	.068	1.080	.087	.067	1.271
$ψ_{β 1}$	$- .011$	.053	1.226	.128	.070	1.126
$ψ_{β 2}$	$- .007$	.049	0.992	.131	.075	1.146
$ψ_{β 3}$	.045	.057	1.110	.091	.082	0.904
$ψ_{β 4}$	.016	.045	1.515	.105	.134	1.346
$ψ_{β 5}$	.033	.071	1.488	.097	.100	1.304

Note. $N = 1, 000$ . RMS = root mean square; SE/SD = standard error/standard deviation.

To examine the robustness of MCMC to its burn-in phase and length, we change the burn-in phase and length from {8,000, 10,000} to {15,000, 15,000} iterations. Specifically, for each of the above 100 data sets with $N = 1, 000$ , we collect 15,000 samples after a burn-in phase of 15,000 iterations to conduct posterior inference. The results are presented in Table 4 (left panel). The parameter estimates are very close to those reported in Table 2.

Table 4.

Estimation Results Under Different Burn-in Phase and Initial Values of the Unknown Parameters in Simulation 2.

Para	Burn-in = 15,000			Initial 2			Para	Burn-in = 15,000			Initial 2
Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD	Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD
$μ_{α 11}$	−.036	.138	0.937	.038	.131	1.014	$μ_{α 12}$	.043	.156	0.973	.034	.139	1.025
$μ_{β 111}$	.028	.068	0.986	.029	.071	0.991	$μ_{β 121}$	−.018	.077	1.005	−.021	.089	1.023
$μ_{β 112}$	.030	.067	0.998	.035	.067	1.048	$μ_{β 122}$	.026	.079	0.956	.038	.071	1.067
$μ_{β 113}$	.017	.073	0.891	−.018	.070	1.043	$μ_{β 123}$	.024	.070	1.005	−.011	.072	1.026
$μ_{β 114}$	−.006	.195	0.949	.023	.194	1.056	$μ_{β 124}$	−.022	.215	1.026	.023	.203	0.927
$μ_{β 115}$	.022	.126	0.949	.017	.119	0.954	$μ_{β 125}$	.005	.145	0.947	.021	.141	1.003
$μ_{α 21}$	.026	.110	1.012	.027	.117	1.081	$μ_{α 22}$	.041	.105	1.019	.059	.106	0.922
$μ_{β 211}$	.022	.058	0.994	.016	.065	0.998	—	—	—	—	—	—	—
$μ_{β 212}$	.018	.063	0.925	.018	.057	0.895	—	—	—	—	—	—	—
$μ_{β 213}$	.023	.080	0.886	.026	.075	0.944	—	—	—	—	—	—	—
$μ_{β 214}$	.057	.145	1.109	.040	.166	1.000	—	—	—	—	—	—	—
$μ_{β 215}$	−.005	.091	1.058	.036	.097	0.977	—	—	—	—	—	—	—
$μ_{α 31}$	.055	.233	1.104	.031	.214	0.885	$μ_{α 41}$	.036	.187	1.081	.042	.172	0.912
$μ_{β 311}$	.026	.105	1.112	−.020	.118	1.071	$μ_{β 411}$	.001	.086	1.078	.018	.073	1.046
$μ_{β 312}$	.035	.105	1.085	.004	.113	1.081	$μ_{β 412}$	.020	.092	1.118	−.036	.088	1.103
$μ_{β 313}$	.027	.114	1.072	.021	.126	0.961	$μ_{β 413}$	.010	.098	1.023	.043	.096	1.074
$μ_{β 314}$	−.015	.278	1.093	.037	.317	1.011	$μ_{β 414}$	.028	.294	0.998	.068	.256	1.077
$μ_{β 315}$	.019	.195	1.125	.011	.116	1.074	$μ_{β 415}$	.036	.165	1.116	.015	.152	0.979
$μ_{α 51}$	.059	.134	1.134	.051	.115	1.097	$μ_{α 52}$	−.015	.143	1.126	.054	.157	0.991
$μ_{β 511}$	.022	.088	1.046	.026	.089	1.001	$μ_{β 521}$	.015	.101	1.007	−.019	.088	1.021
$μ_{β 512}$	.028	.104	0.989	.019	.091	1.063	$μ_{β 522}$	.022	.097	1.060	.028	.097	1.129
$μ_{β 513}$	.020	.093	1.067	.029	.097	1.109	$μ_{β 523}$	.019	.097	1.021	−.031	.088	1.108
$μ_{β 514}$	.059	.243	1.096	.081	.271	1.096	$μ_{β 524}$	.043	.266	1.087	.066	.250	0.996
$μ_{β 515}$	−.019	.156	1.106	.031	.156	0.963	$μ_{β 525}$	−.018	.180	1.016	.058	.151	0.925
$γ_{α 1}$	.046	.086	0.916	.059	.099	0.906	$γ_{α 2}$	.041	.089	0.969	−.037	.074	0.903
$γ_{β 11}$	.028	.055	0.892	.050	.059	1.031	$γ_{β 12}$	.026	.050	1.050	.046	.060	1.008
$γ_{β 21}$	.029	.053	0.925	.039	.045	0.918	$γ_{β 22}$	−.024	.051	1.029	.044	.054	0.960
$γ_{β 31}$	.027	.052	0.864	.046	.049	1.101	$γ_{β 32}$	−.031	.049	1.024	.043	.053	0.865
$γ_{β 41}$	−.018	.126	0.913	.041	.135	0.912	$γ_{β 42}$	.020	.114	1.108	−.017	.133	0.891
$γ_{β 51}$	−.014	.077	0.986	.033	.073	1.047	$γ_{β 52}$	.015	.081	1.032	.048	.072	1.040
$ψ_{α}$	.058	.060	1.106	.059	.067	1.057
$ψ_{β 1}$	.015	.057	1.118	−.001	.048	1.121
$ψ_{β 2}$	−.003	.049	1.179	−.005	.051	1.014
$ψ_{β 3}$	.049	.054	0.956	.031	.053	1.105
$ψ_{β 4}$	.024	.050	1.638	.015	.049	1.588
$ψ_{β 5}$	.025	.069	1.353	.027	.061	1.287

Note. $N = 1, 000$ . RMS = root mean square; SE/SD = standard error/standard deviation.

To evaluate the robustness of MCMC to the initial values of the unknown parameters, we consider two different initial values as follows: Initial 1 (used in simulation 1): the elements of $μ_{r k}$ and $γ$ are 0, and $ψ_{α} = ψ_{β 1} = \dots = ψ_{β 5} = 0.1$ . Initial 2: the elements of $μ_{r k}$ and $γ$ are 2, and $ψ_{α} = ψ_{β 2} = \dots = ψ_{β 5} = 1.5$ . Table 4 (right panel) shows that the estimation results are also close to those reported in Table 2.

Finally, we investigate the robustness of the proposed model to the misspecification of distributional assumptions. We consider a model setup that is the same as simulation 1, except that the distribution of $ε_{i} = (ε_{α}, ε_{β 1}, \dots, ε_{β 5})^{'}$ is specified as follows: (1) $ε_{i j} \sim t_{(3)}$ , (2) $ε_{i j} \sim 0.5 N (- 1, 1) + 0.5 N (1, 1)$ , (3) $ε_{i j} \sim U (- 1, 1)$ , and (4) $ε_{i j} \sim gamma (6, 4) - 1.5$ , where $ε_{i j}$ stands for the j th element of $ε_{i}$ , and $t_{(3)}$ stands for the t distribution with degree of freedom = 3. We simulate 100 data sets with $N = 1, 000$ under each of cases (1)–(4) and conduct parameter estimation using the proposed MCMC which assumes that $ε_{i} \sim N (0, Ψ)$ . The results obtained under the four nonnormal cases are presented in Tables 5 and 6. Most of the parameter estimates (except the SE estimates of the error variances) are robust to this model misspecification. Considering that regression coefficients are of our major interest, we expect that the inaccurate SE estimates of error variances would not distort the inference on the predictor variables in the first- and second-level models.

Table 5.

Estimation Results Under $ε_{i j} \sim t_{(3)}$ and $ε_{i j} \sim 0.5 N (- 1, 1) + 0.5 N (1, 1)$ in Simulation 2.

Para	$t_{(3)}$			0.5N(−1,1) + 0.5N(1,1)			Para	$t_{(3)}$			0.5N(−1,1) + 0.5N(1,1)
Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD	Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD
$μ_{α 11}$	.041	.141	0.890	.050	.133	1.055	$μ_{α 12}$	.022	.141	0.873	.053	.164	0.965
$μ_{β 111}$	−.047	.086	0.979	.007	.075	0.932	$μ_{β 121}$	−.043	.083	1.116	.004	.073	1.042
$μ_{β 112}$	−.046	.095	0.930	−.004	.068	1.019	$μ_{β 122}$	−.038	.101	0.966	−.002	.072	1.100
$μ_{β 113}$	−.018	.097	0.889	−.002	.086	0.864	$μ_{β 123}$	−.001	.121	0.925	.002	.082	0.950
$μ_{β 114}$	−.056	.107	0.923	−.044	.107	0.917	$μ_{β 124}$	−.073	.188	0.904	−.023	.132	0.976
$μ_{β 115}$	−.067	.128	1.003	−.014	.124	0.975	$μ_{β 125}$	−.082	.138	1.092	−.052	.141	0.972
$μ_{α 21}$	−.049	.129	1.061	.079	.116	1.079	$μ_{α 22}$	.060	1.52	1.097	.076	.121	1.035
$μ_{β 211}$	−.047	.098	0.841	−.014	.068	0.944	—	—	—	—	—	—	—
$μ_{β 212}$	−.040	.111	0.863	−.010	.078	0.860	—	—	—	—	—	—	—
$μ_{β 213}$	.009	.091	0.916	.026	.086	0.960	—	—	—	—	—	—	—
$μ_{β 214}$	.029	.083	0.951	−.014	.148	1.061	—	—	—	—	—	—	—
$μ_{β 215}$	−.063	.118	1.059	−.034	.110	0.983	—	—	—	—	—	—	—
$μ_{α 31}$	.020	.195	1.131	.041	.159	1.087	$μ_{α 41}$	.030	.126	0.992	.068	.171	1.075
$μ_{β 311}$	−.030	.133	1.002	.013	.113	1.093	$μ_{β 411}$	−.066	.101	1.107	−.031	.082	1.066
$μ_{β 312}$	−.027	.137	1.081	−.009	.101	1.034	$μ_{β 412}$	−.050	.129	0.966	−.009	.103	1.027
$μ_{β 313}$	−.011	.151	1.054	−.020	.115	1.061	$μ_{β 413}$	−.040	.167	0.879	.010	.099	1.089
$μ_{β 314}$	.023	.137	1.094	−.054	.209	1.129	$μ_{β 414}$	.042	.146	0.953	.059	.210	1.039
$μ_{β 315}$	−.088	.160	1.107	−.043	.208	1.117	$μ_{β 415}$	−.077	.178	1.008	−.009	.151	1.137
$μ_{α 51}$	.034	.148	1.064	.038	.180	1.092	$μ_{α 52}$	−.009	.094	1.109	.052	.108	1.137
$μ_{β 511}$	−.028	.108	1.070	−.018	.083	1.082	$μ_{β 521}$	−.039	.115	1.100	−.007	.082	1.097
$μ_{β 512}$	.005	.122	1.069	−.008	.093	1.105	$μ_{β 522}$	−.050	.133	0.991	.010	.093	1.088
$μ_{β 513}$	−.019	.141	1.010	.012	.108	1.009	$μ_{β 523}$	−.016	.148	0.969	.013	.106	1.029
$μ_{β 514}$	.014	.167	1.057	.043	.218	1.140	$μ_{β 524}$	.046	.132	1.058	.056	.219	1.068
$μ_{β 515}$	−.063	.190	1.125	−.029	.160	1.060	$μ_{β 525}$	−.069	.173	1.140	−.034	.172	1.091
$γ_{α 1}$	−.041	.136	0.926	.042	.101	0.886	$γ_{α 2}$	.011	.084	0.919	.025	.096	0.906
$γ_{β 11}$	−.062	.086	0.915	.005	.063	0.913	$γ_{β 12}$	−.053	.115	0.960	.009	.058	0.940
$γ_{β 21}$	−.063	.087	0.909	−.004	.051	1.095	$γ_{β 22}$	−.065	.133	0.923	.003	.056	1.018
$γ_{β 31}$	.021	.085	0.878	.008	.056	0.944	$γ_{β 32}$	−.055	.148	0.924	.016	.055	0.940
$γ_{β 41}$	−.075	.172	0.917	−.027	.152	0.883	$γ_{β 42}$	−.079	.132	1.105	−.019	.088	0.875
$γ_{β 51}$	−.075	.110	0.963	−.022	.079	1.083	$γ_{β 52}$	−.049	.093	0.934	−.043	.143	0.926
$ψ_{α}$	.570	.165	0.810	.190	.085	0.955
$ψ_{β 1}$	.399	.115	0.720	.104	.073	0.945
$ψ_{β 2}$	.391	.143	0.801	.112	.074	0.958
$ψ_{β 3}$	.665	.135	0.740	.147	.063	1.065
$ψ_{β 4}$	.036	.043	0.839	.041	.052	1.351
$ψ_{β 5}$	.293	.165	0.800	.079	.073	1.203

Note. $N = 1, 000$ . RMS = root mean square; SE/SD = standard error/standard deviation.

Table 6.

Estimation Results Under $ε_{i j} \sim U (- 1, 1)$ and $ε_{i j} \sim gamma (6, 4) - \frac{3}{2}$ in Simulation 2.

Para	$U (- 1, 1)$			$Gamma (6, 4) - \frac{3}{2}$			Para	$U (- 1, 1)$			$Gamma (6, 4) - \frac{3}{2}$
Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD	Para	Bias	RMS	SE/SD	Bias	RMS	SE/SD
$μ_{α 11}$	.049	.169	0.816	.039	.124	0.871	$μ_{α 12}$	.054	.180	0.929	.032	.141	0.869
$μ_{β 111}$	.023	.071	1.043	.014	.072	0.973	$μ_{β 121}$	.028	.082	0.942	.011	.085	0.942
$μ_{β 112}$	.024	.066	1.057	.036	.076	0.955	$μ_{β 122}$	.027	.067	1.110	.004	.079	1.005
$μ_{β 113}$	.026	.068	0.990	.004	.081	0.882	$μ_{β 123}$	.018	.072	0.955	.014	.084	0.929
$μ_{β 114}$	.015	.211	0.922	−.003	.108	0.948	$μ_{β 124}$	.020	.231	0.888	.025	.138	0.9330
$μ_{β 115}$	.006	.123	1.049	−.012	.122	1.080	$μ_{β 125}$	.005	.129	1.055	−.035	.130	1.092
$μ_{α 21}$	.078	.107	1.043	.059	.121	1.012	$μ_{α 22}$	.069	.112	1.034	.071	.116	1.002
$μ_{β 211}$	.017	.058	0.992	−.006	.068	0.929	—	—	—	—	—	—	—
$μ_{β 212}$	.009	.069	0.849	−.001	.069	1.006	—	—	—	—	—	—	—
$μ_{β 213}$	.014	.074	0.937	.005	.088	0.903	—	—	—	—	—	—	—
$μ_{β 214}$	.055	.158	1.007	.022	.158	0.972	—	—	—	—	—	—	—
$μ_{β 215}$	−.028	.094	1.034	−.036	.093	1.074	—	—	—	—	—	—	—
$μ_{α 31}$	.042	.167	0.987	.051	.163	1.110	$μ_{α 41}$	.077	.164	1.127	.076	.200	1.004
$μ_{β 311}$	.023	.095	1.085	.016	.132	0.990	$μ_{β 411}$	.023	.093	1.091	.019	.092	1.100
$μ_{β 312}$	.015	.110	1.118	.004	.122	1.072	$μ_{β 412}$	−.005	.079	1.095	.008	.085	1.112
$μ_{β 313}$	−.002	.109	1.102	.027	.139	0.972	$μ_{β 413}$	.010	.088	1.107	.002	.101	1.010
$μ_{β 314}$	−.043	.230	1.063	.024	.194	1.109	$μ_{β 414}$	.049	.262	1.027	.027	.282	1.044
$μ_{β 315}$	−.011	.222	1.020	−.011	.184	1.021	$μ_{β 415}$	−.025	.149	1.179	.017	.169	1.088
$μ_{α 51}$	.013	.155	1.091	.072	.172	1.091	$μ_{α 52}$	.076	.162	1.097	.071	.159	1.118
$μ_{β 511}$	−.043	.086	1.085	−.005	.093	1.096	$μ_{β 521}$	−.108	.087	1.137	.012	.108	0.947
$μ_{β 512}$	.009	.092	1.107	.036	.094	1.117	$μ_{β 522}$	.016	.098	1.043	.004	.099	1.037
$μ_{β 513}$	.034	.116	0.880	.017	.102	1.002	$μ_{β 523}$	.007	.099	1.105	.015	.103	0.994
$μ_{β 514}$	.047	.208	1.046	.037	.258	1.124	$μ_{β 524}$	.079	.232	1.027	.063	.277	1.042
$μ_{β 515}$	−.034	.171	1.066	−.012	.183	0.989	$μ_{β 525}$	−.023	.164	1.102	−.030	.151	1.106
$γ_{α 1}$	.039	.074	1.061	.045	.129	0.942	$γ_{α 2}$	.028	.086	0.908	.057	.155	1.053
$γ_{β 11}$	.023	.048	1.029	.018	.057	0.930	$γ_{β 12}$	.052	.048	1.024	.029	.085	0.957
$γ_{β 21}$	.029	.057	0.864	.016	.058	0.993	$γ_{β 22}$	.015	.051	0.975	.040	.097	0.958
$γ_{β 31}$	.050	.049	0.913	.021	.155	0.939	$γ_{β 32}$	.048	.054	0.840	.035	.158	0.909
$γ_{β 41}$	.010	.123	0.975	.063	.101	0.891	$γ_{β 42}$	.031	.141	0.875	.017	.152	0.842
$γ_{β 51}$	.053	.079	0.945	−.050	.176	1.141	$γ_{β 52}$	.017	.083	0.907	−.038	.080	1.005
$ψ_{α}$	.051	.059	1.211	.092	.079	1.001
$ψ_{β 1}$	.029	.059	1.052	.031	.065	0.956
$ψ_{β 2}$	.027	.064	0.947	.036	.069	0.921
$ψ_{β 3}$	.061	.054	1.065	.054	.074	0.773
$ψ_{β 4}$	.024	.056	1.466	.032	.065	1.525
$ψ_{β 5}$	.025	.055	1.270	.046	.071	1.249

Note. $N = 1, 000$ . RMS = root mean square; SE/SD = standard error/standard deviation.

A Study on the Partially Ordered Adolescent Smoking Behavior Dataset

In this section, we illustrate the proposed methodology by using a longitudinal data set extracted from NLSY (Center for Human Resource Research 2003). The NLSY is one of the most comprehensive longitudinal studies of youth with a representative sample of adolescent respondents aged from 12 to 18 years at the time of the first interview in 1997 and respondents aged from 28 to 34 years at the 16th round (the last round) of interviews in 2013–2014. We focused on years 1998, 1999, 2000, 2002, and 2004 because the questionnaire design in these years is relatively consistent.

This analysis focused on adolescent smokers, who can be classified according to their different smoking behaviors. Following the medical and sociological literature (Ham and Hope 2003; Zhang and Ip 2012), we applied three related questions, namely, “ever smoked,” “number of cigarettes smoked per day in the last 30 days,” and “number of smoking days in the last 30 days,” to form the categories of our outcome smoking behavior variable, including nonsmoker, former smoker, light-nonfrequent smoker (smoked on $\leq 3$ days in the past 30 days and consumed ≤3 cigarettes on each day), light-frequent smoker (smoked on $> 3$ days in the past 30 days and consumed ≤3 cigarettes on each day), heavy-nonfrequent smoker (smoked on $\leq 3$ days in the past 30 days and consumed $> 3$ cigarettes on each day), and heavy-frequent smoker (smoked on $> 3$ days in the past 30 days and consumed $> 3$ cigarettes on each day). Zhang and Ip (2012) considered a partially ordered structure of smoking categories, namely, {nonsmoker}, {former smoker, light-nonfrequent smoker}, {light-frequent smoker, heavy-nonfrequent smoker}, and {heavy-frequent smoker}, in which former smoker and light-nonfrequent smoker as well as light-frequent smoker and heavy-nonfrequent smoker were treated as incomparable. However, they ignored the longitudinal nature of the smoking behavior variable in NLSY. In this study, we followed the suggestion of an anonymous referee to regard {former smoker} as a lighter smoking category comparing to {light nonfrequent smoker}. The proposed poset structure of smoking categories is presented in Figure 4. Given that different adolescent smoking behaviors may be related to different risk factors, we identified the unique risk factors of each smoking behavior while keeping all ordinal or nominal information among the categories in the longitudinal setting.

Figure 4.

The partially ordered set structure of adolescent smoking categories in the National Longitudinal Surveys of Youth dataset.

The following $s = 5$ time-varying covariates were considered as potential influential factors of adolescent smoking behavior in the first-level model: age ( $x_{i t 1}$ ), enrollment ( $x_{i t 2}$ − education enrollment status), household ( $x_{i t 3}$ − household size), parental help ( $x_{i t 4}$ − whether to turn to parents for help with personal problems, 1 = yes and 0 = no), and residence ( $x_{i t 5}$ − whether to live in a rural area, 1 = urban and 0 = rural). The first three continuous longitudinal covariates were standardized according to the method proposed by Little (2013). We included $q = 2$ individual-level covariates, namely, gender ( $z_{i 1}$ , 1 = female and 0 = male) and race ( $z_{i 2}$ , 1 = white and 0 = nonwhite), in the second-level model to examine the interdependence among the repeated responses and to consider the possible source of the correlation. A detailed description of the responses and covariates is provided in Online Appendix B (which can be found at http://smr.sagepub.com/supplemental/). A total of $N = 2, 563$ youths with $T = 5$ years of complete records of the variables were included in the data set.

We applied the following two-level model to analyze the described longitudinally and partially ordered adolescent smoking behavior data set. For the first-level model, the partially ordered smoking behavior variable only has one disjoint component, which means that $Y^{(1)}$ only takes the value of 0. Therefore, we considered the conditional models of $[Y^{(2)} | Y^{(1)}]$ and $[Y^{(3)} | Y^{(1)}, Y^{(2)}]$ . $[Y^{(2)} | Y^{(1)} = 0]$ is ordinally distributed across categories 1–5 that represent totally weekly ordered antichains {nonsmoker}, {former smoker}, {light-nonfrequent smoker}, {light-frequent smoker, heavy-nonfrequent smoker}, and {heavy-frequent smoker}, respectively; and $[Y^{(3)} | Y^{(1)} = 0, Y^{(2)} = 4]$ is nominally distributed across categories 0 and 1 that denote {light-frequent smoker} and {heavy-nonfrequent smoker}, respectively. Similar to the model formulation presented in Two-level Model for Repeated Partially Ordered Responses subsection, the following hierarchical models were considered.

For $[Y^{(2)} = 1, 2, 3, 4, 5 | Y^{(1)} = 0]$ , we proposed the proportional odds model:

log (\frac{\sum_{j \leq k} π_{i t j}^{(2.0)}}{\sum_{j > k} π_{i t j}^{(2.0)}}) = α_{i 1 k} + β'_{i 11} x_{i t}, k = 1, 2, 3, 4,

where $x_{i t} = (x_{i t 1}, \dots, x_{i t 5})$ = (age, enrollment, household, parental help, residence), and $β_{i 11}$ is invariant to k because of the proportional odds assumption imposed on the cumulative logit model (8).

For $[Y^{(3)} = 0, 1 | Y^{(1)} = 0, Y^{(2)} = 4]$ , where 0 refers to the category {light-frequent smoker} and 1 refers to the category {heavy-nonfrequent smoker}, we proposed the conditional nominal model:

log (\frac{π_{i t 1}^{(3.04)}}{π_{i t 0}^{(3.04)}}) = α_{i 21} + β'_{i 21} x_{i t} .

The second-level model to examine the influence of time-invariant individual characteristics, such as gender and race, on the varying coefficients $α_{i r k}$ and $β_{i r k}$ was proposed as follows: for $r = 1, 2$ ,

α_{i r k} = μ_{α r k} + γ_{α 1} z_{i 1} + γ_{α 2} z_{i 2} + ε_{α i}, k = 1, \dots, K_{r} - 1,

β_{i r 1} = μ_{β r 1} + γ_{β 1} z_{i 1} + γ_{β 2} z_{i 2} + ε_{β i},

where $z_{i} = (z_{i 1}, z_{i 2})^{'} = (g e n d e r, r a c e)'$ ; $μ_{α r k}$ , $γ_{α 1}$ , $γ_{α 2}$ , and $ε_{α i}$ are scalars; $β_{i r 1}$ , $μ_{β r 1}$ , $γ_{β 1}$ , $γ_{β 2}$ , and $ε_{β i}$ are five-dimensional vectors; $K_{1} = 5$ ; and $K_{2} = 2$ .

The proposed Bayesian method was applied to conduct the statistical inference. The hyperparameters in the prior distributions were set the same as prior 1 in simulation 1. The algorithm had converged within 10,000 iterations in the analysis of this data set. Thus, we collected 10,000 posterior observations after discarding the first 10,000 burn-in iterations. The estimation results of the model parameters, together with their SE estimates, are presented in Table 7. We also conducted the analysis using prior 2 specified in simulation 2. The estimation results are similar to those reported in Table 7 and are not reported. All the estimated variances ${\hat{ψ}}_{α}$ and ${\hat{ψ}}_{β j}, j = 1, \dots, s$ are substantially different from 0, which shows the strong interdependency of the repeated measurements nested within the same individual. In the second-level model (10), ${\hat{γ}}_{α 1} = 1.01 (0.20)$ and ${\hat{γ}}_{α 2} = - 0.80 (0.23)$ . Thus, female ( $z_{i 1} = 1$ ) and nonwhite ( $z_{i 2} = 0$ ) teenagers have larger values of $α_{i r k}$ , which in turn result in high log odds in the first-level models (8) and (9). In the second-level model (11), ${\hat{γ}}_{β 11} = 1.16 (0.20)$ and ${\hat{γ}}_{β 21} = - 0.41 (0.20)$ , thereby suggesting that only the effects of age and enrollment status on the log odds in the first-level models (8) and (9) change according to gender. In addition, ${\hat{γ}}_{β 12} = - 1.36 (0.21)$ , ${\hat{γ}}_{22} = - 0.58 (0.11)$ , ${\hat{γ}}_{β 32} = 0.49 (0.07)$ , ${\hat{γ}}_{β 42} = 0.26 (0.04)$ , and ${\hat{γ}}_{β 52} = - 0.42 (0.14)$ , thereby indicating that all the five predictor effects on the log odds in the first-level models (8) and (9) change according to race. The specific interpretations of the estimated coefficients are provided below.

Table 7.

Regression Coefficients in the Analysis of the Adolescent Smoking Behavior Dataset.

Coefficients		Ordinal Model (8)			Nominal Model (9)
Coefficients		Para	Estimation	SE	Para	Estimation	SE
Random intercepts ( $α_{i r k}$ )	Intercept	$μ_{α 11}$	1.86	0.22	$μ_{α 21}$	2.43	0.26
		$μ_{α 12}$	4.65	0.17
		$μ_{α 13}$	5.59	0.11
		$μ_{α 14}$	7.06	0.23
	Gender ( $z_{1}$ : 0 = male)	$γ_{α 1}$	1.01	0.20	$γ_{α 1}$	1.01	0.20
	Race ( $z_{2}$ : 0 = nonwhite)	$γ_{α 2}$	−0.80	0.23	$γ_{α 2}$	−0.80	0.23
	Error ( $ε_{α}$ )	$ψ_{α}$	4.54	1.45	$ψ_{α}$	4.54	1.45
Random coefficients ( $β_{i r 1}$ )
Age ( $x_{1}$ : years)	Intercept	$μ_{β 111}$	−3.44	0.15	$μ_{β 211}$	−3.58	0.30
	Gender ( $z_{1}$ )	$γ_{β 11}$	1.16	0.20	$γ_{β 11}$	1.16	0.20
	Race ( $z_{2}$ )	$γ_{β 12}$	−1.36	0.21	$γ_{β 12}$	−1.36	0.21
	Error ( $ε_{β 1}$ )	$ψ_{β 1}$	4.83	1.03	$ψ_{β 1}$	4.83	1.03
Enrollment ( $x_{2}$ : 1 = lowest to 11 = highest)	Intercept	$μ_{β 112}$	2.73	0.20	$μ_{β 212}$	0.39	0.34
	Gender ( $z_{1}$ )	$γ_{β 21}$	−0.41	0.20	$γ_{β 21}$	−0.41	0.20
	Race ( $z_{2}$ )	$γ_{β 22}$	−0.58	0.11	$γ_{β 22}$	−0.58	0.11
	Error ( $ε_{β 2}$ )	$ψ_{β 2}$	3.86	0.27	$ψ_{β 2}$	3.86	0.27
Household size ( $x_{3}$ : 1 = smallest to 20 = largest)	Intercept	$μ_{β 113}$	0.33	0.32	$μ_{β 213}$	−1.52	0.44
	Gender ( $z_{1}$ )	$γ_{β 31}$	0.15	0.11	$γ_{β 31}$	0.15	0.11
	Race ( $z_{2}$ )	$γ_{β 32}$	0.49	0.07	$γ_{β 32}$	0.49	0.07
	Error ( $ε_{β 3}$ )	$ψ_{β 3}$	0.42	0.18	$ψ_{β 3}$	0.42	0.18
Parental help ( $x_{4}$ : 0 = no)	Intercept	$μ_{β 114}$	0.36	0.10	$μ_{β 214}$	0.59	0.36
	Gender ( $z_{1}$ )	$γ_{β 41}$	−0.05	0.11	$γ_{β 41}$	−0.05	0.11
	Race ( $z_{2}$ )	$γ_{β 42}$	0.26	0.04	$γ_{β 42}$	0.26	0.04
	Error ( $ε_{β 4}$ )	$ψ_{β 4}$	0.28	0.10	$ψ_{β 4}$	0.28	0.10
Residence ( $x_{5}$ : 0 = rural)	Intercept	$μ_{β 115}$	0.84	0.07	$μ_{β 215}$	−0.20	0.29
	Gender ( $z_{1}$ )	$γ_{β 51}$	−0.06	0.12	$γ_{β 51}$	−0.06	0.12
	Race ( $z_{2}$ )	$γ_{β 52}$	−0.42	0.14	$γ_{β 52}$	−0.42	0.14
	Error ( $ε_{β 5}$ )	$ψ_{β 5}$	0.41	0.18	$ψ_{β 5}$	0.41	0.18

We start from the conditional ordinal model (8) of $[Y^{(2)} | Y^{(1)} = 0]$ with five ordinal categories, {nonsmoker}, {former smoker}, {light-nonfrequent smoker}, {light-frequent smoker, heavy-nonfrequent smoker}, and {heavy-frequent smoker}. Given the positive value of ${\hat{γ}}_{α 1}$ and the negative value of ${\hat{γ}}_{α 2}$ , female ( $z_{i 1} = 1$ ) and nonwhite ( $z_{i 2} = 0$ ) teenagers are more likely to be lighter smoker than their male and white counterparts. The results for $β_{i 1 k}$ further show the effects of $x_{i t} = (x_{i t 1}, \dots, x_{i t 5})^{'} = (a g e, e n r o l l m e n t, h o u s e h o l d, p a r e n t a l h e l p, r e s i d e n c e)^{'}$ on adolescent smoking behavior and how these effects vary across gender and/or race. In our model setting, the cohort of male ( $z_{i 1} = 0$ ) and nonwhite ( $z_{i 2} = 0$ ) adolescents serves as the baseline cohort. The estimated mean effects of age ( $x_{i t 1}$ ), education ( $x_{i t 2}$ ), household ( $x_{i t 3}$ ), parental help ( $x_{i t 4}$ ), and residence ( $x_{i t 5}$ ) for the baseline cohort are ${\hat{μ}}_{β 111} = - 3.44 (0.15)$ , ${\hat{μ}}_{β 112} = 2.73 (0.20)$ , ${\hat{μ}}_{β 113} = 0.33 (0.32)$ , ${\hat{μ}}_{β 114} = 0.36 (0.10)$ , and ${\hat{μ}}_{β 115} = 0.84 (0.07)$ , which imply that those male and nonwhite teenagers who are younger, enrolled in higher educational programs, with access to parental help, or came from urban areas tend to have less abusive cigarette smoking behaviors, and that the effect of household size is not substantial. Similar conclusions can be obtained for other three cohorts, such as male ( $z_{i 1} = 0$ ) and white ( $z_{i 2} = 1$ ) cohort, female ( $z_{i 1} = 1$ ) and nonwhite ( $z_{i 2} = 0$ ) cohort, and female ( $z_{i 1} = 1$ ) and white ( $z_{i 2} = 1$ ) cohort. These findings agree with the results reported in the literature (see, e.g., Alexander et al. 2001; Blum et al. 2000; Botvin et al. 1992; Dusenbury et al. 1992; Maccoby and Mnookin 1992; Patten et al. 2003; Plotnikoff, Bercovitz, and Loucaides 2004; Repetti, Taylor, and Seeman 2002b; Tseng et al. 2001; Tyas and Pederson 1998; Wetter et al. 1999; Wills, Vaccaro, and McNamara 1992). Despite similar conclusions for distinct teenager cohorts, the magnitudes of the covariate effects vary according to gender and race. For instance, for male ( $z_{i 1} = 0$ ) and white ( $z_{i 2} = 1$ ) cohort, the estimated mean effects of $x_{i t}$ become ${\hat{μ}}_{β 111} + {\hat{γ}}_{β 12} = - 3.44 - 1.36 = - 4.80$ , ${\hat{μ}}_{β 112} + {\hat{γ}}_{β 22} = 2.73 - 0.58 = 2.15$ , ${\hat{μ}}_{β 113} + {\hat{γ}}_{β 32} = 0.33 + 0.49 = 0.82$ , ${\hat{μ}}_{β 114} + {\hat{γ}}_{β 42} = 0.36 + 0.26 = 0.62$ , and ${\hat{μ}}_{β 115} + {\hat{γ}}_{β 52} = 0.84 - 0.42 = 0.42$ . Apparently, the effect of age on adolescent smoking behavior is much larger for the male and white cohort than for the baseline cohort, whereas the effects of other factors also vary from this cohort to the baseline cohort.

For the conditional nominal model (9) of $[Y^{(3)} | Y^{(1)} = 0, Y^{(2)} = 4]$ with two incomparable categories, {light-frequent smoker} and {heavy-nonfrequent smoker}, the results of the second-level models for $β_{i 21}$ are explained as follows. For the male ( $z_{i 1} = 0$ ) and nonwhite ( $z_{i 2} = 0$ ) cohort, the estimated mean effects of age and household size are ${\hat{μ}}_{β 211} = - 3.58 (0.30)$ and ${\hat{μ}}_{β 213} = - 1.52 (0.44)$ , while the effects of education, parental help, and residence area are insignificant. Those male and nonwhite teenagers who are younger or with a larger household size tend to be light-frequent smokers conditional on that they are either light-frequent smokers or heavy-nonfrequent smokers. A possible reasoning behind this finding is that light-frequent smoking is more of a regular habit, while heavy-nonfrequent smoking seems to be an occasional risky behavior; adolescents who are older or with a larger family tend to have better emotional competence and are thus less likely to go to extremes (Hessler and Katz 2010; Morand 1999). Likewise, the effect of age on adolescent smoking behavior is the largest for the male and white cohort but the smallest for the female and nonwhite cohort, whereas the effect of household size on adolescent smoking behavior is the largest for the male and nonwhite cohort but the smallest for the female and white cohort.

In summary, the above analysis revealed how teenagers’ age, education enrollment status, household size, parental help, and residence area influence their smoking behavior. The results of the second-level model provided additional insights into how these effects vary according to certain individuals’ characteristics, such as gender and race. These insights cannot be revealed by the conventional single-level regression.

Discussion

In this study, we propose a Bayesian two-level model to analyze longitudinal data with a partially ordered structure. The model comprises a first-level random coefficients model to investigate the relationship between time-varying covariates and partially ordered outcomes and a second-level regression model to examine the effects of individual characteristics on the random coefficients. The proposed model accommodates the exact structure of the poset so that all the ordinal information among the categories is kept and the incomparable categories among the same ordinal level can be well separated in the analysis of the partially ordered data. The model also accounts for the interdependency among the repeated partially ordered responses and examines how individual characteristics affect the first-level random coefficients. The simulation results show that the proposed method performs well. Applying the proposed method to the longitudinal adolescent smoking behavior data set provides new insights into the effects of risk factors for adolescent smoking behaviors that are characterized through partially ordered categories. Such data structure cannot be accommodated by the existing methods. The analysis results boost our comprehensive understanding of adolescent smoking behaviors and help us come up with an effective strategy to prevent adolescents from smoking.

Though other methods, such as GEE (generalized estimating equation) and robust SE, can do the same or similar analyses without causing much numerical problems, they may encounter difficulty when the model and data structures become increasingly complicated. For instance, in the presence of missing data, GEE can accommodate only the missing data that are missing completely at random. By contrast, MCMC can easily manage various types of missing data by adding additional posterior sampling step(s). Thus, the proposed MCMC method has a high potential to be extended to address various challenging problems for which other methods are either unfeasible or work poorly.

In this study, the partially ordered structure of the response is constructed based on the relations between categories in the poset. Such a poset structure must be determined beforehand according to the nature and aim of the substantive study and/or experts’ knowledge. This structure may vary according to the objective of the study. For instance, in the study of adolescent smoking behavior, {former smoker} can be regarded as a lighter category than {light-nonfrequent smoker} if our focus is on adolescents’ misconduct behavior. However, if the primary goal of a study is to investigate the effect of smoking on the risk of lung cancer, then no clear evidence can tell how to rank the above two smoking categories because which category is associated with a higher risk of lung cancer is unknown.

Notably, the parameters in the conditional models for partially ordered data only allow for conditional interpretations. Moreover, we assume that the regression coefficients in the second-level model (6) are invariant to r and k. Although extensions to allow r- and/or k-specific regression coefficients are straightforward, such extensions dramatically increase the number of parameters and thus cause difficulties in interpretation.

Supplemental Material

Supplemental Material, Appendix_(1) - Bayesian Two-level Model for Repeated Partially Ordered Responses: Application to Adolescent Smoking Behavior Analysis

Supplemental Material, Appendix_(1) for Bayesian Two-level Model for Repeated Partially Ordered Responses: Application to Adolescent Smoking Behavior Analysis by Xiaoqing Wang, Haotian Wu, Xiangnan Feng, and Xinyuan Song in Sociological Methods & Research

Footnotes

Acknowledgment

The authors thank the editor, the associate editor, and the two reviewers for their valuable comments, which substantially improved this article.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by GRF grants 14301918 and 14303017 from the Research Grant Council of HKSAR, the National Natural Science Foundation of China (11671054, 11871263, 71802166), and the Direct Grant of the Chinese University of Hong Kong.

ORCID iD

Xinyuan Song

Supplemental Material

Supplemental material for this article is available online.

References

Alexander

Piazza

Mekos

Valente

. 2001. “Peers, Schools, and Adolescent Cigarette Smoking.” Journal of Adolescent Health 29:22–30.

Ansari

Jedidi

. 2000. “Bayesian Factor Analysis for Multilevel Binary Observations.” Psychometrika 65:475–96.

Birkhoff

. 1940. Lattice Theory. Vol. 25. Providence, RI: American Mathematical Society.

Blum

R. W.

Beuhring

Shew

M. L.

Bearinger

L. H.

Sieving

R. E.

Resnick

M. D.

. 2000. “The Effects of Race/Ethnicity, Income, and Family Structure on Adolescent Risk Behaviors.” American Journal of Public Health 90:1879.

Botvin

G. J.

Baker

Goldberg

C. J.

Dusenbury

Botvin

E. M.

. 1992. “Correlates and Predictors of Smoking among Black Adolescents.” Addictive Behaviors 17:97–103.

Center for Human Resource Research. 2003. A Guide to the Rounds 1-5 Data, National Longitudinal Survey of Youth, 1997. Columbus, OH: U.S. Department of Labor, Bureau of Labor Statistics.

Congdon

2007. Bayesian Statistical Modelling. Hoboken, NJ: John Wiley & Sons.

Cormen

T. H.

2009. Introduction to Algorithms. Cambridge, England: MIT Press.

Diggle

2002. Analysis of Longitudinal Data. Oxford, England: Oxford University Press.

10.

Dunson

D. B.

2000. “Bayesian Latent Variable Models for Clustered Mixed Outcomes.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 62:355–66.

11.

Dusenbury

Kerner

J. F.

Baker

Botvin

James-Ortiz

Zauber

. 1992. “Predictors of Smoking Prevalence among New York Latino Youth.” American Journal of Public Health 82:55–58.

12.

Dushnik

Miller

E. W.

. 1941. “Partially Ordered Sets.” American Journal of Mathematics 63:600–10.

13.

Feng

X.-N.

H.-T.

Song

X.-Y.

. 2017. “Bayesian Adaptive Lasso for Ordinal Regression with Latent Variables.” Sociological Methods & Research 46:926–53.

14.

Geman

. 1984. “Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images.” IEEE Transactions on Pattern Analysis and Machine Intelligence 6:721.

15.

Ham

L. S.

Hope

D. A.

. 2003. “College Students and Problematic Drinking: A Review of the Literature.” Clinical Psychology Review 23:719–59.

16.

Hastings

W. K.

1970. “Monte Carlo Sampling Methods Using Markov Chains and Their Applications.” Biometrika 57:97–109.

17.

Hedeker

Gibbons

R. D.

. 2006. Longitudinal Data Analysis. Vol. 451. New York, NY: John Wiley.

18.

Hessler

D. M.

Katz

L. F.

. 2010. “Brief Report: Associations between Emotional Competence and Adolescent Risky Behavior.” Journal of Adolescence 33:241–46.

19.

Hox

J. J.

Moerbeek

van de Schoot

. 2010. Multilevel Analysis: Techniques and Applications. New York, NY: Routledge.

20.

E. H.

Chen

S.-H.

Quandt

S. A.

. 2016. “Analysis of Multiple Partially Ordered Responses to Belief Items with Don’t Know Option.” Psychometrika 81:483–505.

21.

E. H.

Zhang

Rejeski

Harris

Kritchevsky

. 2013. “Partially Ordered Mixed Hidden Markov Model for the Disablement Process of Older Adults.” Journal of the American Statistical Association 108:370–84.

22.

Laird

N. M.

Ware

J. H.

. 1982. “Random-effects Models for Longitudinal Data.” Biometrics 38:963–74.

23.

Lee

S.-Y.

Song

X.-Y.

. 2004. “Evaluation of the Bayesian and Maximum Likelihood Approaches in Analyzing Structural Equation Models with Small Sample Sizes.” Multivariate Behavioral Research 39:653–86.

24.

Little

T. D.

2013. Longitudinal Structural Equation Modeling. New York, NY: Guilford Press.

25.

Longford

N. T

. 1995. “Random Coefficient Models.” Pp. 519–70 in Handbook of Statistical Modeling for the Social and Behavioral Sciences, edited by Gerhard

Clifford

Sobel

M. E.

. Heidelberg, Germany: Springer.

26.

Maccoby

E. E.

Mnookin

R. H.

. 1992. Dividing the Child: Social and Legal Dilemmas of Custody. Cambridge, England: Harvard University Press.

27.

Metropolis

Rosenbluth

A. W.

Rosenbluth

M. N.

Teller

A. H.

Teller

. 1953. “Equation of State Calculations by Fast Computing Machines.” The Journal of Chemical Physics 21:1087–92.

28.

Meulders

E. H.

Boeck

. 2005. “Latent Variable Models for Partially Ordered Responses and Trajectory Analysis of Anger-related Feelings.” British Journal of Mathematical and Statistical Psychology 58:117–43.

29.

Morand

D. A.

1999. “Family Size and Intelligence Revisited: The Role of Emotional Intelligence.” Psychological Reports 84:643–49.

30.

Muñoz

M. A.

Chang

F. C.

. 2007. “The Elusive Relationship between Teacher Characteristics and Student Academic Growth: A Longitudinal Multilevel Model for Change.” Journal of Personnel Evaluation in Education 20:147–64.

31.

Patten

C. A.

Offord

K. P.

Ames

S. C.

Decker

P. A.

Croghan

I. T.

Dornelas

E. A.

Pingree

Boberg

E. W.

Gustafson

D. H.

Ahluwalia

J. S.

Wolter

T. D.

Hurt

R. D.

. 2003. “Differences in Adolescent Smoker and Nonsmoker Perceptions of Strategies That Would Help an Adolescent Quit Smoking.” Annals of Behavioral Medicine 26:124–33.

32.

Petraitis

Flay

B. R.

Miller

T. Q.

. 1995. “Reviewing Theories of Adolescent Substance Use: Organizing Pieces in the Puzzle.” Psychological Bulletin 117:67.

33.

Peyhardi

Trottier

Guédon

. 2016. “Partitioned Conditional Generalized Linear Models for Categorical Responses.” Statistical Modelling 16(4):297–321.

34.

Plotnikoff

R. C.

Bercovitz

Loucaides

C. A.

. 2004. “Physical Activity, Smoking, and Obesity among Canadian School Youth: Comparison between Urban and Rural Schools.” Canadian Journal of Public Health/Revue Canadienne de Sante’e Publique 95:413–18.

35.

Raudenbush

S. W.

Bryk

A. S.

. 2002. Hierarchical Linear Models: Applications and Data Analysis Methods. Vol. 1. Thousand Oaks, CA: Sage.

36.

Repetti

R. L.

Taylor

S. E.

Seeman

T. E.

. 2002a. “Risky Families: Family Social Environments and the Mental and Physical Health of Offspring.” Psychological Bulletin 128:330.

37.

Repetti

R. L.

Taylor

S. E.

Seeman

T. E.

. 2002b. “Risky Families: Family Social Environments and the Mental and Physical Health of Offspring.” Psychological Bulletin 128:330.

38.

Schane

R. E.

Ling

P. M.

Glantz

S. A.

. 2010. “Health Effects of Light and Intermittent Smoking.” Circulation 121:1518–22.

39.

Skrondal

Rabe-Hesketh

. 2004. Generalized Latent Variable Modeling: Multilevel, Longitudinal, and Structural Equation Models. Boca Raton, Florida: Chapman and Hall/CRC.

40.

Song

X.-Y.

Z.-H.

Cai

J.-H.

E. H.-S.

. 2013. “A Bayesian Modeling Approach for Generalized Semiparametric Structural Equation Models.” Psychometrika 78:624–47.

41.

Tanner

M. A.

Wong

W. H.

. 1987. “The Calculation of Posterior Distributions by Data Augmentation.” Journal of the American Statistical Association 82:528–40.

42.

Tseng

Yeatts

Millikan

Newman

. 2001. “Area-level Characteristics and Smoking in Women.” American Journal of Public Health 91:1847–50.

43.

Tyas

S. L.

Pederson

L. L.

. 1998. “Psychosocial Factors Related to Adolescent Smoking: A Critical Review of the Literature.” Tobacco Control 7:409–20.

44.

Verbeke

Molenberghs

. 2009. Linear Mixed Models for Longitudinal Data. New York, NY: Springer Science & Business Media.

45.

Weschler

Davenport

Dowdall

Moeykens

Castillo

. 1994. “Health and Behavioral Consequences of Binge Drinking in College.” Journal of the American Medical Association 272:1672–77.

46.

Wetter

D. W.

Kenford

S. L.

Smith

S. S.

Fiore

M. C.

Jorenby

D. E.

Baker

T. B.

. 1999. “Gender Differences in Smoking Cessation.” Journal of Consulting and Clinical Psychology 67:555.

47.

Willms

J. D.

Raudenbush

S. W.

. 1989. “A Longitudinal Hierarchical Linear Model for Estimating School Effects and Their Stability.” Journal of Educational Measurement 26:209–32.

48.

Wills

T. A.

Vaccaro

McNamara

. 1992. “The Role of Life Events, Family Support, and Competence in Adolescent Substance Use: A Test of Vulnerability and Protective Factors.” American Journal of Community Psychology 20:349–74.

49.

Wilson

1992. “The Ordered Artition Model: An Extension of the Partial Credit Model.” Applied Psychological Measurement 16:309–25.

50.

Zhang

E. H.

. 2012. “Generalized Linear Model for Partially Ordered Data.” Statistics in Medicine 31:56–68.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.23 MB