Bayesian analysis of a multivariate spatial ordered probit model

Abstract Many phenomena are correlated and spatially dependent. In addition, some of them are ordinal discrete responses. Thus, a model is needed to capture the interactions of multivariate ordinal outcomes and spatial dependences. Following Smith and LeSage (2004) and Jeliazkov et al. (2008), this study proposes a new algorithm for a multivariate spatial ordered probit (MSOP) model to address this need. In applying this model, the parameters are calculated using the Bayesian inference based on Markov chain Monte Carlo (MCMC) sampling. The validity and accuracy of the MSOP model is verified by simulated datasets, and the model performs very well with the simulated data. In addition, this study illustrates the model by applying it to two response variables, self-rated health, and life satisfaction of elderly people in 18 representative provinces in mainland China. The empirical results show that the spatial dependences are indispensable on the response variables.


Introduction
Many social phenomena of interest are spatially dependent, for example, land use change (Zhou and Kockelman 2008), high-yielding variety rice adoption (Holloway et al. 2002), and even the presidential election (Smith and LeSage 2004). In addition, some response variables are discrete, binary, or polychotomous. Thus, spatial probit models are necessary to meet the demand. Spatial probit models are the extension of probit models that take spatial dependence into account.
The traditional estimation method of probit models is the maximum likelihood (ML) method, which is based on the asymptotic theory. Thus, in the case of a small sample size, accuracy of the ML method is questionable. Albert and Chib (1993) introduced latent continuous variables into probit models and proposed the Bayesian inference for the estimation. Meanwhile, they proved that latent continuous variables can be obtained by using truncated normal distributions. Following Albert and Chib (1993), numerous studies have verified this approach as being effective (Chen and Dey 2000;Hasegawa 2010;Jeliazkov et al. 2008;LeSage 1999).
Spatial probit models were first introduced by McMillen (1992), wherein an expectation-maximization algorithm was used for the estimation. In addition, some researchers examined the ML method (Arbia 2006;LeSage 1999), but it is difficult to estimate (Beron and Vijverberg 2004). Following Albert and Chib (1993), LeSage (2010) proposed the Markov chain Monte Carlo (MCMC) method to estimate the spatial probit models. Smith and LeSage (2004) extended the model by including an additive error specification. Furthermore, Kakamu and Wago (2007) extended that model to the spatial panel probit model by including panel data, and they applied it to the business cycle in Japan. All the response variables in the above spatial probit models are binary, but many of the response data are of multiple categories. Thus, Wang and Kockelman (2009) extended Smith and LeSage (2004)'s model to the dynamic spatial ordered probit model, wherein they also considered time-series dependence.
All of the above studies focused on the univariate response. Regarding multivariate responses, Chen and Dey (2000), Chib and Greenberg (1998), and Jeliazkov et al. (2008) illustrated the multivariate probit models with binary and ordinal responses. Using such a model, Hasegawa (2010) analyzed tourists ' satisfaction, and Lauer (2003) analyzed the impact of family background, cohort, and gender on educational attainment in France and Germany.
Currently, there are few models considering multivariate responses in spatial econometrics. Since many phenomena are correlated and spatially dependent, it is better to study them together, rather than separately. Thus, a model that captures the interactions of multivariate responses and spatial dependence is necessary. In addition, some of the responses are ordered discrete data. Therefore, the model must also deal with ordered discrete response variables instead of continuous response variables.
It is difficult to analyze the multivariate ordered probit model with more than two response variables using the classical non-Bayesian method. Such models can be easily analyzed using the Bayesian method (Hasegawa 2010). However, existing Bayesian multivariate probit models do not consider the spatial effect. Following the probit model with spatial dependences in Smith and LeSage (2004) and models for multivariate ordinal outcomes in Jeliazkov et al. (2008), this study proposes a multivariate spatial ordered probit (MSOP) model to meet the above requirements. The MSOP model is the first to analyze the spatial dependences of multivariate ordered response variables. In applying this model, the parameters are calculated using the Bayesian inference based on MCMC sampling.
The rest of this study is structured as follows. Section 2 introduces the model specification for the MSOP model. Section 3 explains the details of Bayesian inference. Section 4 presents the simulated dataset and its simulation results. Section 5 explains the empirical data, and discusses the empirical results for the two ordinal responses: self-rated health and life satisfaction. Section 6 presents the conclusions.

