A nonparametric random coefficient approach for life expectancy growth using a hierarchical mixture likelihood model with application to regional data from North Rhine-Westphalia (Germany)

Background Life expectancy is of increasing prime interest for a variety of reasons. In many countries, life expectancy is growing linearly, without any indication of reaching a limit. The state of North Rhine–Westphalia (NRW) in Germany with its 54 districts is considered here where the above mentioned growth in life expectancy is occurring as well. However, there is also empirical evidence that life expectancy is not growing linearly at the same level for different regions. Methods To explore this situation further a likelihood-based cluster analysis is suggested and performed. The modelling uses a nonparametric mixture approach for the latent random effect. Maximum likelihood estimates are determined by means of the EM algorithm and the number of components in the mixture model are found on the basis of the Bayesian Information Criterion. Regions are classified into the mixture components (clusters) using the maximum posterior allocation rule. Results For the data analyzed here, 7 components are found with a spatial concentration of lower life expectancy levels in a centre of NRW, formerly an enormous conglomerate of heavy industry, still the most densely populated area with Gelsenkirchen having the lowest level of life expectancy growth for both genders. The paper offers some explanations for this fact including demographic and socio-economic sources. Conclusions This case study shows that life expectancy growth is widely linear, but it might occur on different levels.

http://www.biomedcentral.com/1471-2288/ 13/36 proposed here is focusing on temporal-spatial modelling of life expectancy with the aim of identifying spatial clusters in life expectancy growth, ultimately targeting on constructing a life expectancy growth map of NRW. The approach is less focussed on explaining differentials in life expectancy by other factors, say socio-economic factors (Gallo et al. [3]), in the sense of an ecological study analysis, although we will take up this string in the discussion.
In a nutshell, the approach is as follows. For each of the 54 regions a straight line model Y t = α + βt + is assumed for the life expectancy Y t at year t. Here, α and β are the intercept and slope of the line, respectively. As also Figure 1 suggests, we will model each of the 54 regions with identical slope but potentially different intercepts or levels of growth. The key point is that we will focus with our cluster-analysis on the variation of the different levels of growth and try to identify different components if present. The paper is organized as follows. In the next section, some background information of the life expectancy data and the region they stem from is given. In Section 'Methods' , the nonparametric mixture model used for the cluster analysis is presented in parallel with the associated likelihoods. It discusses the EM algorithm (Dempster,Laird and Rubin [4], McLachlan and Krishnan [5], McLachlan and Peel [6]) for finding maximum likelihood estimates as well as classifying the regions into the mixture components (clusters). Section 'Results' presents the results of the analysis including maps of the estimated cluster structure. The paper closes with a discussion which tries to put the findings into perspective.

Data
NRW is the most populous state of Germany, with four of the country's ten largest cities. Its capital is Düsseldorf. The state consists of five provinces (Regierungsbezirke), until 2010 divided into 31 rural districts (Kreise) and 23 urban districts (kreisfreie Städte), forming the above mentioned total of 54 districts which is the basis of our analysis. The underlying dataset 'LifeexpectancyNRW.xls' consists of two sheets, separately aggregated according to gender, each with N = 1134 observations. They include the following variables: -Region: the Municipality Code Number for each of the 54 districts in North Rhine-Westphalia, e.g. "1" for Düsseldorf, -LifeE: life expectancy in each region, -Year: calendar year from 1990 to 2010, recoded as 1 to 21 for this analysis, -UrbanRural: indicator, whether a region is rural (=0) or urban (=1).
Life expectancy is an important demographic indicator which is computed on basis of the life-table technique. A birth cohort is followed over time and, on the basis of the number of persons that died in every life year, mortality rates are determined which allow the computation of life expectancy. Life expectancy can be calculated conditional upon having reached any given age though it is typically considered from birth as done here. To provide timely life expectancy the current force of mortality (here for NRW) is applied to a hypothetical cohort and provides the data used in this study. Life expectancy computed in this way has to be interpreted as the expected life time for a newborn for the period in which the life table used was computed. For more details on life table techniques see Hinde [7].
The Table 1 shows the names of the regions in association with the numbering used in this analysis. Note that sometimes identical names occur such as Aachen (16) and Aachen (20). The explanation is that these correspond to different areas: the first refers to the city area whereas the second to the ensconcing rural vicinity.
Life expectancy is linearly growing in all regions in NRW as Figure 1 indicates. But growth occurs on different http://www.biomedcentral.com/1471-2288/13/36 levels where the level depends on the area under consideration (see again Figure 1). To explore these regional differences further, a likelihood-based cluster analysis is suggested in the following.

