# Bayesian Analysis for Dynamic Generalized Linear Latent Model with Application to Tree Survival Rate

Authors: Yu-sheng Cheng; Mei-wen Ding; Ye-mao Xia; Wen-fa Zhan
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Logistic regression model is the most popular regression technique, available for modeling categorical data especially for dichotomous variables. Classic logistic regression model is typically used to interpret relationship between response variables and explanatory variables. However, in real applications, most data sets are collected in follow-up, which leads to the temporal correlation among the data. In order to characterize the different variables correlations, a new method about the latent variables is introduced in this study. At the same time, the latent variables about AR (1) model are used to depict time dependence. In the framework of Bayesian analysis, parameters estimates and statistical inferences are carried out via Gibbs sampler with Metropolis-Hastings (MH) algorithm. Model comparison, based on the Bayes factor, and forecasting/smoothing of the survival rate of the tree are established. A simulation study is conducted to assess the performance of the proposed method and a pika data set is analyzed to illustrate the real application. Since Bayes factor approaches vary significantly, efficiency tests have been performed in order to decide which solution provides a better tool for the analysis of real relational data sets.


---

## Body

## 1. Introduction

Logistic regression model, a widely appreciated model for analyzing categorical response data especially for binary/dichotomous data in social science research, marketing, biomedical studies, and ecology, has been received a substantive attention [1, 2]. The primary concern of logistic regression analysis is to build the interrelationships between explanatory variables and responses and explain the variability of the response probability in terms of the variability in the observed covariates. The common analysis for the logistic regression model is based on the independent subjects, and statistical inference is carried out via the maximum likelihood approach (see [3–5]).However, in practice, many studies need to account for correlations among the items within time as well as correlations among the same items between times. Independent assumptions cannot capture such correlations. In order to interpret the correlations among the responses, one popular choice is to introduce the latent variables. Even more, taking into account the temporal correlations, a dynamic system is established which leads to the latent variable models [6–8]. However, this dynamic system is associated with dynamic factor. So, a Bayesian approach is proposed in this paper. With an alternative approach about the ML, the development is based on the generalized logistic model in the well-known linear dynamic factor analysis model. To cope with this complicated model, Gibbs sampler [9] with MH algorithm [10, 11] is implemented to generate a sequence of observations from the appropriate joint posterior distribution. The joint Bayes estimates of the parameters and latent variables scores associated with standard error estimates can be obtained directly based on the simulated observations. Beyond the estimation problem, another main objective of this paper is to introduce the Bayes factor (see [12]) to test various hypotheses about the posited model. In general, the computation of the Bayes factor involving multiple integrals is rather intractable and, more often, are evaluated via various computational methods (see [13–19], and among others). Among easy-to-construct, we adopted the well-known path sampling technique [20] for computing Bayes factor. The nice feature of path sampling lies in its simplicity in implementation and effectiveness over such an important sampling (e.g., [21]) and bridge sampling [22].A basic intent of modeling dynamic model is to represent series of observations generated in time and to predict the future evolution of the series. The well-known Kalman recursion [23] is commonly employed for forecasting and provides the optimal forecasts in the case where the state transitions are linear and dynamic system is Gaussian. However, if Gaussian assumptions, in particular, of the distributions of the measurement errors are violated, Kalman recursion only yields the best linear predictor for linear dynamic system which is substantially different from the optimal forecasts. Many authors have suggested various alternatives to the Kalman filter; see, among others, Kitagawa [24], Meinhold and Singpurwalla [25], Carlin et al. [26], and Smith and West [27]. Motivated by the key idea in Carlin et al. [26], in this paper, the missing data is treated as the future observations and augmented them with parameters of latent variables as observed data. Therefore, it can use Bayesian analysis approach to deal them. Simulated observations generated by Gibbs sampler from the joint posterior distribution can be directly applied to forecasts and/or smoothing.This paper is organized as follows. In Section2, we describe the problem under consideration. Section 3 gives the estimation procedure and explores the Bayesian approach in detail. For model selection, Bayes factor and Bayesian forecasting are introduced in Section 4. Simulation study and a real example are given in Section 5. Some concluding remarks are presented in Section 6.


## 2. Materials and Methods about Model

