# Assessing Heterogeneity for Factor Analysis Model with Continuous and Ordinal Outcomes

Authors: Ye-Mao Xia; Jian-Wei Gou
Journal: Journal of Applied Mathematics
 (2016)

---

## Abstract

Factor analysis models with continuous and ordinal responses are a useful tool for assessing relations between the latent variables and mixed observed responses. These models have been successfully applied to many different fields, including behavioral, educational, and social-psychological sciences. However, within the Bayesian analysis framework, most developments are constrained within parametric families, of which the particular distributions are specified for the parameters of interest. This leads to difficulty in dealing with outliers and/or distribution deviations. In this paper, we propose a Bayesian semiparametric modeling for factor analysis model with continuous and ordinal variables. A truncated stick-breaking prior is used to model the distributions of the intercept and/or covariance structural parameters. Bayesian posterior analysis is carried out through the simulation-based method. Blocked Gibbs sampler is implemented to draw observations from the complicated posterior. For model selection, the logarithm of pseudomarginal likelihood is developed to compare the competing models. Empirical results are presented to illustrate the application of the methodology.


---

## Body

## 1. Introduction

Owing to its wide applications in behavioral and social science researches, analysis of factor analysis models with mixed data structure has received a lot of attention; see [1–6]. However, most of these methods are mainly developed within particular parametric distribution families such as the exponential family or normal scale mixture family, which have a limited role in dealing with the distributional deviations, in particular heterogeneity or multimodality of the data. Though some robust methods are developed to downweight the influence of the outliers [7–12], most of them are still confined to dealing with unimodality and are less effective for the asymmetric and/or multimodal problems.Recently, some authors focused on the Bayesian semiparametric modeling for latent variables model. For multivariate categorical data analysis, Kottas et al. [13] extended the traditional multivariate probit model [14–16] to a flexible underlying prior probability model. The usual single multivariate normal model for the latent variables is replaced by a mixture of normal priors with infinite number of components. And, for the latent variable model with fixed covariates and continuous responses, Lee et al. [17] established the semiparametric Bayesian hierarchal model for the structural equation models (SEMs) by relaxing the common normal distribution of exogenous factors to follow a finite-dimensional Dirichlet process [18]. Song et al. [19] developed a semiparametric Bayesian procedure for analyzing the latent variable model with unordered categorical data. For some recent advances in semiparametric analysis for factor analysis model, see [20–23] among others.In this paper, we developed a Bayesian semiparametric approach for analyzing factor analysis model with mixed continuous and ordinal responses. The methods are twofold. Firstly, we extended Kottas, Müller, and Quintana’s model to a more general multivariate model which contains factor variables. This extension aims to interpret the relationships between measurements and latent variables and explore correlations among the multiple manifest variables. Moreover, we treat the threshold parameters as unknown and estimate them simultaneously with other model parameters, thus providing a more flexible approach to fit the data. Secondly, we introduce the truncated Dirichlet process prior as the prior of the mean vector and variance-covariance parameters of unique errors and latent variables. This facilitates the interpretation of heterogeneity in the mean and/or covariance structure across the subjects.This paper is organized as follows. We first introduce the Bayesian semiparametric modeling framework for factor analysis model with continuous and ordinal variables. We then present the Markov chain Monte Carlo procedure for parameters estimation and model selection. Simulation studies and a real example are provided to illustrate the performance of the proposed procedure. We close with some remarks and concluding comments.


## 2. Model Description

### 2.1. Factor Analysis Model with Continuous and Ordinal Responses

Suppose that ap-dimensional mixed observed vector y
i
=
(
x
i
T
,
z
i
T
)
T contains r continuous variables x
i
=
(
x
i
1
,
…
,
x
i
r
)
T and s
=
p
-
r ordinal variables z
i
=
(
z
i
1
,
…
,
z
i
s
)
T with z
i
j taking an integral value in S
j
=
{
0,1
,
…
,
b
j
} for j
=
1
,
…
,
s, i
=
1
,
…
,
n. We assume that the observed ordinal vector z
i is related to the unobserved continuous vector u
i
=
(
u
i
1
,
…
,
u
i
s
)
T through(1)
z
i
j
=
l
if
τ
j
l
<
u
i
j
≤
τ
j
,
l
+
1
,
l
∈
S
j
,where {
τ
j
l
:
l
=
0
,
…
,
b
j, j
=
1
,
…
,
s
} is a set of unknown threshold parameters that define the categories: -
∞
=
τ
j
0
<
τ
j
1
<
⋯
<
τ
j
b
j
<
τ
j
,
b
j
+
1
=
∞. Hence, for the jth variable z
i
j, there are b
j
+
1 categories.Lety
i
∗
=
(
x
i
T
,
u
i
T
)
T denote the vector of continuous observed measurements and unobserved variables. For subject i, we formulate the dependence among y
i
j
∗’s through the following measurement model:(2)
y
i
∗
=
μ
+
Λ
ω
i
+
ϵ
i
,where μ is a p
×
1 intercept vector, Λ is a p
×
m factor loading matrix, ω
i is an m
×
1 vector of latent variables, and ϵ
i is a p
×
1 vector of measurement errors which is independent of ω
i. In many applications, ω
i may represent the hypothesized factors underlying manifest responses and/or unobserved heterogeneity not explained by covariates.The latent variable model with mixed continuous and ordinal responses defined by (1) and (2) faces two sources of identification problems. The first one is associated with the determinacy of latent variables y
∗ in modeling of categorical variables, and the second one is related to the uniqueness of the factor loadings matrix. To solve the first problem, we use the common method [24] to fix endpoints τ
j
1 and τ
j
b
j(j
=
1
,
…
,
s) at preassigned values. For the second problem, we follow the usual practice in structural equation modeling to identify the covariance matrix of y
i
∗ by fixing appropriate elements in Λ at preassigned values.Letθ be the parametric vector formed by the unknown parameters contained in {
μ
,
Ψ
ϵ
,
Φ
} and let ϑ denote the free parameters contained in factor loading matrices Λ and τ
=
(
τ
1
T
,
…
,
τ
s
T
)
T with τ
j
=
(
τ
j
1
,
…
,
τ
j
b
j
)
T. Based on the assumptions of (2), the conditional distribution of y
i
∗ given (
θ
,
ϑ
) is a normal distribution with mean vector μ and covariance matrix Σ
(
θ
,
ϑ
)
=
Λ
Φ
Λ
T
+
Ψ
ϵ.Note that the latent factors here play an important role in characterizing the associations between the observed variables. It can be seen clearly thatz
i and x
i are dependent when ω
i is integrated out. The marginal density of y
i is given by (3)
p
y
i
∣
ϑ
,
θ
=
∫
p
x
i
∣
ω
i
,
ϑ
,
θ
p
z
i
∣
ω
i
,
ϑ
,
θ
p
ω
i
∣
ϑ
,
θ
d
ω
iwith(4)
p
z
i
∣
ω
i
,
ϑ
,
θ
=
∏
j
=
1
s
Φ
c
τ
j
z
i
j
+
1
-
μ
r
+
j
-
Λ
r
+
j
T
ω
i
ψ
ϵ
r
+
j
-
Φ
c
τ
j
z
i
j
-
μ
r
+
j
-
Λ
r
+
j
T
ω
i
ψ
ϵ
r
+
j
,in which Φ
c
(
·
) is the standard normal cumulative distribution function.


### 2.2. Bayesian Semiparametric Hierarchical Modeling

