A Bayesian semiparametric dynamic two-level structural equation model for analyzing non-normal longitudinal data

Journal of Multivariate Analysis 121 (2013) 87–108

Contents lists available at ScienceDirect

Journal of Multivariate Analysis journal homepage: www.elsevier.com/locate/jmva

A Bayesian semiparametric dynamic two-level structural equation model for analyzing non-normal longitudinal data Song Xin-Yuan a,∗ , Chen Fei b , Lu Zhao-Hua a a

Department of Statistics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong

b

Statistics & Mathematics College, Yunnan University of Finance and Economics, China

article

info

Article history: Received 30 May 2012 Available online 29 June 2013 AMS subject classification: 62H12 Keywords: Dynamic structural equation model Bayesian semiparametric modeling Blocked Gibbs sampler Lν -measure

abstract Analyses of non-normal data and longitudinal data to study changes in variables measured repeatedly over time have received considerable attention in social and psychological research. This paper proposes a dynamic two-level nonlinear structural equation model with covariates for analyzing multivariate longitudinal responses that are mixed continuous and ordered categorical variables. To cope with the non-normal continuous data, the corresponding residual errors at both first-level and second-level models are modeled through a Bayesian semiparametric modeling on the basis of a truncated and centered Dirichlet process with stick-breaking priors. The first-level model is defined for measures taken at each time point nested within individuals for investigating their characteristics that vary with time; while the second level is defined for individuals to assess their characteristics that are invariant with time. An algorithm based on the blocked Gibbs sampler is implemented for estimation of parameters. An efficient model comparison statistic, namely the Lν -measure, is also introduced. Results of a simulation study indicate that the performance of the Bayesian semiparametric estimation is satisfactory. The proposed methodologies are applied to a real longitudinal study concerning cocaine use. © 2013 Elsevier Inc. All rights reserved.

1. Introduction In substantive research, it is common to encounter latent variables that should be assessed through several observed variables. In many substantive research that involves many observed and latent variables, the main objective is in assessing inter-relationships among observed and latent variables. Due to the presence of latent variables, the classical regression models are not appropriate to achieve the goal. Structural equation models (SEMs) are a flexible class of confirmatory models for analyzing inter-relationships among observed and latent variables, and covariates [1,14]. In general, SEMs involve a measurement equation which is a confirmatory factor analysis model for relating the appropriate observed variables with the latent variables, and a regression type structural equation which regresses the outcome latent variables on explanatory latent variables and covariates to examine their inter-relationships. SEMs have been extensively applied to behavioral, educational, biomedical, social and psychological research. The development of longitudinal models with observed variables measured repeatedly over time has drawn considerable attention in the statistic literature, such as the linear mixed models [22] and the generalized linear mixed models [4], among others. For studying the dynamic changes of the relationships among latent and observed variables over time, Song et al. [20] recently developed a two-level SEM on the basis of the Bayesian approach. Like many SEMs, their model was developed under the critical assumption that the distribution of residual errors is normal. As this assumption may

∗

Corresponding author. E-mail address: [email protected] (X.-Y. Song).

0047-259X/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jmva.2013.06.001

88

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

not be true in substantive research, development of robust methods for analyzing non-normal data has received a great deal of attention in the field of SEMs. For instance, see [11,16] for robust methods that were developed on the basis of the tdistribution or the t-distribution with random weights for analyzing some single-level SEMs. However, statistical inferences using the t-distribution are parametric. In many situations, containing inference to a specific form may limit the scope of application. For example, as the t-distribution is symmetric, the above mentioned robust methods are not effective to cope with skewed or mixture distributions. In general, the limitation of the parametric methods motivates the development of Bayesian semiparametric modeling in statistics. Inspired by the wide applications of the Bayesian semiparametric methods to statistical models [13,10,15,21], it is natural to extend the two-level dynamic parametric SEM to a semiparametric SEM for solving the non-normality problem of the error residuals. More specifically, in this article, we will develop a Bayesian semiparametric two-level dynamic SEM for analyzing the dynamic behaviors of observed and latent variables, and covariates over time. The first level model is formulated for accumulating the characteristics that are variant with time, while the second level model is focused on the characteristics that are invariant over time, although time-invariant fixed covariates are also accommodated in the structural equation; see Eq. (5). The residual errors corresponding to the continuous observed variables in models at both levels are formulated by a semiparametric modeling with the truncated and centered Dirichlet process [10]. A block Gibbs sampler [9] is implemented to obtain the Bayesian estimates and a criterion-based statistic is used for model comparison. We illustrate our developed methodology with a longitudinal dataset about cocaine use and the related phenomena. Cocaine use is a major social and health problem in the United States. Cocaine is consistently detected in the urine specimens from approximately one-third of arrestees tested across the nation each year; and it is the illegal drug mentioned most often in emergency department records. There are various latent traits that influence cocaine use; such as the impact of psychiatric problems on cocaine-dependent patients (see [2]) and the social support (see [6,23]). To illustrate the application of the proposed methodologies, we analyze the UCLA longitudinal dataset collected from cocaine-dependent patients at intake, one year, two years, and 12 years after treatment [7]. As the distributions of some of the continuous observed variables are non-normal, the proposed Bayesian semiparametric approach is useful in analyzing this dataset. The contributions of this article include: (i) establishment of a novel Bayesian semiparametric dynamic two-level SEM for analyzing longitudinal non-normal continuous and ordered categorical data; (ii) development of an algorithm to compute the Bayesian estimates of the parameters in the proposed model, with derivations of conditional distributions under the new model settings; and (iii) introduction of a Bayesian statistic to address the complicated issue about model comparison for the semiparametric dynamic two-level SEMs. The rest of the article is organized as follows. Section 2 presents the proposed Bayesian semiparametric two-level SEM for analyzing non-normal data with mixed continuous and ordered categorical variables measured at multiple time points. The Bayesian semiparametric approach for estimation and model comparison is discussed in Section 3, together with some details on the block Gibbs sampler and the model comparison statistic. A simulation study and a real example of the longitudinal study of cocaine use are presented in Sections 4 and 5, respectively. A discussion is presented in Section 6. Some technical details are given in the Appendix. 2. A Bayesian semiparametric two-level dynamic SEM Although the proposed two-level dynamic SEM for analyzing longitudinal data is similar to that given in [20], the distributions of the residual errors are very different. Instead of assuming normal distributions for the residual errors as in their model, we propose a semiparametric modeling for coping with non-normal data. As a result, the development of statistical methods in this paper is very different from that in [20]. To introduce notation, we first present the basic equations before the semiparametric formulation of the model. 2.1. Basic measurement and structural equations of the model We consider a set of observation of a p × 1 random vector (z′gt , x′gt )′ for each individual g, observed at multiple time points t, for g = 1, . . . , G and t = 1, . . . , T , where xgt and zgt = {zgtj ; j = 1, . . . , s} are observation vectors of continuous and ordered categorical variables, respectively. The dataset {(z′gt , x′gt )′ : g = 1, . . . , G, t = 1, . . . , T } can be regarded as an observed sample of hierarchical observations that were measured at different time points and were nested in the individuals. Let ugt = (w′gt , x′gt )′ , where wgt is an s × 1 vector of unobserved continuous variables whose information is given by zgt . The relationship between an ordered categorical variable zgtj and its underlying continuous variable wgtj is defined by zgtj = k + 1 if αk,gtj ≤ wgtj < αk+1,gtj , for k = 0, . . . , mtj − 1, where {−∞ = α0 < α1 < · · · < αm = ∞} is the set of threshold parameters that define the m categories. The following two-level SEM is proposed to model ugt : ugt = yg + vgt

(1)

where yg and vgt are the random vectors that, respectively, model the second-level effect of the g-th individual and the first-level effect with respect to the g-th individual at time t. The second-level model for examining characteristics of individuals that are invariant over time is defined by: For g = 1, . . . , G, yg = A0 cg0 + Λ0 ωg0 + ϵg0 ,

(2)

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

89

where A0 and Λ0 are matrices of unknown coefficients, cg0 is a vector of fixed covariates, ωg0 is a q0 -dimensional vector of latent variables, and ϵg0 = {ϵdg0 , ϵcg0 } is a vector of residual errors, in which ϵdg0 and ϵcg0 correspond to wgt and xgt , respectively. We assume that ωg0 is identically and independently distributed (i.i.d.) as N (0, Φ0 ), and it is independent of ϵg0 . The purpose of this second-level model is for studying the relationships between ugt and ωg0 with respect to different individuals but invariant with time. The measurement model for the more important first-level random vector vgt is defined as follows: For g = 1, . . . , G, vgt = At cgt + Λt ωgt + ϵgt ,

t = 1, . . . , T

(3)

in which the definitions of At , cgt , Λt , ωgt , and ϵgt are similar to those given in Eq. (2); except they are defined at time point t nested within the individual g. Here, the first element of cgt is fixed at 1 at all the time points, so that the first column of ′ At is a vector of intercepts. Moreover, let (η′gt , ξ gt )′ be a partition of ωgt , where ηgt (q1 × 1) and ξ gt (q2 × 1) are outcome and explanatory latent vectors, respectively, and q = q1 + q2 . For g = 1, . . . , G, it follows from (2) and (3) that the random vector ugt conditional on yg has the following structure: ugt = A0 cg0 + Λ0 ωg0 + ϵg0 + At cgt + Λt ωgt + ϵgt ,

t = 1, . . . , T .

(4)

Eq. (4) accounts for dependency among the observed variables measured for the individual g at a given time point t via cgt and ωgt . The relationships of vgt and the corresponding cgt and ωgt that change dynamically over time can be assessed by estimating At and Λt . The following autoregressive nonlinear structural equation is defined for assessing the effects of ξ gt on ηgt : For g = 1, . . . , G, t = 1, . . . , T ,

ηgt = B0 dg0 + Bt dgt + Γ ∗t F∗t (ωg1 , . . . , ωg ,t −1 , ξ gt ) + δgt ,

(5)

where B0 , Bt , and Γ t are matrices of unknown coefficients, dg0 (which is different from the time-invariant covariate cg0 in the second level model) and dgt are vectors of fixed covariates that are invariant and variant over t, respectively; δgt is a vector of residual errors, and F∗t = (ft1 , . . . , ftr ) is a vector-valued function in which ftj is a differentiable function of the independent latent vector ξ gt at the current time point and the latent vectors ωg1 , . . . , ωg ,t −1 at previous time points. The vector of residual errors δgt is independent of ϵgt , ωg1 , . . . , ωg ,t −1 , and ξ gt . As the structural equation involves latent variables rather than observed variables that produce the data, we assume that δgt are independently and identically distributed ′ ′ (i.i.d.) as N (0, Ψ δ t ), where Ψ δ t is a diagonal matrix. Moreover, let ξ g = (ξ g1 , . . . , ξ gT )′ ; the distribution of ξ g is assumed to be N (0, Φ), where Φ contains the variance–covariance matrix Φtt of ξ gt , and the covariance matrix Φti ,tj of ξ gti and ξ gtj at different time points ti and tj . Note that the empirical distributions corresponding to the categorical data are allowed to be skewed and/or bimodal. To cope with the non-normal continuous data, the related residual errors will be modeled by a Bayesian semiparametric approach. ∗

2.2. Bayesian semiparametric modeling of the residual errors In classical structural equation modeling, it is assumed that xgt follows a multivariate normal distribution given the latent vectors ωgt and ωg0 . This assumption may not be true in substantive research. Hence, it is desirable to relax the basic normal assumption of the corresponding residual errors in the measurement equation by a Bayesian semiparametric approach. For g = 1, . . . , G and t = 1, . . . , T , let ϵdgt = (ϵgt1 , . . . , ϵgts ), ϵcgt = (ϵgt ,s+1 , . . . , ϵgt ,p ), ϵdg0 = (ϵg01 , . . . , ϵg0s ), and ϵcg0 = (ϵg0,s+1 , . . . , ϵg0,p ), where p is the dimension of yg or vgt . It is assumed that ϵgt ,s+1 , . . . , ϵgt ,p and ϵg0,s+1 , . . . , ϵg0,p are independent. We model these residual errors as follows. For those associated with the ordered categorical variables and t = 0, 1, . . . , T , we assume that each random error ϵdgt follows N (0, Ψ dt ) where Ψ dt is a diagonal matrix. For those associated with continuous variables, inspired by the key idea in Bayesian nonparametric approach [10,9], we model ϵcgt as follows: For i = s + 1, . . . , p,

[ϵgti |ζgti ] ∼ N (ζgti , ψtii ), i.i.d.

[ζgti |Pit ] ∼ Pit ,

t = 0, . . . , T ,

Pit ∼ Pit ,

t = 0, . . . , T ,

(6) (7)

where (·) are random probability measures, and ζgti ’s are random variables introduced to define the prior distribution of ϵgti ’s. Note that t = 0 corresponds to the second-level model, and t > 0 corresponds to the first-level model. Inspired by the approach given in [24], we consider a truncated and centered Dirichlet process (TCDP) with stick-breaking priors as follows: For i = s + 1, . . . , p and t = 0, . . . , T ,