In this section, the problem and model structure are presented as follows, which will be considered throughout the paper. Fori
=
1
,
…
,
T, j
=
1
,
…
,
K, let y
i
j denote the number of living tree within the jth plot in the ith month. We assume that y
i
j follows the binomial distribution with the unknown survival rate p
i
j
(1)
[
y
i
j
∣
p
i
j
]
~
B
i
n
(
n
i
j
,
p
i
j
)
,
where n
i
j is the total number of tree being investigated at time i in thejth plot. Further, the probability p
i
j is in relation to the covariates via the following link function:
(2)
log
⁡
i
t
(
p
i
j
)
=
x
i
j
T
β
+
λ
j
ξ
i
+
ε
i
j
,
in which “log
⁡
i
t” denotes the log
⁡
i
t function, x
i
j
=
(
1
,
x
i
j
1
,
…
,
x
i
j
r
)
T is a (
r
+
1
) vector of fixed covariates, β is a (
r
+
1
) regression parametric vector, ξ
i is the common factor, and ε
i
=
(
ε
i
1
,
…
,
ε
i
K
)
′ is the unique error with distribution N
K
(
0
,
Ψ
ε
) where Ψ
ε
=
diag
⁡
{
ψ
ε
1
,
…
,
ψ
ε
p
} is the diagonal matrix. The factor loading vector is Λ
=
(
λ
1
,
…
,
λ
K
)
′, which measures the effect of ξ
i on the functions of the response probability. Moreover, it assumes that ξ
i satisfies the following one-order autoregressive model
(3)
ξ
i
=
ϕ
ξ
i
-
1
+
δ
i
,
where |
ϕ
|
<
1 and δ
i is the iid random errors distributed according to the normal distribution N
(
0
,
σ
2
) (σ
2
>
0). For the initial variable ξ
0, we assume that ξ
0
~
N
(
μ
0
,
σ
0
2
) with mean μ
0 and variance σ
0
2. δ
i and ε
i
j are mutually independent and both are independent of ξ
i. In this paper, we treat μ
0 and σ
0
2 to be known and fix them at some preassigned values.Letϑ
i
j
=
x
i
j
T
β
+
λ
j
ξ
i
+
ε
i
j, b
(
ϑ
i
j
)
=
log
⁡
(
1
+
exp
⁡
(
ϑ
i
j
)
). Then, (1), (2), and (3) can be reformulated as follows:
(4)
p
(
y
i
j
∣
ϑ
i
j
)
=
C
n
i
j
y
i
j
exp
⁡
{
y
i
j
ϑ
i
j
-
n
i
j
b
(
ϑ
i
j
)
}
,
ϑ
i
j
=
x
i
j
T
β
+
λ
j
ξ
i
+
ε
i
j
,
ε
i
~
N
(
0
,
Ψ
ε
)
,
ξ
i
=
ϕ
ξ
i
-
1
+
δ
i
,
ξ
0
~
N
(
μ
0
,
σ
0
2
)
,
δ
i
~
N
(
0
,
σ
2
)
,
δ
i
⊥
ε
i
,
ξ
0
⊥
ε
i
,
δ
i
⊥
ξ
0
.
The main aim of introducing the latent variables is to characterize the temporal correlation between the consecutive observations and explain the interrelationships among the response probability. To give more details, note that for i
>
i
′, and j
,
l
=
1
,
…
,
K, the correlation between ϑ
i
′
j and ϑ
i
l is given by
(5)
Corr
(
ϑ
i
′
j
,
ϑ
i
l
)
=
λ
j
λ
l
σ
ξ
2
(
i
)
ϕ
i
-
i
′
λ
j
σ
ξ
2
+
ψ
ε
j
λ
l
2
σ
ξ
2
+
ψ
ε
l
,
in which
(6)
σ
ξ
2
(
i
)
=
σ
0
2
ϕ
2
i
+
σ
2
1
-
ϕ
2
i
1
-
ϕ
2
.
A path diagram of latent variables and manifest variables at time i
-
1 and i is given in Figure 1. Following the conventions of path diagrams, the ellipses represent the latent variables, the rectangles denote the observed measurements, and the arrows represent the direct effect.Figure 1
Path diagram about interaction with the latent variables and manifest variables.For ease of exposition, lety
i
=
(
y
i
1
,
…
,
y
i
K
)
′, ϑ
i
=
(
ϑ
i
1
,
…
,
ϑ
i
K
)
′, and X
i
=
(
x
i
1
,
…
,
x
i
K
)
′, which denote the vector of random observations, the vector of latent variables, and the vector of covariates cross K items at time i, respectively. Let Y
=
{
y
i
:
i
=
1
,
…
,
T
} be the collection of the observed data, let Ω
=
{
ξ
i
:
i
=
0
,
…
,
T
} be the collection of the factor variables, let Ξ
=
{
ϑ
i
:
i
=
1
,
…
,
T
} be the collection of the latent variables, and let θ be the vector of unknown parameters in β, Λ, Ψ
ε, ϕ, and σ
2 involved in the model and then the joint distribution of (Y
,
Ξ
,
Ω) conditioning on θ is given by
(7)
p
(
Y
,
Ξ
,
Ω
θ
)
=
p
(
Y
∣
Ξ
,
θ
)
p
(
Ξ
∣
Ω
,
θ
)
p
(
Ω
∣
θ
)
=
(
∏
i
=
1
T
∏
j
=
1
K
p
(
y
i
j
∣
ϑ
i
j
)
p
(
ϑ
i
j
∣
ξ
i
,
θ
)
)
×
p
(
ξ
0
)
∏
i
=
1
T
p
(
ξ
i
∣
θ
)
.
In the framework of Bayesian analysis, the following priors for unknown parameter vector θ are used to complete Bayesian specifications: p
(
θ
)
=
p
(
β
)
p
(
Λ
,
Ψ
ε
)
p
(
ϕ
)
p
(
σ
2
) and
(8)
β
~
N
r
+
1
(
β
0
,
Σ
0
)
,
ϕ
~
Uniform
(
-
1
,
1
)
,
σ
-
2
~
Gamma
(
α
0
,
τ
0
)
,
(
Λ
∣
Ψ
ε
)
~
∏
k
=
1
K
N
(
Λ
0
k
,
ψ
ε
k
H
0
k
)
,
Ψ
ε
~
∏
k
=
1
K
Gamma
(
α
ε
0
k
,
β
ε
0
k
)
,
in which β
0, Σ
0, Λ
0
k, H
0
k, α
0, τ
0, α
ε
0
k, and β
ε
0
k are the known superparameters.