Letp
(
y
i
∗
∣
θ
,
ϑ
) be the conditional density of y
i
∗ given (
θ
,
ϑ
) and denote by F a prior distribution function of θ. Suppose that F is proper; we define the following mixture density: (5)
p
y
i
∗
∣
F
,
ϑ
=
∫
p
y
i
∗
∣
θ
,
ϑ
F
d
θ
,in which F
(
d
θ
) is the conditional distribution of θ given F. By taking a prior for ϑ and restricting F to be a parametric family of distributions indexed by θ, we complete the Bayesian parametric model specification. However, this restriction severely constrains the estimation of θ and produces estimators that shrink data values toward the same points. A more flexible modeling for y
i
∗ is to treat F as random and assign a prior for it. For this end, we introduce a latent variable vector θ
i
=
{
μ
i
,
Ψ
ϵ
i
,
Φ
i
} and assume that, given θ
i, y
i
∗’s are conditionally independent and drawn from p
(
y
i
∗
∣
θ
i
,
ϑ
). Furthermore, we suppose that θ
i’s are independent and identically distributed (i.i.d.) according to F with a prior P on it. As a result, we break the mixture model p
(
y
i
∗
∣
F
,
ϑ
) into (6)
y
i
∗
∣
θ
i
,
ϑ
~
i
n
d
p
y
i
∗
∣
θ
i
,
ϑ
,
θ
1
,
…
,
θ
n
∣
F
~
i
i
d
F
,
F
~
P
,where “ind” means “independent” and P is a prior of F.We consider the following truncated version of Dirichlet process forF:(7)
P
·
=
P
G
·
=
∑
k
=
1
G
π
k
δ
θ
k
∗
·
,in which δ
θ
k
∗
(
·
) denotes a discrete probability measure concentrated on atom θ
k
∗ and π
k
(
k
=
1
,
…
,
G
), independent of θ
k
∗, are random weights constructed through the following stick-breaking procedure:(8)
π
1
=
V
1
,
π
k
=
1
-
V
1
⋯
1
-
V
k
-
1
V
k
,
k
=
2
,
…
,
G
-
1
,
π
G
=
1
-
V
1
⋯
1
-
V
G
-
1
,with V
k
~
i
.
i
.
d
.
B
e
t
a
(
1
,
α
); θ
k
∗’s are i.i.d with common distribution F
0.Truncated Dirichlet process prior (7) can be considered as a truncation version of Dirichlet process [25–30] in the nonparametric Bayesian analysis. It can be shown that, under (7) and (8), for any Borel set A in R
p,(9)
E
F
A
=
F
0
A
,
V
a
r
F
A
=
F
0
A
1
-
F
0
A
α
+
1
1
+
α
G
.This indicates that F
0 can be served as the starting point or guess of F and α determines the concentration of the prior around F
0. In practice, the value of G is either set to a large, predetermined value (e.g., G
≥
100) or chosen empirically. For instance, Ishwaran and Zarepour [31] suggested that the adequacy of the truncation level, G, can be assessed by evaluating moments of the tail probability. Our simulation results have shown that G
=
100 is more than adequate for the model considered in the present context.Now, we specify the distributionF
0. Recalling that by convention θ
k
∗ is the collection of {
μ
k
∗
,
Ψ
ϵ
k
∗
,
Φ
k
∗
}, hence, we assume that(10)
F
0
μ
k
∗
,
Ψ
ϵ
k
∗
,
Φ
k
∗
∣
ν
,
Σ
ν
,
R
=
N
μ
k
∗
∣
ν
,
Σ
ν
·
∏
j
=
1
m
G
a
m
m
a
-
1
ψ
ϵ
k
j
∗
∣
α
ϵ
0
j
,
β
ϵ
0
j
·
W
i
s
h
a
r
t
-
1
Φ
k
∗
∣
ρ
0
,
R
-
1
,where ν
,
Σ
ν, and R are hyperparameters, Σ
ν
=
diag
⁡
{
σ
ν
1
,
…
,
σ
ν
p
} is a diagonal matrix with the kth diagonal element σ
ν
k, and R is an m
2
×
m
2 positive definite matrix; G
a
m
m
a
-
1
(
α
ϵ
0
j
,
β
ϵ
0
j
) refers to the inverse gamma distribution with shaper parameters α
ϵ
0
j and scale parameters β
ϵ
0
j, respectively, and W
i
s
h
a
r
t
-
1 denotes the inverse Wishart distribution [32].ModelingF in (7) into the random probability measure and incorporating the latent variable ω
i into (5) generate the following hierarchical model: for i
=
1
,
…
,
n,(11)
y
i
∗
∣
ω
i
,
θ
i
,
ϑ
~
i
n
d
N
μ
i
+
Λ
ω
i
,
Ψ
ϵ
i
,
ω
i
∣
θ
i
~
i
n
d
N
0
,
Φ
i
,
θ
i
∣
F
~
i
i
d
F
,
F
~
P
G
·
,where P
G is given by (7) and (8).


## 2.1. Factor Analysis Model with Continuous and Ordinal Responses

Suppose that ap-dimensional mixed observed vector y
i
=
(
x
i
T
,
z
i
T
)
T contains r continuous variables x
i
=
(
x
i
1
,
…
,
x
i
r
)
T and s
=
p
-
r ordinal variables z
i
=
(
z
i
1
,
…
,
z
i
s
)
T with z
i
j taking an integral value in S
j
=
{
0,1
,
…
,
b
j
} for j
=
1
,
…
,
s, i
=
1
,
…
,
n. We assume that the observed ordinal vector z
i is related to the unobserved continuous vector u
i
=
(
u
i
1
,
…
,
u
i
s
)
T through(1)
z
i
j
=
l
if
τ
j
l
<
u
i
j
≤
τ
j
,
l
+
1
,
l
∈
S
j
,where {
τ
j
l
:
l
=
0
,
…
,
b
j, j
=
1
,
…
,
s
} is a set of unknown threshold parameters that define the categories: -
∞
=
τ
j
0
<
τ
j
1
<
⋯
<
τ
j
b
j
<
τ
j
,
b
j
+
1
=
∞. Hence, for the jth variable z
i
j, there are b
j
+
1 categories.Lety
i
∗
=
(
x
i
T
,
u
i
T
)
T denote the vector of continuous observed measurements and unobserved variables. For subject i, we formulate the dependence among y
i
j
∗’s through the following measurement model:(2)
y
i
∗
=
μ
+
Λ
ω
i
+
ϵ
i
,where μ is a p
×
1 intercept vector, Λ is a p
×
m factor loading matrix, ω
i is an m
×
1 vector of latent variables, and ϵ
i is a p
×
1 vector of measurement errors which is independent of ω
i. In many applications, ω
i may represent the hypothesized factors underlying manifest responses and/or unobserved heterogeneity not explained by covariates.The latent variable model with mixed continuous and ordinal responses defined by (1) and (2) faces two sources of identification problems. The first one is associated with the determinacy of latent variables y
∗ in modeling of categorical variables, and the second one is related to the uniqueness of the factor loadings matrix. To solve the first problem, we use the common method [24] to fix endpoints τ
j
1 and τ
j
b
j(j
=
1
,
…
,
s) at preassigned values. For the second problem, we follow the usual practice in structural equation modeling to identify the covariance matrix of y
i
∗ by fixing appropriate elements in Λ at preassigned values.Letθ be the parametric vector formed by the unknown parameters contained in {
μ
,
Ψ
ϵ
,
Φ
} and let ϑ denote the free parameters contained in factor loading matrices Λ and τ
=
(
τ
1
T
,
…
,
τ
s
T
)
T with τ
j
=
(
τ
j
1
,
…
,
τ
j
b
j
)
T. Based on the assumptions of (2), the conditional distribution of y
i
∗ given (
θ
,
ϑ
) is a normal distribution with mean vector μ and covariance matrix Σ
(
θ
,
ϑ
)
=
Λ
Φ
Λ
T
+
Ψ
ϵ.Note that the latent factors here play an important role in characterizing the associations between the observed variables. It can be seen clearly thatz
i and x
i are dependent when ω
i is integrated out. The marginal density of y
i is given by (3)
p
y
i
∣
ϑ
,
θ
=
∫
p
x
i
∣
ω
i
,
ϑ
,
θ
p
z
i
∣
ω
i
,
ϑ
,
θ
p
ω
i
∣
ϑ
,
θ
d
ω
iwith(4)
p
z
i
∣
ω
i
,
ϑ
,
θ
=
∏
j
=
1
s
Φ
c
τ
j
z
i
j
+
1
-
μ
r
+
j
-
Λ
r
+
j
T
ω
i
ψ
ϵ
r
+
j
-
Φ
c
τ
j
z
i
j
-
μ
r
+
j
-
Λ
r
+
j
T
ω
i
ψ
ϵ
r
+
j
,in which Φ
c
(
·
) is the standard normal cumulative distribution function.


