Logistic regression has been considered as the standard approach for the analysis of ordered categorical outcome data in application area such as biomedical and social science studies. For the ordinal outcomes, the most common regression models fall under the set of cumulative logit link models proposed by Walker and Duncan (1967) because of a direct result of its simplicity and ease of interpretation. Later it is called the proportional odds model by McCullagh (1980) since it assumes that the relationship between the cumulative log-odds and the explanatory variables does not depend on the response category. The cumulative logit model guarantees that for C ordered groups, the cumulative probability of belonging to group c is less than or equal to the cumulative probability of belonging to group c + 1, for c = 1, …, C − 1 (McKinley et al., 2015). Various inference processes have been proposed for the cumulative logit link models and comprehensive details about models and the analysis of ordered category data can be found in Agresti (2010).

Another commonly used form of modeling the cumulative logit model is based on an unobserved continuous latent variable. McCullagh and Nelder (1989) showed that the cumulative log-odds model could be derived from an underlying unobserved continuous random variable falling into an interval on the real line. A fully Bayesian method for modeling polychotomous ordinal categories based on the latent variable with probit link was developed in Albert and Chib (1993) by using the data augmentation approach. It might be easier to use a normal distribution with probit link instead of logistic distribution. For example, Yi et al.2004 proposed a Bayesian approach via MCMC algorithms to extend the composite model space approach to detect genomewide multiple interaction quantitative trait loci (QTL) for ordinal traits on the basis of a threshold model. A Bayesian method for the cumulative logit link model has been proposed by Lang (1999) such that because of the complexity of the likelihood function, a random-walk Metropolis-Hastings (MH) algorithm with a multivariate-Gaussian proposal distribution is used to estimate the posterior distribution of parameter. Holmes and Held (2006) suggested Gibbs samplers based on the full posterior distribution for logistic multinomial regression models based on auxiliary variable approach, and O’Brien and Dunson (2004) developed a multivariate logistic regression framework that provides a marginal logistic structure for each of the outcomes via a multivariate t-distribution approximation to the multivariate logistic distribution.

For the mean trend in regression, the linear model is commonly used in various field. However, with many explanatory variables, there might be chance of the inference problem of high-dimensional covariates. Principal component (PC) regression has been known as a tool for data analysis and dimension reduction in applications throughout science and engineering. The standard principal component regression can be represented using the reduced singular-value decomposition (SVD) of a matrix of explanatory variables. Then with orthogonal linear combination of predictors, PC, from SVD, the p-predictor linear regression transforms to the k factor linear regression considering subset of PC with high variance. With relevant subset of PC, PC regression model can be used for the prediction of responses. Because only high variance PC corresponding to the higher eigenvalues is considered in the model, the low variance components are discarded in the model, thus for the purpose of the prediction, other shrinkage method such as the ridge regression estimator is perhaps a better choice compared to the PC regression which has a discrete shrinkage effect (Frank and Friedman, 1993). Even though there are many controversies for the interpretation and usage of PC regression, it is still a useful tool when there exists a multicollinearity problem among explanatory variables in the regression models.

The standard PC regression performs principal component analysis (PCA) on the centered explanatory variables matrix, and selects subset of the PC based on some appropriate criteria. Then, the observed vector of responses is regressed linearly on the selected PC as covariates, and the regression coefficients are estimated via ordinary least square method. The estimated PC regression coefficients are transformed back to the parameters of the actual predictors in the original model using the selected PCA loadings. In these steps, the principal components are obtained from the SVD of the explanatory variables matrix, thus the PC regression estimator obtained from using these principal components may produce unsatisfactory results for the prediction of the responses. Bair et al.2006 suggested a technique called supervised principal components which is similar to conventional PCA except that it uses a subset of the predictors selected based on their association with the outcome. To choose the reduced rank of data, a threshold is considered and cross-validation of the log-likelihood ratio test was considered to estimate the best value of threshold. In Bayesian framework, West (2003) discussed the standard principal component regression with shrinkage prior on the PC regression coefficients allowing for varying degrees of shrinkage in each of the orthogonal PC dimensions.