## 3. Parameters Estimation via Gibbs Sampler with MH Algorithm

For the Bayesian analysis, we are required to evaluate the complicated posterior distributionp
(
θ
∣
Y
) which involves the high-dimensional integrals. Data augmentation [28] technique is used to cope with the posterior analysis in relation to the complicated p
(
θ
∣
Y
). Specifically, the observed data Y are augmented with the latent quantities {
Ξ
,
Ω
} in the posterior analysis. A sequence of random observations will be generated by the Gibbs sampler [9], coupled with the Metropolis-Hastings algorithm [10, 11] from the joint posterior distribution p
(
Ξ
,
Ω
,
θ
∣
Y
), specifically, at the mth iteration with current values (
Ξ
(
m
)
,
Ω
(
m
)
,
θ
(
m
)
)
=
(
Ξ
(
m
)
,
Ω
(
m
)
,
β
(
m
)
,
Λ
(
m
)
,
Ψ
ε
(
m
)
,
ϕ
(
m
)
,
σ
2
(
m
)
). We do the following:(i)
drawΞ
(
m
+
1
) from p
(
Ξ
∣
Ω
(
m
)
,
β
(
m
)
,
Λ
(
m
)
,
Ψ
ε
(
m
)
,
ϕ
(
m
)
,
σ
2
(
m
)
,
Y
);
(ii)
drawΩ
(
m
+
1
) from p
(
Ω
∣
Ξ
(
m
+
1
)
,
β
(
m
)
,
Λ
(
m
)
,
Ψ
ε
(
m
)
,
ϕ
(
m
)
,
σ
2
(
m
)
,
Y
);
(iii)
drawβ
(
m
+
1
) from p
(
β
∣
Ξ
(
m
+
1
)
,
Ω
(
m
+
1
)
,
Λ
(
m
)
,
Ψ
ε
(
m
)
,
ϕ
(
m
)
,
σ
2
(
m
)
,
Y
);
(iv)
draw(
Λ
(
m
+
1
)
,
Ψ
ε
(
m
+
1
)
) from p
(
Λ
,
Ψ
ε
∣
Ξ
(
m
+
1
)
,
Ω
(
m
+
1
)
,
β
(
m
+
1
)
,
ϕ
(
m
)
,
σ
2
(
m
)
,
Y
);
(v)
drawϕ
(
m
+
1
) from p
(
ϕ
∣
Ξ
(
m
+
1
)
,
Ω
(
m
+
1
)
,
β
(
m
+
1
)
,
Λ
(
m
+
1
)
,
Ψ
ε
(
m
+
1
)
,
σ
2
(
m
)
,
Y
);
(vi)
drawσ
2
(
m
+
1
) from p
(
σ
2
∣
Ξ
(
m
+
1
)
,
Ω
(
m
+
1
)
,
β
(
m
+
1
)
,
Λ
(
m
+
1
)
,
Ψ
ε
(
m
+
1
)
,
ϕ
(
m
+
1
)
,
Y
).It has been shown [9, 29] that under mild conditions and, for sufficiently large m, such as J, the joint distribution of (
Ξ
(
m
)
,
Ω
(
m
)
,
θ
(
m
)
) converges at an exponential rate to the desired posterior distribution p
(
Ξ
,
Ω
,
θ
∣
Y
). Hence, p
(
Ξ
,
Ω
,
θ
∣
Y
) can be approximated by the empirical distribution of {
(
Ξ
(
m
)
,
Ω
(
m
)
,
θ
(
m
)
)
:
m
=
J
+
1
,
…
,
J
+
M
}, where M is chosen to give sufficient precision to the empirical distribution. The convergence of the hybrid Gibbs sampler with MH algorithm can be monitored by the “estimated potential scale reduction (EPSR)” values as suggested by Gelman and Rubin [30]. To implement Gibbs sampler with MH algorithm, it requires the full conditional distributions for each parameters and latent quantities. For saving space, these technical details are omitted.Simulated observations obtained from the posterior values can be used for statistical inference via straightforward analysis procedure. For brevity, let{
(
Ξ
(
m
)
,
Ω
(
m
)
,
θ
(
m
)
)
:
m
=
1
,
…
,
M
} be the random observations of (
Ξ
,
Ω
,
θ
) generated by the Gibbs sampler from p
(
Ξ
,
Ω
,
θ
∣
Y
). The joint Bayesian estimate of θ and Ω can be obtained easily via the corresponding sample means of the generated observations as follows:
(9)
θ
^
=
∑
m
=
1
M
θ
(
m
)
M
,
Ω
^
=
∑
m
=
1
M
Ω
(
m
)
M
.
The consistent estimates of covariance matrix of (Ω
,
θ) can be obtained as follows:
(10)
Cov
(
θ
^
)
≈
∑
m
=
1
M
(
θ
(
m
)
-
θ
^
)
(
θ
(
m
)
-
θ
^
)
′
(
M
-
1
)
,
Cov
(
Ω
^
)
≈
∑
m
=
1
M
(
Ω
(
m
)
-
Ω
^
)
(
Ω
(
m
)
-
Ω
^
)
′
(
M
-
1
)
.


