International Journal of Health Geographics a Model for Spatial Variations in Life Expectancy; Mortality in Chinese Regions in 2000

Background: Life expectancy in China has been improving markedly but health gains have been uneven and there is inequality in survival chances between regions and in rural as against urban areas. This paper applies a statistical modelling approach to mortality data collected in conjunction with the 2000 Census to formally assess spatial mortality contrasts in China. The modelling approach provides interpretable summary parameters (e.g. the relative mortality risk in rural as against urban areas) and is more parsimonious in terms of parameters than the conventional life table model.


Background
This paper develops a model for mortality contrasts between 31 administrative divisions of China. Despite major gains in life expectancy in China and improved health and living standards, there is evidence of considerable social and health inequality [1], and poor access to health care in rural and less developed areas [2]. As a result, health and life expectancy improvements have been uneven, typically greater in developed eastern and southern areas of China.
Assessment of age-related inequalities in mortality is therefore important, and life table analysis is often used to this end. However, conventional life table analysis proceeds by independent analysis of each region, often by spreadsheet; for example, see the Population Analysis Spreadsheets developed by the US Census Bureau [3]. This approach to life table analysis takes no account of spatial structure in mortality risks due to socio-economic and environmental factors that affect neighbouring regions similarly. Life table analysis most commonly also adopts (albeit often implicitly) a saturated fixed effects model that does not model correlation in mortality rates between adjacent age groups, and is heavily parameterised.
Here the modelling approach pools strength over ages and regions using random effects methods that recognize correlations between adjacent ages and areas. The number of parameters involved in the model becomes an unknown, often referred to as an effective parameter total, and with a suitably designed model this total may be considerably less than the parameters needed under fixed effects approaches [4].
On the other hand, it is possible that a mortality model oversimplifies such that predictions from it do not accurately reproduce the observed mortality data. Many models in the disease mapping literature assume multiplicative age and area effects, under which area effects on mortality relative risk are constant across age groups [5,6]. The modelling approach here, while relatively parsimonious, explicitly recognizes that age and area effects in observed mortality schedules may not be multiplicative (or additive in the logs of mortality risk), in which case predictive discrepancies will be apparent when a multiplicative model is used. Instead an interaction scheme between age and area dimensions is proposed by adapting the bilinear interaction scheme of Lee and Carter [7], who consider mortality forecasts in age and time dimensions.