## 2.2. Bayesian Semiparametric Hierarchical Modeling

Letp
(
y
i
∗
∣
θ
,
ϑ
) be the conditional density of y
i
∗ given (
θ
,
ϑ
) and denote by F a prior distribution function of θ. Suppose that F is proper; we define the following mixture density: (5)
p
y
i
∗
∣
F
,
ϑ
=
∫
p
y
i
∗
∣
θ
,
ϑ
F
d
θ
,in which F
(
d
θ
) is the conditional distribution of θ given F. By taking a prior for ϑ and restricting F to be a parametric family of distributions indexed by θ, we complete the Bayesian parametric model specification. However, this restriction severely constrains the estimation of θ and produces estimators that shrink data values toward the same points. A more flexible modeling for y
i
∗ is to treat F as random and assign a prior for it. For this end, we introduce a latent variable vector θ
i
=
{
μ
i
,
Ψ
ϵ
i
,
Φ
i
} and assume that, given θ
i, y
i
∗’s are conditionally independent and drawn from p
(
y
i
∗
∣
θ
i
,
ϑ
). Furthermore, we suppose that θ
i’s are independent and identically distributed (i.i.d.) according to F with a prior P on it. As a result, we break the mixture model p
(
y
i
∗
∣
F
,
ϑ
) into (6)
y
i
∗
∣
θ
i
,
ϑ
~
i
n
d
p
y
i
∗
∣
θ
i
,
ϑ
,
θ
1
,
…
,
θ
n
∣
F
~
i
i
d
F
,
F
~
P
,where “ind” means “independent” and P is a prior of F.We consider the following truncated version of Dirichlet process forF:(7)
P
·
=
P
G
·
=
∑
k
=
1
G
π
k
δ
θ
k
∗
·
,in which δ
θ
k
∗
(
·
) denotes a discrete probability measure concentrated on atom θ
k
∗ and π
k
(
k
=
1
,
…
,
G
), independent of θ
k
∗, are random weights constructed through the following stick-breaking procedure:(8)
π
1
=
V
1
,
π
k
=
1
-
V
1
⋯
1
-
V
k
-
1
V
k
,
k
=
2
,
…
,
G
-
1
,
π
G
=
1
-
V
1
⋯
1
-
V
G
-
1
,with V
k
~
i
.
i
.
d
.
B
e
t
a
(
1
,
α
); θ
k
∗’s are i.i.d with common distribution F
0.Truncated Dirichlet process prior (7) can be considered as a truncation version of Dirichlet process [25–30] in the nonparametric Bayesian analysis. It can be shown that, under (7) and (8), for any Borel set A in R
p,(9)
E
F
A
=
F
0
A
,
V
a
r
F
A
=
F
0
A
1
-
F
0
A
α
+
1
1
+
α
G
.This indicates that F
0 can be served as the starting point or guess of F and α determines the concentration of the prior around F
0. In practice, the value of G is either set to a large, predetermined value (e.g., G
≥
100) or chosen empirically. For instance, Ishwaran and Zarepour [31] suggested that the adequacy of the truncation level, G, can be assessed by evaluating moments of the tail probability. Our simulation results have shown that G
=
100 is more than adequate for the model considered in the present context.Now, we specify the distributionF
0. Recalling that by convention θ
k
∗ is the collection of {
μ
k
∗
,
Ψ
ϵ
k
∗
,
Φ
k
∗
}, hence, we assume that(10)
F
0
μ
k
∗
,
Ψ
ϵ
k
∗
,
Φ
k
∗
∣
ν
,
Σ
ν
,
R
=
N
μ
k
∗
∣
ν
,
Σ
ν
·
∏
j
=
1
m
G
a
m
m
a
-
1
ψ
ϵ
k
j
∗
∣
α
ϵ
0
j
,
β
ϵ
0
j
·
W
i
s
h
a
r
t
-
1
Φ
k
∗
∣
ρ
0
,
R
-
1
,where ν
,
Σ
ν, and R are hyperparameters, Σ
ν
=
diag
⁡
{
σ
ν
1
,
…
,
σ
ν
p
} is a diagonal matrix with the kth diagonal element σ
ν
k, and R is an m
2
×
m
2 positive definite matrix; G
a
m
m
a
-
1
(
α
ϵ
0
j
,
β
ϵ
0
j
) refers to the inverse gamma distribution with shaper parameters α
ϵ
0
j and scale parameters β
ϵ
0
j, respectively, and W
i
s
h
a
r
t
-
1 denotes the inverse Wishart distribution [32].ModelingF in (7) into the random probability measure and incorporating the latent variable ω
i into (5) generate the following hierarchical model: for i
=
1
,
…
,
n,(11)
y
i
∗
∣
ω
i
,
θ
i
,
ϑ
~
i
n
d
N
μ
i
+
Λ
ω
i
,
Ψ
ϵ
i
,
ω
i
∣
θ
i
~
i
n
d
N
0
,
Φ
i
,
θ
i
∣
F
~
i
i
d
F
,
F
~
P
G
·
,where P
G is given by (7) and (8).


## 3. Parameters Estimation and Model Selection

### 3.1. Prior Specifications and Estimation via Blocked Gibbs Sampler