## 4. Model Selection and Bayesian Forecasting

Model selection is an important issue for Bayesian generalized logistic analysis since it is of interest to determine whether the latent variables are involved or whether the significant effects exist among the competing covariates. Consider the problem of comparing competing modelsM
1 and M
0 which are nested or nonnested. Let p
(
Y
∣
M
1
) and p
(
Y
∣
M
0
) denote the marginal density of the data Y under M
1 and M
0, respectively. A popular choice for selecting models (e.g., [16, 18]) is achieved via the following Bayes factor:
(11)
BF
10
=
p
(
Y
∣
M
1
)
p
(
Y
∣
M
0
)
=
∫
p
(
Y
,
Ξ
,
Ω
∣
θ
,
M
1
)
p
(
θ
∣
M
1
)
d
θ
d
Ξ
d
Ω
∫
p
(
Y
,
Ξ
,
Ω
∣
θ
,
M
0
)
p
(
θ
∣
M
0
)
d
θ
d
Ξ
d
Ω
,
where p
(
θ
∣
M
k
) and p
(
Y
,
Ξ
,
Ω
∣
θ
,
M
k
) are the prior density of θ and the joint probability density of (
Y
,
Ξ
,
Ω
) given the values of θ under M
k, respectively. Kass and Rafter [18] provide the following categories to furnish appropriate guidelines in Table 1.Table 1
The evidence about Bayes factor againstM
0 [18].
logBF10
BF10
Evidence againstM
0
0 to 1
1 to 3
Not worth more than a bare mention
1 to 3
3 to 20
Positive
3 to 5
20 to 150
Strong
>5
>150
Very strongComputingBF
10 or log BF
10 involves the high-dimensional integrations of likelihood with respect to (
Ξ
,
Ω
,
θ
) which is rather difficult. Various techniques are explored to address this problem (see [7]). Among easy-to-construct, we consider the path sampling method [20]. The core of path sampling is to construct a series of linked models which link up the competing models.Specifically, consider a class of densities indexed by a continuous parametert in [
0
,
1
] such that t
=
0 corresponding to M
0 and t
=
1 to M
1. Note that
(12)
p
(
Ξ
,
Ω
,
θ
∣
Y
,
t
)
=
1
m
(
t
)
p
(
Y
,
Ξ
,
Ω
,
θ
∣
t
)
,
where m
(
t
) is the normalizing constant of p
(
Y
,
Ξ
,
Ω
,
θ
∣
t
) given by
(13)
m
(
t
)
=
m
(
Y
∣
t
)
=
∫
p
(
Y
,
Ξ
,
Ω
,
θ
∣
t
)
d
Ξ
d
Ω
d
θ
.
A path using the parameter t in [
0
,
1
] is constructed so that BF
10
=
m
(
1
)
/
m
(
0
). Taking logarithm and differentiating (13) with respect to t give
(14)
d
log
⁡
m
(
t
)
d
t
=
∫
d
d
t
log
⁡
p
(
Y
,
Ξ
,
Ω
,
θ
∣
t
)
p
(
Ξ
,
Ω
,
θ
∣
Y
,
t
)
d
Ξ
d
Ω
d
θ
=
E
Ξ
,
Ω
,
θ
U
(
Ξ
,
Ω
,
θ
,
Y
,
t
)
,
in which E
Ξ
,
Ω
,
θ
[
U
(
Ξ
,
Ω
,
θ
,
Y
,
t
)
] denotes the expectation with respect to the distribution p
(
Ξ
,
Ω
,
θ
∣
Y
,
t
) and
(15)
U
(
Ξ
,
Ω
,
θ
,
Y
,
t
)
=
d
log
⁡
p
(
Y
,
Ξ
,
Ω
,
θ
∣
t
)
d
t
.
Hence,
(16)
log
⁡
BF
10
=
log
⁡
m
(
1
)
m
(
0
)
=
∫
0
1
E
Ξ
,
Ω
,
θ
U
(
Ξ
,
Ω
,
θ
,
Y
,
t
)
d
t
≅
1
2
∑
s
=
0
S
(
U
-
s
+
1
+
U
-
s
)
Δ
t
s
,
where {
t
s
:
s
=
0
,
…
,
S
+
1
} are the fixed grids such that t
0
=
0
<
t
1
<
⋯
<
t
s
+
1
=
1.0, Δ
t
s
=
t
s
+
1
-
t
s, and
(17)
U
-
s
=
1
M
∑
m
=
1
M
U
(
Ξ
(
m
)
,
Ω
(
m
)
,
θ
(
m
)
,
Y
,
t
s
)
,
in which {
(
Ξ
(
m
)
,
Ω
(
m
)
,
θ
(
m
)
)
:
m
=
1
,
…
,
M
} are observations sampled from p
(
Ξ
,
Ω
,
θ
∣
Y
,
t
). The MCMC method proposed for estimation produced can be directly applied to simulate the above observations for computing BF
10 or 2log BF
10.With all the complete conditionals available for sampling, a predictive or estimated values for a futurey
T
+
1 can be obtained, provided that X
T
+
1 was available, which offers a solution to the so-called forecasting problem. More generally, the predictions of the future observations Y
*
=
{
y
T
+
1
,
…
,
y
T
+
t
0
} can be obtained in terms of the observed data Y and the covariates {
X
1
,
…
,
X
T
+
t
0
}. As mentioned in introduction, we treat Y
* as missing data and augment them with Ω
*
=
{
ξ
0
,
…
,
ξ
T
+
t
0
}, and Ξ
*
=
{
ϑ
1
,
…
,
θ
T
+
t
0
} and data Y in the posterior analysis. Gibbs sampler now is implemented to sample observations from the posterior p
(
Y
*
,
Ξ
*
,
Ω
*
,
θ
∣
Y
). Predictive values for Y
*are given by Y
*
=
E
[
Y
*
∣
Y
] which can be solved by
(18)
Y
^
*
=
E
[
Y
*
∣
Y
]
=
∑
m
=
1
M
Y
*
(
m
)
M
,
where {
Y
*
(
m
)
:
m
=
1
,
…
,
M
} are the random observations generated by the Gibbs sampler with MH algorithm from p
(
Y
*
,
Ξ
*
,
Ω
*
,
θ
∣
Y
). Consistent estimates of covariance matrix of (18) can be obtained via sample covariance matrix.


