Comparison of Bayesian Spatio-Temporal Models for Chronic Diseases

This paper discusses a comprehensive statistical approach that will be useful in answering health-related questions concerning mortality and incidence rates of chronic diseases such as cancer and hypertension. The developed spatio-temporal models will be useful to explain the patterns of mortality rates of chronic disease in terms of environmental changes and social-economic conditions. In addition to age and time effects, models include two components of normally distributed residual effects and spatial effects, one to represent average regional effects and another to represent changes of subgroups within region over time. Numerical analysis is based on male lung cancer mortality data from the state of Missouri. Gibbs sampling is used to obtain the posterior quantities. As a result, all models discussed in this article fit well in stabilizing the mortality rates, especially in the less populated areas. Due to the richness of hierarchical settings, easy interpretation of parameters and ease of implementation, any models proposed in this paper can be applied generally to other sets of data.


Introduction
The objective of this paper is to develop a comprehensive statistical approach that will be useful in answering health-related questions concerning mortality and incidence rates of chronic diseases such as cancer and hypertension.The methods will employ hierarchical models which include demographic variable (i.e., age) and geographic variable to identify regional differences.The proposed models will also include a longitudinal variable that will describe temporal trends in mortality and incidences for different population groups.The developed spatiotemporal models will be useful to explain the patterns of mortality rates of chronic disease in terms of environmental changes and social-economic conditions.
Within the last two decades, many authors have proposed fully (hierarchical) Bayesian or empirical Bayes estimates of mortality rates with spatial correlation, along with Binomial or Poisson sampling processes.This paper is closely related to the general development in this area, which is typified, among others, by Tsutakawa (1985Tsutakawa ( , 1988)), Clayton and Kaldor (1987), Ghosh et al. (1998), Best et al. (1998), Sun et al. (2000), and Kim and Oleson (2008).These studies are based on hierarchical models that use a Poisson model for the first stage and incorporate covariates by various modeling of the Poisson parameters at the second stage.Spatial (Regional) effects are usually treated as random with some distributions whose parameters must be estimated.Hyper-parameters are typically estimated by maximum likelihood methods, and the Poisson parameters and rates are estimated by Bayesian or empirical Bayes method.Tsutakawa (1985) shows that the empirical Bayes method could seriously underestimate the uncertainty in the estimated rates.Ghosh et al. (1998) and Best et al. (1998) consider special cases of pairwise difference prior (cf.Besag, York and Mollié, 1991) for spatial effects without specifying the degree of spatial correlation.Waller et al. (1997) and Knorr-Held (2000) propose spatio-temporal interaction models where the spatial effects are nested within time so that they can examine the evolution of heterogeneity and spatial patterns over time.Sun et al. (2000) and Kim and Oleson (2008) assume normal distributions for residual effects which account for the extra-Poisson variations.However, our specification of simultaneously autoregressive (SAR) and conditional autoregressive (CAR) priors includes parameters to represent the degree of spatial correlation among regional areas while we assume both normal and gamma distributions for residual effects based on the choice of the first-stage model.
In this paper, we will modify the spatio-temporal models proposed by Sun et al. (2000).Their log-linear mixed model of mortality rates includes two components of normally distributed residual effects and spatial effects, one to represent average regional effects and another to represent changes of subgroups within region over time.Our modifications are made in two directions: first, replace the normally distributed residual effects by those which have log-gamma distributions; second, replace the CAR (Besag, 1974) effects by those which are SAR effects (Ord, 1975;Cliff and Ord, 1981).In order to examine the effect of these modifications we also consider making changes in one direction at a time, giving us a total of four possible models.Those four models are compared using Missouri lung cancer data.
The developed methods will not only be applicable to data obtained from other states or countries, but also be useful in analyzing mortality due to other chronic diseases.Moreover, the methods can be applied more generally to measure problems such as crime rates, accident rates, and high school dropouts rates.
In addition, this study motivates more general research in statistics such as model selection, formulation of prior distributions for hierarchical models, and designs for computer simulation.
The subsequent sections are in the following orders.In Section 2 we introduce multiplicative and log-linear mixed models with residual effects for mortality rates where, for a given period of time, age effects are fixed and spatial effects are random.A CAR model by Besag (1974) and an SAR model by Ord (1975) are chosen to describe the random spatial effects.The necessary prior distributions for the parameters and hyper-parameters are given.In Section 3, we derive the available full conditional distributions for the four proposed models.Our choice of SAR and CAR models allows us to sample all regional effects simultaneously without computing the inverses of matrices each time, which improves the efficiency of the computation greatly.Numerical results on age effects, spatial effects, mortality rates, variance components, spatial correlations and residuals are presented in Section 4. The convergence of Gibbs sampling is investigated using Gelman and Rubin's (1992) diagnostics.Also, model selection using cross-validation predictive density and deviance information criterion is done and followed by conclusions in Section 5.