Model specification
Following the probit model with spatial dependences in Smith and LeSage (2004) and models for multivariate ordinal outcomes in Jeliazkov et al. (2008), this study proposes a multivariate spatial ordered probit (MSOP) model, wherein observed response variables y iks are the censored forms of latent response variables z iks : y iks ¼ l s if c s, l s À1 < z iks c s, l s , for s ¼ 1, :::, S, and l s ¼ 1, :::, L s , where i indexes spatial regions, k indexes individuals inside each region (k ¼ 1, :::, n i ), and s indexes individual response variables. Therefore, y iks denotes the sth observed response variable of individual k in region i. l s indexes the classifications of response variable y iks . In addition, the cutpoints c s, l s are specified as in Chen and Dey (2000): The latent response variables z iks are specified as follows: for i ¼ 1, :::, m, k ¼ 1, :::, n i , and s ¼ 1, :::, S, where x iks ¼ ðx iks1 , x iks2 , :::, x iksq s Þ 0 is a q s -dimensional vector of the explanatory variables. Correspondingly, b s ¼ ðb s1 , b s2 , :::, b sq s Þ 0 is a q s -dimensional vector of the regression parameters for x iks : d is denotes the regional effect, which is used to capture all unobserved common effects for observations within region i for the sth outcome. n iks denotes the individualistic effect of individual k in region i for the sth outcome. The regional effects may exhibit spatial autocorrelation. The following spatial autoregressive model can be used when d is is spatially correlated: w ij d js þ is , for i ¼ 1, :::, m, and s ¼ 1, :::, S, where q is the spatial coefficient, which reflects the degree of spatial influence. Weights w ij reflect the closeness between regions i and j. A commonly used definition is the queen contiguity, which defines w 0 ij ¼ 1 for regions sharing a common edge or a common vertex; otherwise, w 0 ij ¼ 0: In addition, w 0 ii ¼ 0: After that, weights w 0 ij are row normalized to w ij as w ij ¼ w 0 ij = P m j¼1 w ij 0 : Note that the "island" situation (no contiguity region) is not suitable in this MSOP model, because it is very different in essence from the data generating process of spatially connected regions. is are assumed to follow a normal distribution; that is, is $ Nð0, r 2 s Þ: In this MSOP model, the response variables are the multivariate ordered variables that an observation has, and the regional effect d is is used to capture all unobserved common effects for observations within region i. Therefore, the spatial coefficient q could be assumed to be constant for all s in the above spatial autoregressive model.
Stacking all S responses, (3) can be written as: where z ik is an S-dimensional vector of the latent response variables, and thus z ik ¼ ðz ik1 , :::, z iks :::, z ikS Þ 0 : The matrix of covariates X ik is organized as: where x iks is a q s -dimensional vector of the explanatory variables, and thus the number of columns of X ik is P S s¼1 q s ¼ q: Correspondingly, b ¼ ðb 0 1 , :::, b 0 s , :::, b 0 S Þ 0 is a q-dimensional vector of the regression parameters. d i ¼ ðd i1 , :::, d is , :::, d iS Þ 0 and n ik ¼ ðn ik1 , :::, n iks , :::, n ikS Þ 0 are S-dimensional vectors of the regional effects and individualistic effects, respectively.
In this model, the unobserved common effects d is may well be correlated between outcome variables. Therefore, is are assumed to be correlated across S response variables. Defining i ¼ ð i1 , :::, is , :::, iS Þ, i is assumed to follow a multivariate normal distribution: where R ½s, s ¼ r 2 s : The individualistic effects n ik are assumed conditionally iid normal variates with zero means and common covariance matrix V; that is, where V is an S Â S symmetric positive definite covariance matrix. Furthermore, the following equations are defined: G iks ¼ c s, l s À1 , c s, l s ð if y iks ¼ l s , for l s ¼ 1, :::, L s , and s ¼ 1, :::, S, G ik ¼ G ik1 Â Á Á Á Â G ikS , for i ¼ 1, :::, m, and k ¼ 1, :::, n i : Defining y ik ¼ y ik1 , :::, y ikS ð Þ 0 as an S-dimensional vector of the response variables, it can be derived that: where 1(Á) is an indicator function.
where n i denotes an n iS -dimensional vector of the individualistic effects for all the individuals in region i (n ik is an S-dimensional vector, and thus the dimension of n i is P n i k¼1 S ¼ n i S). z i denotes an n iS -dimensional vector of the latent response variables for all the individuals in region i, and X i denotes an n i S Â q matrix of covariates for all the individuals in region i (X ik is an S Â q matrix of covariates, and thus the dimension of X i is P n i k¼1 S Â q ¼ n i S Â q). Based on the above definitions, (5) can be rewritten as: where 1 i denotes a n i Â 1 vector of 1's. Moreover, stacking all regions, and then the MSOP model can be summarized as: , I S , I n , and I mS are identity matrices. Defining B q ¼ I ms À q W I S ð Þ , and B q is assumed nonsigular. Then, (9) can be expressed as: : Therefore, the density of d is: : As a result, the density of z is: When we pay attention to the interaction of the multivariate ordered variables an observation has, and want to capture the unobserved common regional effects, the MSOP model in this study could be applied. For example, the MSOP model can be used to examine the interaction of individual education, medical, and resident satisfaction. These satisfactions are ordered and correlated. They are not only affected by the individuals' conditions and supporting investment in the region where they are located, but also affected by the surrounding regions (when the supporting investment in other regions is higher than that in the local region, individual satisfaction in the local region may tend to decrease).