Model and associated likelihoods
We assume that the life expectancy Y it in region i and year t is available for i = 1, · · · , n and t = 1, · · · , T. Note that we use the index t starting from 1 instead of the Christian calendar. We further assume that a (latent) component structure is present within the population of regions which has not been observed directly and that within a component the life expectancy in time follows a simple straight model for j = 1, · · · , J and that within this component j and region i the data follow a multivariate normal distribution with diagonal covariance matrix and a common element describing this diagonal where y i = (y i1 , . . . , y iT ) T is the T-vector of observations of life expectancy in area i. Note that this model allows straight lines with potentially J different levels. Also note that σ 2 is the variance parameter of the meanzero normal random error it .
We should point out that we assume here that repeated observations of life expectancy are independent for the 21 observation years conditional upon component membership j. This assumption is crucial but not untypical for random effects modelling (see also McLachlan and Peel [6]). We mention that covariance structures could be modelled leading to a multivariate normal distribution for y i (and ultimately to mixtures of multivariate normals). However, we prefer to remain in the spirit of random effects modelling where we assume that covariance structures are coped with by the introduction of random effects.
Since we do not observe component membership we only take the marginal distribution as a nonparametric mixture where the p j represents the unknown weights of the components in the population.
Consequently, the observable mixture model likelihood is which needs to be maximized in α j , p j for j = 1, · · · , J and β and σ 2 . Note the special form of the likelihood in its hierarchical structure. Conditional upon component membership it assumes independence in time.
Note that this form of random effects modelling is not uncommon for this situation (Arminger et al. [8]; Goldstein [9]; Aitkin [10], Ng et al. [11], Ram and Grimm [12], Muthén and Asparouhov [13]; Rabe-Hesketh and Skrondal [14]). The central idea is that the random effect copes with the temporal and/or the spatial structure of the data.
Since the observed likelihood function is difficult to maximize in the parameters we consider the unobserved likelihood typical for mixture problems of this kind. Let z ij denote the unobserved indicator informing about component membership. In other word, z ij = 1 if the i-th http://www.biomedcentral.com/1471-2288/13/36 region belongs to component j and 0 otherwise. Then the unobserved complete likelihood is showing that the likelihood can now be separately maximized in α j , for j = 1, · · · , J, β, σ 2 on one hand, and p j for j = 1, · · · , J on the other hand. This is best done with the EM algorithm. The problem is well-posed in the sense that if J is fixed the likelihood is bounded and can be maximized. However, if the likelihood is also maximized with respect to J as suggested in the approach by Laird [15] and Lindsay [16,17], then the likelihood becomes unbounded as σ 2 approaches 0 (see also Wang [18]). To overcome this problem we follow Aitkin [10] and keep J fixed for determining the maximum likelihood estimates and then stepwise vary J. We then select the number of components on the basis of the Bayesian Information Criterion (BIC) to search for the best model.
We use the following definition for the BIC where r is the number of estimated parameters and n is the number of regions (here n = 54). Models are considered suitable with small BIC-value. Another criterion could be where T is the considered number of years (here T = 21). Note that N = nT so that two different penalty terms are possible, namely log(n) and log(nT), respectively. Given the choice of modelling which considers each area as clustered in time, we find (6) the more appropriate selection criterion which uses the number of areas n in the penalty term. This also seems to correspond to common practice (Muthén and Asparouhov [13]). For completeness, we shall compute both.