In this paper, we propose a Bayesian approach for the cumulative log-odds PC regression model to deal with ordinal response based on an approximation of a t-distribution to a logistic distribution. The purpose that we propose a Bayesian cumulative logit PC regression model is to make MCMC implementation easier and to consider multicollinearity problem with interpreting parameters as the cumulative log-odds model. The close correspondence between the logistic and t distribution was noted by Albert and Chib (1993) such that for probabilities between 0.001 and 0.999, logistic quantiles are approximately a linear function of t(8) quantiles. O’Brien and Dunson (2004) instead considered a multivariate t-distribution with σ² = π²(ν − 2)/3ν, a value chosen to make the variance of the univariate t and logistic distribution equal, and ν = 7.3, a value chosen to minimize the integrated squared distance between the univariate t and univariate logistic densities. The reason for considering a t-distribution approximation to the logistic is that if we take back transformation of logistic link, it will facilitate interpretability by explaining the parameter of interest in the form of odd ratios. Even with an approximation, it becomes possible to interpret the effects of predictors to the probabilities of each category directly. In contrast, with a probit link, coefficients for probit models can be interpreted as the difference in standard normal score associated with each one-unit difference in a predictor variable, which is not intuitive to interpret its effect. Thus we prefer to use a logistic model based on the t-distribution approximation. Also, for the mean, we consider PC regression with shrinkage priors of West (2003) on the regression coefficients, and we choose the number of subset of PC based on the Bayesian information criterion (BIC). Usually, BIC seems to be a popular criterion for model selection because of a simple form and an easy application. Raftery (1995) argued that Bayes factors allow the direct comparison of non-nested models, in a simple way and BIC provides a simple and accurate approximation to Bayes factors for the Bayesian approach to hypothesis testing, model selection and accounting for model uncertainty. In the frame of Bayesian inference, it is easy to obtain BIC using posterior mean from Gibbs samples of posterior distribution of parameters. Thus we consider that for our Bayesian inference of the cumulative log-odds model with PC regression, BIC might be a good model selection criterion to choose the number of subset of PC within regression model framework.

In the remainders of the paper, we present the details about the cumulative log-odds model with the PC based regression and subsequently describe about the Bayesian inference with conjugate priors. We then apply the proposed method to real data concerning sprout- and scab-damaged kernels of wheat. We finally discuss about conclusion.

For the Bayesian analysis of ordered categorical outcomes, we express the ordinal logistic regression model via underlying continuous variables. Then with PC regression, we extend to the cumulative log-odds PC regression model.

2.1. Cumulative logistic model

For the polychotomous ordinal outcomes of C categories, we usually assume that

The correspondence between the observed outcome and the latent variable is defined by

where X_i is the i^th observed predictors of length p, β is a p × 1 vector of regression parameter, γ = (γ₀, γ₁, …, γ_C)′ is a vector of cutpoints with −∞ = γ₀ < γ₁ < … < γ_C−1 < γ_C = ∞, and ?(μ,σ²) is a density function of logistic distribution with mean μ and scale parameter σ². With given notation and σ² = 1, the cumulative logit model for the j^th ordinal outcome can be expressed as

for j = 1, 2, …, C − 1. Here, the logit link is monotone increasing, and for the ease of notation, we let P(Y_i ≤ j|X_i, β, γ) ≡ η_{i j}, then the probabilities of each category are computed by

In (2.1), the effect of the predictors is independent of the category involved, and subject to the constraint −∞ = γ₀ < γ₁ < … < γ_C−1 < γ_C = ∞, the cumulative probability will be such that 0 < η_i1 < η_i2 < … < η_i(C−1) < 1, strict stochastic ordering (McCullagh, 1980).

As we discussed in Section 1, the cumulative logistic model is known as the proportional odds model (McCullagh, 1980). It is because the cumulative log-odds ratio for two sets of predictors, X₁ and X₂ is given by