Multiplicative and log-linear mixed models
The Poisson distribution arises naturally when data set takes the form of counts; for instance, a major area of application is epidemiology, the study of incidence of diseases.Let y ijk denote the frequencies of deaths due to a specific cause for the kth time period within the jth age group of the ith county.Each cell has its population size n ijk .We assume that y ijk have independent Poisson distributions with means n ijk p ijk .
In order to model the effects of county, age and time on p ijk , first we consider the multiplicative model (Tsutakawa, 1988;Kim et al., 2002) having the form, where t k are the midpoints of the kth time period and t i are the effects of the ith county, θ * j the effects of the jth age group, µ * j the overall rates of change over time for the jth age group, W * i the rates of change in the ith county over time, and * ijk unexplained residual effects for the ijkth cell.We may express this multiplicative model in the log-linear form, where Z i , θ j , µ j , W i and ijk are natural logarithms of Z * i , θ * j , µ * j , W * i and * ijk , respectively.
In the specification of the model we treat θ * j and µ * j (or equivalently, θ j and µ j ) as fixed effects while Z * i , W * i and * ijk (or equivalently, Z i , W i and ijk ) are treated as random effects.The parameterization will be made unique by placing restrictions on the distributions of the random effects.Depending on the specification of the distributions for the random effects, several cases can be constructed.We focus on two distributions for residual effects ijk and two distributions for the spatial effects (Z i , W i ), resulting in four distinct combinations.These are described below.

Distributions of residual effects
To estimate the rates of chronic diseases such as lung cancer, the extra variation (residual) arises from socio-economic status, air pollution, smoking prevalence, and many other factors not included in the model.These factors are unknown or hard to quantify in most instances.The role of the residual effects is not only to account for the lack of fit, but also to use it as a tool to detect irregular patterns and changes (Kim et al., 2002).The residual effect will sum up all those county, age group, and period effects which are not spatially diffused (i.e., local).Here we have two distributions for residual effects.
• Normal: Suppose the ijk are independent and identically distributed (iid) normal with mean 0 and variance δ 0 .

Distributions of spatial effects
We represent the spatial correlation among the counties by a CAR model introduced by Besag (1974) and an SAR model by Ord (1975).Although CAR and SAR models have been available for more than 30 years, their relative merits have not been fully investigated.
• SAR: The SAR model is defined by the autoregressive relation where τ i are iid N (0, δ u ).We define the adjacency matrix C with entries c ij = 1 if county i and j are neighbors and c ij = 0 otherwise, with the diagonal entries c ii = 0.It can be shown that (I − ρC) is nonsingular when ρ is restricted to some interval depending the eigenvalues of C.Then, it follows that U = (U 1 , U 2 , • • • , U I ) has a multivariate normal distribution with mean 0 and covariance matrix δ u (I − ρC) −2 .
• CAR: The CAR model is defined by the conditional probability density function where δ u > 0 and |ρ| < 1.This model is a modified version of that used by Besag (1974) and Ripley (1981).This modified version was also used by Clayton and Kaldor (1987) and Ripley (1988).Sun, Tsutakawa and Speckman (1999) show that if λ −1 1 < ρ < λ −1 I where λ 1 and λ I are smallest and largest eigenvalues, with where (ρ 1 , ρ 2 ) represent the degrees of spatial correlation between regions and (δ 1 , δ 2 ) the dispersion parameters.Note that the choices of γ are 1 for CAR and 2 for SAR.Contrary to CAR model, it is difficult to write the conditional distribution of (Z i |Z j , j = i) for the SAR model in a higher dimension.Most statisticians prefer to use a CAR model, however, the computational methods for SAR and CAR are almost identical except for the choice of γ.Although it may appear more suitable to work with the gamma distribution under the multiplicative model and the normal distribution under the log-linear model, such pairing is not crucial for Bayesian hierarchical modelings (cf.Kim et al., 2002).