LetΘ
∗
=
{
θ
k
∗
:
k
=
1
,
…
,
G
}. To implement Bayesian analysis, blocked Gibbs sampler is used to simulate observations from the posterior. The key for blocked Gibbs sampler is to recast model (11) completely by introducing the cluster variables L
=
(
L
1
,
…
,
L
n
)
T such that θ
i
=
θ
L
i
∗
. Consequently, the semiparametric hierarchical model (11) can be reformulated as the following framework:(12)
y
i
∗
∣
ω
i
,
θ
i
,
ϑ
~
i
n
d
N
μ
i
+
Λ
ω
i
,
Ψ
ϵ
i
,
ω
i
∣
θ
i
~
i
n
d
N
0
,
Φ
i
,
L
i
=
·
∣
π
~
i
i
d
∑
k
=
1
G
π
k
δ
k
·
,
π
,
Θ
∗
~
p
π
p
Θ
∗
,
ϑ
~
p
ϑ
,
τ
~
p
τ
,where p
(
ϑ
) is a prior of ϑ, p
(
π
) is the stick-breaking prior given by (8) with [
V
i
∣
α
]
~
i
i
d
B
e
t
a
(
1
,
α
), and p
(
Θ
∗
) is the joint distribution of Θ
∗ given by(13)
p
Θ
∗
∣
ν
,
Σ
ν
,
R
=
∏
k
=
1
G
p
θ
k
∗
∣
ν
,
Σ
ν
,
R
=
∏
k
=
1
G
p
μ
k
∗
,
Ψ
ϵ
k
∗
,
Φ
k
∗
∣
ν
,
Σ
ν
,
R
,
μ
k
∗
,
Ψ
ϵ
k
∗
,
Φ
k
∗
∣
ν
,
Σ
ν
,
R
~
i
i
d
F
0in which F
0
(
·
∣
ν
,
Σ
ν
,
R
) is given in (10).For the Bayesian analysis, we need to specify priors for the parameters involved in the model. The whole parameters can be divided into two parts: parametric component part{
ϑ
,
τ
} and nonparametric component part {
ν
,
Σ
ν
,
R
,
α
}. For the parametric components, we assume that p
(
ϑ
,
τ
)
=
p
(
ϑ
)
p
(
τ
) with(14)
p
Λ
k
=
D
N
Λ
0
k
,
H
ϵ
0
k
,
p
τ
=
∏
j
=
1
s
p
τ
j
=
∏
j
=
1
s
p
τ
j
,
2
,
…
,
τ
j
,
b
j
-
1
∝
∏
j
=
1
s
I
τ
j
,
2
<
⋯
<
τ
j
,
b
j
-
1
,where Λ
k is a p
×
1 column vector that contains unknown parameters in the kth row of Λ.For the hyperparameterβ
=
{
ν
,
Σ
ν
,
R
,
α
}, we consider the following conjugate priors:(15)
ν
~
N
μ
0
,
Σ
0
,
R
~
W
i
s
h
a
r
t
-
1
ρ
0
ϕ
,
R
0
ϕ
,
σ
ν
k
~
G
a
m
m
a
-
1
κ
1
,
κ
2
,
α
~
G
a
m
m
a
τ
1
,
τ
2
.The hyperparameters μ
0
,
Σ
0
,
Λ
0
k, H
ϵ
0
k, R
0
ϕ
,
α
ϵ
0
j
,
β
ϵ
0
j
,
ρ
0
,
ρ
0
ϕ, κ
1, κ
2, τ
1, and τ
2 in (10), (14), and (15) are treated as known.LetY
(
n
×
p
)
=
(
y
1
,
…
,
y
n
)
T, Ω
=
(
ω
1
,
…
,
ω
n
)
T, and Y
∗
=
(
y
1
∗
,
…
,
y
n
∗
)
T. Posterior analysis in relation to the complex p
(
ϑ
,
β
∣
Y
) is carried out through the data augmentation technique [33]. Specifically, we treat the latent quantities {
Ω
,
Y
∗
,
π
,
Θ
∗
,
L
} as missing data and augment them with the observed data. A sequence of random observations is generated from the joint posterior distribution p
(
Ω
,
Y
∗
,
ϑ
,
π
,
Θ
∗
,
L
,
β
∣
Y
) by the blocked Gibbs sampler [31, 34], coupled with the Metropolis-Hastings algorithm [35, 36]: given {
Ω
(
l
)
,
Y
∗
(
l
)
,
π
(
l
)
,
Θ
∗
(
l
)
,
L
(
l
)
} at the lth iterationdrawΩ
(
l
+
1
) from p
(
Ω
∣
Y
∗
(
l
), ϑ
(
l
), π
(
l
), Θ
∗
(
l
), L
(
l
), β
(
l
), Y
),
draw(
ϑ
(
l
+
1
)
,
Y
∗
(
l
+
1
)
) from p
(
ϑ
,
Y
∗
∣
Ω
(
l
+
1
), π
(
l
), Θ
∗
(
l
), L
(
l
), β
(
l
), Y
),
draw(
π
(
l
+
1
)
,
Θ
∗
(
l
+
1
)
) from p
(
π
,
Θ
∗
∣
Ω
(
l
+
1
), Y
∗
(
l
+
1
), ϑ
(
l
+
1
), L
(
l
), β
(
l
), Y
),
drawL
(
l
+
1
) from p
(
L
∣
Ω
(
l
+
1
), Y
∗
(
l
+
1
), ϑ
(
l
+
1
), π
(
l
+
1
), Θ
∗
(
l
+
1
), β
(
l
), Y
),
drawβ
(
l
+
1
) from p
(
β
∣
Ω
(
l
+
1
), Y
∗
(
l
+
1
), ϑ
(
l
+
1
), π
(
l
+
1
), Θ
∗
(
l
+
1
), L
(
l
+
1
), Y
),and form {
Ω
(
l
+
1
)
,
Y
∗
(
l
+
1
)
,
π
(
l
+
1
)
,
Θ
∗
(
l
+
1
)
,
L
(
l
+
1
)
}. It can be shown that as l tends to infinity, the empirical distribution of {
Ω
(
l
)
,
Y
∗
(
l
)
,
π
(
l
)
,
Θ
∗
(
l
)
,
L
(
l
)
} converges to p
(
Ω
,
Y
∗
,
ϑ
,
π
,
Θ
∗
,
L
,
β
∣
Y
) at any geometrical rate. The full conditional distributions and the implementation of the above algorithm are given in the Appendix.


### 3.2. Model Selection

Model selection is an important issue in Bayesian semiparametric modeling for latent variable model since it is of practical interest to compare different modelings for factor analytic models. Formal Bayesian model selection is accomplished by comparing the marginal predictive distribution of data across models. Consider the problem of comparing competing modelsM
1 and M
2. Let p
(
Y
∣
M
1
) and p
(
Y
∣
M
2
) denote the marginal density of data Y under M
1 and M
2, respectively. A popular choice for selecting models is achieved via Bayes factor (BF) (e.g., [37–39]). However, in view of the fact that computing BF involves the high-dimensional density which is hard to estimate well, we prefer comparing the following logarithm of pseudomarginal likelihood (LPML) [40, 41]:(16)
L
P
M
L
Y
=
∑
i
=
1
n
log
⁡
C
P
O
i
,where C
P
O
i is known as the conditional predictive ordinate (CPO) defined as(17)
C
P
O
i
=
p
y
i
∣
Y
i
=
∫
1
p
y
i
∣
Y
i
,
ϑ
,
Θ
∗
,
π
p
ϑ
,
Θ
∗
,
π
∣
Y
d
ϑ
d
Θ
∗
d
π
-
1
=
∫
1
p
y
i
∣
ϑ
,
Θ
∗
,
π
p
ϑ
,
Θ
∗
,
π
∣
Y
d
ϑ
d
Θ
∗
d
π
-
1
.Here, Y
(
i
) is the data set Y with y
i removed. Obviously, from (17), we can see that C
P
O
i is the marginal posterior predictive density of y
i given Y
(
i
) and can be interpreted as the height of this marginal density at y
i. Thus, small values of LPML imply that Y does not support the model.Based on MCMC sample{
(
Θ
∗
t
,
π
t
,
L
t
,
ϑ
t
,
β
t
)
:
t
=
1
,
…
,
T
} already available in the estimation, a consistent estimate for LPML can be obtained via ergodic average given by (18)
L
P
M
L
^
Y
=
-
∑
i
=
1
n
log
⁡
1
T
∑
t
=
1
T
1
p
y
i
∣
ϑ
t
,
Θ
∗
t
,
π
t
.It is noted that, under our proposed model,(19)
p
y
i
∣
ϑ
,
Θ
∗
,
π
=
∫
p
y
i
∣
ω
i
,
L
i
,
ϑ
,
Θ
∗
p
ω
i
∣
L
i
,
ϑ
,
Θ
∗
p
L
i
∣
π
d
ω
i
d
L
iwhich is complicated due to the existence of Ω and L. This can be solved by Monte Carlo approximation. Specifically, given the current values {
ϑ
(
l
)
,
Θ
∗
(
l
)
,
π
(
l
)
} at the lth iteration, we draw (i) L
i
l
,
h from p
(
L
i
∣
π
(
l
)
) and (ii) ω
i
l
,
h from p
(
ω
i
∣
L
i
l
,
h
,
ϑ
(
l
)
,
Θ
∗
(
l
)
) for h
=
1
,
…
,
H and then evaluate p
(
y
i
∣
ϑ
(
l
)
,
Θ
∗
(
l
)
,
π
(
l
)
) at the observation y
i through (20)
p
^
y
i
∣
ϑ
l
,
Θ
∗
l
,
π
l
≈
1
H
∑
h
=
1
H
p
y
i
∣
ω
i
l
,
h
,
L
i
l
,
h
,
ϑ
l
,
Θ
∗
l
.Obviously, the distributions involved in (i) and (ii) are standard and sampling is rather straightforward and fast.


## 3.1. Prior Specifications and Estimation via Blocked Gibbs Sampler