Hence, the ratio of corresponding odds is independent of category j and depends only on the difference between the covariate values X₂ − X₁.

2.2. Principal component regression model

In statistics, principal component regression (PCR) is a regression analysis technique that is based on principal component analysis. Massy (1965) introduced PCR to simplify the sample space of explanatory variables by a transformation to PC, for the cases such that the original independent variables are highly collinear with one another or there are a great number of potential explanatory variables. PCR regresses on a small number of regressors accounting for most of the variability of the explanatory variables. Massy (1965) also proposed a method choosing PCs which are most highly correlated with the response to be include in the regression and thus to explain the most variation therein.

Reiss and Ogden (2007) extended PCR to a functional version of PCR with splines, and introduced the reduced singular-value decomposition (SVD) of a matrix of explanatory variables to form a new design matrix of PCR. We consider the SVD of n × p explanatory matrix X, X = U∑V^T. Here, U is a n × n orthonormal matrix and the column vectors are the left singular vectors of X which can be obtained by the spectral decomposition of XX^T. V is a p × p orthonormal matrix and the column vectors are the right singular vector of X from the spectral decomposition of X^TX. Also ∑ is a n × p diagonal matrix of singular values (square root of eigenvalues) of XX^T or X^TX in descending order.

Let U_K, ∑_K, and V_K be the selected K versions of the three SVD matrices, then we assume that

where F = U_K∑_K is a n × K factor matrix and $\mathit{A}={\mathit{V}}_L^T$ is a K × p SVD loading matrix with K < min(n, p). Here, AA^T = I, ${\mathit{F}}^T\mathit{F}={\mathbf{\Sigma }}_K^{\text{2}}$ and ∑_K is the K × K upper left submatrix of ∑ which we assume to be nonsingular. For the linear model, we usually assume that y = Xβ + ε where y is the n-vector of responses and ε ~ N(0, σ²I) is the n-vector error term. Then K component PCR can be expressed as such that

where ε ~ MVN_n(0, σ²I). Here, the regression transforms via Xβ = Fθ where θ = Aβ is the length K vector of regression parameters for the PCR. The original regression parameter β can be obtained by the least-norm inverse β = A^Tθ.

2.3. Cumulative logistic principal component regression model

A useful feature of the logistic distribution is that the logistic density and the t-distribution density closely approximate one other when the degrees of freedom ν and the variance ξ² are chosen appropriately (O’Brien and Dunson, 2004). This fact allows us that an approximation to the posterior distribution can be derived by considering an appropriate t-distribution for the logistic distribution. Albert and Chib (1993) noted that logistic quantiles are approximately a linear function of t(8) quantiles. O’Brien and Dunson (2004) considered a multivariate t-distribution with ξ² = π²(ν − 2)/3ν, a value chosen to make the variance of the univariate t and logistic distribution equal, and ν = 7.3, a value chosen to minimize the integrated squared distance between the univariate t and univariate logistic densities. Thus with properties of cumulative logistic regression that are discussed in Section 2.1, if we combine the t-approximation to the logistic regression and the PCR to consider the multicollinearity among X which is discussed in section 2.2, the cumulative logistic PC regression model with latent variable Z_i can be expressed that for i = 1, …, n and j = 1, …, C,

where ν = 7.3 and ψ = π²(ν − 2)/3ν, a value chosen to make the variance of the univariate t and logistic distribution equal in our model.

For the t-distribution with location μ, scale ξ and degrees of freedom ν, the density function can be expressed in the following scale mixtures of normal form,

where Gamma(.|a, b) is the gamma density function. It means that the t-distribution can be expressed as the following hierarchical form such that

This representation is considered for the model in (2.3), and based on the model in (2.3), we consider the Bayesian methods for the inference of parameters, and the details are in the following section.