Prior distributions of parameters
To complete the Bayesian inference, the following prior distributions of the parameters are introduced: b $ N c, T ð Þ; V $ IW n 0 , Q 0 ð Þ; R $ IW n , Q ð Þ; and q $ U k À1 min , k À1 max À Á : IW n 0 , Q 0 ð Þ and IW n , Q ð Þ denote the Inverse-Wishart distributions with degree of freedom n 0 and n and S Â S scale matrix Q 0 and Q , respectively. k min and k max are the minimum and maximum eigenvalues of W: In addition, because c has the limitation of À1 ¼ c s, 0 < c s, 1 ¼ 0 < c s, 2 < Á Á Á < c s, L s À1 ¼ 1 < c s, L s ¼ 1, p c ð Þ / 1 for c s, 2 < ::: < c s, L s À2 , s ¼ 1, :::, S:

MCMC sampler
The joint posterior density function of the parameters is as follows: where the probability of y iks equals 1 when c s, l s À1 < z iks c s, l s , and 0 otherwise. Given z and c, y ik are conditionally i.i.d, and thus Before giving the parameter estimations, Cowles (1996)'s findings must be explained first. The sample of c can be obtained using a uniform distribution, and almost all the literature on the theme of spatial ordered probit models used such an estimation method (Wang and Kockelman, 2009). Cowles (1996) pointed out that the MCMC sampling of cutpoints c by using such a uniform distribution is not stable, and she suggested joint sampling of latent data z and cutpoints c: Following Cowles (1996), the subsequent studies (Albert and Chib 2001;Chen and Dey 2000) offered a better way to obtain parameters c to remove the ordering constraint using the one-toone map. In this study, the conditional distribution of c s and z s are obtained jointly, and the sampling of c uses transformations as in Chen and Dey (2000): Therefore, the estimation of c s is transformed to the estimation of j s : Defining j s ¼ j s, 2 , :::, j s, LsÀ2 ð Þ 0 , the prior for j s is assumed to follow a normal distribution: The parameters are estimated by using a set of full conditional distributions (FCDs). Supporting materials A presents details of FCDs of all the parameters. Then, the parameters are estimated by MCMC sampling as follows: Let z s ¼ z 11s , :::, z 1n i s , :::, z m1s , :::, z mn m s ð Þ 0 denotes the vector of the sth element z iks from z ik i ¼ 1, :::, m; k ¼ 1, :::, n i ð Þ , and let z Às denotes the vector removing z s from z: Similarly, let z ik, Às denotes the vector removing z iks from z ik : Furthermore, defining c ¼ c 1 , :::, ð c s , :::, c S Þ 0 , c s ¼ c s, 2 , :::, c s, LsÀ2 ð Þ 0 , and thus c Às denotes the vector removing c s from c: Sample c s , z s jb, d, q, R , V, c Às , z Às , y for s ¼ 1, :::, S as follows: , :::, V s, S Þ, v Às, s ¼ V 1, s , :::, V sÀ1, s , V sþ1, s , :::, V S, s ð Þ 0 , and V Às, Às denotes an S À 1 ð ÞÂ S À 1 ð Þ matrix removing the elements v s, s , v Às, s , and v s, Às from V: Defining h ik ¼ X ik b þ d i , and thus h ik, Às denotes the vector removing h iks from h ik : (b) Sample c s, l s by using the transformations as c s, l s ¼ for 2 < l s < L s À 2, where the samples of j s, l s are generated by using a M-H algorithm as explained in Supporting materials A7.