Data sources
As mentioned by Banister & Hill [8] the most complete mortality data for China to date have been gathered from all households in the three most recent nationwide censuses of 1982, 1990, and 2000. Over a yearly period ending at the census date, information about household deaths was collected by single year of age and sex. These are population wide collections of death details parallel to the full population enumeration provided by the census. Here the focus is on mortality data over a calendar year (1.11.1999 to 31.10.2000) collected in conjunction with the 5 th population census of China conducted on November 1, 2000. The mortality data can be disaggregated to a relatively low geographic level [9], but here the framework is the 31 administrative divisions of China. These divisions are a mix of 22 provinces, 5 autonomous regions and 4 direct-controlled municipalities (Beijing, Tianjin, Shanghai, and Chongqing), though with status equal to that of the provinces [ Populations and deaths in towns and cities within each of the 31 divisions are amalgamated to provide a simplified comparison of urban and rural mortality.
While the deaths data collected in conjunction with recent China censuses are the most complete available, there is still under-recording. A full discussion of death undercounting is provided by Bannister & Hill [8] who obtained correction factors using the general growth balance method [11]. In deriving mortality rates and life expectancies in this paper, adjustment factors from the paper by Banister and Hill are applied to correct for mortality under-recording. Specifically, populations at risk are reduced to correspond to recorded deaths; thus the male mortality adjustment factor was 1.113, so recorded deaths were retained as the response variables but analyzed in relation to census populations scaled by 1/1.113 = 0.898; for females the adjustment factor is 1.181 so female populations are scaled by 0.847. As reported below (Table 5) China wide life table estimates from the models used here are close to those from Banister and Hill that correct for death undercount.
So the observed data consists of deaths and (scaled) populations differentiated by age (21 groups from 0-4,5-9 through to 100+), by gender, and by an urban-rural categorisation of populations within each administrative division. For simplicity we refer to the administrative divisions simply as divisions below. To provide some background on the varying socio-economic character of these divisions, Appendix 1 tabulates comparative data on indices of literacy, rurality, income and employment. The wide contrasts in development and living standards are relevant to interpreting mortality differences.
Other analyses of the 2000 census mortality data are provided by Heilig [12]. Lai [13], and Cai [9], while Lai et al [14] provide analysis of corresponding data for 1990. These analyses adopt conventional fixed effects approaches, in contrast to the comprehensive approach to smoothing over areas and ages as obtained by the random effects methods used in this paper. Whereas a fixed effects model of conventional life table analysis would involve 31 × 21 × 2 × 2 = 2604 parameters when applied to the data in this paper, the approach developed here is shown to be considerably more parsimonious, while also reproducing the data accurately.
With repeated sampling from the posterior density of parameters using Monte Carlo Markov Chain (MCMC) methods, the stochastic variation in life table parameters to be readily obtained: for example, one may obtain 95% intervals for division life expectancies or more specific outputs such as the difference in urban as against rural life expectancy within each division.

The life table model
Let i denote division, x denote age, s denote gender (1 = females, 2 = males) and r denote urban vs rural populations in each division (r = 1 for urban, 2 for rural). The data are deaths y risx and scaled populations P risx . Because of the aggregated spatial scale and relatively large death counts involved, the life table model assumes deaths y risx to follow a negative binomial density, namely where ξ risx represents the modelled mortality count. The The initial multiplicative model 1 for mortality counts ξ risx involves a) an overall mortality level parameter by gender s, κ s ; b) parameters η xs to represent the age-sex mortality rates for age × and gender s; these are assumed to follow a first order random walk that reflects correlation in rates between successive ages; c) parameters γ s to represent a gender specific rural mortality differential; d) division effects φ is to represent spatially correlated mortality contrasts that are likely to be gender differentiated.
The age parameters η xs describe the typical age profile of mortality: relatively high child mortality, followed by low mortality for older children and young adults, and then rising at older ages. The age parameters in model 1 apply across all provinces in line with the multiplicative model -the validity of which the analysis here is seeking to assess. For example, it may well be that some degree of regional variation in age effects is in fact present in the Chinese mortality data, and model elaborations explained below allow for this. The final set of parameters φ is reflect unmeasured risk factors for (or influences on) mortality that are themselves spatially patterned [16]. For example, differences in environment, health care, climate, economic development and so on are likely to be spatially correlated.
Then model 1 specifies log(ξ 1isx ) = log(P 1isx ) + κ s + η xs + φ is where P risx are scaled populations from the 2000 China Census. Detailed assumptions about the age and spatial effects (in terms of prior densities) are considered in Appendix 2.
In model 1 the age effects η xs are assumed to be independent of area (division) effects, in line with the widely applied multiplicative model [5]. Multiplicativity refers to the original mortality risk scale, when all effects are exponentiated in (3). The multiplicative assumption leads to a parsimonious model but the actual mortality pattern may not conform to the simplifying assumption of a uniform age gradient through all divisions. For example, it may be that age related discrepancies from the multiplicative model occur, in line with the overall mortality level in a division. So under a multiplicative model, infant and child mortality may be underpredicted in relatively backward high mortality divisions and overpredicted in the more developed lower mortality divisions.
To reflect such possibilities a multiplicative interaction between area and age is introduced in the log relative risk scale. This provides a more general though still relatively parsimonious model (model 2), namely where the θ is represent area mortality contrasts, especially for particular age groups. In fact, for parameter identification the θ is sum to zero, while the φ xs sum to 1, namely . Higher weights for a particular age express the fact that these ages are particularly subject to the spatial variation expressed in θ is . For example, divisions with high θ is may tend to have high child or young adult mortality (this anticipates findings later in the paper), and so the age-sex φ xs parameters will be higher at such ages.
An analogous scheme, but in an age-time rather than agearea context, is the Lee-Carter model used for mortality forecasting [7]. For death counts y xt with negative binomial or Poisson means ξ xt the Lee-Carter model would involve a log-bilinear form where age weights φ x are highest for those age groups showing most improvement according to mortality trend parameters θ t which are typically correlated in time.
The overall parameter combination φ xs θ is in equation (4) provides a relatively parsimonious representation of agesex-province mortality effects involving 104 parameters, that avoids using 31*21*2 = 1302 age-area interaction parameters ψ xis . If model discrepancies remain despite the generalisation in equation 4 (e.g. in predictions matching observations), one might introduce ψ xis , but this is likely to be at the cost of model parsimony if applied across all ages. The option considered here is to selectively add area effects ψ xis for age bands x according to an age level binary indicator δ x .
So with a low prior probability π δ (e.g. π δ = 0.05) that δ x = 1, additional area-gender effects specific to the selected age are added. This might be expected to occur particularly for any ages where the mechanism in model 2 still leaves discrepancies in fit. The potential additional effects ψ xis are assumed to be normally distributed random effects with mean zero and an age specific variance parameter ζ x , namely ψ xis ~ N(0, ζ x ). So for model 3, one has δ x B ern(π δ ) and