Choice of priors
One objective of specifying priors is to stabilize the estimates of mortality rates by 'smoothing' the crude map.In hierarchical models, the selection of the hyper-prior distribution of dispersion parameter is crucial, in general, when the sample size is not large enough.The crude map using the standardized mortality rates (Mason and McKay, 1974), for each (i, j, k) estimates (usually the maximum likelihood estimation) of p ijk based only on y ijk , often features large outlying mortality rates in less (sparsely) populated areas, so that the map is visually dominated by the rates having the highest uncertainty.One problem with the maximum likelihood is that a zero incidence case in any area provides a zero estimate of mortality rate in that area.This approach also fails to account for the possible similarity of spatial effects in adjacent areas.To choose hyper-parameters for the computation, we begin with proper but very flat priors.For both Poisson-Lognormal and Poisson-Gamma Structures, we define p = (p 111 , p 112 , δ 1 , δ 2 ) and ρ = (ρ 1 , ρ 2 ) where δ and ρ are vectors of variance components and spatial correlations, respectively.For the Poisson-Lognormal structure, we assume that for given (Z, θ, µ, W), the first stage prior distributions of p ijk follow equation (2.2), where ijk are iid N (0, δ 0 ).As mentioned earlier, given (Z, θ, µ, W), v ijk = log(p ijk ) have independent normal distributions with means and variance δ 0 .The following conditional independence assumptions are also needed.
Here the hyper-parameters (ξ mj , δ mj ), (ξ sj , δ sj ) and (a l , b l ) are fixed constants.On the other hand, for the Poisson-Gamma structure, we assume that for given (Z, θ, µ, W), the first stage prior distributions of p ijk follow equation (2.1), where * ijk are iid G(ψ −1 , ψ −1 ).Thus, for given (Z, θ, µ, W), p ijk have independent gamma distributions with shape and scale parameters as } −1 , respectively.The conditional independence assumptions are same as those for Poisson-Lognormal structure except replacing δ 0 with ψ.The following (conditional) prior distributions are assumed.
In order to compare Poisson-Lognormal and Poisson-Gamma structures numerically, we decide to use the values of above hyperparameters which give approximately equal first and second moments in their distributions under the two structures.Detailed values will be given in the following Section.

Estimation via MCMC
For the data set considered in this article, we use the Missouri male lung cancer mortality data.This data set contains 115 counties, 4 age groups (45-54, 55-64, 65-74, 75+), and 4 time periods (1973-1977, 1978-1982, 1983-1987, 1988-1992).There are a few large urban areas and many rural communities in Missouri.The state is geologically diverse and has experienced considerable changes in mortality patterns in the recent decades.The total number of parameters (p, θ, µ, Z, W, δ, ρ) and δ 0 or ψ is equal to 2, 083.We use a Gibbs sampling (Gelfand and Smith, 1990) to obtain the posterior moments and densities.All of the full conditional distributions for a Gibbs sampling and generating random samples are listed in the Appendices A through C. Some of them are standard distributions such as normal (univariate or multivariate) and inverse gamma, while others require sampling from log-concave densities.

Specification of prior distributions
When we choose hyper-parameters for the computation, we decide to use the values of hyper-parameters which give approximately equal first and second moments in their distributions under the two different structures.This will enable us to compare the posterior distributions numerically, somewhat free of the impact of prior selections (cf.Kim et al., 2002).In both structures, the values of hyper-parameters for the variance components are (a l , b l ) = (a ψ , b ψ ) = (3.0,0.1).This is equivalent to assuming that the prior of the reciprocal of the variance components has mean 30 and variance 300.Under Poisson-Lognormal structure, the values of hyper-parameters for age and overall slope of age effects are set to (ξ mj , δ mj ) = (ξ sj , δ sj ) = (−0.5,0.75).These values are comparable to those set (a θ j , b θ j ) = (a µ j , b µ j ) = (3.0,2.0) for Poisson-Gamma structure.It appears that the posterior estimates for each of the four different models are quite robust in terms of the choices of hyper-parameters.

Numerical Results
We define PLS for Poisson-Lognormal-SAR, PLC for Poisson-Lognormal-CAR, PGS for Poisson-Gamma-SAR, and PGC for Poisson-Gamma-CAR models.