In Bayesian framework, usually an unknown parameter is treated as a random variable, and uncertainty about the parameters is inferred by its probability distribution, the so-called posterior distribution. The posterior distribution is a combination of evidence collected from data that is formally expressed by the likelihood function and a prior belief that is specified with a prior distribution. To update a posterior distribution for the interest parameter, a Markov chain Monte Carlo (MCMC) algorithm is used, through which it enables to iteratively draw sample from even a complex distribution without analytically solving for the point estimate.

3.1. Specifications of prior distribution

In our work, we consider shrinkage priors on the underlying factor regression parameters θ (West, 2003). West (2003) discussed that a particularly flexible class of such priors assumed independent t-distributions for the elements of θ, so allowing for varying degrees of shrinkage in each of the orthogonal factor dimensions. The shrinkage priors on the PC regression parameters are the following: for k = 1, …, K,

independently, for some r > 0. The tuning parameter r is the degree-of-freedom parameter for the implied t-distribution for θ_k that follows on marginalization over the random precision δ_k. Other prior distributions on the parameters of the model in (2.3) and prior distribution on the hyper-parameter ρ are the following,

where j = 1, …, C − 1 and IG is the inverse gamma distribution. Here, the normal distribution of a₀ is indexed by two hyperparameters, which means we can specify a particular prior distribution by fixing two features of the distribution. However, without loss of generality, we usually set α₀ = 0 and η² with large values such as 100.

3.2. Sampling scheme

With model in (2.3) and priors in (3.1) and (3.2), the Gibbs sampling is straightforward because the most of conditional posterior distributions are in closed form except for the tuning parameter r. The conditional posterior distribution of r is no closed form, but the sampling from the conditional posterior distribution can be processed by a Metropolis-Hastings algorithm. Thus, we iterate sampling from their full conditional posterior distributions until convergence.

With given the number of PC, K, the joint posterior distribution has the form of

where G = diag(ρ/k²δ_k)_{k=1,..., K} and ${\mathit{\theta }}^T{\mathit{G}}^{-1}\mathit{\theta }=1/\rho \hspace{0.17em}{\mathrm{\Sigma }}_{k=1}^K{k}^2{\delta }_k{\theta }_k^2$. Then, the posterior sampling step is carried out using the Gibbs sampling with priors as discussed earlier. For the tuning parameter r, we consider a Metropolis-Hastings algorithm because the conditional distribution is no closed form. But for other parameters, we use the Gibbs because of the conjugacy. Further details on the sampling scheme with the full conditional distributions are provided in the Appendix A.

There is an identifiability problem in intercept a₀, cutpoints γ, and variance σ². Hirk et al.2019 pointed out that in ordinal models absolute location and absolute scale of the underlying latent variable are not identifiable. It means that in our model, only the quantities θ/σ^* and (γ_j − a₀)/σ^*, j = 1, …, C − 1 are identifiable. Here, σ^*2 is the full variance of the error term. Thus, to secure the identifiability of scaled intercepts a₀/σ^* and scaled cutpoints γ_j/σ^*, we consider a typical constraints on the parameters such that γ₁ = 0 and σ² = 1 in our model. It is a common practice in Bayesian logistic regression models for both univariate and multivariate binary outcomes. The fixing of the cutpoints at zero in the process of MCMC is to prevent the identifiability problem of the parameter estimation. Thus, with these constraints, we update the posterior samples of the parameters iteratively.

As we discussed in Section 1, we choose the number of subset of PC, K, based on the Bayesian information criterion (BIC). In Bayesian analysis, Bayes factors is a tool of model selection. To compute the Bayes factor, calculating the integrated likelihood over the parameters space is involved and it can be a high-dimensional and intractable integral. Exact analytic integration might be possible for exponential family distributions with conjugate priors including normal linear models. However, Kass and Raftery (1995) argued that it is more often intractable and thus must be computed by numerical methods, and introduced various asymptotic approximation to the marginal density of the data. BIC (The Schwarz criterion) was one of the asymptotic approximation tools which Kass and Rafter (1995) introduced. Raftery (1995) argued that Bayes factors allow the direct comparison of non-nested models in a simple way and BIC provides a simple and accurate approximation to Bayes factors for the Bayesian approach to hypothesis testing, model selection and accounting for model uncertainty. BIC can be applied as a standard procedure even when the priors of the parameters are hard to set precisely. Another advantage of BIC is that in the frame of Bayesian inference, it is easy to obtain BIC using posterior mean from Gibbs samples of posterior distribution. Thus we consider that for our Bayesian inference of the cumulative log-odds model with PC regression, BIC might be a good model selection criterion to choose the number of subset of PC within regression model framework.