LetΘ
∗
=
{
θ
k
∗
:
k
=
1
,
…
,
G
}. To implement Bayesian analysis, blocked Gibbs sampler is used to simulate observations from the posterior. The key for blocked Gibbs sampler is to recast model (11) completely by introducing the cluster variables L
=
(
L
1
,
…
,
L
n
)
T such that θ
i
=
θ
L
i
∗
. Consequently, the semiparametric hierarchical model (11) can be reformulated as the following framework:(12)
y
i
∗
∣
ω
i
,
θ
i
,
ϑ
~
i
n
d
N
μ
i
+
Λ
ω
i
,
Ψ
ϵ
i
,
ω
i
∣
θ
i
~
i
n
d
N
0
,
Φ
i
,
L
i
=
·
∣
π
~
i
i
d
∑
k
=
1
G
π
k
δ
k
·
,
π
,
Θ
∗
~
p
π
p
Θ
∗
,
ϑ
~
p
ϑ
,
τ
~
p
τ
,where p
(
ϑ
) is a prior of ϑ, p
(
π
) is the stick-breaking prior given by (8) with [
V
i
∣
α
]
~
i
i
d
B
e
t
a
(
1
,
α
), and p
(
Θ
∗
) is the joint distribution of Θ
∗ given by(13)
p
Θ
∗
∣
ν
,
Σ
ν
,
R
=
∏
k
=
1
G
p
θ
k
∗
∣
ν
,
Σ
ν
,
R
=
∏
k
=
1
G
p
μ
k
∗
,
Ψ
ϵ
k
∗
,
Φ
k
∗
∣
ν
,
Σ
ν
,
R
,
μ
k
∗
,
Ψ
ϵ
k
∗
,
Φ
k
∗
∣
ν
,
Σ
ν
,
R
~
i
i
d
F
0in which F
0
(
·
∣
ν
,
Σ
ν
,
R
) is given in (10).For the Bayesian analysis, we need to specify priors for the parameters involved in the model. The whole parameters can be divided into two parts: parametric component part{
ϑ
,
τ
} and nonparametric component part {
ν
,
Σ
ν
,
R
,
α
}. For the parametric components, we assume that p
(
ϑ
,
τ
)
=
p
(
ϑ
)
p
(
τ
) with(14)
p
Λ
k
=
D
N
Λ
0
k
,
H
ϵ
0
k
,
p
τ
=
∏
j
=
1
s
p
τ
j
=
∏
j
=
1
s
p
τ
j
,
2
,
…
,
τ
j
,
b
j
-
1
∝
∏
j
=
1
s
I
τ
j
,
2
<
⋯
<
τ
j
,
b
j
-
1
,where Λ
k is a p
×
1 column vector that contains unknown parameters in the kth row of Λ.For the hyperparameterβ
=
{
ν
,
Σ
ν
,
R
,
α
}, we consider the following conjugate priors:(15)
ν
~
N
μ
0
,
Σ
0
,
R
~
W
i
s
h
a
r
t
-
1
ρ
0
ϕ
,
R
0
ϕ
,
σ
ν
k
~
G
a
m
m
a
-
1
κ
1
,
κ
2
,
α
~
G
a
m
m
a
τ
1
,
τ
2
.The hyperparameters μ
0
,
Σ
0
,
Λ
0
k, H
ϵ
0
k, R
0
ϕ
,
α
ϵ
0
j
,
β
ϵ
0
j
,
ρ
0
,
ρ
0
ϕ, κ
1, κ
2, τ
1, and τ
2 in (10), (14), and (15) are treated as known.LetY
(
n
×
p
)
=
(
y
1
,
…
,
y
n
)
T, Ω
=
(
ω
1
,
…
,
ω
n
)
T, and Y
∗
=
(
y
1
∗
,
…
,
y
n
∗
)
T. Posterior analysis in relation to the complex p
(
ϑ
,
β
∣
Y
) is carried out through the data augmentation technique [33]. Specifically, we treat the latent quantities {
Ω
,
Y
∗
,
π
,
Θ
∗
,
L
} as missing data and augment them with the observed data. A sequence of random observations is generated from the joint posterior distribution p
(
Ω
,
Y
∗
,
ϑ
,
π
,
Θ
∗
,
L
,
β
∣
Y
) by the blocked Gibbs sampler [31, 34], coupled with the Metropolis-Hastings algorithm [35, 36]: given {
Ω
(
l
)
,
Y
∗
(
l
)
,
π
(
l
)
,
Θ
∗
(
l
)
,
L
(
l
)
} at the lth iterationdrawΩ
(
l
+
1
) from p
(
Ω
∣
Y
∗
(
l
), ϑ
(
l
), π
(
l
), Θ
∗
(
l
), L
(
l
), β
(
l
), Y
),
draw(
ϑ
(
l
+
1
)
,
Y
∗
(
l
+
1
)
) from p
(
ϑ
,
Y
∗
∣
Ω
(
l
+
1
), π
(
l
), Θ
∗
(
l
), L
(
l
), β
(
l
), Y
),
draw(
π
(
l
+
1
)
,
Θ
∗
(
l
+
1
)
) from p
(
π
,
Θ
∗
∣
Ω
(
l
+
1
), Y
∗
(
l
+
1
), ϑ
(
l
+
1
), L
(
l
), β
(
l
), Y
),
drawL
(
l
+
1
) from p
(
L
∣
Ω
(
l
+
1
), Y
∗
(
l
+
1
), ϑ
(
l
+
1
), π
(
l
+
1
), Θ
∗
(
l
+
1
), β
(
l
), Y
),
drawβ
(
l
+
1
) from p
(
β
∣
Ω
(
l
+
1
), Y
∗
(
l
+
1
), ϑ
(
l
+
1
), π
(
l
+
1
), Θ
∗
(
l
+
1
), L
(
l
+
1
), Y
),and form {
Ω
(
l
+
1
)
,
Y
∗
(
l
+
1
)
,
π
(
l
+
1
)
,
Θ
∗
(
l
+
1
)
,
L
(
l
+
1
)
}. It can be shown that as l tends to infinity, the empirical distribution of {
Ω
(
l
)
,
Y
∗
(
l
)
,
π
(
l
)
,
Θ
∗
(
l
)
,
L
(
l
)
} converges to p
(
Ω
,
Y
∗
,
ϑ
,
π
,
Θ
∗
,
L
,
β
∣
Y
) at any geometrical rate. The full conditional distributions and the implementation of the above algorithm are given in the Appendix.


## 3.2. Model Selection