Findings
Models are estimated using the WINBUGS program [17]. Methods of assessing model fit and adequacy are discussed in Appendix 3, and model fit is summarised in Table 1. For all models, two chains were run for 5000 iterations with convergence by 1000 using Gelman-Rubin scale reduction factors [18]. The Deviance Information Criterion (DIC) -see Appendix 3-under the multiplicative model 1 is 35860 with 99 estimated parameters, while the log of the pseudo marginal likelihood is -17934 (Table 1). Just over 5% of the cases are not included in the 95% intervals of predictions y rep,risx . So in overall terms the model is producing predictions that are concordant with the data. On the other hand, an examination of the pattern of model discrepancies at the age group level ( Table  2) shows that a relatively high number of death counts for ages 0-4 and the oldest ages are not predicted satisfactorily; note that there are 2604/21 = 124 observations in each age group so under model 1, 38 of the 124 death counts in age group are not satisfactorily predicted.
Application of model 2, as in equation (4), results in considerably improved global fit measures such as the DIC and log(psML). The estimated model dimension d e is actually lower than model 1, possibly indicating a higher degree of pooling strength under this model [15]. There is also improved fit in the lowest age band with the average log(CPO) in the 0-4 age band falling from -9 to -8.3 (Table 2) and the number death counts in 0-4 age group not predicted well falling to 30. On the other hand, deaths at ages over 90 still show discrepant predictions.
The division parameters θ is in equation 4 (see Table 3) are highest in the southwestern divisions of Yunnan, Qinghai, Tibet and Guizhou (relatively undeveloped in economic terms, with below average income per head -see Appendix 1), and lowest in developed divisions such as Beijing, Shanghai and Tianjin. Figure 1 plots shows the male and female age weights φ xs ; these show that child and young adult mortality is most strongly linked to spatial  (Figure 2) show the expected pattern of increase with age.
In an attempt to reduce model discrepancies further, model 3 is applied with Pr(δx = 1) = πδ = 0.05. This model produces a further improvement in the global fit measures and reduces predictive discrepancies at age group level ( Table 2). As can be seen from Table 1, the better fit comes at the cost of a much higher effective parameter count with de = 346.
One may expect from the discrepancy pattern for model 2, that the need for additional age-area interactions might be greater at the youngest and oldest ages. In fact, there is unequivocal evidence that ψ xis are needed at ages 0-4 and for the last three age bands (90-95, 95-99, 100+) where the posterior probability for inclusion is Pr(δ x = 1|y) = 1. That is, age-area interactions are retained at these ages in all MCMC iterations. For the 85-90 age band, Pr(δ x = 1|y) = 0.95. Ten of the remaining probabilities Pr(δ x = 1|y) are under 0.05, and the remainder range between 0.05 and 0.28. So age-specific area effects are only necessary conclusively for a minority of age groups.
In order to summarise age-sex effects from model 3 in tabular form, the 31 administrative divisions are aggregated to three broad zones [19]: the western zone (the six provinces of Shaanxi, Gansu, Qinghai, Sichuan, Yunnan and Guizhou, the three autonomous regions of Ningxia, Xinjiang and Tibet, and the Chongqing municipality); the eastern and coastal zone (Beijing, Tianjin, Hebei, Liaoning, Shanghai, Jiangsu, Zhejiang, Fujian, Shandong, Guangdong, Guangxi, and Hainan), and the middle zone (Shanxi, Inner Mongolia, Jilin, Heilongjiang, Anhui, Jiangxi, Henan, Hubei, Hunan). Table 4 shows the death rates (posterior means) m zsx = ξ zsx /P zsx over the three zones z. Figures 3a and 3b plot the logs of these rates. The latter show the wider contrasts in mortality at lower ages, with excess child and younger adult mortality in the Western zone. Also apparent from Table 4 is higher female than male child mortality (at ages 0-4), as reported elsewhere [20].