## 5. Results of Experiment

### 5.1. Simulation Study

A simulation study is presented to assess the performance of our proposed procedure. The data set is simulated from the model defined by (1), (2), and (3) with T
=
100, K
=
4, and r
=
3. The true population values of the unknown parameters are taken as β
=
(
-
0.6
,
-
0.7
,
0.6
)
′, Λ
=
(
1.0
,
1.2
,
1.0
,
1.8
)
′, ϕ
=
0.5, σ
2
=
1.0, μ
0
*
=
0.0, and σ
0
2
*
=
1.0 in which elements with asterisk are treated as fixed known parameters for model identification (see [6]). For i
=
1
,
…
,
100, j
=
1
,
…
,
4, we first sampled independently the observations from the normal distribution N
(
0
,
I
3
) to form the covariates x
i
j, where I
k is a k
×
k identity matrix, and then generate ξ
i from model (3) with ξ
0
=
1.0. Based on these settings, we generate the observations Y
=
{
y
i
j
:
i
=
1
,
…
,
100, j
=
1
,
…
,
4
} from the binomial distribution Bin (
50
,
p
i
j
) with n
i
j
=
50. Priors inputs for the hyperparameters in the prior distributions are taken as (I) β
0
=
0, Σ
0
=
1000
I
3, Λ
0
k
=
0, H
0
k
=
1000, α
0
=
2.0, τ
0
=
2.0, α
ε
0
k
=
4.0, and β
ε
0
k
=
4.0; (II) β
0 and Λ
0
k are set to be equal to the true values, respectively; Σ
0
=
I
3, H
0
k
=
1.0, α
0
=
4.0, τ
0
=
4.0, α
ε
0
k
=
9.0, and β
ε
0
k
=
8.0. Note that prior (I) approaches the noninformative prior and prior (II) gives the informative prior.At first, four parallel methods are conducted for different starting values as a pilot study to obtain some idea about the number of the Gibbs sampler iterations in getting convergence. For MH algorithm, the acceptance rate is also adjusted in posterior sampling to give about 40% [31]. Hence, in the following analysis, experiment always takes such values for the tuned parameters. Figure 2 gives the estimated potential scaled reduction (EPSR) values [30] of the parameters against the number of the iterations in Gibbs sampler under prior (I).Figure 2
EPSR values change with numbers of iterations under prior (I).In Figure2, it can be seen clearly that all these EPSR values are less than 1.2 in about 2000 iterations. To be conservative, we collect 5000 observations after 5000 “burn-in” in computing the absolute bias (BIAS) and the root mean squares (RMS) of the estimates and the true values in 100 replications. Results obtained from the Gibbs sampler are summarized in Table 2, where SD denotes the standard deviation of the estimates. Based on these results, it can be seen that our proposal is rather accurate and effective. It also shows that there are no significant differences of estimates between prior (I) and prior (II). Therefore, the proposal method is more robust against the choices of values of hyperparameters. An interesting phenomenon is that the estimates of variance parameters in unique errors under prior (I) are more accurate than those under prior (II). This reflects that it provides more information in estimation.Table 2
Statistical results of Gibbs under prior (I) and prior (II) with relevant parameters.
Para.
Priors
(I)
(II)
BIAS
RMS
SD
BIAS
RMS
SD
β
0
0.0127
0.0456
0.0508
0.0228
0.0510
0.0467
β
1
0.0267
0.0526
0.0469
0.0231
0.0533
0.0431
β
2
-
0.0202
0.0472
0.0490
-
0.0246
0.0455
0.0451
λ
2
0.0510
0.0943
0.1114
0.0160
0.0839
0.0964
λ
3
0.0294
0.0971
0.1087
0.0154
0.0880
0.0961
λ
4
0.0178
0.0739
0.0822
-
0.0063
0.0545
0.0714
ψ
ε
1
0.1617
0.1701
0.1017
0.0573
0.0923
0.0953
ψ
ε
2
0.1845
0.1928
0.1110
0.0591
0.0821
0.1015
ψ
ε
3
0.1791
0.1869
0.1096
0.0595
0.0851
0.1020
ψ
ε
4
0.1313
0.1412
0.0909
0.0255
0.0569
0.0829
ϕ
-
0.0115
0.0694
0.0664
-
0.0110
0.0692
0.0662
σ
2
-
0.1013
0.1944
0.1905
-
0.0345
0.1464
0.1725For model comparison, we consider the following competing models to illustrate the performance of the Bayes factor.M
1 is the aforementioned model given by (4). Consider
(19)
M
0