Posterior Quantities
A summary of the posterior quantities of the parameters (θ, µ, δ, ρ), δ 0 and ψ for four different models is presented in Table 1.As seen in Table 1, all posterior means except estimates for the spatial correlations ρ, from four different models based on 50, 000 Gibbs cycles, are about the same while the posterior standard deviations for (θ, ρ) disagree.It is interesting to note that posterior standard deviations, except for ρ 1 , from Poisson-Gamma structure are smaller than or equal to those from the Poisson-Lognormal structure.Such results are likely due to the different tails of the gamma distribution with respect to the lognormal distribution.Through the MCMC simulations we have compared the marginal correlations of Z i between adjacent counties and also between non-adjacent counties.We conclude that marginal correlations based on adjacent counties are, on the average, larger than those based on non-adjacent counties.For a given specific county, we have observed that correlations with adjacent counties are consistently larger than correlations with non-adjacent counties.We also approximate the spatial covariance matrix of Z, δ 1 A −γ 1 , using the posterior estimates of (δ 1 , ρ 1 ).We notice that almost all of diagonals and absolute values of off-diagonals of spatial covariance matrix based on Poisson-Gamma structure are less than those based  1977 1978−1982 1983−1987 1988−1992 Figure 3: Plots of annual mortality rates for four selected areas over time period given age group.Solid line is for raw rates; dotted line is for PLC; and dashed line is for PLS.
on Poisson-Lognormal structure.The generalized variances, which are the determinants of the matrices, of spatial covariance matrix of Z indicate that the SAR model gives less variability in Z i since the generalized variances under SAR models are smaller than those under CAR models.
• Mortality rates p: Disease maps show the patterns of possible changes over time and the concentration of counties of high and low rates.Figures 1 and 2 map the frequency and Bayesian posterior estimates (based on the PGC model for the purpose of illustration since other models show the similar patterns) of the rate p ijk , which have been rescaled to represent the annual rates per 10 5 population.In Figure 2, we note that most of the extreme rates in Figure 1 have disappeared.In particular, the isolated cases of high rates among the two lower age groups have disappeared as well as the isolated lower rates among the two upper age groups.It can clearly be seen how the rates have increased over time for the two upper age groups and how they tend to decrease over time for the youngest group.In spite of the general trends, it is worth to note that in each age group there exist several counties where the trend is reversed or in the opposite direction of the general trend.
Figure 3 shows plots of the annual mortality rates for four selected areas.Based on population sizes, we choose the smallest (Worth, county=113), medium (Perry, county=79), largest (St. Louis, county=95) and a special county (St. Louis city, county=115) for the purpose of demonstration.Mortality rates between the Poisson-Gamma and Poisson-Lognormal structures are close enough to each other that a difference is negligible.On the other hand, frequency estimates vary drastically for less populated counties.The Bayesian posterior estimates for densely populated counties are close to the frequency rates with small standard errors, while less populated counties have posterior means around the overall mean rates with larger standard errors.We observe that mortality rates from the Poisson-Gamma structure for less populated counties are relatively larger (toward the overall mean rates) than those from the Poisson-Lognormal structure because of more shrinkage effects.Table 2 gives the summary of the relative errors.The relative errors are defined as |p ijk − p * ijk |/p ijk , where p * ijk is from one model and pijk from another, based on posterior means or standard deviations of the mortality rates.The relative errors support that the same structure generates the similar values of posterior estimates.• Spatial Correlations ρ: Table 1 and density plots (not shown) of the posterior distributions of the ρ indicate the spatial correlation ρ 1 for Z is clearly positive but ρ 2 for W are centered, roughly about zero.We construct 95% Bayesian credible intervals for each of (ρ 1 , ρ 2 ) and find that the interval for ρ 1 excludes zero while that for ρ 2 does not.The positive values of ρ 1 indicate the similarities between neighboring counties.Based on ρ 1 , Poisson-Lognormal structure gives larger posterior means with smaller posterior standard deviations than the Poisson-Gamma structure does.Also, the CAR model gives larger posterior means than does the SAR.
• Variance Components δ: The variability in Z, W and e may be seen in the posterior estimates of their variance δ 1 , δ 2 and δ 0 or ψ in Table 1.The concentration of the posterior distributions away from zero indicates the importance of Z, W and e in the model, though the estimates are small.
• Residual Effects eps: We have computed the ordinary sums of squares due to errors of the fitted models.About 67.00% (SAR) and 67.04% (CAR) of the total variation in log(p ijk ) are explained by the linear fit under the Poisson-Lognormal structure, while 65.59% (SAR) and 65.70% (CAR) of the total variation in pijk under the Poisson-Gamma structure.In other words, about 32.96% ∼ 34.41% of the total variation remains in the residuals, which is due to other factors that are not included in the model.Therefore, the residual effects are needed.The detailed discussion on including residual effects in hierarchical models is in Kim et al. (2002).
• Convergence Diagnostic: In monitoring the convergence of Gibbs sampling, we have run several parallel MCMC chains with different starting values and used graphical monitoring of the chains for a representative subset of the parameters.It is noted that the estimates are quite stable after about 10, 000 iterations.Among 50, 000 cycles of Gibbs sampling, we have chosen to use the last 40, 000 cycles for an estimation purpose following a burn-in of first 10, 000 cycles.