Findings on expectancy variations
Of particular interest as model outputs are the expectancy gap between urban and rural areas, and the contrasts in life expectancy between divisions at different stages of development. As to urban-rural differentials within divisions, the parameter γ s is estimated in model 3 as 0.48 (with standard deviation 0.02) for females, translating into an average rural mortality relative risk of 1.62 with urban areas at 1. For males the corresponding parameter is 0.37 (with relative risk 1.45). So females are particularly disadvantaged by rural location. The rural health disadvantage in modern China has been attributed to the disproportionate subsidies given to urban hospitals, and to the fact that income inequalities between rural and urban areas have been rapidly followed by health inequalities [21]. China's health system has been ranked by the WHO as the lowest in the world in terms of health equity: urban residents make up only about 20% of China's total population, but enjoy about 80% of health resources [22]. Table 5 shows the profile of life expectancies at birth resulting from model 3. These are by administrative division and urban-rural populations within divisions. Also shown are adjusted estimates of life expectancy presented both by Heilig [12] and Cai [9], that were made by China's National Bureau of Statistics [23]. The China wide life expectancies for males and females according to the NBS are respectively 69.6 and 73.3. Cai [9] mentions that "although the NBS did not provide information on how the adjustments were made, the comparison indicates that adjustments were made for the undercounting of deaths in the census". The China wide life expectancies for males and females from model 3 are respectively 69.7 and 72.7. These are very close to estimates made by Banister & Hill [8], that also adjust for death undercounting, namely 69.7 and 72.8.  urban life expectancy densities found 12 with significant positive kurtosis (and none with significant negative kurtosis). So assessing whether life expectancy at birth E 0i in division i exceeded that in division j might be problematic under classical analysis, but under a Bayesian sampling approach the posterior probability Pr(E 0i > E 0j |y) is easily obtained. Table 5 shows generally higher expectancies for women; men have considerably higher tobacco and alcohol consumption than women in China [24,25], and are involved in more traffic accidents and work-related diseases and deaths than women [26].      -19  20-24  25-29  30-34  35-39  40-44  45-49  50-54  55-59  60-64  65-69  70-74  75-79  80-84  85-89  90-  ing that rural populations within healthier divisions benefit in improved survival terms as well as urban population groups.