:
p
(
y
i
j
∣
ϑ
i
j
)
=
C
n
i
j
y
i
j
exp
⁡
⁡
{
y
i
j
ϑ
i
j
-
n
i
j
b
(
ϑ
i
j
)
}
,
ϑ
i
j
=
x
i
j
T
β
+
λ
j
ξ
i
,
ξ
i
,
iid
N
(
0
,
σ
2
)
.
Note that M
0 is nested to M
1 in the sense that all the parameters in model M
0 are included in model M
1. M
0 indicates that no temporal correlation exists among the random effects. The linked model between M
1 and M
0 is taken as
(20)
M
t


:
p
(
y
i
j
∣
ϑ
i
j
)
=
C
n
i
j
y
i
j
exp
⁡
⁡
{
y
i
j
ϑ
i
j
-
n
i
j
b
(
ϑ
i
j
)
}
,
ϑ
i
j
=
x
i
j
T
β
+
λ
j
ξ
i
,
ξ
i
=
t
ϕ
ξ
i
-
1
+
δ
i
,
in which t
=
1 and t
=
0 correspond to M
1 and M
0, respectively. For computation, we choose S
=
20 in [
0
,
1
] and draw 5000 observations after 5000 burn-ins deleted from the posterior distribution for each t
s to calculate the log BF
10. Figure 3 gives the histogram of the values of log BF
10 across 100 replications. The estimate of log BF
10 is 11.2222 with standard deviation 6.6242. Based on the guidelines given by Kass and Raftery [18], M
1 is selected with positive evidence.Figure 3
Histogram of logBF
10 of posterior distribution for each t
s.To investigate the smoothing and forecasting of our proposal, we regenerate a data set with sample size 100 and treat the last 5 × 4 observations as unobserved. The Gibbs sampler with MH algorithm is used to compute the first ninety-five smoothing values and the last five forecasting values of the latent variablesξ
s. Figure 4 gives the true values and smoothing/forecasting values based on 5000 simulated observations after deleting first 5000 observations. It can be seen that it fits well between the estimated values and true values.Figure 4
Comparison with true values (dashed line) and forecasting values (pointed line).


### 5.2. Pika’s Data

An application is considered in relation to a study of Guo et al., 2009 [32], to illustrate the developed methodology. The data set is collected from Qinshui Forest farm of Jincheng Coal Industry Group sited on the northwest of Qinshui County of Jincheng City in Shanxi Province, 112°19′40′′–112°19′52 east longitude and 35°48′43′′–35°48′57 north latitude. Details about experiment design can be found in Guo et al., 2009 [32]. The original data set is constituted of the number of the living tree and pika within four plots in five months. The primary concern is to assess the relationship between the survival rate of young plantation and the pika population and model the optimal pika population per hectare. The data set is reanalyzed by Xia et al. [33] to investigate the correlations among the observed responses. In these studies the latent variables model is established through AR(
1
) model (3) coupled with model (1) and (2), which interprets the correlations among the outcomes resulting from the time dependence.Firstly, the following hypothesis is considered:H
0


:
ϕ
=
0. If H
0 is true, the model becomes the common generalized logistic model considered by Xia et al. [33]. Path sampling is used to compute log BF
10. The linked model is defined similarly as (20). S
=
20 grids were chosen in [
0
,
1
], and, for each t
s (s
=
1
,
…
,
20), 10000 random observations were drawn from the each posterior distributions and the first 5000 observations were deleted in view of the burn-in phase. The logarithm of Bayes factor is 3.0143. It shows that there exists significant evidence against H
0; hence, our proposal is appropriate.To assess the effect of forecasts of the proposal methodology, the original data set is divided into two parts: one is formed by the previous sixteen observations in four months within four plots and the other is made up of four observations in the following fifth month. The proposal algorithm computes the predictive values of the proportion of the living tree in the fifth month based on the previous 20 observations. Figure5 shows the forecasts and actual values of the proportions of the living tree within four plots. It can be seen that fitted values match the actual values closely.Figure 5
Comparison with the predictions for survival rate (solid line) and the actual values (dashed line).