For the simulation study and data analysis, to analyze the data in the lack of the information of prior, diffuse prior distributions were considered. For this, we set the hyperparameters of prior distribution of intercept, the mean and variance, of a₀ as α₀ = 0 and η² = 100. For the posterior inference, we obtained the posterior median and the 95% credible intervals of the parameter estimates from the proposed MCMC algorithm ran for 30, 000 iterations with a burn-in period of 15, 000 iterations and a thinning by taking every fifth posterior sample. Also, to fit the Bayesian cumulative logistic PC regression on data, we need to decide the optimal number of components. As discussed in Section 3.2, we compute the BIC of the Bayesian cumulative logistic PC regression for the number of components from 2 to a large value L (L can be vary because of the data structure). Based on the BIC plot, we choose the number of components with the smallest BIC value.

For the comparison of the suggested model, we consider the EM based inference with linear mean trend using polr (proportional-odds logistic regression) function in MASS package (Venables and Ripley, 2002) and bayespolr function in arm package (Gelman and Su, 2020) in R (R core team, 2021). The polr function is developed based on the traditional maximum likelihood method using Fisher scoring or the Newton-Raphson algorithm. The bayespolr function is developed based on the weakly informative priors for logistic regression (Gelman et al., 2008, 2015) assuming a Student-t prior distribution with mean 0, degrees-of-freedom (df) and scale with df and scale chosen to provide minimal prior information to constrain the coefficients to lie in a reasonable range. The full Bayesian computation using the usual MCMC methods can be applied, but for a quick iterative calculation to get a point estimate of the regression coefficients and standard errors, an approximate EM algorithm is applied.

4.1. Simulation study

We conduct a simulation study to evaluate how well the proposed method performs with considering a high-dimensional and highly correlated explanatory matrix. For a high-dimensional predictors, we set the number of observations as n = 50 and the number of explanatory variables as p = 70. To generate highly correlated explanatory matrix, we set the correlation coefficients of X_i and X_j as ${\sigma }_{ij}^X={(0.8)}^{?i-j?}$ and the standard deviation of X_i is 0.1. With given covariance matrix and mean zero vector, we generate n×p explanatory matrix from a multivariate normal distribution. For the true value of regression parameter, we set 35 ones and 35 zeros such that β = (1, …, 1, 0, …, 0). Incorporating a regression model into the mean Xβ, an ordered categorical outcome y were determined based on the logistic distribution and the predefined cutpoints (−∞, −1.9, 0, 1.8,∞).

To decide the number of components that are considered in the Bayesian cumulative logistic PC regression, we computed BIC for each number of components and the BIC plot up to 15 number of components is in Figure 1. In Figure 1, we observe that the minimum value of BIC is obtained when the number of principal components is 8, k = 8. Thus, for the remaining analysis and model comparison, we employ 8 principal components of explanatory matrix as regressors and after fitting the Bayesian cumulative logistic PC regression, we re-transfer it to the regression parameter as we discussed in Section 2.2. Also, to check up the convergence of MCMC, we draw trace plots and autocorrelation plots for few selected θs and the plot is in Figure B.1 in Appendix B. The plots how that the proposed Gibbs sampler converged fast and mixed well.

