In demography and epidemiology, the analysis of mortality data referenced over a Lexis structure allows to capture the simultaneous effects of age, period, and cohort (with their potential interactions); extract relevant mortality patterns; detect features of the population dynamics; and predict mortality trends by different ages and causes of death. The notion of Lexis surface was introduced by Arthur and Vaupel (Arthur and Vaupel, 1984) to generalize the concept of Lexis diagram: population data indexed over small time and age intervals can be interpreted as a surface with constant values in each time-age square. When the time and the age intervals are short, the discrete surface can be approximated by a continuous surface which may be interpreted as a density across the time and age domains (Arthur and Vaupel, 1984). One of the first applications can be found in Caselli et al. (1985), in which the authors introduce the contour map technique to represent surfaces of mortality rates.
At first, the approach was descriptive and oriented to develop statistical methods and graphical tools to display the data, useful in comparative studies. For instance in Caselli et al. (1987), the authors investigate differences and similarities between the mortality patterns observed in the surfaces of France and Italy, in the period 1900-1979. Since the publication of the seminal paper by Lee and Carter (1992), interest has been oriented towards the development of statistical models and computational methods to model mortality data in terms of surfaces. The main reason lies in the potential use of Lexis mortality data to forecast mortality dynamics and predict population scenarios by different ages. These aspects are particularly relevant for policy making in social security and public health, and also for risk evaluations by insurance companies. The Lee-Carter (LC) model (Lee and Carter, 1992)
represents the reference model in the literature and it is currently adopted as benchmark by several international institutions including the US Census Bureau and the United Nations. The LC model is a parametric model that can be used to produce mortality and life expectancy forecasts through a two steps procedure: the estimation of parameters by the ordinary least squares method with singular value decomposition and the prediction of mortality rates referring to age and time patterns by the autoregressive integrated moving average method.
Extensions of the LC approach have been proposed by several authors: Renshaw and Haberman modify the LC model to account for mortality reduction factors (Renshaw and Haberman, 2003b)2003a). In Brouhns et al. (2002) a Poisson log-bilinear regression model is proposed to improve the flexibility of the LC approach while Li and Lee (2005) use the LC method to forecast the mortality of a group of populations. Further extensions of the LC model can be found in Lee (2000), Lee and Miller (2001), and Li et al. (2009), while the use of the LC approach in a Bayesian framework is presented in Czado et al. (2005), Pedroza (2006), Katrien et al. (2015), and Wisniowski et al. (2015).
A Bayesian model, different from the LC approach, is presented in Dellaportas et al. (2001) to extend the parametric model introduced by Heligman and Pollard (1980) and account for incomplete life tables. With a different perspective Cairns et al. (2011) present a Bayesian approach to model the joint behaviour in mortality patterns of two dependent populations.
Another approach to the modelling and the estimation of mortality surfaces considers the use of nonparametric smoothing techniques. For instance, in Currie et al. (2004) the application of the P-spline method is presented, in Luoma et al. (2012) smoothing splines are considered in a Bayesian setting, while Li et al. (2016) implement a two dimensional kernel smoother to account for cohort effects. An application of a semiparametric approach with the use of the B-spline technique can be found in Barbi and Camarda (2011) in which the relative importance on the elderly mortality of the cohort effect versus the period effect is discussed.
The problem of the forward mortality surface, that is the prediction of the whole period-age structure of mortality, is addressed in Bauer et al. (2012) by the use of a finite dimensional Brownian motion while Gao and Shang (2017) present a functional data analysis to forecast multiple populations. Comparative studies of different approaches and models are reported in Haberman and Renshaw (2011), Giacometti et al. (2012), Danesi et al. (2015), and Novokreshchenova (2016).
Carstensen (2007) presents and discusses period-age-cohort models with disease data. The role of period-age, age-cohort and period-age-cohort models in epidemiology is reviewed in Clayton and Schifflers (1987a, b), particularly with respect to the analysis of cancer data.
A fundamental reference for the interpretation of mortality surfaces is the seminal paper by Vaupel et al. (1979) in which the notion of frailty is introduced to interpret the extra variation observed in mortality rates.
In the present work we propose a new spatial statistics approach to estimate the force of mortality across a Lexis structure in order to decompose the data into two interpretable components: a primary component accounting for the main smooth effect of period, age, and cohort; and a secondary component representing additional mortality in excess or in defect of that primary pattern. We present an extensive analysis using data from the Human Mortality Database about 37 countries. For several countries we observed a band of extra mortality in the secondary surface across the time domain, in the age interval between 60 and 90 years, with a slightly positive slope. The band is observed in the most populated countries and represents an overdispersion effect due to the presence of heterogeneity in the data.
The paper is organized as follows. Section 2 reviews the methodological background concerning the analysis of mortality data across a Lexis structure, introduces the basic idea of our approach, and discusses the similarities with other models. In Section 3 we present the hierarchical Bayesian model with mathematical details: likelihood, priors, and posterior. Furthermore, the Bayesian estimation of the two components and computational aspects of the Markov Chain Monte Carlo (MCMC) algorithm are discussed in detail. Section 4 describes the datasets we considered whilst Section 5 presents results and comments. In particular, estimates of the Bayesian mortality surface with the primary and secondary components are presented for the cases of Sweden and Italy, countries with different demographic transitions and mortality patterns. Finally, in Section 6 conclusions are drawn and future developments are briefly discussed.
2 A spatial approach for Lexis mortality data
A Lexis diagram is a two dimensional system of coordinates, for example representing calender time and age, fundamental to investigate demographic phenomena such as the mortality events occurring in populations of interest (Vandeschrick, 2001)
. A Lexis diagram for individuals in a population allows to plot the death and life of each individual. Each death event is represented by a point, with continuous time and age coordinates, whilst each individual of the population is represented by a segment with unit slope and connecting the moment of birth on the time coordinate with the point event. In the case of mortality events, the diagonal segments represent the individual life lines and the deaths of a population can be interpreted as the realization of a point process on the time-age coordinates system(Keiding, 1990).
The fundamental relevance of the Lexis diagram consists in the possibility to jointly represent three demographic dimensions: period, age, and cohort. In fact, across the diagonal directions with unit slope, it is possible to interpret the events as cohort related.
In this paper, we consider a new approach of the period-age-cohort structure and it is useful to interpret the diagram as a graph. We consider a finite time domain (e.g. a set of calendar years), a finite age domain (e.g. a set of age classes), and the two dimensional Lexis lattice , with knots. For each knot , we denote by the observed death count, by the exposed population, and by the force or intensity of mortality at the time and the age .
Under the assumption that each death count
is the realization of an independent Poisson random variable with respective mean, the maximum likelihood estimate of the mortality intensity can be obtained in each knot of the Lexis graph as the empirical mortality rate
When mortality is represented in log-scale on the Lexis graph, the empirical rates produce an empirical mortality surface (Arthur and Vaupel, 1984; Vaupel et al., 1997) that shows the main pattern of the force of mortality. On the other hand, if we are also interested to investigate the simultaneous effects of period, age, and cohort or the presence of latent structures due to heterogeneity factors, a statistical modelling approach is necessary (Carstensen, 2007).
In 1992, Lee and Carter introduced their seminal model (Lee and Carter, 1992). At each point of the Lexis graph, the mortality rate is modelled on the log-scale as the additive result of a baseline schedule across the age domain, , and a period-age interaction, defined by a multiplicative term, , that is
where is an independent Gaussian error.
The LC model captures the main pattern of the underlying force of mortality and the interaction period-age but it does not include a direct cohort component, accounting for the effect of potential selection factors (Renshaw and Haberman, 2006). Furthermore, some relevant issues arise on the identification of the parameters and the computation of their estimates so that constraints on the parameter estimates of and are necessary (Renshaw and Haberman, 2003b).
In the literature several extensions of the LC model have been considered, particularly with different representations of the interaction term and with different assumptions about the sampling error (Brouhns et al., 2002; Renshaw and Haberman, 2003b, 2006). Nevertheless, the common approach consists in modelling the components and in terms of factor effects or regression functions.
Our proposal originates from a different point of view. It is very reasonable to assume that the force of mortality at any point of the Lexis lattice is comparable and rather similar to the force of mortality of knots close in time and age, which we call here neighbouring knots. Therefore the mortality pattern should be rather smooth with respect to time, age, and cohort.
The smoothness of functions defined on a lattice system is often an important assumption in order to derive powerful and coherent statistical models. For instance, in the field of geographical epidemiology, the general risk of a certain disease is supposed to be smooth across the spatial domain of interest (Lawson, 2006); in image analysis, the intensity measured at some pixel of a satellite or medical image can be typically assumed similar to the measurements recorded in its nearby pixels (Li, 2009; Winkler, 2003). Spatial statistics is a collection of models and inferential methods to study functions which change smoothly in the space where they are defined (Banerjee et al., 2014; Cressie, 1991; Schabenberger and Gotway, 2005) and we follow this approach in order to model mortality data across the Lexis lattice. In fact, by analogy, a mortality surface can be considered as an image in which the domain of the pixels is defined across time and age instead of latitude and longitude.
In each knot we assume that the mortality intensity results from the effect of two components, and , acting at different regularity scales of the Lexis structure. The first component is a locally smooth term which takes into account for the main effect of time, age, and cohort. In addition to this smooth component, we account for additional mortality in excess or in defect. This second feature, that could be the result of heterogeneous latent factors acting on the population of interest, is modelled by a set of independent shocks .
Through the canonical link of the Poisson likelihood, we assume the following log-linear model
in which is a fixed offset as, for instance, the general baseline risk of the whole surface at the log-scale. We assume that the two components in Equation (1) are random effects: the first one, , is locally structured and similar to the effects in the neighbouring knots while the second one, , does not have any regular structure on the Lexis lattice. Notice that this second component is not a residual component but a term accounting for additional mortality with a non smooth pattern. We need to specify the definition of similarity and nighbourhood and details are given in Section 3.
The idea to represent the mortality intensity by two components is not new. It was introduced by Vaupel, Manton, and Stallard in their seminal paper (Vaupel et al., 1979) to account for the potential heterogeneity introduced by latent selection factors, the frailty by cohort. The authors considered the possibility that individuals, even though in the same cohort, can experience the event of interest (death) with different attitude or susceptibility: the individual frailty. Then, under the assumption that two mortality intensities at the individual level are proportional to each other by the ratio of the respective frailties, for each point , the authors derived the following relation
in which is the force of mortality for a standard individual with unit frailty while
is the expected frailty at the population level, obtained by averaging the individual frailty through a Gamma distribution, seeVaupel et al. (1979) for details. At the log-scale, Equation (2) is similar to Equation (1), but we make different assumptions on the two components which allow us to discover new features in the Lexis surface.
3 A hierarchical Bayesian model for the force of mortality
We adopt a Bayesian approach and formalize a hierarchical model for Equation (1) in which the randomness of the two components, the locally structured effect and the independent shock
, is modelled in terms of their prior probability distributions. These priors, and the additional hyperpriors, are combined with the Poisson likelihood through the Bayes theorem and in order to derive the joint posterior distribution, necessary to make inference on the quantities of interest. A similar approach was introduced in 1991 by Besag, York and Mollie(Besag et al., 1991) in the field of disease mapping and then widely applied, including ecological analysis (Clayton et al., 1993; Richardson and Best, 2003) and geographical demography (Divino et al., 2009).
Given the Lexis lattice , for each knot , we assume that the mortality count is the realization of a conditional independent Poisson random variable given the mean . Therefore, the full likelihood is given by
in which and denote the sets and respectively. Furthermore, we assume that the Poisson mean , at any point , transformed by the canonical logarithmic link, can be represented as follows
where represents the expected mortality count under the constant baseline intensity acting on the whole surface, that is . It means that the expected count , at any knot , is equal to the expected count under a constant intensity, , adjusted by the sum (at the exponential scale) of the two effects and , that is . Furthermore, from Equation (3) it follows
that models the force of mortality in terms of , , and . Notice that Equation (4) corresponds to Equation (1) with . The two terms and represent different features of the mortality intensity: is the locally structured component accounting for the smooth effect of time, age, and cohort whilst is the no structured component accounting for additional (but not residual) mortality in excess or in defect. Next, we define and precisely.
3.2 Priors and hyperpriors
Let us denote by and the sets of the components and respectively, that is and . We assume that the smooth term is apriori distributed as an intrinsic Gaussian Markov random field (Besag, 1974; Cressie, 1991; Li, 2009; Rue and Held, 2005) with pairwise interactions
in which the symbol denotes the symmetric relation of the first order Markov neighbourhood (Li, 2009) among the sites of , and is the positive precision parameter of the field , governing the strength of these pairwise interactions. We consider two knots and as neighbours if they are adjacent along the three direction of the Lexis graph: time (latitudinal), age (longitudinal), and cohort (diagonal); and along the diagonal with negative unit slope to account for further correlation between parallel cohorts. Therefore, the first order Markov neighbourhood system of a specific knot in red on the Lexis lattice in Figure 1 is represented by the eight adjacent knots, in yellow in Figure 1. In other words, in each point of the Lexis graph, given the neighbours, the component
represents a conditional autoregressive term with Gaussian distribution(Cressie, 1991).
< Figure 1 >
As is a set of independent variables, we assume that each is a priori distributed as a Gaussian random variable centered in zero and with positive precision parameter
, so that the joint distribution ofis given by
On top of these two priors we need to define the hyperpriors for the precision parameters. Following Mollie (1999), we assume that and are apriori independently distributed as Gamma variables, that is
are the respective hyperparameters. In this paper we consider vague hyperpriors(Congdon, 2014) with , , , and .
The offset can be estimated consistently by the maximum likelihood estimate
with the assumption that a constant mortality intensity is acting across the whole surface. It is an appropriately weighted average across time of the crude death rate (CDR), where the CDR in a certain year equals all deaths over all exposures that year. It can also be interpreted as the general effect in a log-linear analysis setting. In this paper we compute it as by Equation (5) and then fix it in the analysis. Notice that it could also be considered as a further parameter in the Bayesian modelling although it plays simply the role of scale quantity to stabilize the estimation and the computation of the model, and it does not have any effect on the decomposition of the force of mortality.
3.3 Posterior distribution
The Bayes theorem gives the joint posterior distribution
in which denotes the set of the sites adjacent to with respect to the first order Markov neighbourhood system in Figure 1. To obtain estimates for and , we have to integrate Equation (3.3), but this is not feasible and we need to approximate such integral by Monte Carlo.
3.4 MCMC computation
To produce samples from the marginal posteriors of , , , and we implemented an MCMC algorithm with Metropolis steps for the quantities and , and Gibbs Sampler steps for the parameters and (Liu, 2004; Robert and Casella, 2004).
In particular, for each and for each
, the Metropolis steps use Gaussian proposals centered on the previous value and with variance calibrated to obtain a global acceptance rate approximately betweenand (Roberts et al., 1997; Roberts and Rosenthal, 2001). The Metropolis acceptance probabilities are based on the following local potentials
respectively. Therefore, in the Metropolis step an old value is replaced by a new proposed value with respective probabilities
with obvious notation.
For and , we derive the conditional posteriors as Gamma distributions with conjugate parameters
respectively, and the Gibbs Sampler steps can be performed directly. Our MCMC algorithm can be summarized in the following scheme.
3.5 Bayesian estimates
The MCMC converges to the posterior in Equation (3.3) when the number of iterations increases to infinity. After the burn-in runs, the MCMC samples are used to compute point estimates for each quantity of interest by the ergodic marginal posterior means
where denotes the number of iterations considered after burn-in. Then, a point estimate of the mortality intensity at each knot of the Lexis graph is obtained by
with as in Equation (5).
The estimated smooth component , for varying , allows to interpret the regular mortality pattern that progresses in time, age, and cohort in a locally homogeneous way. Different populations, for example countries, will have different which can be compared to understand population specific smooth mortality structures. The shock component , on the other hand, collects additional but not residual effects of the mortality surface, which are not smooth. Therefore, allows to determine population specific features that can not be explained by structured smoothness.
4 The Human Mortality Database
We considered data at the time September 2016 from the open access Human Mortality Database (HMD, www.mortality.org), a project that since 2002 has involved the Department of Demography at the University of California, Berkeley (US) and the Max Planck Institute for Demographic Research in Rostock (Germany), and supported by other research institutes. When used for international comparisons, HMD data are generally assumed to be of superior quality compared to data from national statistical organizations. The reason is that raw HMD data are processed based on a common state-of-the-art procedure. The adjustments includes handling of death counts with missing ages and estimation of intercensal exposure times.
We considered the death count () and the population exposed to risk (), for female and male separately, and with time age intervals equal to 1-year 1-year, for the following 37 countries: Australia, Austria, Belarus, Bulgaria, Canada, Chile, Czech Republic, Denmark, Estonia, Finland, France, Germany, Greece, Hungary, Iceland, Ireland, Israel, Italy, Japan, Latvia, Lithuania, Luxembourg, Netherlands, New Zealand, Norway, Poland, Portugal, Russia, Slovakia, Slovenia, Spain, Sweden, Switzerland, Taiwan, United Kingdom, Ukraine, and United States. Belgium was omitted because of the presence of missing values in the data. The age domain of the data is the same for all the countries considered: 1-year classes from to , and a final cumulative class with individuals aged more than years. On the other hand, the time domain is not the same for each country: the Chilean data are the shortest series with the timespan of only years (1992-2005) while the Swedish dataset is the longest with years (1771-2014) of observations. The average timespan across all datasets is years. Table 1 reports a summary.
< Table 1 >
5 Results and Discussion
For each country, we estimated the two components and with the respective precisions and , and produced three estimated mortality surfaces: the surface of the marginal posterior means , the estimated primary smooth surface , and the estimated secondary surface representing additional mortality on top of . In addition, we considered also the surface of the empirical rates in log-scale as comparison term. In Figure 2 the case of Canada, male population, is presented as an illustration.
< Figure. 2 >
The Bayesian surface (top right in Figure 2) is the overall estimate based on the raw empirical rates (in log-scale) (top left in Figure 2). These two surfaces are rather similar in the case of Canada male and essentially in all our reconstructions, the purpose of this paper being the decomposition of into and , the smooth component and the additional free component. However, we see a difference between and in Figure 2, namely that in estimates for the most elderly population are given, while empirical estimates are missing. When information is poor, particularly at the elderly ages with small death counts and small populations at risk, the empirical mortality rates are affected by large variations or can not be even calculated at the log-scale, so the white areas in the top border of the surfaces. The presence in the Bayesian model of a spatial smooth component allows to estimate missing values as Bayesian averages of neighbouring rates.
The interesting contribution of our approach is the decomposition represented by (bottom left in Figure 2) and (bottom right in Figure 2). We see that is a smooth surface in the Canadian male case, it is almost identical to , as is in practice around zero everywhere. There is in this case a small exception, except for extra mortality at birth (first row in the bottom of ). Because we have a constant smoothing strength in (represented by a constant ), adding this extra mortality to the smooth component would apparently be problematic, as would be less smooth between the first and second row, compared to the rest of the surface. Hence extra mortality in .
5.1 Sweden and Italy
Our analysis of all 37 countries led to two different typical patterns concerning the secondary surfaces and well represented by the examples of Sweden and Italy. In Figure 3 and Figure 4, the surfaces of Sweden (first column) and Italy (second column), for female and male population respectively, are reported.
< Figure 3 >
< Figure 4 >
In the case of Sweden, the primary smooth surface is the only important component and represents in detail the whole shape of the mortality pattern, as it was for Canada in Figure 2, corresponding to estimated values being close to zero (in the case of the male population some weak effect can be noticed). On the other hand, in the case of Italy, for both genders, the secondary component presents strong and interesting patterns. First, at the birth, particularly at the beginning of the period and across the whole time domain, the surface shows relevant levels of additional mortality which are related to the infant mortality dynamics, and which appear to be stronger than what expected by an assumption of mortality smooth variation across the Lexis surface. Most importantly, we observe the presence of a large band across time, with a slightly positive slope and approximately in the age range 60-90. In this band the force of mortality accelerates with respect to the pattern explained by the primary component and displayed in the surface .
The visual analysis of mortality surfaces is also interesting to detect historical dynamics. For instance, in Figure 3 and Figure 4, the effects of the 1918 Spanish Influenza epidemic (both countries) and of two world wars (Italy only) are clearly visible. Elevated mortality represented by thin vertical lines in 1918 can be seen in all graphs. In addition there is a clear effect for Italian men of World War I (1914-1918) and World War II (1940-1945). For Italian women the effect is very small. Sweden was neutral under both wars, and hence no war effects appear.
< Table 2 >
5.2 The mortality profiles
The Bayesian procedure allows also to decompose estimated mortality profiles for specific cohorts, for specific ages, or for specific time points. In Figure 6, Figure 7, and Figure 8 we show, as examples, the mortality profiles of the cohort 1902, for the infant mortality during the period 1872-2012, and by age in the year 2012 respectively; for Sweden and Italy, female and male population.
< Figure 6 >
< Figure 7 >
< Figure 8 >
The significant presence of the extra mortality band is seen in Figure 6 and Figure 8 in the secondary profiles of Italy, compared to the null effect present in the plots of Sweden. In Figure 7 we can observe the decreasing pattern of the infant mortality during the period 1872-2012, both in Italy and Sweden. It is interesting to notice the different patterns in the secondary profiles. In the case of Italy, female and male population, our decomposition allows to highlight a significant difference in the levels of mortality between the primary component and the secondary component: the infant mortality globally decreases (total and primary profiles) but since the 80s an increasing and then stable level of additional mortality with respect to the primary smooth profile is present in the Italian maps (secondary profiles). The extra mortality is in the range 0.5-1.5 in log-scale (corresponding to 1.6-4.5) compared to the baseline mortality rate in the magnitude of in log-scale ( and in mortality rate scale, for female and male respectively). This is an interesting observation that deserves further investigation.
5.3 Computational details
The convergence of the MCMC chains was inspected by visual checks and by the Gelman-Rubin statistics (Gelman and Rubin, 1992). In Figure 9 we show the MCMC trajectories and histograms of the two precision parameters in the case of the Norwegian data (male population), across 1,000,000 iterations. The second chain converges to a posterior marginal which is flatter than the posterior marginal of the first chain, but convergence is attained quickly for both parameters. In our applications, we used 100,000 total iterations: 70,000 of burn-in with the subsequent 30,000 MCMC samples used to compute the marginal posterior mean for each quantity of interest.
< Figure 9 >
5.4 The precision ratio
The potential presence of extra mortality in the secondary surface can be detected also from the estimates of the two precision parameters and . These parameters summarize the mortality variation across the whole Lexis graph. In fact, in each point , the a priori conditional variance of , given potential values in its neighbours , is equal to (Mollie, 1999), with denoting the number of knots in . When the surface is not relevant, the estimated values show small variation with a high value of the parameter . On the other hand, when in the secondary surface there is a relevant presence of additional mortality, the variation among the values increases, reducing the level of precision. As in Mollie (1999), we introduce a measure to assess the respective roles played by each component and define the following precision ratio
that summarizes the information included in and . In fact, small values of reflect the presence of non smooth variation and a secondary surface with relevant patterns whilst large values of denote that the variation of the primary smooth component dominates, with the surface being not important.
In Table 2 the estimates of the precision parameters and the precision ratio for each country, female and male population, are reported. There are ten countries (France, Germany, Italy, Japan, Poland, Russia, Spain, United Kingdom, Ukraine, and United States) with particularly small values of the precision ratio. They correspond to the cases in which we detected the presence of extra mortality in the secondary surface.
The information represented by the precision ratio can be also displayed graphically as in Figure 5 in which we plot the precision ratio at the log-scale and by gender.
< Figure 5 >
Three well defined clusters can be noticed. The first cluster (highlighted by a green circle) includes the set of countries mentioned above which present extra mortality in the surface , in the second cluster (yellow circle) there are Canada and Netherlands with a moderate presence of additional mortality, while all the other countries for which the mortality pattern is well explained by the primary component belong to the last cluster (red circle). In general, from the bottom left corner of the diagram and along the diagonal, the countries are displayed by increasing levels of the precision ratio, corresponding to decreasing levels of extra mortality.
Interestingly, the green cluster includes the most populated countries in our datasets, Canada and Netherlands with a mid size of population belong to the yellow cluster, while the countries in the red cluster have all small populations, see for example United Nations, Department of Economic and Social Affair - Population Division (2015). Therefore, it seems that there is a potential positive association between the presence of extra mortality in the secondary surface and population size of the respective country.
5.5 The overdispersion band
As we reported above, the countries in the green cluster in Figure 5 (France, Germany, Italy, Japan, Poland, Russia, Spain, United Kingdom, Ukraine, and United States) present in the secondary surface significant levels of extra mortality across the time domain, particularly at the birth and in a band with a slightly positive slope in the age interval 60-90.
In order to explain the relevance of these patterns, we need to further clarify the respective roles played by the two terms and in the estimation procedure. Both terms represent effects which are random in the Bayesian setting with appropriate priors: a smooth Gaussian Markov random field and a set of independent Gaussian variables respectively. Therefore, why would the estimation procedure need to assign additional mortality to some areas of the secondary surface? It is obvious that the use of the independent shocks , which by definition fit the data with a larger level of flexibility, would be necessary in addition to the smooth components only when the marginal variation of the data in those areas is larger than in the rest of the surface. In particular, when this variation is not compatible with the smoothing level of the rest of the surface. Hence, the mortality patterns potentially observable in the secondary surface represent patterns of overdispersion (Cox, 1983; Xekalaki, 2014), beyond what the Gaussian Markov random field can accommodate. Therefore, is not a residual surface but it accounts for additional mortality patterns with a level of variation larger than the level in .
Overdispersion, especially with count data, is a well known phenomenon (Breslow, 1984; Lindsey, 1995). It arises when the data exhibit a larger variation than expected under a reference model and may be generally determined by several factors: latent dependence among the elements of the population, contagion and change of the individual behaviours during the survey, clustering in the structure of the population and heterogeneity inside the population (Xekalaki, 2014). In the particular case of spatial count data, the overdispersion may be determined by two relevant causes as pointed out in Clayton et al. (1993): the clustering and the heterogeneity acting on the response counts across the spatial domain of interest, that in our case is the Lexis graph. Therefore, following Clayton et al. (1993), we can interpret the model in Equation (4) as a log-linear model with constant average and with two random terms, and , accounting for extra variation due to the clustering and extra variation due to the heterogeneity respectively. Hence, while the primary surface shows the smooth pattern in terms of overdispersion effects by clustering, the patterns potentially present in the secondary surface represent overdispersion effects due to the heterogeneity in the counts and display specific features which are locally non smooth as the infant mortality dynamics or the excess of mortality during the period of the Spanish Influenza epidemic for instance.
In this setting the overdispersion band in the age interval 60-90 can be interpreted as a heterogeneity effect due to large counts of deaths. In fact, we notice that the bands occur in large countries in ages 60-80 to begin with, and in ages 70-90 in recent years. These are the age intervals in which many individuals die (see also the plots in Figure 6 and Figure 8), with higher marginal variation and because of many different causes of death.
Furthermore, the slope agrees with an increasing mean age at death (i.e. observed deaths). Over the years, the bands show a tendency to become more concentrated around the mean age. This is in agreement with what we know about the age distribution of mortality in many countries: not only an increase in life expectancy/mean age at death, but also a compression of mortality around the mean.
The bands are not visible in countries with small populations because overdispersion is difficult to detect when the counts are small as pointed ou by Davison (2008): “large amounts of data will be needed to detect overdispersion when the counts are small”. Indeed, when data of small countries are aggregated with respect to a shared time domain, as we did with Nordic countries (Denmark, Finland, Iceland, Norway, and Sweden) across the period 1878-2013 for both female and male populations, then significant levels of extra mortality become clearly visible in the secondary surface, as displayed in Figure 10. Each country alone does not present a relevant secondary pattern but when the data are aggregated, the increasing in size makes the overdispersion band visible. Therefore, the bands do not appear in countries with small populations, because potential overdispersion is difficult to detect when the population exposed to risk is small.
< Figure 10 >
5.6 The Vaupel model
As mentioned in Section 2, the idea to decompose the force of mortality into two components was introduced by Vaupel et al. (1979). Indeed, when represented in log-scale, the model in Equation (2) can be rewritten as
which is very similar with Equation (4). In particular, is the force of mortality of a standard individual with unit frailty and, for varying and , represents the main standard pattern of mortality across the Lexis graph, as the primary surface previously introduced. Therefore, in each point of the Lexis graph the term and the quantity play similar roles as they account for the main trend over time and age.
The second term in the Vaupel model represents the average frailty at the population level for the surviving individuals in the cohort . Vaupel et al. (1979) consider that each individual in the cohort , at any point of the Lexis graph, has frailty which is a random quantity and Gamma distributed with parameters . These parameters are fixed to initial values at the time of birth and then updated by the recursive use of the cumulative hazard function, see Vaupel et al. (1979) for details (in the notation of their paper the parameters are considered depending only on the age as they are the same for every cohort and independent of the time point ). Under this assumption the term is derived as the expectation
where denotes the operator of expectation at the Lexis point with the respective Gamma distribution. As Vaupel et al. (1979)
assert, “the Gamma distribution was chosen because it is analytically tractable and readily computable”, but its shape is not very different from the Log-Normal distribution that was also considered by the authors as a potential choice to model the individual frailty. Therefore the termand the secondary component (at the exponential scale) originate from similar distributions, Gamma vs. Log-Normal, and play similar roles as they both account for additional mortality and “acknowledging the known heterogeneity in populations” (Vaupel et al., 1979).
Thus, if we consider the following similitude
in terms of roles played in the respective frameworks, it is clear that our model is a simple extension of the Vaupel model in the direction of assuming that the term is a Markov random field smooth across the Lexis graph. Therefore, we also see that the secondary component can be interpreted as frailty at the population level (in log-scale).
The analysis of mortality data in terms of surface over a Lexis structure allows to capture the simultaneous effects of period, age, and cohort. Following an approach originally developed in spatial statistics, we proposed a model which allows to decompose the force of mortality into two components: a primary smooth component that gives the main features of age-specific mortality trends over time and a secondary component without any Lexis structure that represents additional mortality patterns not captured by the main trend.
The model was estimated in a Bayesian framework. Therefore the two components were specified in terms of prior distributions: a set of conditional autoregressive terms jointly distributed as an intrinsic Gaussian Markov random field and a set of independent Gaussian variates respectively. The smoothing of the primary component and the flexibility of the secondary component were controlled by the respective precisions and , included into the model as hyperparameters with vague Gamma priors. In order to compute Bayesian estimates we used MCMC computation with Metropolis and Gibbs sampler steps.
It is important to emphasize that a Bayesian approach can be sensitive to the choice of priors and hyperpriors. In our model, a relevant question concerns the distribution of the precision parameters and . Indeed, they modulate the trade-off between smoothing and no smoothing across the whole Lexis surface and, although they are sensitive to the choice of the respective hyperparameters, the effect is not on the estimation of the force of mortality but mainly on the variance of this estimation (Penttinen et al., 2003). Thus, the mortality surfaces based on the Bayesian posterior estimates are at most weakly affected by prior choices on the hyperparameters.
We presented an extensive application to data from the Human Mortality Database, where 37 countries are considered. For each country, three mortality surfaces were computed: the Bayesian surface, estimating the whole force of mortality; the primary surface, smooth across period, age, and cohort; and the secondary surface of potential additional mortality. We interpreted the potential extra mortality in the secondary surface as an overdispersion effect. A significant amount of extra mortality in the secondary component is necessary only when the smooth component alone is not able to account for the variation of the data.
Particularly interesting are two relevant patterns over the time domain: an excess of mortality at the birth related to the infant mortality dynamics and a band of overdispersion in the age interval 60-90 with a slightly positive slope in time in agreement with an increasing mean age at death and a compression of mortality around the mean. The patterns were detectable only in countries with large population size. In countries with small population size those patterns are not absent, they are simply hidden and the lack of detection is due to the small amount of statistical information in terms of small counts.
There can be several reasons behind the excess mortality represented by the bands: of course one possible hypothesis is that better health systems delay the death of a part of the elderly population, which therefore has to die “in excess” in the later years. The slight positive trend of the band might indicate that these excess deaths are further delayed. We notice that these effects are in addition to the general smooth trends in population aging, which are present in the primary smooth component.
The issue is also related to the similitude between our Bayesian decomposition and the model proposed by Vaupel et al. (1979), similitude that allows to interpret the secondary surface in terms of frailty patterns at the population level. These relevant questions deserve further investigations.
- Arthur and Vaupel (1984) Arthur, W. B. and Vaupel, J. (1984). Some general relationships in population dynamics. Population Index, 50(2), 214–228.
- Banerjee et al. (2014) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2014). Hierarchical Modeling and Analysis for Spatial Data. Chapman & Hall, London.
- Barbi and Camarda (2011) Barbi, E. and Camarda, C. G. (2011). Period and cohort effect on elderly mortality: a new relational model for smoothing mortality surfaces. Statistica, 71(1), 51–69.
- Bauer et al. (2012) Bauer, D., Benth, F. E., and Kiesel, R. (2012). Modeling the forward surface of mortality. SIAM Journal on Financial Mathematics, 3(1), 639–666.
- Besag (1974) Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society. Series B, 36(2), 192–236.
- Besag et al. (1991) Besag, J., York, J., and Mollie, A. (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics, 43(1), 1–21.
- Breslow (1984) Breslow, N. E. (1984). Extra-Poisson variation in log-linear models. Journal of the Royal Statistical Society. Series C (Applied Statistics), 33(1), 38–44.
- Brouhns et al. (2002) Brouhns, N., Denuit, M., and Vermunt, J. K. (2002). A Poisson log-bilinear regression approach to the construction of projected lifetables. Insurance: Mathematics and Economics, 31, 373–393.
- Cairns et al. (2011) Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., and Khalaf-Allah, M. (2011). Bayesian stochastic mortality modelling for two populations. ASTIN Bulletin, 41, 29–59.
- Carstensen (2007) Carstensen, B. (2007). Age-period-cohort models for the Lexis diagram. Statistics in Medicine, 26, 3018–3045.
- Caselli et al. (1985) Caselli, G., Vaupel, J. W., and Yashin, A. (1985). Mortality in Italy: contours of a century of evolution. Genus, XLI(1-2), 39–55.
- Caselli et al. (1987) Caselli, G., Vallin, J., Vaupel, J. W., and Yashin, A. (1987). Age-specific mortality trends in France and Italy since 1900: period and cohort effects. European Journal of Population, 3, 33–60.
- Clayton and Schifflers (1987a) Clayton, D. and Schifflers, E. (1987a). Models for temporal variation in cancer rates. I: age-period and age-cohort models. Statistics in Medicine, 6, 449–467.
- Clayton and Schifflers (1987b) Clayton, D. and Schifflers, E. (1987b). Models for temporal variation in cancer rates. II: age-period-cohort models. Statistics in Medicine, 6, 469–481.
- Clayton et al. (1993) Clayton, D., Bernardinelli, L., and Montomoli, C. (1993). Spatial correlation in ecological analysis. International Journal of Epidemiology, 22(6), 1193–1202.
- Congdon (2014) Congdon, P. (2014). Applied Bayesian Modelling. John Wiley & Sons, New York.
- Cox (1983) Cox, D. R. (1983). Some remarks on overdispersion. Biometrika, 1(7), 269–274.
- Cressie (1991) Cressie, N. (1991). Statistics for Spatial Data. John Wiley & Sons, New York.
- Currie et al. (2004) Currie, I. D., Durban, M., and Eilers, P. H. C. (2004). Smoothing and forecasting mortality rates. Statistical Modelling, 4, 279–298.
- Czado et al. (2005) Czado, C., Delwarde, A., and Denuit, M. (2005). Bayesian Poisson log-bilinear mortality projections. Insurance: Mathematics and Economics, 36, 260–284.
- Danesi et al. (2015) Danesi, I. L., Haberman, S., and Millossovich, P. (2015). Forecasting mortality in subpopulations using Lee-Carter type models: a comparison. Insurance: Mathematics and Economics, 62, 151–161.
- Davison (2008) Davison, A. C. (2008). Statistical Models. Cambridge University Press, London.
- Dellaportas et al. (2001) Dellaportas, P., Smith, A. F. M., and Stavropoulos, P. (2001). Bayesian analysis of mortality data. Journal of the Royal Statistical Society. Series A, 164, 275–291.
- Divino et al. (2009) Divino, F., Egidi, V., and Salvatore, A. S. (2009). Geographical mortality patterns in Italy: a Bayesian analysis. Demographic Research, 20, 435–466.
- Gao and Shang (2017) Gao, Y. and Shang, H. L. (2017). Multivariate functional time series forecasting: application to age-specific mortality rates. Risk, 5(2), 21.
- Gelman and Rubin (1992) Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7, 457–511.
- Giacometti et al. (2012) Giacometti, M., Bertocchi, M., Rachev, S. T., and Fabozzi, F. J. (2012). A comparison of the Lee–Carter model and AR–ARCH model for forecasting mortality rates. Insurance: Mathematics and Economics, 50(1), 85–93.
- Haberman and Renshaw (2011) Haberman, S. and Renshaw, A. E. (2011). A comparative study of parametric mortality projection models. Insurance: Mathematics and Economics, 48, 35–55.
- Heligman and Pollard (1980) Heligman, L. and Pollard, J. (1980). The age pattern of mortality. Journal of the Institute of Actuaries, 107(1), 49–80.
- Katrien et al. (2015) Katrien, A., Bardoutsos, A., and Ouburg, W. (2015). Bayesian Poisson log-bilinear models for mortality projections with multiple populations. European Actuarial Journal, 5(2), 245–281.
- Keiding (1990) Keiding, N. (1990). Statisitcal inference in the Lexis diagram. Philosophical Transactions of the Royal Society of London, 332(1627), 487–509.
- Lawson (2006) Lawson, A. B. (2006). Statistical Methods in Spatial Epidemiology. John Wiley & Sons, New York.
- Lee (2000) Lee, R. D. (2000). The Lee-Carter method for forecasting mortality, with various extensions and applications. North American Actuarial Journal, 4(1), 80–91.
- Lee and Carter (1992) Lee, R. D. and Carter, L. R. (1992). Modeling and forecasting U.S. mortality. Journal of the American Statistical Association, 87(419), 659–671.
- Lee and Miller (2001) Lee, R. D. and Miller, T. (2001). Evaluating the performance of the Lee-Carter method for forecasting mortality. Demography, 38(4), 537–549.
- Li et al. (2016) Li, H., O’Hare, C., and Vahid, F. (2016). Two-dimensional kernel smoothing of mortality surface: an evaluation of cohort strength. Journal of Forecasting, 35(6), 553–563.
- Li et al. (2009) Li, J.-H., Hardy, M. R., and Tan, K. S. (2009). Uncertainty in mortality forecasting an extension to the classical Lee Carter approach. Astin Bulletin, 39(1), 137–164.
- Li and Lee (2005) Li, N. and Lee, R. D. (2005). Coherent mortality forecasts for a group of populations: an extension of the Lee–Carter method. Demography, 43(3), 575–594.
- Li (2009) Li, S. Z. (2009). Markov Random Field Modeling in Image Analysis. Springer, New York.
- Lindsey (1995) Lindsey, J. K. (1995). Modelling Frequency and Count Data. Oxford University Press, Oxford.
- Liu (2004) Liu, J. S. (2004). Monte Carlo Strategies in Scientific Computing. Springer, New York.
- Luoma et al. (2012) Luoma, A., Puustelli, A., and Koskinen, L. (2012). A Bayesian smoothing spline method for mortality modelling. Annals of Actuarial Science, 6, 284–306.
- Mollie (1999) Mollie, A. (1999). Bayesian and empirical Bayes approaches to disease mapping. In A. Lawson, A. Biggeri, D. Böhning, E. Lesaffre, J.-F. Viel, and R. Bertollini, editors, Disease Mapping and Risk Assessment for Public Health, pages 14–29. John Wiley & Sons, New York.
- Novokreshchenova (2016) Novokreshchenova, A. (2016). Predicting human mortality: quantitative evaluation of four stochastic models. Risk, 4(4), 48.
- Pedroza (2006) Pedroza, C. (2006). A Bayesian forecasting model: predicting U.S. male mortality. Biostatistics, 7(4), 530–550.
- Penttinen et al. (2003) Penttinen, A., Divino, F., and Riiali, A. (2003). Spatial hierarchical Bayesian models in ecological applications. In P. J. Green, N. L. Hjort, and S. Richardson, editors, Highly Structured Stochastic Systems, pages 271–288. Oxford University Press, London.
- Renshaw and Haberman (2003a) Renshaw, A. E. and Haberman, S. (2003a). Lee–Carter mortality forecasting: a parallel generalized linear modelling approach for england and wales mortality projections. Appl. Stat., 52, 119–137.
- Renshaw and Haberman (2003b) Renshaw, A. E. and Haberman, S. (2003b). On the forecasting of mortality reduction factors. Insurance: Mathematics and Economics, 32, 379–401.
- Renshaw and Haberman (2006) Renshaw, A. E. and Haberman, S. (2006). A cohort-based extension to the Lee-Carter model for mortality reduction factors. Insurance: Mathematics and Economics, 38, 556–570.
- Richardson and Best (2003) Richardson, S. and Best, N. (2003). Bayesian hierarchical models in ecological studies of health–environment effects. Environmetrics, 14, 129–147.
- Robert and Casella (2004) Robert, C. P. and Casella, G. (2004). Monte Carlo Statistical Methods. Springer, New York.
- Roberts and Rosenthal (2001) Roberts, G. O. and Rosenthal, J. S. (2001). Optimal scaling for various Metropolis-Hastings algorithms. Statistical Sciences, 16(4), 351–367.
- Roberts et al. (1997) Roberts, G. O., Gelman, A., and Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. The Annals of Applied Probability, 7(1), 110–120.
- Rue and Held (2005) Rue, H. and Held, L. (2005). Gaussian Markov Random Fields. Chapman & Hall, London.
- Schabenberger and Gotway (2005) Schabenberger, O. and Gotway, C. A. (2005). Statistical Methods for Spatial Data Analysis. Chapman & Hall, London.
- United Nations, Department of Economic and Social Affair - Population Division (2015) United Nations, Department of Economic and Social Affair - Population Division (2015). World Population Prospects. Retrieved from https://esa.un.org/unpd/wpp/.
- Vandeschrick (2001) Vandeschrick, C. (2001). The Lexis diagram, a misnomer. Demographic Research, 4, 97–124.
- Vaupel et al. (1979) Vaupel, J. W., Manton, K. G., and Stallard, E. (1979). The impact of heterogeneity in individual frailty on the dynamics of mortality. Demography, 16(3), 439–454.
- Vaupel et al. (1997) Vaupel, J. W., Wang, Z., Andreev, K. F., and Yashin, A. I. (1997). Population Data at a Glance. Odense University Press, Odense.
- Winkler (2003) Winkler, G. (2003). Image Analysis, Random Fields and Markov Chain Monte Carlo Methods. Springer, New York.
- Wisniowski et al. (2015) Wisniowski, A., Smith, P. W. F., Bijak, J., Raymer, J., and Forster, J. J. (2015). Bayesian population forecasting: extending the Lee-Carter method. Demography, 52(3), 1035–1059.
- Xekalaki (2014) Xekalaki, E. (2014). On the distribution theory of over-dispersion. Journal of Statistical Distributions and Applications, 1(1), 1–19.