Simulated dataset
A simulation study is implemented to demonstrate the validity and accuracy of the model. The number of response variables is set as S ¼ 2 (y 1 and y 2 ) for ease of understanding and explanation.
Each response variable is assumed to have q s ¼ 2 explanatory variables. The explanatory variables could be set the same for different responses. In an empirical study, different response variables may need different explanatory variables. Additionally, the model becomes more complex when the explanatory variables are different. Therefore, such variables are set to be different for the two response variables. They are generated using a standard uniform distribution [U 0, 1 ð Þ]. Note that the simulation does not include the intercept.
A key element of a Monte Carlo simulation exercise is the relative contribution of Xb, d, and n to the variation in z: A frequent error is allowing Xb to be too dominant and d and n to be rather unimportant. In such a case, b s is always well estimated. To avoid this phenomenon, the regression parameters are set as b d is also a key element to the variation in z: The vector of regional effects d is generated using Spatial coefficient q is a key parameter in generating d of the MSOP model, and thus q is set as 0.3, 0.5, and 0.7 in different experiments to test the validity and accuracy. R is set the same in different experiments as R ¼ 0:02 0:005 0:005 0:01 : Moreover, the weight matrix W is essential to generate d. In the simulated dataset, the number of regions m is set as 6, 10, 15, and 20, and that of individuals in each region n i is set as 100, 200, and 300 in different experiments to test the estimation effect of different sample sizes. The regions are assumed to be located as shown in Figure 1. The weight matrix is generated based on queen contiguity. For example, region 15 is contiguous with regions 8, 9, 10, 14, 16, and 20. For the cases of m ¼ 6, 10, and 15, the regions are assumed to be located as the first 6, 10, and 15 regions in Figure 1.
In empirical datasets, some ordered discrete responses are concentrated at the midpoints, and some are more highly concentrated at the endpoints. The former is known as balanced (bellshaped) frequency distribution, and the latter is unbalanced frequency distribution. Since there are two responses, the setting of covariance matrix V is extremely important: the greater is the covariance, the stronger is the interaction of the response. Additionally, V contributes greatly to the variation in z because individualistic effects n is generated using the multivariate normal distribution N 0, I n V ð Þ : To verify the performance of the model under these two frequency distributions, V is set as V ¼ 0:1 0:09 0:09 0:1 and V ¼ 1 0:9 0:9 1 in different experiments, where weak interaction 0.09 and small variance 0.1 is set to obtain balanced (bell-shaped) frequency tables; and strong interaction 0.9 and big variance 1 is set to obtain unbalanced frequency tables. In addition, the cutpoints are set as the same in different experiments as follows: Based on the above settings, the vector of latent responses z is generated using Eq. (8), and then the observed variables y iks are generated using (1). There are two settings for regression parameters b, three settings for spatial coefficient q, four settings for the number of regions m, and three settings for the individuals in each region n i , and thus there are 72 experiments for each type of frequency distribution. Tables A1 and A2 in Supporting materials provide the complete balanced and unbalanced frequency tables, respectively.