## 5.1. Simulation Study

A simulation study is presented to assess the performance of our proposed procedure. The data set is simulated from the model defined by (1), (2), and (3) with T
=
100, K
=
4, and r
=
3. The true population values of the unknown parameters are taken as β
=
(
-
0.6
,
-
0.7
,
0.6
)
′, Λ
=
(
1.0
,
1.2
,
1.0
,
1.8
)
′, ϕ
=
0.5, σ
2
=
1.0, μ
0
*
=
0.0, and σ
0
2
*
=
1.0 in which elements with asterisk are treated as fixed known parameters for model identification (see [6]). For i
=
1
,
…
,
100, j
=
1
,
…
,
4, we first sampled independently the observations from the normal distribution N
(
0
,
I
3
) to form the covariates x
i
j, where I
k is a k
×
k identity matrix, and then generate ξ
i from model (3) with ξ
0
=
1.0. Based on these settings, we generate the observations Y
=
{
y
i
j
:
i
=
1
,
…
,
100, j
=
1
,
…
,
4
} from the binomial distribution Bin (
50
,
p
i
j
) with n
i
j
=
50. Priors inputs for the hyperparameters in the prior distributions are taken as (I) β
0
=
0, Σ
0
=
1000
I
3, Λ
0
k
=
0, H
0
k
=
1000, α
0
=
2.0, τ
0
=
2.0, α
ε
0
k
=
4.0, and β
ε
0
k
=
4.0; (II) β
0 and Λ
0
k are set to be equal to the true values, respectively; Σ
0
=
I
3, H
0
k
=
1.0, α
0
=
4.0, τ
0
=
4.0, α
ε
0
k
=
9.0, and β
ε
0
k
=
8.0. Note that prior (I) approaches the noninformative prior and prior (II) gives the informative prior.At first, four parallel methods are conducted for different starting values as a pilot study to obtain some idea about the number of the Gibbs sampler iterations in getting convergence. For MH algorithm, the acceptance rate is also adjusted in posterior sampling to give about 40% [31]. Hence, in the following analysis, experiment always takes such values for the tuned parameters. Figure 2 gives the estimated potential scaled reduction (EPSR) values [30] of the parameters against the number of the iterations in Gibbs sampler under prior (I).Figure 2
EPSR values change with numbers of iterations under prior (I).In Figure2, it can be seen clearly that all these EPSR values are less than 1.2 in about 2000 iterations. To be conservative, we collect 5000 observations after 5000 “burn-in” in computing the absolute bias (BIAS) and the root mean squares (RMS) of the estimates and the true values in 100 replications. Results obtained from the Gibbs sampler are summarized in Table 2, where SD denotes the standard deviation of the estimates. Based on these results, it can be seen that our proposal is rather accurate and effective. It also shows that there are no significant differences of estimates between prior (I) and prior (II). Therefore, the proposal method is more robust against the choices of values of hyperparameters. An interesting phenomenon is that the estimates of variance parameters in unique errors under prior (I) are more accurate than those under prior (II). This reflects that it provides more information in estimation.Table 2
Statistical results of Gibbs under prior (I) and prior (II) with relevant parameters.
Para.
Priors
(I)
(II)
BIAS
RMS
SD
BIAS
RMS
SD
β
0
0.0127
0.0456
0.0508
0.0228
0.0510
0.0467
β
1
0.0267
0.0526
0.0469
0.0231
0.0533
0.0431
β
2
-
0.0202
0.0472
0.0490
-
0.0246
0.0455
0.0451
λ
2
0.0510
0.0943
0.1114
0.0160
0.0839
0.0964
λ
3
0.0294
0.0971
0.1087
0.0154
0.0880
0.0961
λ
4
0.0178
0.0739
0.0822
-
0.0063
0.0545
0.0714
ψ
ε
1
0.1617
0.1701
0.1017
0.0573
0.0923
0.0953
ψ
ε
2
0.1845
0.1928
0.1110
0.0591
0.0821
0.1015
ψ
ε
3
0.1791
0.1869
0.1096
0.0595
0.0851
0.1020
ψ
ε
4
0.1313
0.1412
0.0909
0.0255
0.0569
0.0829
ϕ
-
0.0115
0.0694
0.0664
-
0.0110
0.0692
0.0662
σ
2
-
0.1013
0.1944
0.1905
-
0.0345
0.1464
0.1725For model comparison, we consider the following competing models to illustrate the performance of the Bayes factor.M
1 is the aforementioned model given by (4). Consider
(19)
M
0


:
p
(
y
i
j
∣
ϑ
i
j
)
=
C
n
i
j
y
i
j
exp
⁡
⁡
{
y
i
j
ϑ
i
j
-
n
i
j
b
(
ϑ
i
j
)
}
,
ϑ
i
j
=
x
i
j
T
β
+
λ
j
ξ
i
,
ξ
i
,
iid
N
(
0
,
σ
2
)
.
Note that M
0 is nested to M
1 in the sense that all the parameters in model M
0 are included in model M
1. M
0 indicates that no temporal correlation exists among the random effects. The linked model between M
1 and M
0 is taken as
(20)
M
t