Pit

t

Pit

=

Ni 

πikt δZ∗˜t (·), ik

k=1

˜ =˜

− µP t ,

i.i.d. Zikt ∗

t z˜ i

Zikt

˜

Zikt ∗

∗

i

∼ N (µ , ψz˜ti ),

k = 1, . . . , Nit ,

(8)

90

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

where δ ∗˜t (·) denotes a discrete measure concentrated at Z˜ikt , πikt are the random weights independent of Z˜ikt ∗ , such that Zik

0 ≤ π ≤ 1 and t ik

Nit

k=1

πikt = 1, and

t

µP t = ∗

i

Ni 

πikt Z˜ikt ∗ .

(9)

k =1

The random weights are defined through the following stick-breaking procedure:

πi1t = νi1t ,

and πikt = νikt

k−1 

(1 − νijt ),

k = 2, . . . , Nit − 1,

(10)

j=1 i.i.d

in which νikt ∼ Beta(1, α˜ it ), k = 1, . . . , Nit − 1 and α˜ it ∼ Gamma(at1i , at2i ). Roughly, the construction of the stick-breaking priors can be thought of as a stick-breaking procedure, where at each stage we randomly and independently break what is left of a stick of unit length and assign the length of this break to the current πikt value. Here the random weight νikt is modeled by Beta(1, α˜ it ), and α˜ it is considered as a random hyper-parameter with a prior distribution Gamma(at1i , at2i ). t Moreover, the prior distributions of µtz˜ i and (ψz˜ti )−1 are taken to be N (0, ψµt i ) and Gamma(c1i , c2it ), respectively. The hypert parameters ψµt i , c1i , c2it , at1i and at2i are preassigned given values. The TCDP (8) is a typical truncated Dirichlet process (TDP) except for the centralization of the atoms Z˜ikt , which constrains the random distribution with mean 0. Under this TCDP, it can be shown that E (ϵgti |Pit ) = 0. Hence, the proposed Bayesian semiparametric model is identified by this centralization of the atoms. The proposed model is reformulated through some classification random variables for achieving efficient Markov chain Monte Carlo (MCMC) sampling. For i = s + 1, . . . , p and t = 0, . . . , T , let Lti = (Lti1 , . . . , LtiG )′ , where Ltig are conditionally ∗ ∗ ∗ independent classification variables that identify the Z˜ikt ∗ associated with ζgti , where ζgti = Z˜iLt ∗t and ζgti = ζgti − µ∗P t . Let ig

i

t∗ πti = (πi1t , . . . , πiNt t ) and Z˜ ti ∗ = (Z˜i1t ∗ , . . . , Z˜iN t ). The semiparametric hierarchical model defined in (6)–(10) can be rewritten i

as

i

ϵgti |Z˜ ti ∗ , Lti ∼ N (ζgti , ψtii ),

∗ ζgti = ζgti − µ∗P t , i

∗ ζgti = Z˜iLt ∗t , ig

Nit

Ltig |π ti ∼



πikt δk (·),

(πti , Z˜ ti ∗ ) ∼ gi1t (πti )gi2t (Z˜ ti ∗ ),

(11)

k=1 t where gi1 (πti ) and gi2t (Z˜ ti ∗ ) are given by (10) and (8), respectively. The prior distributions for the unknown parameters in the model are given as follows. For t = 0, 1, . . . , T , let ψtii be the variance of ϵgti for i = 1, . . . , s, and A′tk and Λ′tk be the kth row of At and Λt , respectively, for k = 1, . . . , p. For t = 1, . . . , T , let ψδ tll be the lth diagonal elements of Ψ δ t and Λ′δ tl be the lth row of Λδ t for l = 1, . . . , q1 , where Λδ t = (Bt , Γ ∗t ). For l = 1, . . . , q1 , let B′0l be the lth row of B0 . The following common conjugate prior distributions in Bayesian analysis of SEM are used: For t = 0, 1, . . . , T ,

−1 π (ψtkk ) ∼ Gamma(α˜ ϵt k , β˜ ϵt k ),

π (Λtk |ψtkk ) ∼ N (Λtϵ k , ψtkk Htϵ k ),

π (Atk ) ∼ N (Atϵ k , HtAk ),

(12)

and for t = 1, . . . , T ,

π (ψδ−tll1 ) ∼ Gamma(α˜ δt l , β˜ δt l ), π (B0l ) ∼ N (B0δl , H0Bl ),

π (Λδtl |ψδtll ) ∼ N (Λtδl , ψδtll Htδl ),

1 Φ− 0 ∼ Wishartq0 (R0 , ρ0 ),

Φ−1 ∼ Wishartq2 ×T (Rw , ρw ),

(13)

where α˜ ϵt k , β˜ ϵt k , α˜ δt l , β˜ δt l , Λtϵ k , Atϵ k , Λtδ l , B0δ l , ρ0 , ρw , and the positive definite matrices Htϵ k , HtAk , Htδ l , H0Bl , R0 , Rw are hyperparameters whose values are assumed to be given from the prior information. Moreover, we use the common noninformative prior distribution for the unknown thresholds of the ordered categorical variables [14]. Given that the proposed model is very complicated, it is necessary to pay special attention to model identification. A common practice in SEMs is to fix some parameters to identify the model. Specifically, the following steps can be employed: (i) we fix some elements in loading matrices Λt , t = 0, 1, . . . , T . Usually, each latent variable has its own manifest indicators and the indicators of different latent variables are not overlapped. Then in each column of Λt , the factor loadings unrelated to the latent variable corresponding to this column should be fixed at zero. Moreover, one of the nonzero factor loadings for each latent variable is fixed at a preassigned value, say one, to introduce a scale to this latent variable. An alternative method for determining the scale of the latent variables is to fix the diagonal elements of Φt at 1.0 but free the nonzero factor

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

91

Fig. 1. Plots of the estimated densities of the random errors in the simulation under the setting (I). The dotted, dash-dot, and solid lines stand for the 2.5% quantile, 97.5% quantile, and average of the estimated densities obtained in 100 replications, respectively; and the dashed lines stand for the true density of the random errors.

loadings [14]. Using this non-overlapping structure of loading matrices, we can identify the measurement equations and obtain a simple interpretation of latent variables [17]. (ii) We fix the threshold parameters α1 and αm−1 for each categorical manifest indicator in order to solve the model indeterminacy caused by the ordinal data. Another important issue is about the identification of ϵg0 and ϵgt in Eq. (4). When the number of time points T equals 1, ϵg0 and ϵgt cannot be identified because they are in an additive form. In that case, ϵg0 should be removed. However, for longitudinal data such as the ones considered in Sections 4 and 5, T is usually larger than 1. In those cases, ϵg0 and ϵgt can be identified as long as the sample size is not too small. In fact, if the distributions of either ϵg0 or ϵgt , t = 1, . . . , T is changed, the joint distribution of ϵg0 + ϵg1 , . . . , ϵg0 + ϵgT will be changed because the correlation among ϵg0 +ϵg1 , . . . , ϵg0 +ϵgT is only from ϵg0 while the distribution of each ϵg0 +ϵgt is determined by ϵgt when the distribution of ϵg0 is given. Like those in many hierarchical models, ϵg0 , g = 1, . . . , G are set here to describe the random errors among different individuals (g = 1, . . . , G). Given that ϵg0 is common for all the time points t = 1, . . . , T and stands for a common part of ugt for t = 1, . . . , T which cannot be explained by A0 + cg0 + Λ0 ω0 , its role cannot be played by the random error ϵgt because ϵgt s are set to be independent of each other. The simulation in Section 4 shows that ϵg0 and ϵgt can be identified and their distributions can be estimated (see Figs. 1 and 2).

3. Parameter estimation and model comparison Let α be the parameter vector of the unknown thresholds, and θ be the parameter vector that contains all the unknown distinct structural parameters in {ψtkk , At , Λt ; k = 1, . . . , p, t = 0, 1, . . . , T } and {Ψ δ t , Bt , Γ ∗t , B0 , Φ0 , Φ; t = 1, . . . , T }, which are involved in Eqs. (2), (3) and (5). Recall that zgt = {zgtj ; j = 1, . . . , s} is a vector of ordered categorical observations, and ugt = (w′gt , x′gt )′ , where xgt is a vector of continuous measurements and wgt is the latent continuous measurement corresponding to zgt . Let ug = {ugt ; t = 1, . . . , T }, Z = {zgt ; g = 1, . . . , G, t = 1, . . . , T }, W = {wgt ; g = 1, . . . , G, t = ′ ′ 1, . . . , T }, and X = {xgt ; g = 1, . . . , G, t = 1, . . . , T }. Moreover, let ωg = (η′g1 , . . . , η′gT , ξ g1 , . . . , ξ gT )′ , Ω1 = {ωg ; g = 1, . . . , G}, Ω0 = {ωg0 ; g = 1, . . . , G}, and Y = {yg ; g = 1, . . . , G} be matrices of various kinds of latent vectors at the first and second levels. Finally, for i = s + 1, . . . , p, t = 0, . . . , T , let L = {Lti }, Z˜ ∗ = {Z˜ ti ∗ }, π = {π ti }, µz˜ = {µtz˜ i }, Ψ z˜ = {ψz˜ti }, α˜ = {α˜ it }, and Ξ = {ζgti ; g = 1, . . . , G}. In these sets, the indices of i and t are ranged from s + 1 to p, and 0 to T , respectively.

92

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

Fig. 2. Plots of the estimated densities of the random errors in the simulation under the setting (II). The dotted, dash-dot, and solid lines stand for the 2.5% quantile, 97.5% quantile, and average of the estimated densities obtained in 100 replications, respectively; and the dashed lines stand for the true density of the random errors.

3.1. Bayesian estimation via the blocked Gibbs sampler Based on the proposed Bayesian semiparametric model and the prior distributions proposed above, a blocked Gibbs sampler [9] is implemented with the following updates in each cycle of the algorithm:

[Ξ |Y, Ω0 , Ω1 , W, θ, α, X], [Y, Ω0 , Ω1 |Ξ , θ, α, X, W],

(14) and

(15)

[θ, α, W|Y, Ω0 , Ω1 , Ξ , X, Z].

(16)