Expectation-maximization (EM) algorithm
To estimate the parameters by maximum likelihood we will use the EM algorithm (Dempster, Laird, and Rubin [4]; McLachlan and Krishnan [5]; McLachlan and Peel [6]). The EM algorithm consists of two steps: the E-Step and the M-Step. The algorithm cycles between these two steps back and forth.

E-step
In the E-step the unobserved indicator variables Z ij are replaced by their expected values conditional upon the current parameter estimates and the data y it These expected values can be easily computed using Bayes theorem as e ij = f ij p j J k=1 f ik p k and can be interpreted as the posterior probability that region i belongs to component j (note e ij ≥ 0 and j e ij = 1). Here

M-step
It is easy to see that the likelihood (5) is maximized for p j asp For the remainder we concentrate on Setting the partial derivative ∂ ∂α j log L = 0 leads tô Similarly, setting all other partial derivatives to 0 we achievê Here e ij , p j , α j , β and σ 2 refer to the values of these parameters in the previous cycle of the EM algorithm. The EM algorithm toggles between E-and M-step until convergence, say until is less than where is a small number such as 0.0001. refers to the absolute difference in each of the parameters between two consecutive cycles s + 1 and s, for example = |α s+1 j − α s j | if we consider the intercept.

Initial values
We need to compute initial values for the variables e ij , p j , α j , β and σ 2 . Only for this purpose we fit the following linear model: for each region i = 1, . . . , n, separately leading to n estimates of a i , b i and σ 2 i denoted asâ i ,b i andσ 2 i . Now we http://www.biomedcentral.com/1471-2288/13/36 use these estimates to get our starting values for the EM algorithm: Then we run the EM algorithm for various values of J starting with the homogeneity case J = 1 to get estimates of e ij , p j , α j , β and σ 2 . With these values we compute the (maximized) likelihood (4). Table 2 shows the results of the EM algorithm for men for J = 1, . . . , 10. As we can see the values for BIC decrease with growing index J until J = 7 just to increase again. The values for BIC 2 show the same behaviour but have the minimum at J = 6. We also run the EM algorithm for the female data for J = 1, . . . , 10. Again the BIC is decreasing but now we find the minimum at J = 8. The optimal J lies between 6 and 8 regardless of the different selection criteria. Overall, it seems appropriate to take J = 7 which splits the data into 7 different category groups. Table 2 provides also estimates for the slope β and the variance σ 2 . The choice of J = 7 appears also justified when we consider the value of σ 2 in dependence of J which becomes stable at J = 7. Note that the slope estimate is stable independent of the choice of J. Details of the full estimation results of the mixing distribution are provided in Table 3.

Maximum posteriori classification Men
Since each e ij describes the probability that region i belongs to component j, we can easily identify to which component each region belongs to according to the maximum posterior probability rule (MAP). The MAP classifies region i into component j where The classification tables are given in Table 4. This classification is the second column of Table 4 where one can also see the matrix e ij (rounded to 2 digits after the decimal point). Note that in all cases the classification is unique in the sense that there is a high classification probability for a particular component. Now we are able to construct a graph wherein the datapoints are coloured by the different components where they belong to (Figure 2).
In addition to the data points we have included in Figure 2 the regression lines for each component with the parameters from Table 2.

Women
In Table 5 we consider women. It again consists of the Region, the corresponding component for each region and the entire matrix e ij for J = 7. With this information we construct the plot of life expectancy for women whereas the data points are colored by

Explaining the cluster structure
The data contain also a variable characterizing each area as rural (= 0) or urban (= 1). The results can be summarized into the following cross-classified tables ( Table 6,  Table 7). A simple chi-square test investigates the relation between these two categorical variables: classification using the performed cluster analysis and the binary variable rural/urban. We find for men: χ 2 = 18.4645 by 6 DF and p-value = 0.00517, which is highly significant. For women we find: χ 2 = 15.3361 by 6 DF and p-value = 0.00178, clearly significant.
We conclude the results section with a final analysis as follows. We have done a separate cluster analysis for men and women. For men, a particular region will be classified into a component, but for women this region might be classified into a different component. To provide a consistent analysis both classifications should be correlated. This is what the last graphic is about. Figure 6 shows the connection between the components for each region for men and women. There we can identify for which regions the life expectancy for both men and women is high or low. For example, region 30 (Münster) and 17 (Bonn) are for men and women in the highest level of life expectancy growth, whereas region 29 (Gelsenkirchen) is in lowest for both gender groups.