Model selection is an important issue in Bayesian semiparametric modeling for latent variable model since it is of practical interest to compare different modelings for factor analytic models. Formal Bayesian model selection is accomplished by comparing the marginal predictive distribution of data across models. Consider the problem of comparing competing modelsM
1 and M
2. Let p
(
Y
∣
M
1
) and p
(
Y
∣
M
2
) denote the marginal density of data Y under M
1 and M
2, respectively. A popular choice for selecting models is achieved via Bayes factor (BF) (e.g., [37–39]). However, in view of the fact that computing BF involves the high-dimensional density which is hard to estimate well, we prefer comparing the following logarithm of pseudomarginal likelihood (LPML) [40, 41]:(16)
L
P
M
L
Y
=
∑
i
=
1
n
log
⁡
C
P
O
i
,where C
P
O
i is known as the conditional predictive ordinate (CPO) defined as(17)
C
P
O
i
=
p
y
i
∣
Y
i
=
∫
1
p
y
i
∣
Y
i
,
ϑ
,
Θ
∗
,
π
p
ϑ
,
Θ
∗
,
π
∣
Y
d
ϑ
d
Θ
∗
d
π
-
1
=
∫
1
p
y
i
∣
ϑ
,
Θ
∗
,
π
p
ϑ
,
Θ
∗
,
π
∣
Y
d
ϑ
d
Θ
∗
d
π
-
1
.Here, Y
(
i
) is the data set Y with y
i removed. Obviously, from (17), we can see that C
P
O
i is the marginal posterior predictive density of y
i given Y
(
i
) and can be interpreted as the height of this marginal density at y
i. Thus, small values of LPML imply that Y does not support the model.Based on MCMC sample{
(
Θ
∗
t
,
π
t
,
L
t
,
ϑ
t
,
β
t
)
:
t
=
1
,
…
,
T
} already available in the estimation, a consistent estimate for LPML can be obtained via ergodic average given by (18)
L
P
M
L
^
Y
=
-
∑
i
=
1
n
log
⁡
1
T
∑
t
=
1
T
1
p
y
i
∣
ϑ
t
,
Θ
∗
t
,
π
t
.It is noted that, under our proposed model,(19)
p
y
i
∣
ϑ
,
Θ
∗
,
π
=
∫
p
y
i
∣
ω
i
,
L
i
,
ϑ
,
Θ
∗
p
ω
i
∣
L
i
,
ϑ
,
Θ
∗
p
L
i
∣
π
d
ω
i
d
L
iwhich is complicated due to the existence of Ω and L. This can be solved by Monte Carlo approximation. Specifically, given the current values {
ϑ
(
l
)
,
Θ
∗
(
l
)
,
π
(
l
)
} at the lth iteration, we draw (i) L
i
l
,
h from p
(
L
i
∣
π
(
l
)
) and (ii) ω
i
l
,
h from p
(
ω
i
∣
L
i
l
,
h
,
ϑ
(
l
)
,
Θ
∗
(
l
)
) for h
=
1
,
…
,
H and then evaluate p
(
y
i
∣
ϑ
(
l
)
,
Θ
∗
(
l
)
,
π
(
l
)
) at the observation y
i through (20)
p
^
y
i
∣
ϑ
l
,
Θ
∗
l
,
π
l
≈
1
H
∑
h
=
1
H
p
y
i
∣
ω
i
l
,
h
,
L
i
l
,
h
,
ϑ
l
,
Θ
∗
l
.Obviously, the distributions involved in (i) and (ii) are standard and sampling is rather straightforward and fast.


## 4. A Simulation Study

In this section, a simulation study to evaluate the performance of the proposed procedure is conducted. The goal is to assess the accuracy of estimates under parametric, partly exchangeable, and semiparametric modelings when data take on the multimodality or heterogeneity. We consider the situation in which each observed vector consists of three-dimensional continuous vector and three-dimensional ordinal vector with threshold valuesτ
j
=
(
-
1
.
0
∗
,
-
0.6,0.3,1
.
0
∗
)
T
(
j
=
1,2
,
3
). We generate Y by first generating Y
∗ with y
i
∗
=
(
x
i
T
,
u
i
T
)
T from the mixture of two factor analytic models with weights 0.45 and 0.55 and then transforming u
i into z
i  (s
=
3) via (1) to create the ordinal observations, where x
i represents a 6
×
1 observed continuous random vector and u
i is a 3
×
1 latent continuous random vector. Each component in the mixture model is specified through the following measurement model: for m
=
1,2,(21)
y
i
∗
=
μ
m
+
Λ
ω
i
m
+
ϵ
i
m
,
ϵ
i
m
~
N
0
,
Ψ
ϵ
m
.The parameters involved in the components of mixture model are taken as μ
(
1
)
=
-
1.5
×
1
6, μ
(
2
)
=
1.0
×
1
6, Ψ
ϵ
(
1
)
=
0.36
I
6, Ψ
ϵ
(
2
)
=
I
6,(22)
Λ
T
=
1
∗
0.8
0.8
0
∗
0
∗
0
∗
0
∗
0
∗
0
∗
1
∗
0.8
0.8
,
Φ
1
=
1
-
0.3
-
0.3
1.0
,
Φ
2
=
1
0.6
0.6
1.0
,in which 1
6 is a 6
×
1 vector with all elements equal to one and I
6 is a 6
×
6 identity matrix. The elements with asterisks involved in loading matrix Λ and threshold parameters {
τ
j
}
j
=
1
3 are treated as fixed for identifying model (see Section 2.1). Based on these settings, random sample with size 500 is generated and 100 replications are completed for each combination.Prior inputs in the prior distributions involved in the parametric components (see (14)) are as follows: H
ϵ
0
k and H
ζ
0 are diagonal matrices with the diagonal elements 1.0, and elements in {
Λ
0
k
,
Π
0
} are equal to the true values, while prior inputs in the prior distribution of superparameter β (see (15)) are μ
0
=
0
9, Σ
0
=
100
I
9, κ
1
=
κ
2
=
0.001, R
0
ϕ
=
0.01
I
2, ρ
0
=
ρ
0
ϕ
=
10, α
ϵ
0
k
=
β
ϵ
0
k
=
2.0, and τ
1
=
τ
2
=
2.0. Note that these values ensure approaching noninformative priors.A few test runs are conducted to explore the effect of truncated levels on the estimates of unknown parameters and the convergence of the blocked Gibbs sampler. We takeG
=
50, 60, 70, 80, 90, 100, 200, and 300 and calculate the total sum of the root mean square (RMS) of estimates (see below for details). The resulting values are 1.9830, 1.7382, 1.6582, 1.5548, 1.4194, 1.4128, 1.4108, and 1.4101, respectively. It can be seen that the total sum of the root mean square (RMS) becomes rather stable when G
≥
80. In the following analysis, we set G
=
100 in our data analysis. For the threshold parameters {
τ
j
k
:
j
=
1,2
,
3
,
k
=
2,3
}, we choose σ
M
H
j
k
2
=
0.002 (see Appendix) in MH algorithm to produce the acceptance rate about 0.40. Figure 1 gives the plots of EPSR (estimated potential scale reduction [42]) values of unknown free parameters in Λ, τ, and α against iterations for three groups of different starting values. It can be seen that the estimates converge in less than 1000 iterations. To be conservative, in the following analysis, we collect 3000 observations after 2000 “burn-in”s deleted to take posterior analysis. We first consider the performance of the proposed LMPL in model comparison. We compare the proposed model with the parametric model (denoted by PARA) and the partly exchangeable model (denoted by PAEX), which approximately correspond to α
=
+
∞ and α
=
0 under our proposal, respectively. The parametric model is defined by (23)
y
i
∗
∣
ω
i
,
θ
~
i
n
d
N
μ
+
Λ
ω
i
,
Ψ
ϵ
,
ω
i
∣
Φ
~
i
i
d
N
0
,
Φ
.The priors of the parameters are given by p
(
μ
)
=
D
N
(
μ
0
,
Σ
0
), p
(
Λ
,
Ψ
ϵ
)
=
D
∏
k
=
1
9
N
(
Λ
0
k
,
H
0
k
)
·
G
a
m
m
a
-
1
(
α
ϵ
0
k
,
β
ϵ
0
k
), and Φ
~
W
i
s
h
a
r
t
-
1
(
10,7.0
I
2
)
.Figure 1
Plot of the values of EPSR of unknown parameters against the number of iterations under different starting values for the simulated data.The partly exchangeable model is given by(24)
y
i
∗
∣
ω
i
,
θ
i
,
ϑ
~
i
n
d
N
μ
i
+
Λ
ω
i
,
Ψ
ϵ
i
,
ω
i
∣
θ
i
~
i
n
d
N
0
,
Φ
i
,where θ
i
=
{
μ
i
,
Ψ
ϵ
i
,
Φ
i
} are i.i.d. with distribution F
0
(
·
∣
ν
,
Σ
ν
,
R
) given in (5); the priors for the unknown parameter vector Λ and hyperparametric vector {
ν
,
Σ
ν
,
R
} are, respectively, given in (14) and (15).Under the foregoing settings for the hyperparameters, observations obtained through the blocked Gibbs sampler are used to compute the values of LPML for each scenario across 100 replications. For the parametric and partly exchangeable model, computing values of LPML is very straightforward and standard. For the semiparametric model, we draw 50 observations for approximatingp
(
y
i
∣
ϑ
,
Θ
,
π
). The values of LPML under parametric model, semiparametric model, and partly exchangeable model are, respectively, −6684.740, −6255.553, and −8487.259 with standard deviations 62.509, 151.480, and 147.742. Based on the LPML criteria, semiparametric model is selected, which is consistent with the fact that the true model takes on the multimodes. Moreover, according to our empirical results, the correct rates of LPML selecting the true model across 100 replications are about 0.93.Table1 gives the biases (BIAS), root mean squares (RMS), and standard deviations (SD) of estimates of unknown parameters across 100 replications under semiparametric models and parametric and partly exchangeable model, respectively. The measures BIAS, RMS, and SD are given as (25)
B
I
A
S
ϑ
^
j
=
ϑ
¯
j
-
ϑ
j
0
,
ϑ
¯
j
=
1
S
∑
r
=
1
S
ϑ
^
j
r
,
R
M
S
ϑ
^
j
=
1
S
∑
r
=
1
S
ϑ
^
j
r
-
ϑ
j
0
2
,
S
E
ϑ
^
j
=
1
S
∑
r
=
1
S
ϑ
^
j
r
-
ϑ
¯
j
2
,where S is the number of replications. It can be seen that estimates obtained through the proposed approach are reasonably accurate. The values of λ
^
j
k under our approach are smaller than those under parametric and exchangeable modelings in terms of the absolute values of BIAS and RMS. The results show that ignoring heterogeneity among the data may lead to biased estimates and incorrect interpretation of the analyzed phenomena. This also reflects that the factor loadings λ
j
k are not robust against the distributional deviations of inceptor, variance of unique errors, and covariances of latent factors.Table 1
Summary of the estimates under the parametric, partly exchangeable, and semiparametric approaches in analyzing simulated data: mixture data.
Para.
PARA
PAEX
SEMI
BIAS
RMS
SD
BIAS
RMS
SD
BIAS
RMS
SD
λ
21
0.130
0.136
0.043
0.216
0.171
0.047
−0.022
0.066
0.070
λ
31
0.131
0.138
0.043
0.212
0.181
0.047
−0.026
0.080
0.070
λ
52
0.140
0.147
0.046
0.263
0.164
0.072
−0.015
0.093
0.092
λ
62
0.144
0.150
0.046
0.266
0.152
0.071
−0.012
0.098
0.092
τ
12
0.060
0.091
0.040
−0.110
0.014
0.052
0.007
0.062
0.038
τ
13
0.070
0.085
0.033
−0.216
0.050
0.060
0.007
0.067
0.040
τ
22
0.070
0.098
0.040
−0.113
0.016
0.051
0.008
0.056
0.038
τ
23
0.077
0.095
0.032
−0.186
0.038
0.059
0.013
0.062
0.042
τ
32
0.062
0.092
0.040
−0.118
0.016
0.051
0.013
0.069
0.038
τ
33
0.083
0.098
0.032
−0.184
0.037
0.059
0.018
0.070
0.042Further simulation study is conducted to assess the performance of the proposed model and parametric model as well as the partly exchangeable model when data are generated from a single normal distribution. The population values of parameters are taken asμ
=
0
6, Ψ
ϵ
=
I
6, and (26)
Φ
=
1
0.3
0.3
1.0
.The values of factor loadings and threshold points are the same as those in previous mixture model. As usual, we take G
=
100 for truncated levels. The sample size is set to 62 which is analogous to the real example. The inputs for superparameters involved in priors are set the same as that in mixture model. The results based on 100 replications are summarized in Table 2.Table 2
Summary of the estimates under the parametric, partly exchangeable, and semiparametric approaches in analyzing simulated data: normal data.
Para.
PARA
PAEX
SEMI
BIAS
RMS
SD
BIAS
RMS
SD
BIAS
RMS
SD
λ
21
0.109
0.103
0.068
0.216
0.201
0.147
0.152
0.125
0.130
λ
31
−0.114
0.108
0.066
0.212
0.201
0.207
0.138
0.137
0.120
λ
52
−0.122
0.116
0.101
0.253
0.134
0.172
0.146
0.139
0.141
λ
62
−0.123
0.114
0.102
0.366
0.136
0.271
0.144
0.118
0.176
τ
12
−0.012
0.004
0.052
−0.151
0.036
0.126
0.005
0.002
0.052
τ
13
−0.018
0.002
0.056
−0.153
0.041
0.126
−0.013
0.003
0.060
τ
22
−0.001
0.003
0.046
−0.110
0.014
0.052
−0.035
0.004
0.045
τ
23
0.016
0.003
0.051
−0.216
0.050
0.060
−0.022
0.004
0.055
τ
32
−0.001
0.003
0.046
−0.113
0.016
0.051
−0.001
0.002
0.048
τ
33
−0.018
0.002
0.051
−0.186
0.038
0.059
−0.012
0.003
0.058Based on Table2, it can be found that the results obtained from our proposal are rather reasonable when compared to normal model, while partly exchangeable model gives serious biases. Moreover, we consider different inputs of superparameters in priors and find that the estimates are rather robust.