Simulation results
In the simulation, the following diffuse priors are used: , k min and k max are the minimum and maximum eigenvalues for different W; The posterior results are calculated using MCMC sampling. Tables 1 and 2 give partial simulation results of posterior means and true values as well as the mean squared error (MSE) of b and q for different cases. The Supporting materials Table A3 and A4 provide the complete results of 72 simulation experiments for the balanced and unbalanced cases, respectively.
The following mean squared error (MSE) is used to compare the simulation effect, which is calculated as follows: where b t ð Þ and q t ð Þ represent the final tth simulation results. When using the MCMC method, the number of iterations should be set as high as possible to obtain the convergent value of the parameters. There is an obvious tradeoff between accuracy and computation time, and this tradeoff becomes more pressing when the sample size and number of choice variables increase. For example, when the sample size is 6,000 (m ¼ 20 and n i ¼ 300), the computation time of one experiment for 11,000 iterations is 90,389 s and 133,473 s for the balanced and unbalanced cases, respectively. The computation time of 11,000 iterations is 3.7 times 3,000 iterations, but there are no obvious improvements in simulation results between these two settings, as shown in Tables 1 and 2. For example, MSE (sum of MSE b and MSE q ) only decreases from 0.0583 to 0.0533 for the balanced case (m ¼ 20, n i ¼ 300, and q ¼ 0:5). Finally, each MCMC simulation is run for 3,000 iterations for tradeoff, and the first 1,000 samples are discarded as the burn-in period, because the sampled parameter distributions appear to stabilize after 1,000 iterations. Thus, the parameters are calculated based on the final 2,000 runs. From Tables 1 and 2, for both balanced and unbalanced cases, the MSE are close for the two settings of b true : For example, when q true ¼ 0:3, the average MSE are 0.0797 and 0.0733 for the " ÃÃ " and " Ã " denote that zero is not included in the 95% and 90% credible intervals, respectively. b MSE ¼ MSE b þ MSE q : balanced cases b true ¼ 0:55, 0:35, 0:41, 0:56 ð Þ 0 and b true ¼ 0:40, 0:49, 0:45, 0:52 ð Þ 0 , respectively. Therefore, the MSOP model is stationary to capture different sizes of b: When q true increases, MSE tends to decrease for both balanced and unbalanced cases. For example, when b true ¼ 0:55, 0:35, 0:41, 0:56 ð Þ 0 , the average MSE are 0.1518, 0.1077, and 0.0867 for the unbalanced cases q true ¼ 0:3, q true ¼ 0:5, and q true ¼ 0:7, respectively. Additionally, when q true ¼ 0:3, zero is included in the 95% and 90% credible intervals for most cases; that is, posterior means of q are not evident. On the contrary, when q true ¼ 0:5 or 0.7, the credible intervals of most posterior means of q do not include zero and are close to the true values. Therefore, the MSOP model in this study captures higher spatial dependencies more accurately.
Furthermore, the MSE are the smallest for most cases of m ¼ 20. In addition, for the same number of regions m, the MSE tend to be smaller with the increase in observations n i . That is, the larger the sample size, the higher the accuracy of the estimate.
Additionally, no matter how large q true is, the MSE of balanced frequency distributions are smaller than the MSE of unbalanced frequency distributions for most cases, and the posterior means of q are far from the true values for most unbalanced cases. Therefore, the MSOP model herein is more suitable for ordered discrete responses with bell-shaped frequency distributions. However, the MSOP model is still suitable for unbalanced frequency distributions when the observations are sufficiently large (n i > 300).
Overall, considering the magnitudes of the parameter values, the average estimation results are quite close to the true values. Even for a small number of regions, when the observations in each region are over 300, the MSOP model performs well with the simulated data, especially for balanced cases. It satisfactorily explains the spatial dependences and interaction of multivariate response variables as well as the effects of the explanatory variables.