Discussion and conclusion
The normal density (2) is frequently used as mixture kernel and appropriate for our application. However, if necessary it allows easy extensions either in the mean structure or the variance-covariance structure. For one, one could allow component-specific variances leading to f (y i |α j , β, σ 2 j ). In addition, one could also think of giving up independence within area i leading to a multivariate normal distribution with either structured or completely unstructured variance-covariance matrix . Furthermore, one could think of modelling componentspecific variance-covariance matrices j . For two, instead of using a common slope model this could be generalized to component-specific slopes . The E-step of the EM algorithm has to be changed appropriately, in the case of a common variance parameter σ 2 leading tô T i e ij , as before, and However, for our data constellation the proposed model (2) is appropriate as also Figure 1 suggests. A more rigorous analysis for this assumption requires fitting the model with random intercepts only and also the model allowing the slopes to be random as well. This has be done using a mixed model approach with a normal random effects assumption. The BIC-values associated with the model fits support the random-intercept only assumption. Also, we have looked at the potential for curvature. This would correspond to an asymptotic change in life expectancy growth and relates to the discussion in Oeppen and Vaupel [1]. For males, the log-likelihood for the model with a quadratic term for year is −2 = 2919.0, whereas the log-likelihood for the model without the quadratic terms is −2 = 2920.2, leading to a likelihood ratio test statistic value of 1.2 with p-value 0.27, clearly not significant. For females, we have similar results.
A qualitatively different approach follows a conditional autoregressive model (CAR) which was originally suggested by Clayton and Kaldor [19] and more recently modified by Rasmussen [20]. In principle, the idea could be also utilized for spatial-temporal modelling as in this case and tries to implement a smoothing element by utilizing spatial information. The key element of CAR models is to model mean and covariance structure of the random effect (here the intercept in the temporal straight line model) by neighboring information. The ultimate goal is to reach a smooth map of the measure of interest (here level of life expectancy growth). This approach is quite meaningful, in particular, if the underlying process is likely to have a smooth characteristic. In our case, however, we were more interested in identifying a potential clustered structure in life expectancy growth for which we thought the likelihood based cluster approach is more appropriate.
Hence we have not followed up on CAR models in this case.
In North Rhine-Westphalia (NRW), there is an apparent continuous rise in life expectancy at birth in men and women within the last twenty years. However, this pattern needs to be contemplated differentially. Our analysis shows that in North Rhine-Westphalia, life expectancy is predominantly higher in rural than in urban districts and differs considerably by region. Within the observed period from 1990 to 2010, levels of growths of life expectancy ranged from 70.3 to 73.7 years in men and from 77.3 to 80.2 in women. Life expectancy in the 54 districts was influenced by a latent categorical variable, which consists of seven categories or clusters. Each of the 54 districts is allocated into one of the seven clusters. This latent variable might be a surrogate variable for socio-economic factors. Life expectancy, as well as its counterpart mortality, strongly depends on factors like education, income, occupational status in addition to the factors sex and age. Most recent analyses of the European Prospective Investigation into Cancer and Nutrition (EPIC) showed that total mortality among men with highest education level is reduced by 43% compared to men with the lowest (hazard ratio (HR): 0.57, 95% confidence interval (CI) 0.52 -0.61). Among women, the reduction was 29% (HR 0.71, 95% CI 0.64 -0.78). In men, social inequalities were highly statistically significant for all causes of death examined. In women, the authors found a less strong, but statistically significant association with social inequalities for all causes of death except for cancer mortality and injuries (Gallo et al. [3]). For the region 29 (Gelsenkirchen), we found the lowest life expectancy for both, men and women. Socio-economic factors (see also Health reporting unit at the NRW Centre for Health [21]) support this finding and point to possible underlying