## 5. A Real Example

To illustrate the proposed procedure with a real example, a political-economic risk data set [43] was analyzed, which is adopted from Henisz’s [44] political constraint index data set (POLCON), Marshall et al. [45] state failure problem sets (PITF), and Alvarez et al.’s [46] ACLP Political and Economic Database (ACLP). The data set is formed by the two economic indicators and three political variables from 62 countries. The first index is the log black market premium (BMP). This is a continuous variable which is usually used as a proxy for illegal economic activity. The second index is log real gross domestic product (GDP). It is used to measure the productivity of a country. The third variable is a measure of independence of the national judiciary. This is a binary variable: it takes 1 if the judiciary is judged to be independent and 0 otherwise. The next measurement, measuring the level of lack of expropriation risk threat (LE), is an ordered categorical variable coded with 0, 1, 2, 3, 4, and 5. The last variable is an expert judgment of measuring lack of corruption (LC). It is also an ordered categorical variable scaled with 0 to 5. The total sample size is 62 and the frequencies of each category occurring are equal to {
34,28
}, {
2,6
,
7,19,14,14
}, and {
5,11,18,11,8
,
9
}, respectively. To unify scales of the continuous variables, the corresponding raw data were standardized.Lety
T = (log
⁡
B
M
P, log
⁡
G
D
P, IJ, LE, LC) be the vector of the observed variables. Based on the objective of this example, it is natural to group (i) {
log
⁡
B
M
P
,
log
⁡
G
D
P
} to an endogenous latent variable that can be interpreted as “economic factor, ξ” and (ii) {
I
J
,
L
E
,
L
C
} to an exogenous genotype latent variable that can be interpreted as “political factor, η.” Hence, the following loading matrix Λ in the measurement equation with ω
i
=
(
η
i
,
ξ
i
)
T is considered: (27)
Λ
T
=
0
∗
0
∗
1
∗
λ
41
λ
51
1
∗
λ
22
0
∗
0
∗
0
∗in which the ones and zeros are treated as known. Although other structures of Λ could be used, here we consider a nonoverlapped structure for clear interpretations of the latent variables: λ
j
k measures the effect of ω
k on the observed variable y
j. Since the third variable is binary and the last two variables are measured on a six-point scale with each involving six thresholds, for model identification, we fix ψ
ϵ
3
=
1 and endpoints of thresholds τ
31
,
τ
35
,
τ
41, and τ
45 at −1.8486, 0.7527, −1.4007, and 1.0574, respectively. These fixed threshold values were chosen via τ
j
k
=
Φ
-
1
(
p
^
j
k
), where p
^
j
k are observed marginal proportions of the categories with z
j
<
k.By primary data analysis, we find that the skewness and kurtosis of the first two variables are{
-
0.1340
,
-
0.4319
} and {
2.0892,2.0958
}, respectively. We also evaluate the predictive density function for continuous variables. Figure 2 gives the contours of posterior predictive density of pair (
y
1
,
y
2
) under parametric model and semiparametric model M
ϵ (see below) based on 60
×
60 grids. It can be seen that the data for pair (
y
1
,
y
2
) are heavy-tailed and the predictive density under semiparametric model captures the high frequency region successfully while parametric model fails. For model comparison, we consider the following competing models: (28)
M
μ
:
y
i
∗
=
μ
i
+
Λ
ω
i
+
ϵ
i
,
ϵ
i
~
N
0
,
Ψ
ϵ
,
ω
i
~
N
0
,
ϕ
;
M
ϵ
:
y
i
∗
=
μ
+
Λ
ω
i
+
ϵ
i
,
ϵ
i
~
N
0
,
Ψ
ϵ
i
,
ω
i
~
N
0
,
ϕ
;
M
μ
ϵ
:
y
i
∗
=
μ
i
+
Λ
ω
i
+
ϵ
i
,
ϵ
i
~
N
0
,
Ψ
ϵ
i
,
ω
i
~
N
0
,
ϕ
;
M
μ
ϕ
:
y
i
∗
=
μ
i
+
Λ
ω
i
+
ϵ
i
,
ϵ
i
~
N
0
,
Ψ
ϵ
,
ω
i
~
N
0
,
ϕ
i
;
M
ϵ
ϕ
:
y
i
∗
=
μ
+
Λ
ω
i
+
ϵ
i
,
ϵ
i
~
N
0
,
Ψ
ϵ
i
,
ω
i
~
N
0
,
ϕ
i
;
M
μ
ϵ
ϕ
:
y
i
∗
=
μ
i
+
Λ
ω
i
+
ϵ
i
,
ϵ
i
~
N
0
,
Ψ
ϵ
i
,
ω
i
~
N
0
,
ϕ
i
.The following two types of prior inputs are, respectively, used for the hyperparameters involved in the parametric components and semiparametric components: (I) Λ
ϵ
0
k
=
Λ
~
ϵ
k, H
ϵ
0
k
=
I
2, γ
0
=
γ
~, H
ζ
0
=
1, μ
0
=
μ
~, Σ
0
=
diag
⁡
{
S
}, α
ϵ
0
k
=
9.0, β
ϵ
0
k
=
(
α
ϵ
0
k
-
1
)
ψ
~
ϵ
k, ρ
0
=
ρ
0
ϕ
=
20, R
0
ϕ
-
1
=
(
ρ
0
-
2
)
ϕ
~, κ
1
=
κ
2
=
8.0, and τ
1
=
τ
2
=
8.0, where θ
~ denotes the maximum likelihood estimates of θ under parametric model from analysis of a “control-group” sample and S is the polychoric correlation matrix obtained on the basis of single confirmatory factor analysis model; (II) λ
ϵ
0
j
k
=
0, H
ϵ
0
k
=
0.01
I
2, γ
0
=
0, H
ζ
0
=
0.01, μ
0
=
0
5
×
1, Σ
0
=
0.01
I
5, α
ϵ
0
k
=
β
ϵ
0
k
=
2.0, ρ
0
=
ρ
0
ϕ
=
10, R
0
-
1
=
ρ
0
-
2, κ
1
=
κ
2
=
0.01, and τ
1
=
τ
2
=
0.01
. Note that prior (I) gives more information than prior (II) since it partly takes advantage of information from sample.Figure 2
Plot of contours of posterior predictive density of pair(
y
1
,
y
2
) under parametric model and semiparametric model M
ϵ: (a) corresponds to parametric model and (b) corresponds to semiparametric model with G
=
50.
(a)
(b)The proposed Bayesian semiparametric approach withG
=
50 was applied to calculate the values of CPO and LPML. We draw 100,000 effective observations from the corresponding posteriors via the blocked Gibbs sampler and divide them into 100 batches equally. Table 3 gives the means and standard deviations of LPML under priors (I) and (II). The following facts can be found. (i) The values of LPML under prior (I) are larger than those under prior (II). This indicates that the LPML tends to choose the model with informative prior. (ii) M
ϵ give the largest value. Among the posited models, M
ϵ is selected. We also compute the values of LPML for parametric model. They are -
287.4262 and -
288.6033 under priors (I) and (II) with standard deviations 6.287 and 5.065, respectively. Therefore, the data support the semiparametric model instead of parametric model.Table 3
Mean and standard deviation (SD) of LPML in the political and economic risk data.
Model
LPML
BAY I
BAY II
Mean
SD
Mean
SD
M
μ
−253.0577
60.7693
−274.9704
14.7176
M
ϵ
−176.5289
110.0227
−188.8361
175.3741
M
μ
ϵ
−193.2270
81.3209
−209.0063
77.0664
M
μ
ϕ
−235.7465
15.3588
−235.7465
11.6911
M
ϵ
ϕ
−267.4536
10.2084
−271.2875
10.4206
M
μ
ϵ
ϕ
−255.2335
27.2208
−258.7124
15.8416Table4 presents the estimates of factor loading λ
j
k as well as their standard deviations with semiparametric and parametric model under prior (I). The factor loading estimates λ
^
k
j in the measurement equation can be interpreted according to a standard confirmatory factor analysis model. The difference between the two approaches is obvious: the estimates of λ
41 and λ
51 under parametric model are only half of those under semiparametric model. Moreover, the standard deviations of estimates with parametric method are uniformly larger than that of semiparametric model. Since we identify illegal economic activity log
⁡
P
C
R with economic factor ξ (λ
12
∗
=
1) and independent of judiciary with political factor η (λ
31
∗
=
1), respectively, the level of economic factor has a negative effect on real gross domestic product, while the level of political factor has positive effect on lack of expropriation risk threat and lack of corruption. The estimate λ
^
22
=
-
0.123 indicates that a one-unit increase in the level of economic factor leads to 0.123-unit decrease in the magnitude of gross domestic product. The interpretation of λ
^
41 and λ
^
51 is similar. The differences of estimates between parametric and semiparametric methods illustrate the effects of heavy tails of the data on the estimates.Table 4
Estimates and standard errors estimates of the parameters in analysis of political and economic risk data.
Parameter
Parametric model
M
ϵ model
Est.
SD
Est.
SD
λ
22
−0.155
0.083
−0.123
0.077
λ
41
0.418
0.104
0.846
0.088
λ
51
0.367
0.090
0.754
0.066
τ
42
−1.336
0.157
−1.340
0.055
τ
43
−0.905
0.167
−0.898
0.061
τ
44
−0.008
0.149
0.004
0.055
τ
52
−0.794
0.142
−0.787
0.046
τ
53
0.001
0.155
0.002
0.053
τ
54
0.566
0.136
0.567
0.035


## 6. Concluding Remarks

Parametric modeling for latent variable model with mixed data structure has long dominated Bayesian inference work, typically developed within the standard exponential family. Such modeling is often confused with handling the multimodal and unknown heterogeneous problems. In dealing with multimodality or increased heterogeneity in data, one naturally resorts to the finite mixture model [47, 48] which is more flexible and feasible to implement due to advances in simulation-based model fitting.Rather than handling the large number of parameters resulting from the finite mixture models with a large number of components, we consider, in this paper, the finite-dimensional Dirichlet process mixture model for latent variable model with continuous and ordinal responses. The core of our proposal is to model the mean vector and/or variance-covariance parameters of unique errors and latent variables into the finite-dimensional stick-breaking priors. This will help to reveal the local dependence structure such as classification groups and clustering among the data. The blocked Gibbs sampler developed by Ishwaran and Zarepour [31], which takes advantage of the block updating and accelerates mixing in Gibbs sampling, is adapted here to cope with the posterior inference.The proposed methodologies in this paper can be applied to more general latent variable models that include the multilevel SEMs [49] and longitudinal latent trait models [5] with discrete variables.



---

*Source: 1648462-2016-03-30.xml*