Subset of estimated values and 95% confidence intervals of selected regression parameters in EM based proportional-odds logistic regression (polr) and in EM based Bayesian proportional-odds logistic regression (bayespolr), and estimated values (posterior median) and 95% credible intervals of regression parameters in Bayesian cumulative logistic PC regression (bayesPC(K=3)) for simulated data are in Table B.1 in Appendix B. Also, estimated values and 95% confidence intervals of re-transformed regression parameters of the suggested model can be found in Figure 2. Because the explanatory matrix is not full rank, EM based proportional-odds logistic regression model does not use full predictors, and it only estimated the regression parameter of 49 explanatory variables and cutpoints. However, all the estimated regression parameters and cutpoints were not significant because of very largely computed variance of estimators. Thus, it is not enough to use this model for the prediction and analysis in this highly correlated high-dimensional data. Bayesian proportional-odds logistic regression shows almost same results of EM based proportional-odds logistic regression model. However, because of dimension reduction effects on PC regression and the considered prior is a generalized shrinkage prior, the propose model provides reasonable parameter estimates. Even though somewhat amplified estimates and largely valued variances were observed, the proposed model detects the true zero values of the last 35 regression parameters in the setting of highly correlated high-dimensional data. Also, the model fit based on the mean squared error could be computed and the computed MSE was 0.06 with given explanatory matrix X which is used to fit the model. Thus, we argue that the proposed model works better for the highly correlated high-dimensional data with providing parameter estimates and good prediction.

4.2. A real data example

A Bayesian approach to cumulative logistic regression model with orthogonal principal components based on singular value decomposition was proposed to consider the multicollinearity among predictors. The suggested method has advantage such as considering dimension reduction and parameter estimation simultaneously. However, if the aim of the analysis is the prediction of the response variable, the suggested model has limitation, because the suggested method uses components with largely valued variance and ignores components with small valued variance. Also, the considered prior is a generalized shrinkage prior and West (2003) discussed the standard principal component regression with shrinkage prior on the PC regression coefficients allowing for varying degrees of shrinkage in each of the orthogonal PC dimensions. Thus, the suggested model might not be good to be considered for the prediction. However, if the multicollinearity problem is serious among explanatory variables, there might be the inference problem of high-dimensional covariates.

We conducted a simulation study to show how well the proposed method performs with considering a high-dimensional and highly correlated predictors. Because the predictor matrix is not full rank, compared to EM based proportional-odds logistic regression (polr) and to EM based Bayesian proportional-odds logistic regression (bayespolr), the proposed Bayesian cumulative logistic PC regression works better for the highly correlated high-dimensional data with providing parameter estimates and good prediction. Also, we applied the proposed method to real data concerning sprout- and scab-damaged kernels of wheat. Because of the shrinkage priors on the components, the estimated regression parameters were small valued estimates. Compared to EM based proportional-odds logistic regression (polr) and to EM based Bayesian proportional-odds logistic regression (bayespolr), the proposed Bayesian cumulative logistic PC regression showed limitation on the prediction. Also, the estimated parameters of density, and the estimates of the first and second cutpoints of the suggested method were very different from other methods numerically. We argue here that the scale of density variable much larger than other explanatory variables and it affects the position of cutpoints in somewhat ways.

Methodologically, the proposed method can be further extended to accommodate latent factor regression models for partitioning variation in the predictor variables into multiple component that reflect common patterns and separating these from noise variation. Although the suggested model considered the orthogonal components, there was a limitation on prediction. With latent factors, there might not be the orthogonal property, but we might be able to find out common structure in explanatory matrix isolating idiosyncratic variation for dimension reduction in predictor space.

The sampling scheme is as follows. With ν = 7.3 and ψ = π²(ν − 2)/3ν

• update Z: for i = 1, …, n and j = 1, …, C

  $$Z_i?\mathit{\theta },{\phi }_i,a_0,{\sigma }^2,\mathit{\delta },r,\rho ,\mathit{\gamma },\mathit{X},\mathit{Y}~N\hspace{0.17em}\left(a_0+{\mathit{F}}_i\mathit{\theta },\hspace{0.17em}\psi {\phi }_i^{-1}{\sigma }^2\right),$$

  with Z_i is truncated at the left and right by γ_j−1 and γ_j if Y_i = j.