To draw observations from (14), we treat TDP as a parameter-expanded version of TCDP. That is, we draw samples from ∗ ˜ Y, Ω0 , Ω1 , W, θ, α, X], and get ζgti [π, Z˜ ∗ , L, µz˜ , Ψ z˜ , α| = Z˜iLt ∗t , g = 1, . . . , G. Then we obtain the observations of Ξ from ig

N t ∗ − µ∗P t , where µ∗P t = k=i 1 πikt Z˜ikt ∗ . [Ξ |Y, Ω0 , Ω1 , W, θ, α, X] by the centralizing transformation ζgti = ζgti i i Conditional on the given Y, Ω0 , Ω1 , W, θ , α and X, ϵgti , i = s + 1, . . . , p, t = 0, 1, . . . , T , are given. Hence, we have ˜ Y, Ω0 , Ω1 , W, θ, α, X) = π (π, Z˜ ∗ , L, µz˜ , Ψ z˜ , α|

p T  

π (πti , Z˜ ti ∗ , Lti , µtz˜ i , ψz˜ti , α˜ it |ψtii , ϵgti , g = 1, . . . , G).

(17)

t =0 i=s+1

To draw samples from the posterior distribution

π (πti , Z˜ ti ∗ , Lti , µtz˜ i , ψz˜ti , α˜ it |ψtii , ϵgti , g = 1, . . . , G),

(18)

for i = s + 1, . . . , p and t = 0, . . . , T , we complete the following five updates in each cycle of the sampler through Eqs. (2) and (4):

[πti , Z˜ ti ∗ |Lti , µtz˜ i , ψz˜ti , α˜ it , ψtii , ϵgti , g = 1, . . . , G],

(19)

[ |π , ˜ , ψtii , ϵgti , g = 1, . . . , G],

(20)

Lti

t i

Zti ∗

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

93

[µtz˜ i |Z˜ ti ∗ , ψz˜ti ],

(21)

[ψz˜ti |Z˜ ti ∗ , µtz˜ i ],

(22)

[α˜ |π ].

(23)

t i

t i

Conditional distributions in relation to (19)–(23) are presented as follows. Derivations are presented in the Appendix. (i) Sampling from [π ti , Z˜ ti ∗ |Lti , µtz˜ i , ψz˜ti , α˜ it , ψtii , ϵgti , g = 1, . . . , G]: Note that

π (πti , Z˜ ti ∗ |Lti , µtz˜ i , ψz˜ti , α˜ it , ψtii , ϵgti , g = 1, . . . , G) = π (πti |Lti , α˜ it )π (Z˜ ti ∗ |Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G).

(24)

It can be shown that the conditional distribution [π ti |Lti , α˜ it ] is a generalized Dirichlet distribution Dir (ati1∗ , bti1∗ , . . . , atiN∗t −1 , btiN∗t −1 ), where atik∗ = 1 + mtik , btik∗ = α˜ it + i

i

Nit

j=k+1

mtij , k = 1, . . . , Nit − 1 and mtik is the number of Ltig whose

value equals to k. Hence, observations can be sampled as follows:

πi1t = νi1t ∗ ,

and πikt = νikt ∗

k−1 

(1 − νijt ∗ ),

k = 2, . . . , Nit − 1,

(25)

j =1

where νikt ∗ ∼ Beta(atik∗ , btik∗ ). t ∗[L] correspond to those values in Z˜ ti ∗ Let {Lti1∗ , . . . , Ltim∗ t } be the unique set of Ltig values, Z˜ ti ∗L = (Z˜iLt ∗t ∗ , . . . , Z˜iLt ∗t ∗ ), and Z˜ i i

which exclude ˜

Zti ∗L ,

imt i

i1

we have

π (Z˜ ti ∗ |Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G) = π (Z˜ ti ∗[L] |µtz˜ i , ψz˜ti )π (Z˜ ti ∗L |Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G).

(26)

Given Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G, it can be shown that Z˜ikt ∗ , k = 1, . . . , Nit are conditionally independent, where the t ∗[L]

components of Z˜ i

follows N (µtz˜ i , ψz˜ti ), and

[Z˜iLt ∗t ∗ |Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G] ∼ N (µtz˜∗ij , ψz˜tij∗ ), ij

in which µtz˜∗ij = ψz˜tij∗ (µtz˜ i /ψz˜ti + in

{k:Ltik =Ltij∗ } ϵkti /ψtii ),



j = 1, . . . , mti ,

(27)

ψz˜tij∗ = (ntij /ψtii + 1/ψz˜ti )−1 , and ntij is the number of times Ltij∗ occurs

Lti .

(ii) Sampling from [Lti |π ti , Z˜ ti ∗ , ψtii , ϵgti , g = 1, . . . , G]: It can be shown that Ltig , g = 1, . . . , G are conditionally independent integers, and t

[

Ltig

|π , ˜ , ψtii , ϵgti , g = 1, . . . , G] ∼ t i

Zti ∗

Ni 

t∗ πigk δk (·),

g = 1, . . . , G,

(28)

k=1

where t∗ t∗ t t ˜ t∗ 2 ˜ t∗ 2 (πig1 , . . . , πigN t ) ∝ (πi1 exp{−(ϵgti − Zi1 ) /(2ψtii )}, . . . , π t exp{−(ϵgti − Z t ) /(2ψtii )}). iN iN i

i

(29)

i

(iii) Sampling from [µtz˜ i |Z˜ ti ∗ , ψz˜ti ]: It can be shown that [µtz˜ i |Z˜ ti ∗ , ψz˜ti ] ∼ N (µtz˜∗i , ψµt ∗i ), where µtz˜∗i = ψµt ∗i ψµt ∗i = 1/(Nit /ψz˜ti + 1/ψµt i ). t (iv) Sampling from [ψz˜ti |Z˜ ti ∗ , µtz˜ i ]: We have [(ψz˜ti )−1 |Z˜ ti ∗ , µtz˜ i ] ∼ Gamma(c1i + Nit /2, c2it +

Nit

k=1

Nit

k =1

Z˜ikt ∗ /ψz˜ti and

(Z˜ikt ∗ − µtz˜ i )2 /2).

 N t −1

t∗ i (v) Sampling from [α˜ it |π ti ]: It can be shown that [α˜ it |π ti ] ∼ Gamma(Nit + at1i − 1, at2i − k= 1 log(1 − νik )). A hybrid algorithm that combines the Gibbs sampler and the Metropolis–Hastings (MH) algorithm can be used to sample from the conditional distribution (15); see [20]. The posterior sample from (16) can be drawn using Gibbs sampler, and the ∗ expressions for the related conditional distributions are provided in the Appendix. Note that, as ζgti = ζgti − µ∗P t , and the i

i-th intercept ati1 at the time point t in the first-level model is additive to ζgti , it should be correspondingly transformed to ati1 = a∗ti1 + µ∗P t + µ∗P 0 , i

i

for i = s + 1, . . . , p and t = 1, . . . , T , where a∗ti1 is the posterior sample of the intercept drawn from (16).

(30)

94

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

3.2. Model comparison Model comparison is an important issue in SEMs. Due to the complexity of the proposed semiparametric model, the commonly used Bayesian model comparison statistics, such as the Bayes factor or the marginal likelihood, are difficult to obtain. For example, to obtain Bayes factor, we need to calculate the ratio of marginal density of observations under two competing models. However, under the proposed semiparametric model, it is difficult to calculate the joint density function of observations and the unknown parameters because the distributions of ϵg0 and ϵgt are described by a semiparametric method and their distribution forms are not specified. Hence, we need to consider a model comparison approach that does not depend on the distribution forms of ϵg0 and ϵgt . Inspired by the predictive approach for model comparison [3], we propose the Lν -measure as the Bayesian model comparison statistic. Lν -measure is a measurement of difference between the observed response and the conditional mean of a future response in an imaginary replicated experiment under the candidate model and the conditional variance of the future response, given the observed data. The estimate of Lν -measure can be obtained via the samples from the posterior distributions of ϵg0 , ϵgt , and the unknown parameters, hence it avoids dependence on the distribution forms of ϵg0 and ϵgt and its computational burden is light. However, the Lν -measure cannot be directly applied to ordered categorical response data. Similar to [3], we first dichotomize the ordered categorical responses as follows: For t = 1, . . . , T , j = 1, . . . , s, g = 1, . . . , G,

 zgtj,k =

1, 0,

if k = zgtj , if k ̸= zgtj . ∗rep′

∗rep′

rep

Let z∗gtj = (zgtj1 , . . . , zgtj,mtj )′ , where mtj is the number of categories for zgtj ; and let (zgt1 , . . . , zgts )′ and xgt be the future dichotomized vector of an ordered categorical response zgt and the future vector of continuous response xgt of an imaginary replicated experiment, respectively. For some 0 ≤ ν < 1, we consider the following multivariate version of Lν -measure: Lν (X, Z) =

G  T 

rep rep ′ [tr {Cov(xrep gt |X, Z)} + ν(E (xgt |X, Z) − xgt ) (E (xgt |X, Z) − xgt )]

g =1 t =1

+

G  T  s 

[tr {Cov(z∗gtjrep |X, Z)} + ν(E (z∗gtjrep |X, Z) − z∗gtj )′ (E (z∗gtjrep |X, Z) − z∗gtj )],

(31)

g =1 t =1 j =1 rep

∗rep

where the expectation is taken with respect to the posterior predictive distribution of p(xgt , zgtj |X, Z). It can be seen from (31) that both the first and the second terms of Lν (X, Z) are a sum of two components. The first component is related to the variability of the prediction, whereas the second component is related to the discrepancy between the prediction and the observed data. Hence, small values of the Lν -measure indicate predictions close to the observed values with low variability. Naturally, the model with the smallest Lν -measure is selected from a collection of competing models. Model comparison based on Lν -measure does not depend on the assumption that either model is ‘‘true’’. Like some other methods (e.g. Bayes factor), Lν -measure is used just to find out the best one among all the competing models for the dataset under study. In the definition of Lν -measure for a competing model, the future observations are from this competing model, not from an assumed ‘‘true’’ model, and the conditional mean of these future observations given the observed data is viewed as the prediction of the observed data based on the competing model. Given that the assumption of a ‘‘true’’ model is not made here, the use of Lν -measure may not be faced with the problem of model mis-specification. rep ∗rep Let Xrep = {xgt ; t = 1, . . . , T , g = 1, . . . , G} and Z∗rep = {zgtj ; j = 1, . . . , s, t = 1, . . . , T , g = 1, . . . , G}. A Monte Carlo sample {(Xrep(l) , Z∗rep(l) ); l = 1, . . . , M } can be simulated from the distribution p(Xrep(l) , Z∗rep(l) |π (l) , Z˜ ∗(l) , L(l) , θ (l) , α(l) ) based on the MCMC sample {(π (l) , Z˜ ∗(l) , L(l) , θ (l) , α(l) ); l = 1, . . . , M } available in the estimation. A consistent estimate of Lν (X, Z) can be obtained via ergodic average of the simulated sample {(Xrep(l) , Z∗rep(l) ); l = 1, . . . , M }. As {(π (l) , Z˜ ∗(l) , L(l) , θ (l) , α(l) ); l = 1, . . . , M } are available in the estimation and simulating observations from the distribution p(Xrep(l) , Z∗rep(l) |π (l) , Z˜ ∗(l) , L(l) , θ (l) , α(l) ) is straightforward, the computational burden for obtaining Lν (X, Z) is light. It has been shown that Lν -measure with ν = 0.5 has attractive theoretical properties [8]; thus, this value of ν will be used in our numerical studies. 4. A simulation study The main objective of this simulation study is to compare the performances of the proposed Bayesian semiparametric estimation with the performances of the Bayesian parametric approach. A dynamic two-level Bayesian semiparametric SEM defined in Section 2 with ugt consisting of nine observed variables at three time points is considered; see Eqs. (2) and (4). The 5th, 6th, and 9th entries in every ugt were transformed to ordered categorical using thresholds (−1.2∗ , −0.6, 0.6, 1.2∗ ). For i = 1, 2, 3, 4, 7, 8, t = 1–3, and both the first and second level models, the following two types of non-normal distributions are considered for residual errors ϵgti that correspond to the continuous observed variables: (I) ϵgti are drawn ∗ ∗ from the mixture of normal distributions, 0.6N (1.0, 0.36) + 0.4N (−1.5, 0.36); and (II) ϵgti = ϵgti − 2, where ϵgti are

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

95

drawn from Gamma(2, 1). The distribution of ϵgti , under (I) and (II) are bimodal and skewed, respectively. As a result, for i = 1, 2, 3, 4, 7, 8, the i-th component of ugt is non-normal. For the second-level model associated with individuals, cg0 is a scalar covariate that is simulated from a standard normal distribution. The true population values of the parameters are given as follows: A′0 = 0.5



1.0 0.3 0.3

 Φ0 =

0.6 0.0∗ 0.0∗

 Λ0 = ′

0.5

0.5

0.3 1.0 0.3

0.5

0.5

0.5

0.5

0.5

0.5 ,



0.3 0.3 1.0



1.0∗ 0.0∗ 0.0∗

0.6 0.0∗ 0.0∗

0.0∗ 0.6 0.0∗

0.0∗ 1.0∗ 0.0∗

0.0∗ 0.6 0.0∗

0.0∗ 0.0∗ 0.6

0.0∗ 0.0∗ 1.0∗

0.0∗ 0.0∗ , 0.6



and ψ0,55 = ψ0,66 = ψ0,99 = 0.2. In this article, parameters with an asterisk are fixed at the preassigned values to identify the model. For the measurement equation of the first-level model (see (4)), cgt is a scalar covariate with value equal to 1, for all g and t. The true values of parameters in At and Λt are given by: A′t = 1.0



 Λt = ′

1.0∗ 0.0∗ 0.0∗

1.0

1.0

0.8 0.0∗ 0.0∗

1.0

0.8 0.0∗ 0.0∗

1.0 0.0∗ 1.0∗ 0.0∗

1.0 0.0∗ 0.8 0.0∗

1.0 0.0∗ 0.8 0.0∗

1.0 0.0∗ 0.0∗ 1.0∗

1.0



0.0∗ 0.0∗ 0.8

0.0∗ 0.0∗ , 0.8



and ψt55 = ψt66 = ψt99 = 0.36 for t = 1–3. Let ωgt = (ηgt , ξ1gt , ξ2gt )′ , the following structural equation is considered:

ηgt = b10 d1g0 + b20 d2g0 + bt dgt + βt −1 ηg ,t −1 + γ1t ξ1gt + γ2t ξ2gt + γ3t ξ1gt ξ2gt + δgt ,

(32)

where the two components in dg0 = (d1g0 , d2g0 )′ are independently simulated from the standard bivariate normal distribution, and the fixed covariate dgt (t = 1–3) is generated from Bernoulli distribution with probability of success 0.6. The true values of the parameters involved in (32) are: b10 = 1.0, b20 = −1.0; and for all t = 1–3, bt = 0.7, βt −1 = 0.6, γ1t = γ2t = 0.6, γ31 = −0.4, γ32 = −0.6, γ33 = −0.8, Ψ δt = 0.5, (φt ,11 , φt ,12 , φt ,22 ) in Φtt equal to (1.0, 0.3, 1.0), and φtj tl ,ik , j ̸= l, covariances of ξigtj and ξkgtl are all 0.3. The total number of parameters is 121. The sample size was 1000, and 100 replications were considered. In this simulation study, Nit = 30 (see (8)). Prior inputs for the hyper-parameters in the prior distributions related to the t parameters involved in the nonparametric component (see Section 2.2) are: ψµt i = 100, c1i = c2it = 0.01, and at1i = at2i = 2. Prior inputs in the prior distributions corresponding to the parameters involved in the parametric component (see (12) and (13)) are: α˜ ϵt k = α˜ ϵ0k = α˜ δt l = 30, H0ϵ k , H0Ak , Htϵ k , HtAk , Htδ l and H0Bl are all identity matrices; R0 = 13 I3 , Rw = 31 I6 , ρ0 = 10, and

ρw = 7; the elements in A0ϵ k and Atϵ k are 0.45 and 0.8, respectively; and the elements in Λ0ϵ k , Λtϵ k , Λtδl and B0δl are equal to the true values. The hyper-parameters β˜ ϵt k and β˜ δt l are treated as random for obtaining more accurate parameter estimates. We take β˜ ϵ0k = β˜ ϵ0 , k = 1, . . . , p, β˜ ϵt k = β˜ ϵ , k = 1, . . . , p, t = 1, . . . , T , and β˜ δt l = β˜ δ , l = 1, . . . , q1 , t = 1, . . . , T , where the prior distributions of β˜ ϵt k and β˜ δt l are all Gamma(2, 0.5). The initial values of the random hyper-parameters β˜ ϵt k and β˜ δt l are taken to be β˜ ϵt k = β˜ δt l = 9 and β˜ ϵ0k = 5.0.

We collect 10,000 observations after 10,000 ‘burn-ins’ in computing the absolute bias (Bias) and the root mean squares (Rms) of the estimates in 100 replications. Results obtained by the Bayesian semiparametric method are reported in Tables 1 and 2. From these tables, it can be seen that the performance of the Bayesian semiparametric method is satisfactory with rather small ‘Bias’ and ‘Rms’ values. To provide more information about the Bayesian semiparametric method on density estimation of the random errors, we have plotted the 2.5% quantile, 97.5% quantile, and average of the estimated densities associated with the residual errors obtained in 100 replications. To save space, only those corresponding to ϵg01 , ϵg13 , ϵg24 , and ϵg38 under the settings (I) and (II) are presented together with true density functions in Figs. 1 and 2, respectively. The behaviors of plots corresponding to others are similar. These figures show that the density estimates obtained by the Bayesian semiparametric method are close to the true density functions under settings (I) and (II). To save space, results on the nuisance threshold parameters are not presented. For comparison, the same datasets under settings (I) and (II) in the 100 replications were analyzed via the Bayesian parametric approach under the normal assumption. The ‘Bias’ and ‘Rms’ values obtained are presented in Tables 3 and 4. Compared the results given in these tables with those presented in Tables 1 and 2, we see that the ‘Bias’ and ‘Rms’ values obtained via the Bayesian parametric approach are generally larger than those obtained via the Bayesian semiparametric approach. In particular, we observe some very large ‘Bias’ and ‘Rms’ values corresponding to some parameters in the secondlevel model. The simulation results show that the performance of the Bayesian semiparametric approach is better than the Bayesian parametric approach. Moreover, these results indicate that even with a pretty large sample size of 1000, the performances of the parametric approach are not acceptable. We have conducted more simulations with smaller simple sizes (750, 500, 250). We found many divergent cases when n = 250, and as expected, larger sample sizes giving better performance for both approaches. Generally, the performances of the Bayesian semiparametric approach are acceptable; and performances of the Bayesian parametric are not satisfactory. To save space, details are not presented.

96

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

Table 1 Simulation results obtained via the Bayesian semiparametric approach under setting (I).

First level

Par

Bias

Rms

Par

λ1,21 λ1,31 λ1,52 λ1,62 λ1,83 λ1,93 λ2,21 λ2,31 λ2,52 λ2,62 λ2,83 λ2,93 λ3,21 λ3,31 λ3,52

0.008 −0.000 0.038 0.049 0.013 0.034 0.001 0.001 0.022 0.026 0.019 0.046 0.003 −0.002 0.023

0.039 0.037 0.081 0.089 0.075 0.071 0.026 0.022 0.069 0.065 0.072 0.091 0.020 0.020 0.071

λ3,62 λ3,83 λ3,93

b10 b20 b1 b2 b3

−0.010

0.040 0.040 0.116 0.118 0.102 0.042 0.029 0.095 0.100 0.123 0.073 0.093

a0,11 a0,12 a0,13 a0,14 a0,15 a0,16 a0,17

−0.007 −0.015 −0.008 −0.001 −0.008

0.055 0.062 0.055 0.047 0.040 0.032 0.059

β1 β2 γ11 γ21 γ31 γ12 γ22

Second level

0.009 −0.019 −0.041 −0.014 0.008 0.009 0.012 0.016 −0.045 0.008 0.038

0.001

−0.008

Rms

Par

Bias

Rms

0.030 0.020 0.029 0.026 0.018 0.011 0.005 0.005 0.017 0.008 0.005 0.012 0.044 0.036 0.028

0.072 0.067 0.077 0.129 0.115 0.110 0.064 0.050 0.047 0.066 0.078 0.053 0.146 0.134 0.128

a2,14 a2,15 a2,16 a2,17 a2,18 a2,19 a3,11 a3,12 a3,13 a3,14 a3,15 a3,16 a3,17 a3,18 a3,19

0.011 0.009 0.009 0.006 0.002 0.010 0.025 0.017 0.019 0.011 0.003 0.005 −0.005 −0.004 0.008

0.075 0.053 0.050 0.064 0.071 0.051 0.152 0.131 0.122 0.069 0.055 0.052 0.069 0.076 0.047

γ32 γ13 γ23 γ33 φ11 φ12 φ13 φ14 φ15 φ16 φ22 φ23

−0.051

0.117 0.089 0.095 0.116 0.137 0.074 0.082 0.078 0.081 0.073 0.150 0.067

φ24 φ25 φ26 φ33 φ34 φ35 φ36 φ44 φ45 φ46 φ55 φ56 φ66

−0.038 −0.020 −0.041 −0.052 −0.021 −0.030 −0.013 −0.089 −0.019 −0.046 −0.049 −0.010 −0.065

0.104 0.067 0.101 0.136 0.075 0.091 0.071 0.170 0.071 0.105 0.142 0.072 0.144

a0,18 a0,19

−0.005 −0.000

0.060 0.035 0.138 0.120 0.089 0.050 0.095

λ0,93 φ0,11 φ0,12 φ0,13 φ0,22 φ0,23 φ0,33

0.017

−0.113 −0.008 −0.013 −0.045 −0.009 −0.007

0.061 0.242 0.100 0.107 0.111 0.069 0.145

a1,11 a1,12 a1,13 a1,14 a1,15 a1,16 a1,17 a1,18 a1,19 a2,11 a2,12 a2,13

λ0,11 λ0,31 λ0,42 λ0,62 λ0,73

Bias

0.019 0.000 −0.047 −0.078 −0.024 −0.040 −0.030 −0.035 −0.025 −0.068 −0.022

0.050 0.067 0.057 0.035 0.049

5. An application of the proposed approach: the longitudinal study of cocaine use We apply the proposed Bayesian semiparametric approach to a dataset obtained from a longitudinal study about cocaine use conducted at the UCLA Center for Advancing Longitudinal Drug Abuse Research. A variety of measures were collected from patients admitted to the West Los Angeles Veterans Affair Medical Center in 1988–1989 and who met the DSM III-R criteria for cocaine dependence [12]. These patients were assessed at baseline, one year after treatment, two years after treatment, and 12 years after treatment (t = t1 , t2 , t3 , t4 ) in 2002–2003. For each time point, the sample size is 227 with some data points containing missing entries. The effect of missing data was taken care by the data augmentation and imputation based on the appropriate conditional distributions. The following observed variables at each t are involved in the current analysis: (i) cocaine use (CC), an ordered categorical variable with codings 1–5 to denote days of cocaine use per month that are fewer than 2 days, between 2 and 7 days, between 8 and 14 days, between 15 and 25 days, and more than 25 days, respectively; (ii) Beck inventory (BI), a continuous variable; (iii) depression (DEP), a continuous variable based on the Hopkins Symptom Checklist-58 scores; (iv) number of friend (NF), an ordered categorical variable with codings 1–5 to denote no friend, 1 friend, 2–4 friends, 5–8 friends, more than 9 friends; (v) ‘have someone to talk to about problem’ (TP); (vi) ‘currently employed’ (EMP); and (vii) ‘alcohol dependence (AD) at baseline’. The last three variables are ordered binary variables with {0,1} for {No, Yes}. In this study, observed variables (CC, BI, DEP, NF, TP) and the fixed covariates EMP measured at four time points were considered for each individual. More specifically, for the gth individual, observed variables (CC, BI, DEP, NF, TP) and EMP were measured at t = t1 , . . . , t4 times points. This gives the observed dataset {ugt ; g = 1, . . . , 226, t = t1 , . . . , t4 } and a set of fixed covariates {dgt (=EMPgt ); g = 1, . . . , 227, t = t1 , . . . , t4 }. Moreover, for the gth individual, the ADg0 at baseline was also considered as an invariant fixed covariate; and this gives the set of fixed covariates {dg0 (=ADg0 ); g = 1, . . . , 227}. These datasets were analyzed by the two-level dynamic SEM defined by Eqs. (1)–(3) and (5), through the proposed Bayesian semiparametric approach. In this two-level SEM, we need to define a model, at the second level related to yg (see Eqs. (1) and (2)) for assessing the invariant characteristics of (CC, BI, DEP, NF, TP) corresponding to the gth individual. Here, the following confirmatory factor analysis model with two factors was used: For every g, yg = Λ0 ωg0 + ϵg0 ,

(33)

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

97

Table 2 Simulation results obtained via the Bayesian semiparametric approach under setting (II).

First level

Par

Bias

Rms

Par

Bias

Rms

Par

Bias

Rms

λ1,21 λ1,31 λ1,52 λ1,62 λ1,83 λ1,93 λ2,21 λ2,31 λ2,52 λ2,62 λ2,83 λ2,93 λ3,21 λ3,31 λ3,52

0.010 0.008 0.021 0.009 0.015 −0.002 0.006 −0.000 0.023 0.025 0.029 0.019 0.001 0.003 0.016

0.039 0.036 0.086 0.080 0.093 0.071 0.030 0.025 0.079 0.074 0.094 0.071 0.024 0.021 0.066

λ3,62 λ3,83 λ3,93

0.012 0.025 0.007 0.013 0.018 0.011 −0.006 −0.002 0.005 −0.004 −0.009 −0.001 0.039 0.031 0.035

0.072 0.098 0.069 0.111 0.099 0.099 0.070 0.056 0.052 0.084 0.086 0.059 0.141 0.118 0.116

a2,14 a2,15 a2,16 a2,17 a2,18 a2,19 a3,11 a3,12 a3,13 a3,14 a3,15 a3,16 a3,17 a3,18 a3,19

−0.006 0.004 0.007 0.001 0.000 −0.006 0.035 0.039 0.030 −0.001 0.001 0.006 0.000 −0.003 0.004

0.065 0.068 0.047 0.071 0.076 0.051 0.142 0.132 0.125 0.072 0.061 0.051 0.070 0.075 0.051

b10 b20 b1 b2 b3

−0.004

0.037 0.043 0.101 0.103 0.108 0.110 0.098 0.042 0.100 0.098 0.120 0.089

γ32 γ13 γ23 γ33 φ11 φ12 φ13 φ14 φ15 φ16 φ22 φ23

−0.042

0.028 0.090 0.085 0.100 0.174 0.077 0.090 0.071 0.096 0.073 0.153 0.065

φ24 φ25 φ26 φ33 φ34 φ35 φ36 φ44 φ45 φ46 φ55 φ56 φ66

−0.039 −0.012 −0.021 −0.068 −0.006 −0.040 −0.002 −0.062 −0.011 −0.021 −0.031 −0.006 −0.002

0.101 0.064 0.095 0.155 0.074 0.087 0.068 0.169 0.072 0.088 0.136 0.070 0.145

a0,11 a0,12 a0,13 a0,14 a0,15 a0,16 a0,17

−0.005 −0.005 −0.000 −0.008 −0.005

0.060 0.069 0.057 0.052 0.047 0.041 0.048

a0,18 a0,19

−0.010 −0.000

0.053 0.034 0.137 0.150 0.081 0.057 0.105

λ0,93 φ0,11 φ0,12 φ0,13 φ0,22 φ0,23 φ0,33

0.019

−0.101 −0.004 −0.020 −0.017 −0.024 −0.003

0.062 0.289 0.107 0.124 0.100 0.080 0.147

0.015 −0.014 −0.048 −0.029 0.010 0.010 −0.002 0.029 −0.023 0.027 0.036

β1 β2 γ11 γ21 γ31 γ12 γ22

Second level

0.003

−0.010

a1,11 a1,12 a1,13 a1,14 a1,15 a1,16 a1,17 a1,18 a1,19 a2,11 a2,12 a2,13

λ0,11 λ0,31 λ0,42 λ0,62 λ0,73

0.017 −0.008 −0.012 −0.049 0.003 −0.041 0.003 −0.036 −0.001 −0.039 −0.001

0.036 0.039 0.048 0.043 0.045

where yg corresponds to (CC, BI, DEP, NF, TP) was used to assess the characteristics of the gth individual invariant over time,

ωg0 = (ωg01 , ωg02 )′ , Ψ 0 = diag(ψ0,11 , . . . , ψ0,55 ),  ∗   ∗ 1.0 0∗ 0∗ 0∗ 0∗ 1.0 Λ′0 = ∗ and Φ0 = 0 λ0,22 λ0,32 λ0,42 λ0,52 φ0,21

 1.0

∗

.

As yg1 was fully defined by CC, λ0,11 was fixed at 1.0. Moreover, as yg1 = ωg01 + ϵg01 , φ0,11 was fixed at 1.0 for model identification. To identify the scale of ωg02 , we may fix either φ0,22 or λ0,22 . In this study, we fixed φ0,22 instead of λ0,22 at 1.0 so that φ0,21 was a correlation coefficient between ωg01 and ωg02 . The latent variable ωg01 represents the invariant portion of cocaine use and ωg02 can be interpreted as an invariant general latent factor over time. The correlation between ωg01 and ωg02 is assessed by φ0,21 . To complete the formulation of the proposed two-level dynamic SEM, we also need to define a first-level model related to vgt (see Eqs. (1), (3) and (5)) for assessing various characteristics that are dynamically changed over time. To achieve the objective, we considered a measurement model vgt = At cgt + Λt ωgt + ϵgt ,

t = t1 , t2 , t3 , t4 ,

(34)

where ωgt = (ηgt , ξ1gt , ξ2gt ) , cgt is fixed at 1.0 so that At is a vector of intercepts, and ′

1.0∗ 0∗ 0∗

 Λ′t =

0∗ 1.0∗ 0∗

0∗

λt ,32 0∗

0∗ 0∗ 1.0∗

0∗ 0∗

λt ,53

 .

For model identification, the thresholds α1 and αm−1 of categorical variables at each time point were fixed at some preassigned values. Specifically, we fixed α1 = Φ ∗−1 (f1∗ ) and αm−1 = Φ ∗−1 (fm∗ ), where Φ ∗ (·) is the distribution function of N [0, 1], f1∗ and fm∗ are the frequency of the first category, and the cumulative frequency of the category with zk < m − 1, respectively. To identify the model with the ordered binary variable TP, ψt ,55 were all fixed to 1.0 for t = t1 , t2 , t3 , t4 . This measurement equation can also be regarded as a confirmatory component of the model with a specific non-overlapping loading matrix Λt for achieving clear interpretation of the latent variables. Hence, latent variables ηgt , ξ1gt , and ξ2gt can

98

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

Table 3 Simulation results obtained via the Bayesian parametric approach under setting (I).

First level

Par

Bias

Rms

Par

Bias

Rms

Par

Bias

Rms

λ1,21 λ1,31 λ1,52 λ1,62 λ1,83 λ1,93 λ2,21 λ2,31 λ2,52 λ2,62 λ2,83 λ2,93 λ3,21 λ3,31 λ3,52

0.002 −0.003 −0.021 −0.025 −0.003 −0.095 −0.002 −0.001 −0.032 −0.037 0.011 −0.070 0.002 −0.002 −0.033

0.042 0.044 0.100 0.094 0.118 0.117 0.028 0.025 0.079 0.082 0.108 0.111 0.023 0.023 0.082

λ3,62 λ3,83 λ3,93

−0.031

a1,11 a1,12 a1,13 a1,14 a1,15 a1,16 a1,17 a1,18 a1,19 a2,11 a2,12 a2,13

0.005 −0.077 0.012 0.008 0.000 −0.003 0.032 0.042 −0.001 −0.002 0.028 0.028 0.025 0.016

0.077 0.098 0.101 0.114 0.103 0.097 0.064 0.060 0.062 0.066 0.077 0.059 0.121 0.108 0.107

a2,14 a2,15 a2,16 a2,17 a2,18 a2,19 a3,11 a3,12 a3,13 a3,14 a3,15 a3,16 a3,17 a3,18 a3,19

0.003 0.037 0.036 −0.002 −0.005 0.028 0.031 0.022 0.024 0.004 0.031 0.032 −0.013 −0.010 0.027

0.074 0.066 0.064 0.065 0.068 0.058 0.136 0.121 0.110 0.068 0.066 0.064 0.070 0.075 0.054

b10 b20 b1 b2 b3

−0.005

0.046 0.043 0.115 0.112 0.106 0.051 0.031 0.140 0.129 0.101 0.091 0.110

γ32 γ13 γ23 γ33 φ11 φ12 φ13 φ14 φ15 φ16 φ22 φ23

0.042 0.007 −0.026 0.065 0.034 0.056 −0.025 0.051 −0.017 0.056 0.127 0.067

0.117 0.101 0.112 0.134 0.191 0.114 0.098 0.103 0.098 0.108 0.258 0.115

φ24 φ25 φ26 φ33 φ34 φ35 φ36 φ44 φ45 φ46 φ55 φ56 φ66

0.032 0.075 0.030 0.054 0.055 −0.013 0.065 0.076 0.067 0.015 0.060 0.071 0.105

0.134 0.120 0.130 0.172 0.112 0.110 0.117 0.218 0.113 0.127 0.186 0.119 0.201

a0,11 a0,12 a0,13 a0,14 a0,15 a0,16 a0,17

−0.004 −0.014 −0.005

0.062 0.070 0.071 0.059 0.044 0.047 0.065

a0,18 a0,19

−0.007

0.068 0.044 0.448 0.467 0.163 0.105 0.211

λ0,93 φ0,11 φ0,12 φ0,13 φ0,22 φ0,23 φ0,33

−0.166

0.193 0.602 0.130 0.151 0.142 0.091 0.656

β1 β2 γ11 γ21 γ31 γ12 γ22

Second level

0.003 0.002 −0.022 −0.012 0.002 0.003 0.035 0.025 0.019 0.009 0.007

0.001 0.016 0.032 −0.009

λ0,11 λ0,31 λ0,42 λ0,62 λ0,73

0.023 −0.154 0.023 −0.135 0.099 −0.098

0.178 −0.044 −0.038 −0.087 −0.002 0.546

be interpreted as ‘cocaine use (CC)’, ‘psychiatric problems’, and ‘social support’, respectively. Finally, we need to choose an autoregressive structural equation for defining Eq. (5). To illustrate the application of L0.5 -measure for model comparison, we considered the competing models with the same measurement equation defined above but with different structural equations defined as follows: M1 : ηg1 = b0 dg0 + b1 dg1 + γ11 ξ1g1 + γ21 ξ2g1 + γ31 ξ1g1 ξ2g1 + δg1 ,

t = t1

ηgt = b0 dg0 + bt dgt + βt −1 ηg ,t −1 + γ1t ξ1gt + γ2t ξ2gt + γ3t ξ1gt ξ2gt + δg1 , t = t2 , t3 , t4 , M2 : ηg1 = b0 dg0 + b1 dg1 + γ11 ξ1g1 + γ21 ξ2g1 + δg1 , t = t1 ηgt = b0 dg0 + bt dgt + βt −1 ηg ,t −1 + γ1t ξ1gt + γ2t ξ2gt + δg1 , t = t2 , t3 , t4 , M3 : ηg1 = γ11 ξ1g1 + γ21 ξ2g1 + γ31 ξ1g1 ξ2g1 + δg1 , t = t1 ηgt = βt −1 ηg ,t −1 + γ1t ξ1gt + γ2t ξ2gt + γ3t ξ1gt ξ2gt + δg1 , t = t2 , t3 , t4 , where dg0 (AD) and dgt (EMP) are fixed covariates that are invariant and variant over time, respectively. In these equations, cocaine use is treated as the outcome latent variables, and it is regressed on the explanatory latent variables related to psychiatric problems and social support, together with their interactions and fixed covariates. In M1 , at t = t1 , cocaine use is first regressed on fixed covariates related to alcohol dependence (dg0 = ADg0 ) and current employment (dg1 = EMPg1 ); and explanatory latent variables related to psychiatric problems (ξ1g1 ) and social support (ξ2g1 ) together with their interaction (ξ1g1 ξ2g1 ). At t = t2 , t3 , t4 , cocaine use at time t −1 (ηg ,t −1 ) is also included in the structural equation. Note that no interaction terms of psychiatric problems and social support are included in M2 , while no fixed covariate about alcohol dependent and current employment are included in M3 . We observed from their histograms that the distributions of BI and DEP are not normal. Hence, we analyzed the dataset through the Bayesian semiparametric approach by modeling the residual errors corresponding to BI and DEP with TCDP. In the TCDP modeling, Nit = 30 is taken (see (8)). Prior inputs for the hyper-parameters in the prior distributions related to t the parameters in the nonparametric component (see Section 2.2) are: ψµt i = 100, c1i = c2it = 2.5, and at1i = at2i = 2.5. Prior inputs of the hyper-parameters in the prior distributions corresponding to the parameters involved in the parametric component (see (12) and (13)) are set as follows: α˜ δt l = 10 and β˜ δt l = 5; α˜ ϵt k = α˜ ϵ0k = 10 and β˜ ϵt k = β˜ ϵ0k = 2, k = 1–5; H0ϵ k ,

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

99

Table 4 Simulation results obtained via the Bayesian parametric approach under setting (II).

First level

Par

Bias

Rms

Par

Bias

Rms

Par

Bias

Rms

λ1,21 λ1,31 λ1,52 λ1,62 λ1,83 λ1,93 λ2,21 λ2,31 λ2,52 λ2,62 λ2,83 λ2,93 λ3,21 λ3,31 λ3,52

0.011 0.015 −0.024 −0.046 −0.022 −0.110 0.008 0.006 −0.020 −0.029 −0.003 −0.083 0.003 0.006 −0.027

0.044 0.048 0.104 0.099 0.129 0.134 0.034 0.028 0.093 0.087 0.122 0.109 0.026 0.028 0.081

λ3,62 λ3,83 λ3,93

−0.034

0.087 0.149 0.113 0.108 0.094 0.102 0.073 0.064 0.062 0.083 0.084 0.061 0.130 0.106 0.102

a2,14 a2,15 a2,16 a2,17 a2,18 a2,19 a3,11 a3,12 a3,13 a3,14 a3,15 a3,16 a3,17 a3,18 a3,19

−0.018 0.033 0.034 −0.012 −0.011 0.012 0.034 0.037 0.026 −0.013 0.029 0.032 −0.011 −0.013 0.023

0.067 0.078 0.060 0.070 0.074 0.053 0.138 0.132 0.120 0.073 0.069 0.063 0.070 0.073 0.057

b10 b20 b1 b2 b3

−0.007

0.044 0.047 0.116 0.102 0.095 0.050 0.033 0.139 0.135 0.119 0.112 0.110

γ32 γ13 γ23 γ33 φ11 φ12 φ13 φ14 φ15 φ16 φ22 φ23

0.139 0.103 0.117 0.166 0.228 0.139 0.112 0.130 0.115 0.127 0.313 0.129

φ24 φ25 φ26 φ33 φ34 φ35 φ36 φ44 φ45 φ46 φ55 φ56 φ66

0.065 0.066 0.087 0.044 0.091 −0.027 0.086 0.164 0.072 0.087 0.057 0.074 0.216

0.162 0.115 0.176 0.203 0.143 0.108 0.132 0.313 0.133 0.162 0.195 0.131 0.337

a0,11 a0,12 a0,13 a0,14 a0,15 a0,16 a0,17

−0.005 −0.012 −0.004 −0.013

0.069 0.074 0.062 0.062 0.052 0.054 0.059

a0,18 a0,19

0.064 0.041 0.503 0.424 0.156 0.107 0.264

λ0,93 φ0,11 φ0,12 φ0,13 φ0,22 φ0,23 φ0,33

−0.199

0.227 0.648 0.123 0.169 0.130 0.106 0.823

β1 β2 γ11 γ21 γ31 γ12 γ22

Second level

0.015 0.011 −0.038 −0.028 0.006 0.009 −0.028 −0.015 0.057 −0.012 −0.019

0.018 0.032 −0.006

a1,11 a1,12 a1,13 a1,14 a1,15 a1,16 a1,17 a1,18 a1,19 a2,11 a2,12 a2,13

λ0,11 λ0,31 λ0,42 λ0,62 λ0,73

0.004

−0.081 −0.015 −0.005 −0.014 −0.019 0.025 0.030 −0.016 −0.020 0.016 0.024 0.017 0.019 0.056

−0.001 −0.064 0.095 0.050 0.092 −0.023 0.091 −0.027 0.084 0.188 0.088

−0.015 0.021 −0.001 0.011 0.125 0.101 −0.151

0.190 −0.014 −0.033 −0.073 −0.026 0.683

H0Ak , Htϵ k , HtAk , Htδ l and H0Bl are all identity matrices; R0 = I2 , Rw = I8 , ρ0 = 3, and ρw = 9; atϵ k = 0.02, t = t1 , t2 , t3 , and atϵ k = 0.25, t = t4 ; λ0ϵ k2 = 0.5, k = 2, 3; λ0ϵ k2 = 0.2, k = 4, 5; λtϵ k = 0.5, k = 5, t = t1 ; and λtϵ k = λ1ϵ 3 = 0.9, k = 3, 5, t = t2 , t3 , t4 . Moreover, the hyper-parameters in Λtδ l corresponding to bt , t = t1 , t2 , t3 , t4 , βt −1 , t = t2 , t3 , t4 , and γkt , k = 1–3, t = t1 , t2 , t3 , t4 are −0.2, −0.3, −0.2, −0.6, 0.0, 0.5, 0.0, 0.5, 0.3, −0.1, 0.3, 0.0, −0.1, 0.0, −0.3, −0.1, 0.4, −0.5 and −0.1, respectively; and the hyper-parameter corresponding to b0 in B0δ l is 0.6. The convergence of the blocked Gibbs sampler is monitored by the ‘estimated potential scale reduction (EPSR)’ values obtained from several parallel sequences of all the parameters generated with different starting values. As suggested by Gelman [5], convergence of these sequences has been achieved if the EPSR values are all less than 1.2. Fig. 3 presents the plots of the EPSR values against the iteration numbers for all the parameters. We observe from this figure that the sequences converged in less than 35,000 iterations. Hence, a total of 40,000 observations are collected after 40,000 ‘burn-in’ iterations via blocked Gibbs sampler for computing parameter estimates. The L0.5 -measure corresponding to M1 , M2 , and M3 are equal to 3753, 3785, and 3789, respectively. Hence, M1 is selected. Moreover, we have calculated the L0.5 -measure of Bayesian parametric model with essentially the same formulation as M1 except that the residual errors are assumed to be normally distributed, where the prior inputs in the prior distributions are the same as those given in the Bayesian semiparametric analysis. For this Bayesian parametric model, L0.5 -measure equals 4088, supporting the choice of semiparametric model. Estimated densities of the residual errors in relation to BI and DEP at the time points t = t1 , t2 , t3 , t4 are plotted in Figs. 4 and 5, respectively. Estimates of the important parameters, such as the factor loadings, and the regression coefficients of the fixed covariates and explanatory latent variables in the structural equation are presented in Fig. 6(a) and (b). Estimates of the parameters in Φ are presented in Table 5, while the remaining parameter estimates corresponding to At , Ψ 0 , Ψ t , and Ψ δ are presented in Table 6. To save space, results on the nuisance threshold parameters are not presented. The following findings can be achieved from Fig. 6(b): (i) The estimates of factor loadings for BI, DEP, and TP are quite large, which indicates substantial associations between the observed variables and the general latent factor. (ii) The general latent factor and cocaine use is positively correlated. From Figs. 6(a) and 7, and Tables 5 and 6, we observe the following phenomenon related to the longitudinal patterns of some important parameters: (i) The estimates of λt ,32 are close to 1.0 at all time points (t = t1 , t2 , t3 , t4 ). As λt ,22 is fixed to 1.0, the loadings of BI and DEP to the latent variable ‘psychiatric problems’ (ξ1gt ) are approximately equal and

100

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

Fig. 3. Convergence of the blocked Gibbs sampler: plots of EPSR values versus iteration numbers for all the estimated parameters. Table 5 Estimates of the parameters and their standard deviations in Φ = [Φti tj ]. t = t1

t = t2

t = t3

t = t4

0.247 (0.065)

−0.058 (0.055) 0.414 (0.141)

0.007 (0.038) −0.022 (0.058) 0.246 (0.065)

0.037 (0.054) 0.166 (0.090) −0.099 (0.065) 0.344 (0.121)

−0.029 (0.035) −0.103 (0.068) 0.046 (0.041) −0.076 (0.064) 0.231 (0.066)

0.033 (0.056) 0.223 (0.083) −0.062 (0.053) 0.149 (0.076) −0.140 (0.058) 0.386 (0.106)

−0.079 (0.047) 0.077 (0.062) −0.104 (0.040) 0.068 (0.060) 0.003 (0.045) 0.093 (0.055) 0.426 (0.081)

0.134 (0.054) 0.009 (0.070) 0.025 (0.056) 0.086 (0.067) −0.026 (0.053) 0.033 (0.067) −0.107 (0.055) 0.370 (0.085)

rather stable over time. This means that the relationship between ‘psychiatric problems’ and {BI, DEP} are stable over time. (ii) The estimates of λt ,53 at later time points are larger than those at former ones. This indicates that the association of TP with ‘social support’ increases at the later follow-up points. (iii) The effect of the invariant fixed covariate AD on ‘cocaine use’ (ηt ) is very minor. (iv) The fixed covariate EMP has a negative effect on ‘cocaine use’ at all time points, which indicates that having a job reduces cocaine use. (v) At baseline, the effect of ‘social support’ on cocaine use is very minor (γˆ21 = 0.062). However, at 12 years after intake, the negative effect of ‘social support’ (γˆ24 = −1.138) on cocaine use becomes very strong. This indicates that in the long run, ‘social support’ is very important on reducing cocaine use. (vi) ‘psychiatric problems’ have a substantial positive effect on cocaine use at the baseline (γˆ11 = 0.506). This effect decreases in the later two years, but becomes more important at 12 years after intake (γˆ14 = 0.779). (vii) There is a substantive interactive effect of ‘psychiatric problem’ and ‘social support’ on cocaine use. Based on the signs of γ1t and γ2t , it seems that there is a suppression between ‘social support’ and ‘psychiatric problems’ on cocaine use. (viii) The effect of cocaine use from t = t1 to t = t2 is negative

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

101

Table 6 Estimates of parameters and their standard deviations in At , Ψ 0 , Ψ t , and Ψ δ . Par

Est.

Second-level

Par

ψ0,11

0.147 (0.034)

ψ0,44

0.210 (0.046)

ψ0,55

0.206 (0.064)

t = t1

a1,11

0.263 (0.144) 0.224 (0.074) 0.050 (0.142) 0.011 (0.065) 0.131 (0.151) −0.006 (0.065) −0.494 (0.251) −0.151 (0.067)

a1,13

0.344 (0.077) −0.013 (0.096) −0.007 (0.065) −0.000 (0.090) −0.111 (0.065) 0.028 (0.086) −0.150 (0.067) 0.030 (0.085)

a1,15

0.242 (0.148) 0.226 (0.079) 0.527 (0.249) 0.188 (0.052) 0.931 (0.391) 0.187 (0.051) 1.044 (0.341) 0.256 (0.098)

a1,12 t = t2

a2,11 a2,12

t = t3

a3,11 a3,12

t = t4

a4,11 a4,12

Est.

Par

a1,14 a2,13 a2,14 a3,13 a3,14 a4,13 a4,14

Est.

ψ1,11 a2,15

ψ2,11 a3,15

ψ3,11 a4,15

ψ4,11

Par

Est.

ψ1,44

0.292 (0.091) 0.460 (0.115) 0.389 (0.108) 0.406 (0.113) 0.310 (0.085) 0.370 (0.085) 0.392 (0.091) 0.675 (0.202)

ψδ 1 ψ2,44 ψδ 2 ψ3,44 ψδ 3 ψ4,44 ψδ 4

Fig. 4. Plots of the estimated densities of the residual errors in relation to BI at the time points t = t1 , t2 , t3 , t4 in the cocaine use example.

(βˆ 1 = −0.371). The impact of cocaine use from 1 year after intake to two years after intake becomes positive and strong (βˆ 2 = 0.592), while the effect of cocaine use from two years after intake to 12 years after intake is very minor (βˆ 3 = 0.003). (ix) The variances φt ,11 , φt ,22 , and φt ,12 are basically constant over time; see Table 5. The above findings coincide with prior literature revealing a strong stability of cocaine use over time. The present results show a rather strong positive association of psychiatric problems with cocaine use at all the time points and a negative association of social support with cocaine use at the last two follow-up points. Basically, both psychiatric problems and social support become stronger factors affecting cocaine use when the treatment effects dissipate at the longer follow-up point.

102

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

Fig. 5. Plots of the estimated densities of the residual errors in relation to DEP at the time points t = t1 , t2 , t3 , t4 in the cocaine use example.

To check the sensitivity of the estimation results to the choice of hyper-parameters, we re-conducted the analysis under other settings of hyper-parameters. For instance, we considered a different setting of hyper-parameters: ψµt i = 10, α˜ δt l = 5,

α˜ ϵt k = α˜ ϵ0k = 5, ρw = 5, and the other hyper-parameters are invariant. The Bayesian estimates of the model parameters

under this prior setting are similar to those reported in Fig. 6 and Tables 5 and 6. For comparison, we also used the Bayesian parametric method to analyze the cocaine use dataset with essentially the same formulation as M1 except that the residual errors are assumed to be normally distributed. Bayesian estimates were obtained with the same prior inputs in the prior distributions as those given in the Bayesian semiparametric analysis. We observe some large differences in some parameter estimates obtained from the different approaches. For instance, the estimates aˆ 3,11 (=0.225), aˆ 4,11 (= −0.294), γˆ31 (= −0.231), γˆ21 (=0.164), γˆ13 (=0.088), γˆ33 (=0.071), γˆ14 (=0.640),

βˆ 1 (= −0.240), βˆ 3 (=0.110), bˆ 3 (= −0.347), bˆ 4 (= −0.244), λˆ 0,22 (=0.638), and λˆ 0,52 (= −0.386) obtained by the parametric

approach are very different from those obtained by the semiparametric approach. To save space, more details about the inferior parametric approach are not presented. Finally, it should be noted that in this application the sample size is rather small and the number of parameters is large; hence the interpretation of the results should be treated with great caution. 6. Discussion Structural equation modeling is a powerful method for assessing inter-relationships among observed and latent variables. In recent years, the growth of SEM is very rapid. In particular, a dynamic two-level SEM for assessing mixed continuous and ordered categorical longitudinal data has been developed by Song et al. [20]. This dynamic SEM is much more general than the basic linear SEM in at least the following two aspects: (i) It is a two-level model that is different from even the twolevel SEM in the literature, and (ii) its structural equation is a general autoregressive nonlinear structural equation that can accommodate effects of independent latent variables at various time points. However, the developments of the dynamic SEM and the related Bayesian methodologies were based on the assumption that the distribution of the residual errors is normal. As many data of in social and psychological research are not normal, it is important to relax the normal assumption. Indeed, developing methods for dealing with the normal assumption has received a great deal of attention for the basic SEM, see [11,16,15,24], among others. However, the contributions provided in these papers were achieved on the basis of standard SEM, and hence cannot be applied to the more general dynamic two-level SEM mentioned above.

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

103

Fig. 6. (a) Path diagram and estimates of some parameters in the first-level model, where observed variables are represented by rectangles, latent variables are represented by ellipse. Here, ηt , ξ1t , and ξ2t , respectively, denote ‘cocaine use’ (CC), ‘psychiatric problems’ and ‘social support’ at t = t1 , t2 , t3 , t4 . (b) Path diagram and estimates of some parameters in the second-level model.

104

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

Fig. 7. (a) Plots of estimates of λ32 , λ53 , b, γ1 , γ2 and γ3 at t = t1 , t2 , t3 , and t4 and (b) plots of estimates of φ11 , φ12 , and φ22 at t = t1 , t2 , t3 , and t4 .

The main objective of this article is to develop a Bayesian semiparametric approach to relax the normal assumption related to the residual errors in both the first level and second level models in the dynamic two-level SEM proposed by [20]. In modeling, the residual errors are modeled through a Bayesian semiparametric approach with a truncated and centered Dirichlet process (TCDP) with stick-breaking priors for relaxing the normal assumption. In the Bayesian analysis, a blocked Gibbs sampler is implemented with the MH algorithm for estimation, and a Lν -measure is computed for model comparison. The above mentioned statistical methods (Dirichlet process, Gibbs sampler, Lν -measure) have been appeared in the literature (for example, [15,24,20]); however, the application and implementation of these statistical methods to the more general dynamic two-level SEM are novel. In particular, the new contributions include (i) the dynamic two-level SEM with semiparametric modeling of residual errors (see details in Section 2), (ii) the conditional distributions derived and presented in Section 3.1 and Appendix, and (iii) the application of Lν -measure to the current dynamic two-level SEM with mixed continuous and ordered categorical longitudinal data. This article also includes a novel application of the developed Bayesian semiparametric methodologies to a cocaine use study. Through this application, we compare our results with those obtained in [20] that use a parametric model with the normal assumption. For clear presentation of the method, the Bayesian semiparametric modeling is only applied to the residual errors that are associated with the continuous variables. However, the key ideas in applying the TCDP with stick-breaking prior and the blocked Gibbs sampler can be extended in developing Bayesian semiparametric methods for dealing with residual errors corresponding to categorical variables. Similarly, these key ideas can also be applied to consider a Bayesian semiparametric modeling of the explanatory latent variables. More research is needed to develop efficient statistical methods for handling other complex data that are frequently encountered in medical research and/or longitudinal studies, such as unordered categorical data that are associated with genotype variables, and non-ignorable missing data. A different approach for handling the violation of the normality assumption is to consider transformation models, which transform response variables so that the resulting model can justify the model assumptions. For instance, [19] considered a single level semiparametric transformation SEM as follows: f(yi ) = Aci + Λωi + ϵi ,

(35)

where f(yi ) = (f1 (yi1 ), . . . , fp (yip )) , fj (·) are unspecified smooth transformation functions, A, ci , Λ, ωi , and ϵi are defined similarly as in Eq. (2). The unknown functions fj (·) are assumed to be strictly monotonically increasing. They proposed using a Bayesian splines approach to estimate the nonparametric transformation functions and unknown parameters, and applied the developed methodology to a study of the dropout behavior on a treatment program preventing polydrug use. Although the idea of the semiparametric SEMs with Bayesian P-splines in [18,19] may be applied to the current two-level SEM for coping with non-normal longitudinal data, the underlying development will be complicated and completely different. As suggested by an anonymous referee, we may consider reducing the number of unknown parameters by modeling the time tendencies of the parameters via some formula. This approach is very valuable if the number of time points, T , is large and the parameters at each time point is of a small amount because it makes the model much more parsimonious and thus reduce the variability. However, it may bring a risk from possible mis-specification of the time tendency formula for T

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

105

the parameters. When T is large but the number of parameters at each point is small, this risk can be reduced by model comparison. In this article, we dealt with the case that T is small and the parameters at each time point is of a large amount. In this situation, modeling the time tendencies of the parameters will be faced with a dilemma. Given that we have many parameters at each time point and we need to choose a formula of time tendency for each one, there will be a large number of candidate models so that model selection will afford a very heavy computational burden. On the other hand, if we reduce the number of candidate models without strong support of prior information, the risk of mis-specification will increase a lot. Based on the above discussion and considering that modeling parameters can only reduce a limited number of parameters as T is small, we did not consider modeling time tendencies of parameters in this article. However, it is an attractive method worth of further study for obtaining a parsimonious model when T is large in a longitudinal study. Acknowledgments This research was supported by grants (446609 and 404711) from the Research Grant Council of Hong Kong Special Administration Region, the grant (11261064) from National Natural Science Foundation of China and the grant (2011FZ151) from Applied and Basic Research Plan of Yunnan Province. The authors are thankful to Dr. Hser from Integrated Substance Abuse Programs, University of California, Los Angeles, USA, for providing the data in the real example. Appendix A. Conditional distribution related to the nonparametric component The conditional distributions related to the nonparametric component are derived as follows. (i) Conditional distribution π (π ti |Lti , α˜ it ): We first note that

π (π | , α˜ ) ∝ π ( |π )π (π |α˜ ) = t i

Lti

t i

Lti

t i

t i

t i

 G 

 π ( |π ) π (πti |α˜ it ). Lti

t i

(A.1)

g =1

Moreover, π ti defined in (10) has a generalized Dirichlet distribution Dir (1, α˜ it , . . . , 1, α˜ it ). It follows that its density has the following form:

π (πti |α˜ it ) = (α˜ it )Ni −1 (πiNt t )α˜ i −1 (1 − Πi1t )−1 · · · (1 − Πit,N t −2 )−1 , t

t

i

(A.2)

i

where Πikt = πi1t + · · · + πikt . From (A.1) and (11), we have Nit   t t t π (πti |Lti , α˜ it ) ∝ (α˜ it )Ni −1 (πiNt t )α˜ i −1 (1 − Πi1t )−1 · · · (1 − Πit,N t −2 )−1 (πikt )mik i

∝ (πi1t )

mti1

mt t i,N −1 i

· · · (πit,N t −1 ) i

i

(πiNt t )

α˜ it +mt t −1 iN i

i

k=1

(1 − Πi1t )−1 · · · (1 − Πit,N t −2 )−1 , i

where mtik is the number of Ltig whose value equals to k. Hence, the conditional distribution π (π ti |Lti , α˜ it ) is a generalized Dirichlet distribution Dir (ati1∗ , bti1∗ , . . . , ati,∗N t −1 , bti,∗N t −1 ), where atik∗ = 1 + mtik and btik∗ = α˜ it + i

Nit

j=k+1

i

mtij .

(ii) Conditional distribution π (Z˜ ti ∗ |Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G): Let Sit be the set of values in {1, 2, . . . , Nit } which exclude {Lti1∗ , . . . , Ltim∗ t }. Based on the semiparametric model defined in Section 2.2, we have i

π (Z˜ ti ∗ |Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G) ∝ π (Z˜ ti ∗ |µtz˜ i , ψz˜ti )π (ϵ1ti , . . . , ϵGti |Z˜ ti ∗ , Lti , ψtii )   t  mi   ∝ π (Z˜ikt ∗ |µtz˜ i , ψz˜ti )  π (Z˜iLt ∗t ∗ |Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G) , {k∈Sit }

j =1

ij

where

π (Z˜iLt ∗t ∗ |Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G) ∝ π (Z˜iLt ∗t ∗ |µtz˜ i , ψz˜ti ) ij

ij



∗ π (ϵgti |ζgti = Z˜iLt ∗t ∗ , ψtii )

{g :Ltig =Ltij∗ }

     1 1 exp − (ϵgti − Z˜iLt ∗t ∗ )2 ∝ exp − t (Z˜iLt ∗t ∗ − µtz˜ i )2 ij ij 2ψz˜ i 2ψtii {g :Ltig =Ltij∗ }        1  ∝ exp − (ntij /ψtii + 1/ψz˜ti )(Z˜iLt ∗t ∗ )2 − 2 µtz˜ i /ψz˜ti + ϵgti /ψtii  Z˜iLt ∗t ∗   2  ij ij t t∗ {g :Lig =Lij }

 ∝ exp −

1 2ψz˜tij∗

 (Z˜iLt ∗t ∗ − µtz˜∗ij )2 , ij

ij

106

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

in which µtz˜∗ij = ψz˜tij∗ (µtz˜ i /ψz˜ti +

ψz˜tij∗ = (ntij /ψtii + 1/ψz˜ti )−1 , and ntij is the number of times Ltij∗ occurs in Lti . Hence, given Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G, we have Z˜ikt ∗ , k = 1, . . . , Nit are conditionally independent, where the t ∗[L] components of Z˜ i follows N (µtz˜ i , ψz˜ti ), and [Z˜iLt ∗t ∗ |Lti , µtz˜ i , ψz˜ti , ψtii , ϵgti , g = 1, . . . , G] ∼ N (µtz˜∗ij , ψz˜tij∗ ), j = 1, . . . , mti . 

{k:Ltik =Ltij∗ } ϵkti /ψtii ),

ij

(iii) Conditional distribution [ |π , Z˜ ti ∗ , ψtii , ϵgti , g = 1, . . . , G]: Note that this conditional distribution is proportional to G t ˜ t∗ |πti )]. It follows that Lti1 , . . . , LtiG are independent, given πti , Z˜ ti ∗ , ψtii , and ϵgti , g = 1, . . . , G. g =1 [π (ϵgti |Lig , Zi , ψtii )π ( t t∗ t∗ From (11), we have Lig follows a multinomial distribution with parameters (πig1 , . . . , πigN t ) which are proportional to Lti t Lig

t i

i

t∗ 2 (πi1t exp{−(ϵgti − Z˜i1t ∗ )2 /(2ψtii )}, . . . , πiNt t exp{−(ϵgti − Z˜iN t ) /(2ψtii )}). i

i

(iv) Conditional distribution [µtz˜ i |Z˜ ti ∗ , ψz˜ti ]: It follows from (8) that this conditional distribution is proportional to:

   Nit   1 1 1 t ∗ t 2 t 2 t t ∗ 2 (Z˜ − µz˜ i ) − exp − (µ ) ∝ exp − t ∗ (µz˜ i − µz˜ i ) ,  2ψz˜ti k=1 ik 2ψµt i z˜ i  2ψµi  

where µtz˜∗i = ψµt ∗i

Nit

k=1

Z˜ikt ∗ /ψz˜ti and ψµt ∗i = 1/(Nit /ψz˜ti + 1/ψµt i ). Hence, [µtz˜ i |Z˜ ti ∗ , ψz˜ti ] ∼ N (µtz˜∗i , ψµt ∗i ).

(v) Conditional distribution [ψz˜ti |Z˜ ti ∗ , µtz˜ i ]: Similarly, this conditional distribution is proportional to

 

 Nit   1 t t (ψz˜ti )−Ni /2 exp − t (Z˜ikt ∗ − µtz˜ i )2 (ψz˜ti )−(c1i +1) exp{−c2it /ψz˜ti }  2ψz˜ i k=1      Nit    1 t t ∝ (ψz˜ti )−(Ni /2+c1i +1) exp −  (Z˜ikt ∗ − µtz˜ i )2 + c2it  (ψz˜ti )−1 .   2 k=1 t Hence, [(ψz˜ti )−1 |Z˜ ti ∗ , µtz˜ i ] ∼ Gamma(c1i + Nit /2, c2it +

Nit

k=1 at1i at2i

(vi) Conditional distribution [α˜ it |π ti ]: As α˜ it ∼ Gamma(

,

(Z˜ikt ∗ − µtz˜ i )2 /2). ), it follows from (A.2) that

π (α˜ |π ) ∝ π (π |α˜ )π (α˜ ) t i

t i

t i

t i

t i

∝ [(α˜ it )Ni −1 (πiNt t )α˜ i −1 (1 − Πi1t )−1 · · · (1 − Πit,N t −2 )−1 ](α˜ it )a1i −1 exp(−at2i α˜ it ) t

t

t

i

i

∝ (α˜ it )(Ni +a1i −1)−1 exp{−(at2i − log(πiNt t ))α˜ it } i     Nit −1    t t log(1 − νikt ∗ ) α˜ it . ∝ (α˜ it )(Ni +a1i −1)−1 exp − at2i −   k=1 t

t

Hence, we have [α˜ it |π ti ] ∼ Gamma(Nit + at1i − 1, at2i −

Nit −1 k=1

log(1 − νikt ∗ )).

Appendix B. Sampling from [θ, α, W |Y , Ω0 , Ω1 , Ξ , X , Z ] Observations can be simulated from [θ, α, W|Y, Ω0 , Ω1 , Ξ , X, Z] by completing the following updates in each cycle of the sampler. (i) Sampling from [Λt , At , ψtkk , k = 1, . . . , p, t = 1, . . . , T |Y, Ω1 , W, X, Ξ ]:

˜ gt = ugt −yg −At cgt −ζ gt , u˜ Agt = ugt −yg − Let ζ gt = (0, . . . , 0, ζgt ,s+1 , . . . , ζgt ,p )′ , Ω1t = (ω1t , . . . , ωGt ), Ct = (c1t , . . . , cGt ), u ˜ t = (u˜ 1t , . . . , u˜ Gt ), U˜ At = (u˜ A1t , . . . , u˜ AGt ), and U˜ ′tk and U˜ Atk′ be the kth rows of U˜ t and U˜ At , respectively. Moreover, Λt ωgt − ζ gt , U −1 −1 −1 + Ω1t Ω′1t )−1 , Λtϵ∗k = Htϵ∗k (Htϵ k Λtϵ k + Ω1t U˜ tk ), β˜ ϵt ∗k = β˜ ϵt k + 2−1 (U˜ ′tk U˜ tk − Λtϵ∗′k Htϵ∗k Λtϵ∗k + Λtϵ′k Htϵ k Λtϵ k ), −1 −1 −1 ˜ A −1 = (HtAk + ψtkk Ct C′t )−1 , and Atϵ∗k = HtAk∗ (HtAk Atϵ k + ψtkk Ct Utk ). It can be shown that, for k = 1, . . . , p and t = 1, . . . , T ,

Htϵ∗k = (Htϵ k

HtAk∗

−1

−1 [ψtkk |Y, Ω1 , W, X, Ξ , At ] ∼ Gamma(G/2 + α˜ ϵt k , β˜ ϵt ∗k ), −1 [Λtk |Y, Ω1 , W, X, Ξ , At , ψtkk ] ∼ N (Λtϵ∗k , ψtkk Htϵ∗k ), −1 [Atk |Y, Ω1 , W, X, Ξ , Λt , ψtkk ] ∼ N (Atϵ∗k , HtAk∗ ).

(ii) Sampling from [Λ0 , A0 , ψ0kk , k = 1, . . . , p|Y, Ω0 , Ξ ]: Let C0 = (c10 , . . . , cG0 ), ζ g0 = (0, . . . , 0, ζg0,s+1 , . . . , ζg0,p )′ ,

˜ = (˜y1 , . . . , y˜ G ), Y˜ A = (˜yA1 , . . . , y˜ AG ), and Y˜ ′k and Y˜ Ak ′ be the kth rows y˜ g = yg − A0 cg0 − ζ g0 , y˜ Ag = yg − Λ0 ωg0 − ζ g0 , Y

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

−1

107

−1

˜ and Y˜ A , respectively. Moreover, H0ϵ k∗ = (H0ϵ k + Ω0 Ω′0 )−1 , Λ0ϵ k∗ = H0ϵ k∗ (H0ϵ k Λ0ϵ k + Ω0 Y˜ k ), β˜ ϵ0k∗ = β˜ ϵ0k + 2−1 (Y˜ ′k Y˜ k − of Y −1 −1 −1 −1 −1 −1 ˜ Ak ). It can be shown that, Λ0ϵ k∗′ H0ϵ k∗ Λ0ϵ k∗ + Λ0ϵ k′ H0ϵ k Λ0ϵ k ), H0Ak∗ = (H0Ak + ψ0kk C0 C′0 )−1 , and A0ϵ k∗ = H0Ak∗ (H0Ak A0ϵ k + ψ0kk C0 Y for k = 1, . . . , p, −1 [ψ0kk |Y, Ω0 , Ξ , A0 ] ∼ Gamma(G/2 + α˜ ϵ0k , β˜ ϵ0k∗ ), −1 [Λ0k |Y, Ω0 , Ξ , A0 , ψ0kk ] ∼ N (Λ0ϵ k∗ , ψ0kk H0ϵ k∗ ),

−1 [A0k |Y, Ω0 , Ξ , Λ0 , ψ0kk ] ∼ N (A0ϵ k∗ , H0Ak∗ ).

(iii) Sampling from [B0 , Φ, Λδ t , Ψ δ t |Ω1 ]: Let D0 = (d10 , . . . , dG0 ), Dt = (d1t , . . . , dGt ), F∗gt = F∗t (ωg1 , . . . , ωg ,t −1 , ξ gt ), ′ ˜ gt = ηgt − B0 dg0 , η˜ ∗gt = ηgt − Bt dgt − Γ ∗t F∗gt , E˜ t = (˜η1t , . . . , η˜ Gt ), E˜ ′tl be the lth F˜ ∗t = (F∗1t , . . . , F∗Gt ), F˜ t = (D′t , F˜ ∗′ t ), η

row of E˜ t . ωξ g = (ξ g1 , . . . , ξ gT )′ , and Ωξ = (ωξ 1 , . . . , ωξ G ). Moreover, Htδ∗l = (Htδ l ′

′

0 −1

−1

−1 + F˜ t F˜ ′t )−1 , Λtδ∗l = Htδ∗l (Htδl Λtδl + F˜ t E˜ tl ),

−1 β˜ t ∗ = β˜ δt l + 2−1 (E˜ ′tl E˜ tl − + H0Bl∗ = [HBl + ( and let B0δ l∗ = H0Bl∗ [H0Bl B0δ l + δlG  T ∗ ˜ gtl /ψδtll )dg0 ], where η˜ gtl is the lth element of η˜ ∗gt . It can be shown that, for l = 1, . . . , q1 and t = 1, . . . , T , g =1 ( t =1 η t ∗ −1 t ∗ Λtδ∗′ Λδl l Hδ l ∗

−1 Λtδ′l Htδl Λtδl ),

[ψδ−tll1 |Ω1 , B0 ] ∼ Gamma(G/2 + α˜ δt l , β˜ δt l∗ ), [B0l |Ω1 , Λδt , Ψ δt ] ∼ N (B0δl∗ , H0Bl∗ ),

−1 ′ −1 , t =1 ψδ tll )D0 D0 ]

T

[Λδtl |Ω1 , B0 , ψδ−tll1 ] ∼ N (Λtδ∗l , ψδtll Htδ∗l ),

1 −1 [Φ−1 |Ω1 ] ∼ Wishartq2 ×T [(Ωξ Ω′ξ + R− w ) , G + ρw ].

(iv) Sampling from [Φ0 |Ω0 ]: It can be shown that 1 −1 −1 ′ Φ− , G + ρ0 ]. 0 |Ω0 ∼ Wishartq0 [(Ω0 Ω0 + R0 )

(v) Sampling from [α, W|Y, Ω1 , Z, θ]: For k = 1, . . . , s and t = 1, . . . , T , let αtk be the parameter vector of the unknown thresholds corresponding to zgtk , Wt = (w1t , . . . , wGt ), Zt = (z1t , . . . , zGt ), Wtk and Ztk be the kth rows of Wt and Zt , respectively. It is natural to use the following non-informative prior distribution for the unknown thresholds in αtk : π(αtk,2 , . . . , αtk,btk −1 ) ∝ c for αtk,2 < · · · < αtk,btk −1 , where c is a constant. Given Y and Ω1 , the ordered categorical data Z and the thresholds corresponding to different rows are conditionally independent, and

π (αtk , Wtk |Y, Ω1 , Ztk , θ) = π (αtk |Y, Ω1 , Ztk , θ)π (Wtk |αtk , Y, Ω1 , Ztk , θ),

(B.1)

with

π (αtk |Y, Ω1 , Ztk , θ) ∝

G  

−1/2

Φ ∗ [ψtkk (αtk,zgtk +1 − ygk − A′tk cgt − Λ′tk ωgt )]

g =1

 −1/2 − Φ ∗ [ψtkk (αtk,zgtk − ygk − A′tk cgt − Λ′tk ωgt )] ,

(B.2)

and π(Wtk |αtk , Y, Ω1 , Ztk , θ) is a product of π (wgtk |αtk , Y, Ω1 , Ztk , θ), where

π (wgtk |αtk , Y, Ω1 , Ztk , θ) ∼ N (ygk + A′tk cgt + Λ′tk ωgt , ψtkk )I(αtk,zgtk ,αtk,zgtk +1 ] (wgtk ),

(B.3)

in which ygk and wgtk are the kth elements of yg and wgt , respectively, IA (w) is an index function which takes 1 if w ∈ A and 0 otherwise, and Φ ∗ (·) denotes the standard normal cumulative distribution function. We consider a MH algorithm for generating observations from the target distribution π (αtk , Wtk |Y, Ω1 , Ztk , θ). The joint proposal density of αtk and Wtk in the MH algorithm is given by (B.1). However, in the proposal distribution, αtk is not simulated from Eq. (B.2). Instead, at the lth MH iteration, a vector of thresholds (αtk,2 , . . . , αtk,btk −1 ) is generated from the univariate truncated normal distribution

αtk,z ∼ N (αtk(l),z , σα2tk )I (αtk,z −1 , αtk(1,)z +1 ](αtk,z ) for z = 2, . . . , btk − 1, where αtk(l),z is the current value of αtk,z at the lth iteration of the Gibbs sampler, and σα2tk is a preassigned constant. The acceptance probability for (αtk , Wtk ) as a new observation is min(1, Rtk ), where Rtk =

π (αtk , Wtk |Y, Ω1 , Ztk , θ)π (α(tkl) , W(tkl) |αtk , Wtk , Y, Ω1 , Ztk , θ) π (α(tkl) , W(tkl) |Y, Ω1 , Ztk , θ)π (αtk , Wtk |α(tkl) , W(tkl) , Y, Ω1 , Ztk , θ)

.

Let w ˜ gtk = ygk + A′tk cgt + Λ′tk ωgt , and it can be shown that

 Φ ∗ [(αtk(l),z +1 − αtk(l),z )/σαtk ] − Φ ∗ [(αtk,z −1 − αtk(l),z )/σαtk ]

btk −1

Rtk =

z =2

×

(l)

Φ ∗ [(αtk,z +1 − αtk,z )/σαtk ] − Φ ∗ [(αtk,z −1 − αtk,z )/σαtk ]

−1/2 −1/2 G  Φ ∗ [ψtkk (αtk,zgtk +1 − w ˜ gtk )] − Φ ∗ [ψtkk (αtk,zgtk − w ˜ gtk )] −1/2

g =1

(l)

−1/2

(l)

Φ ∗ [ψtkk (αtk,zgtk +1 − w ˜ gtk )] − Φ ∗ [ψtkk (αtk,zgtk − w ˜ gtk )]

.

As Rtk does not depend on Wtk , it is unnecessary to generate a new Wtk in any iteration in which the new value of αtk is not accepted. For an accepted αtk , a new Wtk is simulated from Eq. (B.3).

108

X.-Y. Song et al. / Journal of Multivariate Analysis 121 (2013) 87–108

References [1] K.A. Bollen, Structural Equation with Latent Variables, Wiley, New York, 1989. [2] R.A. Browne, P.M. Monti, M.G. Myers, R.A. Martin, T. Rivinus, M.E.T. Dubreuil, D.J. Rohsenow, Depression among cocaine abusers in treatment: relation to cocaine and alcohol use and treatment outcome, American Journal of Psychiatry 155 (1998) 220–225. [3] M.H. Chen, D.K. Dey, J.G. Ibrahim, Bayesian criterion based model assessment for categorical data, Biometrika 91 (2004) 45–63. [4] D.J. Diggle, P. Heagerty, K.Y. Liang, S.L. Zeger, Analysis of Longitudinal Data, second ed., Oxford University Press, Oxford, U.K, 2002. [5] A. Gelman, Inference and monitoring convergence, in: W.R. Gilks, S. Richardson, Spiegelhalter (Eds.), Markov Chain Monte Carlo in Practice, Chapman and Hall, London, 1996, pp. 131–144. [6] Y. Hser, Prediction long-term stable recovery from heroin addiction: findings based on a 33-year follow-up study, Journal of Addiction Diseases 26 (2007) 51–60. [7] Y. Hser, M.E. Stark, A. Paredes, D. Huang, M.D. Anglin, R. Rawson, A 12-year follow-up of a treated cocaine-dependent sample, Journal of Substance Abuse Treatment 30 (2006) 219–226. [8] J.G. Ibrahim, M.H. Chen, D. Sinha, Bayesian Survival Analysis, Springer, New York, 2001. [9] H. Ishwaran, L.F. James, Gibbs sampling methods for stick-breaking priors, Journal of the American Statistical Association 96 (2001) 161–173. [10] H. Ishwaran, M. Zarepour, Markov chain Monte Carlo in approximate Dirichlet and Beta two-parameter process hierarchical models, Biometrika 87 (2000) 371–390. [11] Y. Kano, M. Berkane, P.M. Bentler, Statistical inference based on pseudo-maximum likelihood estimators in elliptical populations, Journal of the American Statistical Association 88 (1993) 135–143. [12] D.N. Kasarabada, M.D. Anglin, E. Khalsa-Denison, A. Paredes, Differential effects of treatment modality on psychosocial functioning of cocainedependent men, Journal of Clinical Psychology 55 (1999) 257–274. [13] A. Kottas, A.E. Gelfand, Bayesian semiparametric median regression modeling, Journal of the American Statistical Association 96 (2001) 1458–1468. [14] S.Y. Lee, Structural Equation Modelling: A Bayesian Approach, Wiley, New York, 2007. [15] S.Y. Lee, B. Lu, X.Y. Song, Semiparametric Bayesian analysis of structural equation models with fixed covariates, Statistics in Medicine 27 (2008) 2341–2360. [16] S.Y. Lee, Y.M. Xia, Maximum likelihood methods in treating outliers and symmetrically heavy-tailed distributions for nonlinear structural equation models with missing data, Psychometrika 71 (2006) 565–585. [17] X.Y. Song, S.Y. Lee, Basic and Advanced Bayesian Structural Equation Modeling: With Applications in the Medical and Behavioral Sciences, Wiley, London, 2012. [18] X.Y. Song, Z.H. Lu, Semiparametric latent variable models with Bayesian P-splines, Journal of Computational and Graphical Statistics 19 (2010) 590–608. [19] X.Y. Song, Z.H. Lu, Semiparametric transformation models with Bayesian P-splines, Statistics and Computing 22 (2012) 1085–1098. [20] X.Y. Song, Z.H. Lu, Y.I. Hser, S.Y. Lee, A Bayesian approach for analyzing longitudinal structural equation models, Structural Equation Modeling. A Multidisciplinary Journal 18 (2011) 183–194. [21] X.Y. Song, Y.M. Xia, S.Y. Lee, Bayesian semiparametric analysis of structural equation models with mixed continuous and unordered categorical variables, Statistics in Medicine 28 (2009) 2253–2276. [22] G. Verbeke, G. Molenberghs, Linear Mixed Models for Longitudinal Data, Springer, New York, 2000. [23] C. Weisner, G.T. Ray, J.R. Mertens, D.D. Satre, C. Moore, Short-term alcohol and drug treatment outcomes predict long-term outcomes, Drug and Alcohol Dependence 71 (2003) 281–294. [24] M. Yang, D.B. Dunson, D.D. Baird, Semiparametric Bayes hierarchical models with mean and variance constraints, Computational Statistics and Data Analysis 54 (2010) 2172–2186.