:
p
(
y
i
j
∣
ϑ
i
j
)
=
C
n
i
j
y
i
j
exp
⁡
⁡
{
y
i
j
ϑ
i
j
-
n
i
j
b
(
ϑ
i
j
)
}
,
ϑ
i
j
=
x
i
j
T
β
+
λ
j
ξ
i
,
ξ
i
=
t
ϕ
ξ
i
-
1
+
δ
i
,
in which t
=
1 and t
=
0 correspond to M
1 and M
0, respectively. For computation, we choose S
=
20 in [
0
,
1
] and draw 5000 observations after 5000 burn-ins deleted from the posterior distribution for each t
s to calculate the log BF
10. Figure 3 gives the histogram of the values of log BF
10 across 100 replications. The estimate of log BF
10 is 11.2222 with standard deviation 6.6242. Based on the guidelines given by Kass and Raftery [18], M
1 is selected with positive evidence.Figure 3
Histogram of logBF
10 of posterior distribution for each t
s.To investigate the smoothing and forecasting of our proposal, we regenerate a data set with sample size 100 and treat the last 5 × 4 observations as unobserved. The Gibbs sampler with MH algorithm is used to compute the first ninety-five smoothing values and the last five forecasting values of the latent variablesξ
s. Figure 4 gives the true values and smoothing/forecasting values based on 5000 simulated observations after deleting first 5000 observations. It can be seen that it fits well between the estimated values and true values.Figure 4
Comparison with true values (dashed line) and forecasting values (pointed line).


## 5.2. Pika’s Data

An application is considered in relation to a study of Guo et al., 2009 [32], to illustrate the developed methodology. The data set is collected from Qinshui Forest farm of Jincheng Coal Industry Group sited on the northwest of Qinshui County of Jincheng City in Shanxi Province, 112°19′40′′–112°19′52 east longitude and 35°48′43′′–35°48′57 north latitude. Details about experiment design can be found in Guo et al., 2009 [32]. The original data set is constituted of the number of the living tree and pika within four plots in five months. The primary concern is to assess the relationship between the survival rate of young plantation and the pika population and model the optimal pika population per hectare. The data set is reanalyzed by Xia et al. [33] to investigate the correlations among the observed responses. In these studies the latent variables model is established through AR(
1
) model (3) coupled with model (1) and (2), which interprets the correlations among the outcomes resulting from the time dependence.Firstly, the following hypothesis is considered:H
0


:
ϕ
=
0. If H
0 is true, the model becomes the common generalized logistic model considered by Xia et al. [33]. Path sampling is used to compute log BF
10. The linked model is defined similarly as (20). S
=
20 grids were chosen in [
0
,
1
], and, for each t
s (s
=
1
,
…
,
20), 10000 random observations were drawn from the each posterior distributions and the first 5000 observations were deleted in view of the burn-in phase. The logarithm of Bayes factor is 3.0143. It shows that there exists significant evidence against H
0; hence, our proposal is appropriate.To assess the effect of forecasts of the proposal methodology, the original data set is divided into two parts: one is formed by the previous sixteen observations in four months within four plots and the other is made up of four observations in the following fifth month. The proposal algorithm computes the predictive values of the proportion of the living tree in the fifth month based on the previous 20 observations. Figure5 shows the forecasts and actual values of the proportions of the living tree within four plots. It can be seen that fitted values match the actual values closely.Figure 5
Comparison with the predictions for survival rate (solid line) and the actual values (dashed line).


## 6. Concluding Remarks

Logistic regression model is the most popular model used to interpret the relationship between the responses and covariates and to assess the effect of covariates on the responsible probability. Bayesian analysis of logistic model has received a lot of attention recently; see Johnson and Albert [34] and Congdon [35]. For example, Choi et al. [36] discussed a Bayesian statistical inference about missing information on the basis of the logistic regression model. They also used Gibbs sampler to get the estimates and the posterior analysis of the posited model. Since their model is not an AR(
1
) model, the underlying development is less complicated. Recently, McCormic et al. [37] proposed an online binary classification procedure based on the dynamic logistic regression and dynamic model averaging. However, their developments are restricted within a single response y
t, which do not need to explore the relationships among the observed variables. Though Xia et al. [33] establish the Bayesian analysis for generalized logistic model, the essential difference lies in that their latent variables are identically and independently distributed according to normal distribution and, hence, can not capture the temporal correlation.In the present paper, dynamic factor model is established to characterize the temporal correlation among responses. One contribution in this paper is the development of a feasible estimation procedure for obtaining the Bayesian estimates of the parameters and latent factor scores. The hybrid Gibbs sampler with MH algorithm was implemented to provide a convenient mechanism for implementing our method. Another contribution is the development of test statistics on the Bayes factor for testing hypothesis of the model. Computation of the Bayes factor in the current complex model is on the basis of path sampling. Further, Bayesian forecasting procedure which takes advantages of the full conditional distribution is presented. Results from empirical studies indicate that these procedures can be usefully applied to real studies.



---

*Source: 783494-2014-07-07.xml*