• update γ: for j = 1, …, C − 1

  $$\begin{array}{l}{\gamma }_j?{\mathit{\gamma }}_{-j},\mathit{\theta },\mathit{\phi },a_0,{\sigma }^2,\mathit{\delta },r,\rho ,\mathit{Z}?\mathit{X},\mathit{Y}\\ ~U\hspace{0.17em}\left(\text{max\hspace{0.17em}}\left\{\text{max\hspace{0.17em}}(Z_i:Y_i=j),{\gamma }_{j-1}\right\},\text{min\hspace{0.17em}}\left\{\text{min\hspace{0.17em}}(Z_i:Y_i=j+1),{\gamma }_{j+1}\right\}\right)\hspace{0.17em}\text{I\hspace{0.17em}}\left({\gamma }_{j-1}<{\gamma }_j<{\gamma }_{j+1}\right),\end{array}$$

  where γ_−j = (γ₁, …, γ_j−1, γ_j+1, …, γ_C).

• update θ

  $$\begin{array}{l}\mathit{\theta }?\mathit{\phi },a_0,{\sigma }^2,\mathit{\delta },r,\rho ,\mathit{\gamma },\mathit{Z},\mathit{X},\mathit{Y}\\ ~{\text{MVN}}_K\hspace{0.17em}\left({\left(\frac{1}{\psi {\sigma }^2}\sum _{i=1}^n{\phi }_i{\mathit{F}}_i^T{\mathit{F}}_i+{\mathit{G}}^{-1}\right)}^{-1}\hspace{0.17em}\left(\frac{1}{\psi {\sigma }^2}\sum _{i=1}^n{\phi }_i{\mathit{F}}_i^T\hspace{0.17em}(Z_i-a_0)\right),\hspace{0.17em}{\left(\frac{1}{\psi {\sigma }^2}\sum _{i=1}^n{\phi }_i{\mathit{F}}_i^T{\mathit{F}}_i+{\mathit{G}}^{-1}\right)}^{-1}\right),\end{array}$$

  where MVN_K is the K-dimensional multivariate normal distribution.

• update φ: for i = 1, …, n

  $${\phi }_i?\mathit{\theta },a_0,{\sigma }^2,\mathit{\delta },r,\rho ,\mathit{\gamma },\mathit{Z},\mathit{X},\mathit{Y}~\text{Gamma\hspace{0.17em}}\left(\frac{\nu }{2}+\frac{1}{2},\frac{\nu }{2}+\frac{1}{2\psi {\sigma }^2}{(Z_i-a_0-{\mathit{F}}_i\mathit{\theta })}^T\hspace{0.17em}(Z_i-a_0-{\mathit{F}}_i\mathit{\theta })\right),$$

• update a₀

  $$\begin{array}{l}{a}_0?\mathit{\theta },\mathit{\phi },{\sigma }^2,\mathit{\delta },r,\rho ,\mathit{\gamma },\mathit{Z},\mathit{X},\mathit{Y}\\ ~N\hspace{0.17em}\left({\left(\frac{1}{\psi {\sigma }^2}\sum _{i=1}^n{\phi }_i+\frac{1}{{\eta }^2}\right)}^{-1}\hspace{0.17em}\left(\frac{1}{\psi {\sigma }^2}\sum _{i=1}^n{\phi }_i(Z_i-{\mathit{F}}_i\mathit{\theta })+\frac{{\alpha }_0}{{\eta }_2}\right),\hspace{0.17em}{\left(\frac{1}{\psi {\sigma }^2}\sum _{i=1}^n{\phi }_i+\frac{1}{{\eta }^2}\right)}^{-1}\right),\end{array}$$

• update σ²

  $${\sigma }^2?\mathit{\theta },\mathit{\phi },a_0,\mathit{\delta },r,\rho ,\mathit{\gamma },\mathit{Z},\mathit{X},\mathit{Y}~\text{IG\hspace{0.17em}}\left(a+\frac{n}{2},\hspace{0.17em}b+\sum _{i=1}^n{\phi }_i\hspace{0.17em}{(Z_i-a_0-{\mathit{F}}_i\mathit{\theta })}^T\hspace{0.17em}(Z_i-a_0-{\mathit{F}}_i\mathit{\theta })\right),$$