Model selection
In order to compare two competing models M 1 and M 2 that have equal prior probabilities of being selected, the Bayes factor BF 12 , the relative weight of evidence for model M 1 compared with model M 2 , is defined as where f (y|M i ) denotes the marginal density of the data y under M i , i=1, 2 (Dey, Chen and Chang, 1997).Generally, the Bayes factor is the ratio of the posterior odds for M 1 versus M 2 to the prior odds for M 1 versus M 2 .Hence, model selection requires no more than the Bayes factor (Jeffreys, 1961;Kass and Raftery, 1995).Unfortunately in complex hierarchical models, in which parameters generally outnumber observations, the Bayes factor cannot be directly applied (Gelfand and Dey, 1994).Instead we use the alternative methods such as the cross-validation predictive density (Dey et al., 1997) and the deviance information criterion (Spiegelhalter et al., 2002).
A Monte Carlo integration of the predictive density yields the estimated conditional predictive ordinate (CPO).Between competing models M 1 and M 2 , each observation favors a specific model of which the CPO is larger.In Table 3, we present the summaries of model selections based on the estimated CPOs.The values in the third and fourth columns indicate that the number of observations (out of 1, 840 total observations) which favor the models M 1 and M 2 , respectively.The last column shows that the percentage of supporting the model M 1 .The value of 50% implies the data set supports two competing models equally.Poisson-Lognormal structure seems to be the ones as being selected.However, the percentages of supporting (or favoring) M 1 over M 2 are not very far away from 50% mark, so no model is strictly superior to others.The second method we use is based on the Deviance Information Criterion (DIC).Spiegelhalter et al. (2002) compare two or more competing models based on the posterior distributions of the deviance where they identify 'fit' as the posterior mean of the deviance, and 'complexity' (i.e. the effective number of parameters, p D ) as the difference between the posterior mean of the deviance and the deviance based on the posterior means of the parameters.The 'fit' and 'complexity' are then added to form a DIC which may be used for model comparison.These quantities can be obtained easily from an MCMC chain.Detailed discussion can be found in Spiegelhalter et al. (2002) and Best et al. (1998).For our Poisson model, the deviance D(p) is defined (McCullagh and Nelder, 1989) as ] .4, we present the summaries of deviance statistic for competing models.The values of p D show that Poisson-Gamma structure has fewer effective parameters (total number of second and third stage parameters is 243), so it may be regarded as more parsimonious.As a whole, all models are virtually indistinguishable in terms of fittings since the differences in DICs, Ds, and p D s are not large across competing models.

Conclusions
Four spatio-temporal models discussed fit well in stabilizing the mortality rates though differences exist such that the Poisson-Gamma has more emphasis on providing shrinkage.All four procedures lead to provide more smooth and accurate maps of rates, especially in less populated areas.As discussed, the spatial effects for the same Poisson-Gamma or Poisson-Lognormal structure are very similar regardless of the choice of a CAR or an SAR model.On the other hand, those for the different structures under the same CAR or SAR model are different.That is, error structures cause the differences in the spatial effects.Poisson-Gamma structure has more tendency to shrink the posterior estimates toward the overall mean rates than does the Poisson-Lognormal structure.As a whole, there exist more differences due to error structures (Lognormal vs. Gamma) than due to spatial correlations (SAR vs. CAR).During MCMC chains, the computational methods for SAR and CAR are almost identical except the values of γ in the covariance matrix of spatial effects.Our choice of CAR model improves the efficiency of the computation greatly over other forms of CAR models appearing in the literature.All posterior means except estimates for spatial correlations ρ, based on four different models, are about the same.It is worthwhile to note that posterior standard deviations from Poisson-Gamma structure are smaller than those from Poisson-Lognormal structure.
In conclusion, due to the richness of hierarchical settings, easy interpretation of parameters and ease of implementation, any models proposed in this paper can be applied generally to other sets of data.We may also use one model to validate the performance based on other model.