Empirical study
The economic development level and the supporting infrastructure (e.g., medical, education, health environment, sanitation, etc.) are quite different among provinces in China. However, these conditions are common in the same province and observations in the same province therefore have some common effects. In this empirical study, the response variables are elderly people's self-rated health and life satisfaction. The responses of the observations in one province may have common regional effects, which may be affected by neighboring provinces. Therefore, the MSOP model is suitable for this empirical study.

Data
The data are taken from the 2010 Chinese Urban and Rural Elderly Population Surveys, conducted by the China Research Center on Aging of the National Committee on Aging (Gao and Li 2016). The observations are people aged 60 and over. The survey covered 20 representative provinces, municipalities, and autonomous regions in mainland China, as shown in Figure 2. Heilongjiang and Xinjiang provinces are dropped in this application as this MSOP model is not suitable for islands. After deleting samples with missing data, there are 17,646 valid survey samples. Since the total number of observations in the MSOP model is a limitation, a 1/3 subsample of 5,882 elderly people are drawn from the overall dataset. The subsamples are randomly selected according to the distribution of original samples in different provinces. In addition, the weights w ij are quantified based on queen contiguity. From Figure 2, the normalized weight matrix (row sum is 1) can be easily obtained. Table 3 explains all the variables used in the MSOP model. In the model, the response variables are self-rated health and life satisfaction of the elderly people aged 60 and over. The explanatory variables include gender, age, education (education level), marriage (marital status), census register (urban and rural), chronic disease, and need for care. Here, age is a continuous variable, which has been changed as [age-min(age)]/[max(age)-min(age)], and thus the range of age is from 0 to 1. The other explanatory variables are discrete variables. Table 4 provides the summary statistics of these variables, and shows the frequency tables are bell-shaped.

Posterior results of the parameters
In the empirical study, the priors are assigned the same as in the simulation part. MCMC simulations are run for 3,000 iterations, and the first 1,000 samples are discarded as the burn-in period. Table 5 gives the posterior means, posterior standard deviations, and posterior medians for the two response variables. From Table 5, gender has an evident effect only on self-rated health, and marriage has an evident effect only on life satisfaction. The other explanatory variables have evident effects on both self-rated health and life satisfaction.
Specifically, for self-rated health, the regression parameters of gender and education are positive. Therefore, self-rated health in elderly males and people with high education levels is more likely to be higher than that in elderly females and people with low education levels, respectively. The regression parameters of age, census register, chronic disease, and need for care are negative. Therefore, self-rated health in younger elderly people, urban elderly people, people without chronic disease, and people without need for care is more likely to be higher than that in older elderly people, rural elderly people, people with chronic disease, and people with need for care, respectively. Table 3. Variables in the model.

Variable
Description Value Dependent variables Self-rated health Self-evaluation of the health condition. Very poor is 1; poor is 2; fair is 3; good is 4; and very good is 5. Life satisfaction Satisfaction with life. Very dissatisfied is 1; dissatisfied is 2; fair is 3; satisfied is 4; and very satisfied is 5.