Conclusion
This paper has sought to show the gains of using a model based approach to life tables for a set of regions. It has considered mortality in 31 Chinese administrative divisions and shown how the major aspects of mortality variation can be summarised in parsimonious models with interpretable parameters. In terms of geographic contrasts the main results are the wide inter-division and urbanrural contrasts in mortality risks and life expectancy.
The modelling approach allows the significance of such contrasts to be readily assessed -and indeed the full posterior sampling density of life The model schemes used have shown that area and age effects do not conform to the multiplicative model considered by Hoem [5] and others -the multiplicative model is equivalent to assuming effects of age and area are additive in the log mortality risk scale. Instead, it is found that spatial contrasts in Chinese division mortality are especially apparent in younger age groups and model fit (for models 2 and 3) is greatly improved by allowing for such non-proportionality.
Much of the disease mapping literature, especially that adopting a Bayesian modelling approach, has assumed additive effects in the log relative risk scale, or has eliminated consideration of age (or other demographic stratifiers) using indirect standardisation, often without prior checks on the validity of such simplifications. The indirect standardisation approach [6] applies a standard age schedule (e.g. national age mortality rates) across all areas and leads to models where event counts y i for area i are compared to expected deaths E i which are treated as an offset.
The present research emphasizes the importance of evaluating such basic assumptions in making epidemiological or health inferences, since relative risk contrasts between areas may not be the same when age-area interactions are allowed for. If mortality and health inequalities between divisions or urban-rural areas are strongly age related then this may be important in strategic health interventions to counteract health imbalances.
While having implications for health policy, the type of analysis undertaken in this paper is descriptive in the sense that one is not seeking to explain life expectancy contrasts in terms of a range of potential risk factors: health behaviours, environmental or climatic factors, income levels or income inequality, ethnicity and/or nationality, education and literacy levels, female status, or economic structure [28,29]. Rather a statistical modelling approach accounting for age and area effects can be seen as a preliminary to, but also providing a distinct perspective to, a causal model that would typically just examine total life expectancy as a response.

Appendix 2 methods: details on priors for models
A fully Bayesian strategy is adopted using Monte Carlo Markov Chain (MCMC) estimation. For bivariate spatial effects, namely λ i = (λ i1 , λ i2 ) in model 1 and θ i = (θ i1 , θ i2 ) in models 2 & 3, a bivariate version of the CAR normal prior [30] is assumed. Thus with n = 31 provinces, and with i and j denoting different provinces, the pairwise difference prior for λ i = (λ i1 , λ i2 ) has precision matrix Ψ λ and joint density A contiguity assumption is made for the geographic interactions so that w ij = 1 for adjacent provinces and zero otherwise. A Wishart prior with identity scale matrix and 2 degrees of freedom is adopted on Ψ λ . The spatial effects are centred at each MCMC iteration (to sum to zero over all provinces).
The age effects η xs in all three models are modelled in terms of normal autoregressive random walks with variances ω s specific to sex, so This prior does not set a level and so male and female age effects are centred at each MCMC iteration (using the car.normal prior in WINBUGS).
A uniform U(0,1000) prior is assumed for the negative binomial α, and N(0,1000) priors on the fixed effects {k s , γ s }. For the precisions ω s (in all models) and ζ x (in model 3), a uniform U(0,1000) prior is assumed. The age weight parameters φ xs in models 2 and 3 are generated using gamma Ga(l,l) priors on parameters f xs , and then scaling the f xs to sum to 1 within genders s. So This is equivalent to a Dirichlet prior on the φ xs 17 .

Appendix 3 methods: model assessment
Comparisons of model fit in the analysis use two criteria: the deviance information criterion (DIC) of Spiegelhalter et al [31], and the pseudo marginal likelihood (psML) based on Monte Carlo estimates of conditional predictive ordinates, or CPOs [32,33]. The DIC is obtained as the posterior mean deviance plus a measure of complexity d e that can also be regarded as estimating the effective dimension of the model, but also reflects features such as parameter precision.
The CPO is a cross-validatory measure of predictive fit, namely  where θ are the model parameters, y h indicates deaths in a particular one of the 2604 subdivision-area-gender-age strata, and y [h] denotes the data in the remaining 2603 strata. The conditional predictive ordinates are indicators of model adequacy at the level of individual observations; low CPOs (or highly negative log CPOs) indicate poorly fitted points. The overall total of log(CPO) provides a log(psML) which will be higher for a model providing a better overall fit.
Model adequacy may also be assessed using predictions from the model -more specifically, replicate observations y rep,risx sampled from the posterior predictive density p(y rep |y). Following Gelfand [34], a summary of how well these predictions match the actual data, y risx , involves a tally of how many actual observations are located within the 95% interval of the corresponding model prediction y rep,risx . If 95% or more of the N = 2 × 31 × 2 × 21 = 2604 death counts y risx are within the 95% intervals of the pre-dictions then the model is judged to be reproducing the observations satisfactorily.