Appendix A. Available full conditional distributions for Poisson-lognormal structure
We now give the full conditional distributions for our hierarchical structures.Note that the Poisson-Lognormal structure shares the same full conditional distributions except for the choice of γ in the covariance structures of spatial effects described in earlier section.This is also true for the Poisson-Gamma structure.
(i) For any (i, j, k), the conditional posterior distribution of v ijk = log(p ijk ) given others (i.e., other parameters) is proportional to where ) .( 5.3) (v) The conditional density of (ρ 1 | others) is proportional to (ix) The conditional density of (ρ 2 | others) is proportional to , where a ijk is given by (5.2).

Appendix B. Available full conditional distributions for Poisson-gamma structure
(i) For any (i, j, k), the conditional distribution of ( } −1 .
(ii) For any j, the conditional distribution of (θ ).Since θ j =log(θ * j ), the Jacobian transformation gives the desired conditional posterior distribution of θ j .
(iii) For any j, the conditional density of ( (vi) Given the eigenvectors of adjacency matrix C, write .

Appendix C. Generating random samples
It is easy to generate random samples from normal (univariate or multivariate) as well as gamma and inverse gamma distributions.If we use the parameter p ijk in the computation of Poisson-Lognormal structure, we can show that the conditional density of p ijk given other parameters is not of a closed form, not even log-concave.Thus, sampling from the conditional density of p ijk given other parameters is not trivial (cf.rejection/acceptance sampling, Metropolis algorithm).However, as shown below, in Poisson-Lognormal structure, v ijk = log(p ijk ) has a log-concave density since v ijk is a monotone (one-to-one) transformation of p ijk .We will also see that the conditional densities of ρ 1 and ρ 2 in all models are log-concave.Additionally, µ j , Z i , W i and η=1/ψ in Poisson-Gamma structure are log-concave.The adaptive rejection sampling method proposed by Gilks and Wild (1992) was used to generate random samples from these log-concave densities.Fact 1 The conditional density of v ijk = log(p ijk ) given in Poisson-Lognormal structure is log-concave.
Proof.Let h(v ijk ) be the exponent in (5.1).The second derivative of h(v ijk ) with respect to v ijk is −n ijk exp(v ijk ) − δ −1 0 , which is negative.
Let Γ be an I × I orthogonal matrix such that the adjacency matrix C can be decomposed as Fact 2 (a) The conditional density of ρ 1 , given in (5.5), is equivalent to (5.17) (b) The conditional density of ρ 2 given in (5.8) is log-concave.
Proof.Part (a) follows from (5.16).For Part (b), the first and second derivatives of h(ρ 1 ) are (5.18) Here, when γ=1 or 2, the second derivative of h(ρ 1 ) is negative.In general, the conditional density of ρ 1 is log-concave if γ ≥ 1.The proof of part (c) follows directly from that of part (b).This proves the results.fact e −µ j (t k − t) − b µ j e −µ j < 0.
For part (b), For part (c), where S ψ and T ψ are defined in equation (5.14).From Bowman and Shenton (1988), we have Hence, For given positive integers (I, J, K) and moderate size of a ψ , the second derivative is negative.This proves the results.

Figure 1 :
Figure 1: Frequency estimates of annual mortality rates per 10 5 population by age group, county and time period.

Figure 2 :
Figure 2: Bayesian posterior estimates of annual mortality rates per 10 5 population by age group, county and time period.

Forp
Poisson-Gamma structure, (a) The conditional density of µ j is log-concave.(b) The conditional density of Z i is log-concave.(c) The conditional density of W i is log-concave.(d) The conditional density of η=1/ψ is log-concave.Proof.For part (a), ijk (t k − t)

Table 1 :
Posterior estimates based on 50, 000 Gibbs The posterior means of θ j increase rapidly with age, except for the two upper groups which have similar values.The posterior means of µ j show a gradual increase with age.The negative value for the youngest group indicates an overall decrease in mortality over the 20 year period.

Table 2 :
Relative errors (RE) between models based on posterior means and standard deviations of mortality rates.

Table 3 :
Summary of model selections based on the estimated CPOs.

Table 4 :
Summary of Deviance Statistic based on 8 different models