Explanatory variables Gender
Male is 1; female is 0. Age Re-ranged age. Change as: (age-min(age)) /(max(age)-min(age)). Education One's education level. Junior high school and above is 1; else is 0. Marriage Marital status. Live with a spouse is 1; else is 0. Census register Census register with categories 'urban' and 'rural'. Rural is 1; urban l is 0. Chronic disease Whether one has a chronic disease or not. Having a chronic disease is 1, else is 0. Need for care Whether one needs daily care from others or not. Need daily care is 1; do not need is 0. For life satisfaction, the regression parameters of age, education, and marriage are positive. Thus, life satisfaction in older elderly people is more likely to be higher than that in younger elderly people. Life satisfaction in people with high education level and people living with a spouse is more likely to be higher than that in people with low education level and people in other marital statuses, respectively. The regression parameters of census register, chronic disease, and need for care are negative. Thus, life satisfaction in urban elderly people, people without chronic disease, and people without the need for care is more likely to be higher than that in rural elderly people, people with chronic disease, and people with the need for care, respectively.
In this empirical study, d is a 36 m ¼ 18, S ¼ 2 ð Þ -dimensional vector of the regional effects. Subvectors d i1 and d i2 represent all unobserved common regional effects for the observations within province i i ¼ 1, :::, 18 ð Þ , that is, for self-rated health and life satisfaction, respectively. In Bayesian analysis, the importance of a parameter is not reflected in the variation, but its credible interval. From Table 5, although the values do not vary significantly, the credible intervals of regional effects d is do not include zero; that is, the regional effects are indispensable on the response variables. Note that the posterior means of d is are not 0 and vary across i (province) because they depend not only on the common q and R , but also on individual's x iks , b s , and z iks .
The posterior mean of spatial coefficient q is 0.8875. It is large, and the credible interval does not include zero. Therefore, the spatial dependences in this study are critical. Furthermore, one province's regional effect has a positive evident effect on the other provinces' regional effects.
The covariance of the individualistic effects for self-rated health n ik1 and life satisfaction n ik2 is captured by V 1, 2 or V 2, 1 , the credible intervals of which do not include zero; that is, the covariance is evident. Because the covariance is positive, the individualistic effects for self-rated health and life satisfaction are positively related.
c 1 and c 2 represent the vectors of cutpoints for self-rated health and life satisfaction, respectively. For example, if 0 < z ik1 0:2804, self-rated health is "poor;" and if z ik2 > 1, life satisfaction is "very satisfied."

Gain of the multivariate spatial ordered probit model
In this empirical study, the parameters could also be estimated individually and separately for each response using a univariate spatial ordered probit (USOP) model as explained in Supporting materials B1. The parameters in the USOP model for self-rated health and life satisfaction are estimated separately by MCMC sampling based on the FCDs as explained in Supporting materials B2.
The estimated results of parameters for the two responses in the USOP models are listed in Table 6. Although the the posterior means of the USOP models and the MSOP model are quite close, USOP is not a special case of MSOP.
(a) q s is allowed to differ between outcomes in the USOP models but not in the MSOP model.
As it turns out, the two q s values are very similar in the USOP models, at 0.8804 and 0.9052. (b) The USOP is not able to estimate r 1, 2 : The MSOP estimate equals 0.0058, implying a correlation coefficient of 0.8744, which is humongous. (c) The USOP is not able to estimate V 1, 2 : The MSOP estimate equals 0.0207, implying a correlation coefficient of 0.2568, which is not humongous but not at all trivial either. (d) If q s was allowed to vary across s, MSOP would nest USOP by setting r 1, 2 and V 12 equal to 0.
Although the parameters are close, accuracy of prediction is different for these two types of models. Here, the correctly predicted percentage is applied to explain the gain of the MSOP model compared with USOP models. The calculation steps for the MSOP model are as follows: 5259 ÃÃ a CR, CD, and NC are the abbreviation of self-rated health, Census Register, Chronic Disease, and Need for Care, respectively. b " ÃÃ " and " Ã " denote that zero is not included in the 95% and 90% credible intervals, respectively. Mean, SD, and Median are posterior mean, posterior standard deviation, and posterior median, respectively.
(1) Obtain the fitted values of c z iks using the posterior means of parameters through (3).
Learning from the CPRs, the accuracy of prediction for the endpoint in the USOP model is very poor. Therefore, the MSOP model is more suitable for correlated ordered responses.