Mülheim
Oberhsn. receiving unemployment benefits (Bonn: 5,738.5/100,000; Münster: 5,090.3/100,000). Large universities are based in Bonn and Münster with thousands of students as city inhabitants. Therefore, the disposable income per inhabitant is above NRW average, but other regions show higher income rates. The proportion of smokers is relatively low in both cities (Bonn: 22.9%; Münster: 23.7%). In 2009, only in five rural districts the proportion of smokers was lower. Results for men and women differ slightly, as was reported for social inequalities in the EPIC cohort, too (Gallo et al. [3]). In men, besides Gelsenkirchen the cities Duisburg (2) and Oberhausen (7) are classified as regions with the lowest life expectancy of NRW. In women, it is only Gelsenkirchen. These findings support results of a socio-spatial cluster analysis conducted in 2007 by Strohmeier et al. [2] which was mentioned already in the introduction. Based on social indicators six clusters were established for NRW,  which classified the 54 districts into six types which were dubbed as follows: poverty pole, family zone, cities dominated by administrative and service units, rising regions / suburban counties, heterogeneous cities, heterogeneous rural districts. Like in our analysis, the poverty pole (representing areas which are in several ways socially deprived) included the cities Gelsenkirchen, Duisburg, and Oberhausen, but also the cities Dortmund and Herne which are all located in the Ruhr area.

Wuppertal
In relation to the NRW health indicators the authors found a significantly lower male and female average life expectancy in the poverty pole. In our analyses also more cities, especially of the Ruhr area, are categorized into the clusters of lower life expectancy. The Ruhr area is an urbanized, high density area comprising 11 cities and 4 counties with about 5 million inhabitants, formerly characterized by heavy industry and now undergoing a structural change towards e.g. information technology and health care industry. An additional underlying cause for lower life expectancy in this area might still be environmental. The Heinz Nixdorf RECALL study (Fuks et al. [22]), which included 4,291 participants from the Ruhr cities Bochum (43), Essen (3) and Mülheim a.d. Ruhr (6), recently confirmed that residential proximity to high road traffic (≤ 50m) and road traffic noise exposure (24h mean noise (Lden) > 65 dB) have a tendency toward higher blood pressure and an elevated prevalence of hypertension. Data from this study also showed that a reduction in distance between residence and major roads by half was associated with a 7.0% (95% CI 0.1 -14.4%) higher coronary artery calcification (CAC) (Hoffmann et al. [23]).
In a subgroup of the RECALL study population, participants residing in Essen (n=1,641) and Mülheim (n=1,742) for which digitized information on inner city roads was available, prevalence of coronary heart disease at high traffic exposure showed significantly elevated OR=1.85 (95% CI 1.21 -2.84, adjusted for cardiovascular risk factors and background air pollution) (Hoffmann et al. [24]). Further analysis showed a stronger effect for men (OR=2.33, 95% CI 1.44 -3.78), which might account for the difference among men and women in our cluster analysis. Another analysis of the RECALL data investigated if the association of road traffic exposure and subclinical cardiovascular disease might be modified by socioeconomic characteristics of individuals or neighborhoods. Participants with low socio-economic status (SES) and simultaneous exposure to high road traffic had highest levels of CAC (Dragano et al. [25]). The prevalence of high CAC was 23.9% in higher-educated men with low traffic exposure but 37.7% in lower-educated men with high road traffic exposure (women: 22.0% vs. 28.1%).
The cluster analysis of life expectancy once more stresses the differences between urban and rural regions in North Rhine-Westphalia. The latent component categorizing the 54 districts into seven categories can be interpreted as a surrogate comprising several underlying factors. The results point to districts where an accumulation of problems has negative impact on health. For males, only three cities are classified into the lowest cluster category, with 5.4% of the total NRW population living there. For women, only Gelsenkirchen is classified into this cluster. Given the emerging insight into possible underlying causes, chances for these cities to improve their outcome may come into closer reach.