==== Front Int J Health Geogr Int J Health Geogr International Journal of Health Geographics 1476-072X BioMed Central London 251 10.1186/s12942-020-00251-z Research Beyond standardized mortality ratios; some uses of smoothed age-specific mortality rates on small areas studies http://orcid.org/0000-0003-4439-2117Perez-Panades Jordi perez_jorpan@gva.es 1 Botella-Rocamora Paloma botella_pal@gva.es 1 Martinez-Beneito Miguel Angel miguel.a.martinez@uv.es 2 1 grid.417564.5Direcció General de Salut Pública i Addiccions, Conselleria de Sanitat Universal i Salut Pública, Avda/Cataluña, 21, 46020 Valencia, Spain 2 grid.5338.d0000 0001 2173 938XDepartament d’Estadística i Investigació Operativa, Universitat de València, C/Dr. Moliner, 50, 46100 Burjassot, Valencia Spain 4 12 2020 4 12 2020 2020 19 5416 8 2020 19 11 2020 © The Author(s) 2020Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.Background Most epidemiological risk indicators strongly depend on the age composition of populations, which makes the direct comparison of raw (unstandardized) indicators misleading because of the different age structures of the spatial units of study. Age-standardized rates (ASR) are a common solution for overcoming this confusing effect. The main drawback of ASRs is that they depend on age-specific rates which, when working with small areas, are often based on very few, or no, observed cases for most age groups. A similar effect occurs with life expectancy at birth and many more epidemiological indicators, which makes standardized mortality ratios (SMR) the omnipresent risk indicator for small areas epidemiologic studies. Methods To deal with this issue, a multivariate smoothing model, the M-model, is proposed in order to fit the age-specific probabilities of death (PoDs) for each spatial unit, which assumes dependence between closer age groups and spatial units. This age–space dependence structure enables information to be transferred between neighboring consecutive age groups and neighboring areas, at the same time, providing more reliable age-specific PoDs estimates. Results Three case studies are presented to illustrate the wide range of applications that smoothed age specific PoDs have in practice . The first case study shows the application of the model to a geographical study of lung cancer mortality in women. This study illustrates the convenience of considering age–space interactions in geographical studies and to explore the different spatial risk patterns shown by the different age groups. Second, the model is also applied to the study of ischaemic heart disease mortality in women in two cities at the census tract level. Smoothed age-standardized rates are derived and compared for the census tracts of both cities, illustrating some advantages of this mortality indicator over traditional SMRs. In the latest case study, the model is applied to estimate smoothed life expectancy (LE), which is the most widely used synthetic indicator for characterizing overall mortality differences when (not so small) spatial units are considered. Conclusion Our age–space model is an appropriate and flexible proposal that provides more reliable estimates of the probabilities of death, which allow the calculation of enhanced epidemiological indicators (smoothed ASR, smoothed LE), thus providing alternatives to traditional SMR-based studies of small areas. Keywords Disease mappingAge-specific smoothingLife expectancyAge-standardized ratesInstituto de Salud Carlos III (ES)PI16/01004Martinez-Beneito Miguel Angel issue-copyright-statement© The Author(s) 2020 ==== Body Background Spatial epidemiology deals with the description and analysis of geographically indexed health data with respect to demographic, environmental, behavioral, socioeconomic, genetic and infectious risk factors [1]. Most of these risk factors, especially those related to mortality, strongly depend on age [2]; therefore, the comparison of raw (unstandardized) rates may be misleading due to differences in the population structure of the units of study, which may be confused with real differences in risk. Age-adjusted (standardized) rates were proposed in the mid-nineteenth century [3] as an answer to this problem, although they are not a completely satisfactory solution. Standardization seeks to remove the effect of having different age distributions for the populations being compared. It produces a single summary index per population, easier to compare than a full set of age-specific rates. The most popular standardization method is the so-called indirect method, which would yield standardized mortality ratios (SMRs) if mortality was the event of interest. SMRs are calculated as the ratio of the number of deaths observed over a specific time interval in a population group, in our case and from now on “spatial units”, to the expected deaths in that population assuming that it had the same age-specific death rates as a reference population. This indicator is commonly used to compare the mortality for different geographical areas to that of the reference population, and its main advantage is that SMRs do not depend on the age-specific rates of each spatial unit, which may be completely unreliable when working with small areas. This property has led to SMRs being, by far, the most commonly used epidemiological indicator for small areas spatial studies. Nevertheless, SMRs do exhibit some particular problems. SMRs of different spatial units are not directly comparable, that is, even if the age-specific rates of two populations were all proportional, this would not mean that the same proportionality holds between their SMRs [4]. In particular, areas with identical age-specific rates do not necessarily yield the same SMRs. Therefore, strictly speaking, SMRs are only valid for comparing the units of study against the standard population of the study, but not against one another. Alternatively, a second age-standardization method for rates, the direct method, produces age-standardized rates (ASRs), a weighted average of the age-specific rates of each unit of study also using a reference population to determine that common set of weights. When mortality is the outcome of study, ASRs have a direct interpretation as the hypothetical death rate that would have occurred in the corresponding unit of study if its age composition were the same as that of the reference population [5]. Therefore, the standardization process of ASRs removes the effect that the age composition of the different population groups could have on raw rates. One advantage of direct over indirect standardization is that spatial units with proportional age-specific rates necessarily yield ASRs keeping that same proportion [4], and so they are also comparable to each other. Nevertheless, ASRs are not as popular as SMRs for small areas geographical analyses since ASRs are only reliable when trustworthy age-specific rates are available for each spatial unit. This usually requires having a considerable amount of observed health events, from now on “deaths”, for every combination of age group and spatial unit, which small areas do not usually have. As a consequence, when ASRs are pursued on small areas, statistical modeling becomes necessary. There have been various approaches to the smooth estimation of ratios, such as mortality rates. In this regard, Kafadar [6] proposed a linear smoother for geographically-defined data that are in form of ratios. Although it shows that the smoothed rate has less variability than the age-adjusted unsmoothed rate, it uses parametric weights and does not consider the dependence between age groups. Ezzati et al. [7] used information on the number of deaths and county-level sociodemographic characteristics to estimate mortality rates and probabilities of death, but merged the counties, pooled the data on deaths and population over 5 years and also did not consider the dependence between age groups. Many more epidemiologic measures show similar problems produced by the age composition of the units of study. For example, life expectancy at birth (LE) is an alternative mortality summary which shows similar problems. This indicator has the advantage of being very intuitive and easier to understand for the majority of people, therefore its use is increasingly demanded from health authorities. Its calculation also depends on the age distribution of mortality for each unit of study. Thus, although directly comparable between spatial units, it also poses the small area estimation problems arising from the age-specific risk estimates for each spatial unit [8]. Several Bayesian hierarchical approaches have been previously proposed to estimate smoothed LEs in small areas studies. Congdon [9] considers the smoothing of LEs but his model does not contain a specific age–space interaction term, so age–space risk variations are just a simple sum of separate age and space components. In other works [10–14] model area and age specific mortality rates in order to subsequently calculate LEs; nevertheless, their approaches do not include a fully structured (dependent) age–space interaction term jointly considering the variation of these two sources of dependence. These works do not use the dependence of contiguous age groups and locations for the age–space interaction term. As a consequence, these models fit the underlying overdispersion but could possibly fail to appropriately describe the space–age group structured interaction of the data. Consequently, the technical difficulties of structuring the dependence structure in age–space terms for LEs have not yet been properly addressed in the literature. In this paper, an autoregressive proposal, based on the original model developed by Botella-Rocamora et al. [15], is made in order to model the age-specific probabilities of death (PoDs) for each spatial unit, which assumes stronger dependence between closer age groups and spatial units. This age–space dependence structure enables information to be transferred between neighboring (consecutive) age groups and neighboring areas, at the same time, providing more reliable age-specific PoDs estimates. Moreover, the model proposed can be easily generalized to the analysis of age–space–time datasets producing spatio-temporal age-specific rates, for example, or include any additional factor, such as sex, which could be of interest to be studied. The smoothed age-specific PoDs from this model could be later used to derive enhanced indicators, such as smoothed ASRs or LEs, of high epidemiological value. This article is organized as follows. "Methods" section shows our modeling proposal for estimating age-specific PoDs in small areas geographical studies and its application to the construction of additional smoothed epidemiological indicators such as ASRs and LEs. The following section shows three case studies illustrating some real uses that smoothed age-specific PoDs could have, enhancing therefore traditional SMRs-based small areas studies. The first case study shows the application of the model introduced in "Methods" section to a real setting, a geographical study of lung cancer mortality in women. Its goal is to illustrate the convenience of considering age–space interaction in geographical studies. In the second case, the model is applied to study mortality from ischaemic heart disease for women in two cities at the census tract level. Smoothed ASRs are derived and compared for the census tracts of both cities, illustrating some advantages of this mortality indicator. In the last case study, the model is applied to estimate smoothed LEs, which is the most widely used synthetic indicator for characterizing overall mortality differences throughout areas. Finally, "Discussion" section summarizes the main contributions of our work. Methods Spatial smoothing of age-specific PoDs Martinez-Beneito [16] proposed a unifying framework for multivariate disease mapping, i.e., for jointly mapping several diseases taking into account the spatial correlation that they might show. A second reformulation of this proposal, the M-model, was proposed in Botella-Rocamora et al. [15] which enables the joint smoothing of tens of diseases in a single study. This approach, with important computational benefits, allows researchers to determine and take advantage of the correlations among diseases, revealing the existence of common risk factors or unknown associations among some diseases. In this model, dependence between spatial units is introduced in a structured way, while the relationship between diseases is introduced in an unstructured way, no prior dependence structure is assumed between diseases. However, this model allows considering also structured relationships between diseases, or geographical patterns in general, (such as an autoregressive structure) if this was considered convenient. Thus, the relationship between geographical patterns of different age groups should be possibly considered as a structured factor, reproducing its ordinal character. This modeling proposal could also be used in our case to model several risk patterns, one per age group, instead of several diseases, although assuming stronger correlations for contiguous age groups. This could be achieved, for example, by assuming spatially dependent first-order autoregressive dependency structures for age within each spatial unit. Our model proposes to model PoDs instead of age-specific rates. Age-specific rates are the relation between observed deaths and the person-years at risk available. Those person-years are not so easily available for small areas studies, which is why we propose to model PoDs instead of rates. Nevertheless, if mortality rates were really our goal, the person-years in the study could be reasonably approximated by the population at risk for each area multiplied by the years in the period of study, as we will see later. Additionally, LEs depend directly on PoDs instead of rates so the modeling of PoDs will be also suitable for deriving these indicators. The autoregressive adaptation of the M-model for PoDs was made as follows. We changed the typical Poisson likelihood made for most disease mapping studies by a binomial distribution since the population of each (small) spatial unit is further divided into several age groups, which makes the denominator of many age specific PoDs quite low. As a consequence, the binomial assumption seems substantially more appropriate in this context as, in contrast to the Poisson case, it implicitly assumes that the number of deaths for each spatial unit and age group may not be higher than the corresponding population. Thus, let Ysa, Nsa and Psa∗ be the number of observed deaths, population and (unknown) PoDs, respectively, for the s-th spatial unit and a-th age group, s=1,…,S; a=1,…,A. We use the asterisk in Psa∗ in order to distinguish the smoothed age-specific PoDs that we are modeling from the raw age-specific PoDs, Psa=Ysa/Nsa. We assume that Ysa follows the binomial distribution: Ysa∼Binomial(Nsa,Psa∗),s=1,…,S,a=1,…,A, and model the logit of the PoDs as: 1 logit(Psa∗)=μa+θsa where μa represents the intercept of the a-th age group and θsa models the PoDs variability between age groups and spatial units. Note that Psa∗ in this model, and in general the PoDs that we are referring to, are the conditional probabilities of death at the age interval a given that you have survived the previous age intervals. By modeling θsa we will induce dependence between age groups and spatial units at the same time. The ideas of the M-model will be used for this end. Specifically, Botella-Rocamora et al. [15] show that it is possible to induce both spatial and multivariate dependence on Θ=(θsa) as the product of matrices: 2 Θ=ΦM where Φ is an S×A matrix whose a-th column Φ·a follows a spatially correlated distribution. In the Φ matrix, the columns are independent of each other, and therefore they do not consider dependence between groups. The second term in this expression M, is an A×A matrix which induces in Θ the dependence between age groups that we pursue. The j-th column of Θ, say Θ·j, contains the spatially referenced logit probabilities for the j-th age group. One can, therefore, interpret that Eq. (2) defines the spatial patterns in Θ as a linear combination of underlying latent variables whose coefficients correspond to the j-th column of M. Thus, 3 Θ·j=Φ·1m1j+…+Φ·AmAj where mij is the (i, j)-th entry in M. Matrix M combines, through linear combinations, the spatial patterns of Φ yielding therefore dependent spatial patterns. In particular, we model the set of columns as independent Proper Conditional Auto-Regressive (PCAR) distributions, that is: Φ·a∼N(0,σ2(D-γaW)-1), where W and D are the adjacency and diagonal matrices typically summarizing the geographical structure in PCAR distributions (see [17] for example). According to Botella-Rocamora et al. [15], for the separable case with γ1=⋯=γA, the variance covariance matrix between age groups Σ is equal to MTM. Therefore, if we want the correlation between age groups to decrease as a function of their distance, M should be modelled in order to reproduce that effect. We will assume for simplicity a first-order autoregressive dependence structure between age groups. In that case, the covariance matrix between age groups would take the form: Σij=ρ|i-j|, where ρ is an autoregressive dependence parameter controlling the strength of the dependence between age groups. Note that the variance parameter of Σ has been removed from this latest expression since the overall variance of Θ=ΦM is already controlled by the term Φ by means of the variance parameter σ2 of its columns. An M matrix inducing the mentioned autoregressive dependence between age groups could be, for instance, the upper-triangular Cholesky matrix of Σ, which has this simple expression: M=1ρρ2⋯ρA-10(1-ρ2)1/2ρ(1-ρ2)1/2⋯ρA-2(1-ρ2)1/200(1-ρ2)1/2⋯ρA-3(1-ρ2)1/2⋮⋮⋮⋱⋮000⋯(1-ρ2)1/2. To complete the model specification, an improper flat uniform prior is proposed for μa and a uniform prior distribution on the interval [-1,1] for ρ. A vague uniform prior on [0, C] for a high enough (non-informative) value of C is also proposed for the standard deviation σ (see Gelman [18] or pages 164–171 of Martinez-Beneito and Botella-Rocamora [17] for more details). Finally, a uniform prior distribution on ]λ(1)-1,λ(I)-1[, where λ(1) and λ(I) are the lowest and highest eigenvalues of D-1/2WD-1/2 respectively, is proposed for γ1,…,γA, as suggested, for example, by Sun et al. [19] or Martinez-Beneito and Botella-Rocamora [17]. Note that the proposed model is equivalent to the spatio-temporal autoregressive model of Martinez-Beneito et al. [20]. This proposal would model the matrix Θ as a first order autoregressive process as: Θsa∼N(ρΘs(a-1),σΘ2). In this case, we would be using the autoregressive component for modeling dependence between age groups instead of temporal dependence. Nevertheless, we have posed our proposal as a M-model for several reasons. First, the M-model can be easily generalized to non-separable dependence structures; in fact, the model just proposed is already inseparable because of the different spatial correlation parameters (γ1,…,γA) considered for the different spatial terms in the model. In a similar manner, different variance parameters (σ12,…,σA2) could be considered for these patterns, yielding therefore more flexible, heteroscedastic, covariance structures [21]. Second, M-models are computationally convenient for modeling multivariate spatial patterns, which makes them an appropriate choice for fitting the data in regular Bayesian inference packages such as WinBUGS, OpenBUGS, Nimble… Finally, and in our opinion most importantly, M-models have been generalized to multidimensional models [22] where several factors, besides the spatial component, are considered. This would make it possible to consider additional spatial patterns that could be correlated with the age-specific PoDs that we are modeling in order to enhance the estimates in the model. For example, mortality data from other sexes or causes of death could be additionally considered in order to yield improved age-specific PoDs in our data set. Alternatively, multidimensional modeling could be used to disaggregate the data, for example by subperiods, race groups… while maintaining a reasonable risk estimate quality by considering dependence between these new groups. In summary, in our opinion, all this makes the M-models an appropiate and flexible proposal for modeling the age-specific PoDs, as proposed above. We should mention that some other previous proposals in the literature have dealt with the joint analysis of several dependence sources, such as age–space–time models. For example, Goicoa et al. [24] proposed an age–space–time separable model where age and time are modelled as first order random walks and spatial dependence is modelled by means of a Leroux et al. [23] dependence structure. Our proposal yields, in principle, more flexible dependence structures since different spatial correlations, for example, for each underlying spatial pattern are considered. Moreover, the autoregressive dependence structure for age groups allows the strength of this dependence source to be adapted. These modelling features could be implemented in INLA (the software package used in [24]), which in principle could take advantage of the sparse structure of the spatial and temporal precision matrices to speed up computations [25]. Nevertheless, this would not be so advisable, since those modelling features would increase the number of parameters to be integrated out by INLA, what would make its fit substantially slower. Similarly, Goicoa et al. [26] propose alternative age–space–time models, but using P-splines for modeling dependence for age and time. Once again, considering spatial random effects of different dependence parameters would seem problematic for this work as well. Thus, in this regard we find our proposal somewhat more suitable. Moreover, although some other Markov chain Monte Carlo (MCMC) based approaches have been proposed for the modeling of age–space–time interactions [27] we find our proposal particularly convenient as it does not require a specific coding of the MCMC algorithm and, in contrast, it can be fitted with regular Bayesian software packages. Epidemiological applications of the smoothed age-specific PoDs The proposal above produces age-specific PoDs, Psa∗, smoothed over space and age groups. It is important to note that any statistical indicator built on P∗=(Psa∗) will inherit the smooth character of P∗, thereby producing smooth, and in principle more reliable, estimates of other more complex elaborate epidemiological indicators. For example, ASRs are weighted averages of the age-specific death rates. The ASRs weights would be given by the age-composition of a reference population, so all the ASRs would represent the observed risks for a common (ideal) population; in this way the age composition of the units of study would in principle no longer be a confounding factor. Traditionally, the age-specific mortality rate for age group a and spatial unit s is given by: 4 Rsa=YsaPersonYearssa≈YsaT·Nsa=PsaT. where T is the number of years in the period of study. If the period of study was just one year then T=1 and then Rsa=Psa. Let poba be the population in age group a in the reference population and let the standard weights be given by: Wa=105·poba∑apoba. Then the ASR for the s-th spatial unit is usually defined as: 5 ASRs=∑aWa·Rsa. The 105 in Wa allows us to undertand the ASRs as the number of deaths that we would expect per 100,000 people from the reference population. Obviously, the unreliability of Psa when the spatial units are small will be transferred to the ASRs, which will make them unreliable indicators. However, we could alternatively define smoothed Age-Standardized Rates (sASRs) by simply replacing Rsa in the previous expression by Rsa∗=Psa∗/T, that is: 6 sASRsa=∑aWa·Rsa∗ which will inherit the smooth character of Psa∗ and will therefore solve the small area estimation problems that traditional ASRs show. Following this same idea, many more epidemiological indicators relying on age-specific PoDs or mortality rates could be smoothed. For example, LEs are usually calculated from life tables as a direct function of the age-specific PoDs, the width of each age interval and the fraction of the last age interval survived for those people dying at each age interval [28]. Additionally, particular care is required for the last age group since this is a right-opened interval, with no upper limit, so it deserves special attention. Nevertheless, if Dsa is the marginal probability of death (which is a combination of the conditional PoDs) for group age a and spatial unit s, and agea is the mean age of death for the people dying at age group a, LEs could be alternatively defined as the expected value of the observed ages of death in the following way: 7 LEs=∑a=1AageaDsa. In a similar manner to that made with life tables, agea could be calculated as the sum of the initial age of each interval and of the proportion of interval survived by those dying during that age interval. Usually, that proportion survived is assumed to be one half of each interval, except for the youngest age group where it is known that perinatal deaths make this proportion shorter [29]. As a consequence, agea is usually taken as the central value for each age interval, which is equivalent to assuming the deaths to be uniformly distributed over those intervals. Regarding ageA, it is usual to assume a constant death rate, or equivalently an exponential distribution for the ages of death in the latest age group [30]. According to the available data, that death rate could be estimated as the age-specific death rate for this group. Therefore, according to the exponential assumption for the ages of death for this interval, the average number of years lived for any person reaching the oldest age group would be the inverse of its death rate, that is: T/PsA. Therefore, ageA would be the sum of this quantity and the initial age for this age group. A detailed description of the life expectancy calculation from a life table can be found, for example, in the Public Health England template [31] produced for this purpose based on the methodology described in [29]. Once the agea values have been calculated, Expression (7) can be formulated as a function of the conditional probabilities of death: 8 LEs=∑a=1A-1agea∏i100) used in many traditional disease mapping analyses since the mentioned sASR statistic takes into account the two extremes of the sASRs distribution instead of comparing the risk of each census tract with that of the reference population. Additionally, the sASRs calculated could be further compared to other additional ASRs corresponding to other areas of study. Thus, our sASR estimates could be compared for example to the ASR for the overall Valencian Region, Spain or even ASRs corresponding to other countries. In this manner, the use and scope of the generated results would be much wider than for SMRs. As an example, the bottom-right plot of Fig. 4 shows the distributions of the sASRs (their posterior means) for the Valencia and Alicante census tracts. As shown in that figure, most of the Valencia census tracts have sASR values smaller than any of the Alicante census tracts. As a matter of percentiles, the 95% percentile of the sASR distribution in Valencia (178.9 deaths per 105 women) is lower than the 5% percentile for Alicante (203.5). In other words, the census tracts with the highest IHD sASRs in Valencia have comparable mortality risks to those census tracts with the lowest mortality by IHD in Alicante. This should prompt the Valencian health authorities to take a close look at this issue, which would surely go unnoticed with a traditional sSMRs-based small area study. Case study 3: smoothed life expectancies Finally, in this third case study, age–space PoDs have been applied to the study of LEs, for each sex, once again for the whole Comunitat Valenciana (CV) for the period 2014–2017. In this case, mortality for all causes has been considered. Mortality data were provided by the mortality registry of the CV. Population data was obtained from the Spanish National Statistics Institute. The spatial aggregation considered was the municipality and 19 age groups: [0, 1), [1, 4], [5, 9), [10, 14),…, [80, 85) and [85,…) were considered for this analysis. Figure 5 shows the municipal LEs for men and women (separately) for the Valencian municipalities. The two maps on the left side of the figure show the raw LEs, directly calculated from the raw age-specific PoDs. The two maps on the right, show the smoothed counterparts of the maps on the left with the corresponding smoothed LEs. For comparability, cuts for these choropleth maps have been defined according to the smoothed LEs septiles for each sex. The legends of the raw LE plots, highlight the high variability of these indicators. Thus, raw LEs for men range from 48.7 to 104.0 and from 50.8 to 110.0 for women, which seems hard to believe. Moreover, 22 and 16 municipalities, respectively for men and women, showed infinite LEs since no deaths were observed in them during the period of study in the final age interval [8]. Smoothed LEs, however, show much more reasonable (moderate) variability than the raw LEs. This result clearly shows the need for smoothed age-specific PoDs for calculating LEs when dealing with small areas. Much clearer spatial patterns can be observed in the smoothed LE maps of Fig. 5, compared to their unsmoothed counterparts. Thus, for men, a clear cluster of low LE levels becomes evident in the eastern-central part of CV. This was not so evident according to the raw LEs. Moreover, some areas where the raw LEs showed high heterogeneity (see for example the north-western part of CV), where most of the municipalities are sparsely populated, now show a much more smoothed (spatial) behaviour which seems considerably more reasonable. For women, something similar happens. The smoothed LEs show a clearer spatial pattern and heterogeneity seems to have disappeared for the sparsely populated areas.Fig. 5 Raw and smoothed life expectancy at birth, for men and women, period 2014–2017 Additionally, the smoothed LEs in Fig. 5 allow interesting conclusions to be drawn, such as, for example, that geographical inequalities in men in terms of LEs (12.7 years between the most extreme municipalities) are substantially higher than those in women (5.5 years). According to the raw LEs these differences are reversed, 55.3 years for men and 59.2 years for women, although evidently these values are clearly untenable. The use of smooth LEs allows us to draw some meaningful conclusion in this regard. Finally, we would like to stress that small areas LE analyses provide a complementary view to that provided by sSMR based analyses. LEs take into account not just the mortality at each spatial unit but also the ages when those deaths take place, so LEs analyses complement the results that regular sSMR small areas analyses typically yield. Discussion Although spatio-temporal and multivariate disease mapping analyses on small areas have become quite popular, the same cannot be said for age–space small areas analyses. Nevertheless, as illustrated, the same methods that are used in both spatial and multivariate contexts could also be regularly used for age–space studies where age groups are disaggregated and separately studied, although considering dependence between them. In this work, the M-model is used in order to model the age-specific PoDs for each spatial unit, which assumes stronger dependence between closer age groups and spatial units. As shown, this age–space dependence structure enables information to be transferred between contiguous age groups and neighboring areas, at the same time, providing therefore more reliable (smoothed) age-specific PoDs estimates. Obviously, there are other alternative possibilities for the age-dependency structure used in the article. However, the one used seems appropriate given that it makes it possible to determine, within the same model, whether or not there is spatial and inter-age group dependency in the data, and even to quantify the strength of these sources of dependency. The small number of age groups in our case studies, and in most of the possible applications that we foresee, could make the discussion and selection of alternative dependency structures for age groups very difficult. There are multiple reasons to use the M-model for age-specific PoDs. The first of them is that it induces both spatial and multivariate dependence by simply multiplying two matrices. Secondly, the M-model can be easily generalized to non-separable dependence structures, resulting in more flexible, heteroscedastic, covariance structures [21]. In addition, these models are computationally convenient for modeling multivariate spatial patterns, which makes them an appropriate choice for fitting the data in regular Bayesian inference packages such as WinBUGS, OpenBUGS, Nimble… Moreover, M-models have already been implemented for the joint analysis of several correlated spatial patterns in INLA, either for a separable model for a large number of spatial patterns [41], or for an inseparable joint analysis of three spatial patterns [42]. Thus INLA should be also born in mind as a potential tool for this kind of analyses. Finally and the most important in our opinion, is that our proposal facilitates the future incorporation of other factors such as time or sex in a straightforward manner, following the proposal of multidimensional models by Martinez-Beneito et al. [22]. This allows us to consider additional spatial patterns that could be also correlated with the age-specific PoDs that we are modeling in order to enhance the estimates in the model current estimates. Beyond sASRs and sLEs calculations, the simple study of PoDs for visualizing the different age-specific spatial patterns is also interesting. Moreover, the study of those PoDs is particularly interesting when an inseparable model provides some fit improvement. For example, in another additional case study not presented here, we also applied our proposal for modeling all-cause mortality in the city of Valencia. In that case, applying the M-model at space–age group level, the smoothed age specific PoDs showed a pattern with stronger spatial dependence for the youngest age groups (higher values for the corresponding γa) than for the oldest age groups. The corresponding PoDs maps showed as if geographical inequalities existed for the youngest age groups, which pointed towards risk excesses for the deprived neighborhoods, but this effect diluted for older age groups since basically old people die all around while young people die, in general, just in these regions. Therefore, the modeling of non-separable age–space structures seems another additional application of our model, although we have not illustrated with such a detail for questions of space. In addition, any statistical indicator built on the smoothes age-specific PoDs inherit that smoothed character, yielding therefore more realiable estimates. Thus, in two of the three case studies we have two of the most relevant epidemiological indicators: the age standardized rate and the life expectancy. We have explored their use as an alternative to sSMR for small areas geographical analyses. On one hand, deriving an ASR version of potential use in small areas studies, with the direct comparation posibilities between areas that this brings, should be considered as an additional benefit of our proposal. On the other hand, smoothed life expectancies provides once again a version of this synthetic health indicator of potential use in small areas studies. Studies based on this indicator that take also into account the age of death could bring new possibilities to spatial epidemiological studies. Nevertheless, additional smoothed epidemiological indicators depending on PoDs may also be derived from the smoothed age–space PoDs, such as the number of potential years of life lost. Our study makes possible the study of those indicators also for small areas. Supplementary information Additional file 1. This document reproduces the analysis made for all three case studies of the article. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Jordi Perez-Panades, Paloma Botella-Rocamora and Miguel Angel Martinez-Beneito contributed equally to this article Supplementary information Supplementary information accompanies this paper at 10.1186/s12942-020-00251-z. Acknowledgements The authors acknowledge the support of the research Grant PI16/01004 (co-funded with FEDER grants) of Instituto de Salud Carlos III. Authors' contributions JP, PB and MAM contributed to the conception of the work, the modelization and the analysis. All authors discussed the results and contributed to the final manuscript. All authors read and approved the final manuscript. Competing interests The authors declare that they have no competing interests. ==== Refs References 1. Richardson S Thomson A Best N Elliot P Interpreting posterior relative risk estimates in disease-mapping studies Environ Health Perspect 2004 112 9 1016 1025 10.1289/ehp.6740 15198922 2. Ahmad OB Boschi-Pinto C Lopez AD Murray CJL Lozano R Inoue M Age standardisation of reports: a new who standard (technical report) 2001 Geneva World Health Organization 3. Neison FGP On a method recently proposed for conducting inquiries into the comparative sanatory condition of various districts J R Stat Soc Lond 1844 7 40 68 10.2307/2337745 4. Fleiss LBPMC Joseph L Statistical methods for rates and proportions 2003 3 New York Wiley 5. Curtin LR, Klein RJ. Direct standardization (age-adjusted death rates), vol. 6. US Department of Health and Human Services, Public Health Service, Centers for Disease Control and Prevention, National Center for Health Statistics. 1995. 6. Kafadar K Smoothing geographical data, particularly rates of disease Stat Med 1996 15 2539 2560 10.1002/(SICI)1097-0258(19961215)15:23<2539::AID-SIM379>3.0.CO;2-B 8961462 7. Ezzati M Friedman AB Kulkarni SC Murray CJL The reversal of fortunes: trends in county mortality and cross-county mortality disparities in the United States PLoS Med 2008 5 0557 0568 8. Eayres D Williams ES Evaluation of methodologies for small area life expectancy estimation J Epidemiol Community Health 2004 58 3 243 249 10.1136/jech.2003.009654 14966240 9. Congdon P A life table approach to small area health need profiling Stat Model 2002 2 63 88 10.1191/1471082x02st026oa 10. Jonker MF van Lenthe FJ Congdon PD Donkers B Burdorf A Mackenbach JP Comparison of Bayesian random-effects and traditional life expectancy estimations in small-area applications Am J Epidemiol 2012 176 10 929 937 10.1093/aje/kws152 23136165 11. Congdon P Estimating life expectancies for us small areas: a regression framework J Geogr Syst 2014 16 1 1 18 10.1007/s10109-013-0177-4 12. Congdon P Area variations in multiple morbidity using a life table methodology Health Serv Outcomes Res Methodol 2016 16 58 74 10.1007/s10742-015-0142-4 27257403 13. Dwyer-Lindgren L Bertozzi-Villa A Stubbs RW Morozoff C Mackenbach JP van Lenthe FJ Mokdad AH Murray CJL Inequalities in life expectancy among us counties, 1980 to 2014: temporal trends and key drivers JAMA Intern Med 2017 177 7 1003 1011 10.1001/jamainternmed.2017.0918 28492829 14. Alexander M Zagheni E Barbieri M A flexible Bayesian model for estimating subnational mortality Demography 2017 54 6 2025 2041 10.1007/s13524-017-0618-7 29019084 15. Botella-Rocamora P Martinez-Beneito MA Banerjee S A unifying modeling framework for highly multivariate disease mapping Stat Med 2015 34 9 1548 1559 10.1002/sim.6423 25645551 16. Martinez-Beneito MA A general modelling framework for multivariate disease mapping Biometrika 2013 100 3 539 553 10.1093/biomet/ast023 17. Martinez-Beneito MA Botella Rocamora P Disease mapping from foundations to multidimensional modeling 2019 Boca Raton CRC Press 18. Gelman A Prior distributions for variance parameters in hierarchical models Bayesian Anal 2006 1 515 534 10.1214/06-BA117A 19. Sun D Tsutakawa RK Speckman PL Posterior distribution of hierarchical models using CAR(1) distributions Biometrika 1999 86 2 341 350 10.1093/biomet/86.2.341 20. Martinez-Beneito MA López-Quílez A Botella-Rocamora P An autoregressive approach to spatio-temporal disease mapping Stat Med 2008 27 2874 2889 10.1002/sim.3103 17979141 21. Corpas-Burgos F Botella-Rocamora P Martinez-Beneito MA On the convenience of heteroscedasticity in highly multivariate disease mapping TEST 2019 28 1229 1250 10.1007/s11749-019-00628-8 22. Martinez-Beneito MA Botella-Rocamora P Banerjee S Towards a multidimensional approach to Bayesian disease mapping Bayesian Anal 2017 12 239 259 10.1214/16-BA995 29707101 23. Leroux BG Lei X Breslow N Halloran ME Berry D Estimation of disease rates in small areas: a new mixed model for spatial dependence Statistical models in epidemiology, the environment and clinical trials 1999 Berlin Springer 24. Goicoa T Ugarte MD Etxebarria J Militino AF age–space–time CAR models in Bayesian disease mapping Stat Med 2016 35 2391 2405 10.1002/sim.6873 26814019 25. Rue H Martino S Chopin N Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations J R Stat Soc B (Stat Methodol) 2009 71 2 319 392 10.1111/j.1467-9868.2008.00700.x 26. Goicoa T Adin A Etxeberria J Militino AF Ugarte MD Flexible Bayesian P-splines for smoothing age-specific spatio-temporal mortality patterns Stat Methods Med Res 2017 28 2 384 403 10.1177/0962280217726802 28847210 27. Quick H Waller LA Casper M Multivariate spatiotemporal modeling of age–space stroke mortality Ann Appl Stat 2017 11 4 2165 2177 10.1214/17-AOAS1068 28. Chiang CL The life table and its applications 1984 Malabar Krieger Publishing 29. Williams E, Dinsdale H, Eayres D, Tahzib F. Calculating life expectancy in small areas. Technical report, South East Public Health Observatory. 2005. 30. Silcocks PBS Jenner DA Reza R Life expectancy as a summary of mortality in a population: statistical considerations and suitability for use by health authorities J Epidemiol Community Health 2001 55 1 38 43 10.1136/jech.55.1.38 11112949 31. PHE: life expectancy template. https://fingertips.phe.org.uk/documents/PHELifeExpectancyCalculator.xlsm. 32. Lunn D Thomas A Best N Spiegelhalter D WinBUGS—a Bayesian modelling framework: concepts, structure, and extensibility Stat Comput 2000 10 325 337 10.1023/A:1008929526011 33. Vergara C, Martinez-Beneito MA. Pbugs. https://github.com/fisabio/pbugs. 34. Brooks SP Gelman A General methods for monitoring convergence of iterative simulations J Comput Graph Stat 1998 7 434 455 35. Carlin BP Gelman A Neal RM Markov chain Monte Carlo in practice: a roundtable discussion Am Stat 1998 52 2 93 100 36. Spiegelhalter DJ Best NG Carlin BP Van Der Linde A Bayesian measures of model complexity and fit (with discussion) J R Stat Soc Ser B (Stat Methodol) 2002 64 583 641 10.1111/1467-9868.00353 37. Duncan EW Mengersen KL Comparing Bayesian spatial models: goodness-of-smoothing criteria for assessing under- and over-smoothing PLoS ONE 2020 15 0233019 10.1371/journal.pone.0233019 38. Adín A Martinez-Beneito MA Botella-Rocamora P Goicoa T Ugarte MD Smoothing and high risk areas detection in space–time disease mapping: a comparison of P-splines, autoregressive and moving average models Stoch Environ Res Risk Assess 2017 31 403 415 10.1007/s00477-016-1269-8 39. Zurriaga O Vanaclocha H Martínez-Beneito MA Botella Rocamora P Spatio-temporal evolution of female lung cancer mortality in a region of Spain: is it worth taking migration into account? BMC Cancer 2008 8 35 1 10.1186/1471-2407-8-35 18173856 40. Besag J York J Mollié A Bayesian image restoration, with two applications in spatial statistics Ann Inst Stat Math 1991 43 1 21 10.1007/BF00116466 41. Vicente G Goicoa T Ugarte MD Bayesian inference in multivariate spatio-temporal areal models using inla: analysis of gender-based violence in small areas Stoch Environ Res Risk Assess 2020 34 1421 1440 10.1007/s00477-020-01808-x 42. Palmi-Perales F, Gomez-Rubio V, Martinez-Beneito MA. Bayesian multivariate spatial models for lattice data with INLA. 2020. arXiv:1909.10804.