• update δ: for k = 1, …, K

  $${\delta }_k?\mathit{\theta },\mathit{\phi },a_0,{\sigma }^2,r,\rho ,\mathit{\gamma },\mathit{Z},\mathit{X},\mathit{Y}~\text{Gamma\hspace{0.17em}}\left(\frac{r}{2}+\frac{1}{2},\frac{r}{2}+\frac{1}{2\rho }{k}^2{\theta }_k^2\right),$$

• update r

  $$\begin{array}{l}r?\mathit{\theta },\mathit{\phi },a_0,{\sigma }^2,\mathit{\delta },\rho ,\mathit{\gamma },\mathit{Z}?\mathit{X},\mathit{Y}\\ ~\pi \hspace{0.17em}\left(r?\mathit{\theta },\mathit{\phi },a_0,{\sigma }^2,\mathit{\delta },\rho ,\mathit{\gamma },\mathit{Z}?\mathit{X},\mathit{Y}\right)\propto \frac{{\left(\frac{r}{2}\right)}^{\frac{r}{2}K}}{{\left\{\mathrm{\Gamma }\hspace{0.17em}\left(\frac{r}{2}\right)\right\}}^K}\hspace{0.17em}\text{exp\hspace{0.17em}}\left\{-r\left(\frac{1}{2}\sum _{k=1}^K{\delta }_k-\frac{1}{2}\sum _{k=1}^K\text{ln\hspace{0.17em}}{\delta }_k\right)\right\}.\end{array}$$

  For the ease of notation, let the target distribution as $P\hspace{0.17em}(r)=\{{(r/2)}^{(r/2)K}\}/\{{\{\mathrm{\Gamma }(r/2)\}}^K\}\hspace{0.17em}\text{exp}\{-r((1/2){\mathrm{\Sigma }}_{k=1}^K{\delta }_k-(1/2){\mathrm{\Sigma }}_{k=1}^K\text{ln\hspace{0.17em}}{\delta }_k)\}$. Consider an exponential distribution with parameter $\lambda =(1/2){\mathrm{\Sigma }}_{k=1}^K{\delta }_k-(1/2){\mathrm{\Sigma }}_{k=1}^K\text{ln\hspace{0.17em}}{\delta }_k$ as a candidate distribution

  $$g^{*}\hspace{0.17em}(r?r_t)=\text{Exponential\hspace{0.17em}}\left(r|\lambda =\frac{1}{2}\sum _{k=1}^K{\delta }_k-\frac{1}{2}\sum _{k=1}^K\text{ln\hspace{0.17em}}{\delta }_k\right),$$

  then for each iteration t,

  • (a) generate a random candidate r′ from the candidate distribution r′ ~ g^* (r|r_t)

  • (b) calculate the acceptance probability

    $$A\hspace{0.17em}(r^{\prime },r_t)=\text{min\hspace{0.17em}}\left(1,\frac{P\hspace{0.17em}(r^{\prime })}{P\hspace{0.17em}(r_t)}\frac{g^{*}\hspace{0.17em}(r_t?r^{*})}{g^{*}\hspace{0.17em}(r^{*}?r_t)}\right)$$

  • (c) r_t+1 = r′ if u ≤ A (r′, r_t) or r_t+1 = r_t if u > A (r′, r_t), where u ~ U(0, 1)

• update ρ

  $$\rho ?\mathit{\theta },\mathit{\phi },a_0,{\sigma }^2,\mathit{\delta },r,\mathit{\gamma },\mathit{Z}?\mathit{X},\mathit{Y}~\text{IG\hspace{0.17em}}\left(\frac{K}{2}-1,\frac{1}{2}\sum _{k=1}^K{k}^2{\delta }_k{\theta }_k^2\right).$$