Posterior results of marginal effects
If x iksc c ¼ 1, :::, q s ð Þ is a continuous variable, the marginal effect is the average change in @Pr y iks ¼l s jb, d, q, R , V, c s , z Às ð Þ @x iksc : If x iksc c ¼ 1, :::, q s ð Þis a discrete variable, the marginal effect is the average change in Pr y iks ¼ l s jb, d, q, R , V, c s , z Às À Á when x iksc of each separate individual goes from 0 to 1. Supporting materials C presents the detailed derivation of the marginal effects of the CR, CD, and NC are the abbreviation of self-rated health, Census Register, Chronic Disease, and Need for Care, respectively. b " ÃÃ " and " Ã " denote that zero is not included in the 95% and 90% credible intervals, respectively. Mean, SD, and Median are posterior mean, posterior standard deviation, and posterior median, respectively.
continuous and discrete variables, and the marginal effects of all the parameters for the two response variables are listed in Table A5. Table 7 presents the marginal effects of the median responses of the explanatory variables. For example, when the individual goes from "do not have a chronic disease" to "have a chronic disease," the probability of life satisfaction being "satisfied" diminishes 1.58%. Note that the actual marginal effect of age should be divided by 40 because it is rescaled from 0 to 1, while the actual age is from 60 to 100.
The median of the two responses are "self-rated health ¼ 3" and "life satisfaction ¼ 4." That is, for most elderly people, their self-rated health is "fair," and life satisfaction is "satisfied." For people's self-rated health to be "fair," the marginal effect of "need for care" is the greatest. Specifically, the probability of self-rated health being "fair" in elderly people without the need for care is 14.10% greater than for those with the need for care. For people's life satisfaction to be "satisfied," the marginal effect of "need for care" is also the greatest. Specifically, the probability of life satisfaction being "satisfied" in elderly people without the need for care is 6.62% greater than for those with the need for care.

Conclusions
This study proposed a new algorithm for a MSOP model. In applying this model, the parameters are calculated using the Bayesian inference based on MCMC sampling. The parameters are estimated using a set of full conditional distributions. To obtain the full conditional posterior distributions, this study assumed diffuse prior distributions for all the parameters. The posterior distributions of all the other parameters are similar to those in former studies on Bayesian spatial econometrics (Kakamu and Wago 2007;Smith and LeSage 2004;Wang and Kockelman 2009) except for cutpoints c: The posterior distributions of c in this study are referred to the method in Chen and Dey (2000). This method removed the ordering constraint using the one-to-one map, and it has been proven to be more stable than using the uniform distribution.
The validity and accuracy of the model are verified by the simulated dataset. The estimation results are quite close to the true values. Therefore, the MSOP model performs very well with the simulated data. In addition, the study illustrated the model by applying it to the survey data of elderly people in 18 representative provinces, municipalities, and autonomous regions in mainland China. The important empirical findings are as follows: (1) The spatial coefficient q is large, and the credible interval does not include zero, and thus the spatial dependences in this study are critical. Furthermore, one province's regional effect has a positive effect on the other provinces' regional effects. (2) Since the covariance V 1, 2 is evident and positive, the individualistic effects for self-rated health and life satisfaction are positively related. (3) The five explanatory variables (age, education, census register, chronic disease, and need for care) have evident effects on both response variables; gender only has an evident effect on self-rated health, and marriage only has an evident effect on life satisfaction. (4) The MSOP model is more suitable than separate USOP for Care, respectively. b " ÃÃ " and " Ã " denote that zero is not included in the 95% and 90% credible intervals, respectively. Mean, SD, and Median are posterior mean, posterior standard deviation, and posterior median, respectively. models for correlated ordered responses. (5) The marginal effect of "need for care" is the greatest for the most elderly people whose self-rated health is "fair" and life satisfaction is "satisfied." In this study, the regional effects are captured for multivariate ordinal responses. Nevertheless, there are still some remaining issues and future research avenues. The observations in region i share a common regional effect d i that is spatially connected with neighboring regions. However, intra-regional interaction also exists. The following model may be used to capture the inter-and intra-regional interactions simultaneously for the multivariate ordinal responses: where q 1 is the intra-regional spatial coefficient, weights w 1 kr reflect the closeness between individuals k and r inside region i; q 2 is the inter-regional spatial coefficient, weights w 2 ij reflect the closeness between regions i and j. Moreover, when the idiosyncrasy should be distinguished for multivariate responses, q 1 and q 2 should be outcome-specific and thus restated as q 1s and q 2s , respectively.