Multilevel survival analysis of health inequalities in life expectancy

The health status of individuals is determined by the interactive effects of multiple genetic, environmental, and socio-demographic factors operating at different levels, such as individuals' micro level socio-economic status and lifestyles, together with the macro level characteristics of the public health system and their residing geographical areas. As a result, health outcomes such as mortality rates and life expectancy are different in different population subgroups. Measuring health inequalities plays a fundamental role in the assessment of the performance or achievement of a health system, as part of the decision-making process that leads to amenable intervention. Although health inequalities have often been interpreted as a variation in health status measured by a specific health outcome across subgroups within a population [1], to reflect the sources of health determinants, a good health system should be measured against two objectives: the best attainable average level of health in the population it serves and the smallest feasible differences among individuals or groups within that population [2]. Most health inequality studies have focused on average differences across groups of people; addressing health inequalities caused by economic [3], social class [4, 5], demographic [6] and educational [7] diversity, but additional understanding about how variation in health is distributed at different levels may provide relevant information for public health intervention [8].

Given the multi-level influences on health, analyses at only the individual-level may ignore the potentially important effects of higher level factors. However, analyses of data aggregated to geographic levels may suffer from ecological bias [9] and has been criticised for not adequately addressing health inequality [1]. Analyses that can account for effects at different levels are key.

The World Health Organisation (WHO) proposed assessing the performance of national health systems using individual level health inequalities alongside the average level of health in a population [3, 10, 11]. In an analysis of individuals' health inequality based on mortality, Gakidou [12] proposed a measure of total health inequality derived from the beta-binomial regression model, which unified treatment of various measures including the Gini coefficient [13] and other estimates of inequalities. Their method has been used to estimate the risk of death in children under two and the between-individuals variance in the risk, based on a two-level data structure of children nested within families. The model provided absolute and relative inequality measures compared across 50 countries. Their model was adequate for mortality outcome based on either cross-sectional survey or cohort data with death event observed for all individuals in the cohort. For time to death or survival time as an outcome with censored time for some individuals based on cohort design, the beta-binomial model will not work. These researchers instead applied the frailty survival model to estimate the risk of death distribution in adults [14] where 'frailty' refers to excess systematic variability or random effects due to unmeasured factors at the community or geographic level which were not explained by individual-level covariates in the model. Having controlled for community-level effects, some inequality indices were also derived and compared between gender and age groups. These marginal models are flexible enough to model distribution of health outcomes with adjustment of individual level confounding, and can provide efficient estimates of health outcomes by accounting for overall variation at higher levels. However, these models do not separate the total variance explicitly by the source of variation, and cannot fit data with more than two levels of hierarchy [15, 16].

In practice, a model with a hierarchical structure of two or more levels is often plausible such as having individuals nested within households within residential areas. Inequality estimates which do not take account of this hierarchical structure may not identify variation at relevant levels and thus fail to pinpoint potential targets for interventions [8].

In recent years, multilevel models (MLM) have been available to further improve analysis of health inequalities over the marginal models [15, 17]. By explicitly modelling separate random effects at each level of the data, the MLMs can provide estimates of the outcome distribution at each level, quantify and distinguish contextual effects (at higher levels) from individual effects and separate variances of social groups from the total variation within a population [15]. The importance of MLM analysis is well recognized in social sciences, public health, health care [18] and epidemiological research [19]. The use of MLM in health inequalities research is increasing [8, 18, 20], but most previous reports are on multilevel logistic models for mortality outcome based on cross-sectional data. There is a lack of reports on multilevel models directly for modelling life expectancy (LE) based on cohort data with censored information on survival time, although multilevel survival models [15, 21] built in the software MLwiN [17] are already available.

In this study we apply multilevel survival models to investigate distributions in life expectancy (from 1970 onwards) among individuals from different socio-economic sub groups and multiple sources of variation of life expectancy in a cohort of elderly Swedish individuals followed up for 30 years. Like any multilevel models, multilevel survival time models consist of two parts, fixed effects and random effects. The model estimates of the fixed effects measure the mean differences between the socio-economic sub groups in life expectancy (so called socio-economic inequality), and estimates of the random effects quantify the variance or stochastic elements in life expectancy that cannot be accounted for by the mean differences. The total variance in life expectancy is disentangled into components for variation between municipalities, between parishes, between households and between individuals. The unexplained heterogeneity in the term of model residuals then is used to calculate some established and meaningful inequality indices for the health outcome of LE at different levels in the data hierarchy.

Our study is mainly aimed at introducing the new method for health inequality research. We illustrate the methods using a real example analysis that explores socio-economic and geographical inequalities in life expectancy in the elderly Swedish cohort. The analysis has four specific aims: (i) to identify components of the total variation in individual life expectancies by the geographic level of their residence, namely, individual, household, parish, and municipality; (ii) to explore the extent to which socio-economic and demographic factors defined at different levels can account for the variation in life expectancy; (iii) to compare the variation in life expectancy within socio-economic and demographic groups within the population; (iv) to explore absolute and relative indices of health inequalities in life expectancy in different population subgroups.

To obtain model based age and gender specific LE for model validation, Model C (Additional file 1, Table S1) which includes age and gender was considered. Adding the three random part estimates to the fixed part estimates in Model C, we obtained the model based median life expectancy of males for the age cohorts 65, 66, 67, 68 in 1970 as 14.4, 13.8, 13.2, 12.2 and 11.8 respectively, and of females as 19.1, 18.2, 17.4, 16.3 and 15.5 accordingly. They were fairly close to the K-M estimators, suggesting acceptable fitting to the data.

Variances in life expectancy by gender and by household poverty

For variability of life remaining or survival time, the raw coefficient of variance (CV) demonstrates greater variation of survival time in men than in women (Table 2), and a graded pattern of reduced variation of the LE followed by the decreasing of household poverty level of individuals. The adjusted variance based on the model E in Additional file 1, Table S1 has presented comparable estimates of variation for men and women (0.569 and 0.446 respectively), and similarly graded changing pattern for subgroups by the household socio-economic level.

	Measures Based on Raw Data			Measures Based on Model E
	K-M LE (SE)	Gini Coef.	CVa	Variance (SE)	II[1,1]	√IM[2,2]
Household SES
Poorest	11.0 (0.50)	0.419	0.655	0.856 (.062)	0.168	0.394
Poor	14.0 (0.27)	0.362	0.562	0.620 (.033)	0.141	0.296
Middle	15.0 (0.07)	0.320	0.511	0.502 (.031)	0.125	0.252
Top	16.0 (0.10)	0.304	0.490	0.479 (.031)	0.117	0.239
Gender
Men	13.0 (0.07)	0.318	0.552	0.569 (.004)	0.237	0.376
Women	18.0 (0.07)	0.317	0.454	0.446 (.004)	0.217	0.289

^a

Derived from the model estimates, Figure 1 presents density functions or the distribution of life remaining by gender, and Figure 2, by household socio-economic groups. The vertical lines in the figures represent the median years of life remaining for different categories of each variable.

This study attempts to improve on both the individual-level approach and the aggregated analysis to assessing health inequality by using multilevel survival models to measure health inequalities in the life expectancy of an elderly Swedish cohort after 30 years of follow-up. The analysis on the large Swedish cohort data has shown that a multilevel log-duration model can be an improvement on existing methodologies in the study of health inequalities of life expectancy among individuals, within areas or across areas in a population, in several aspects. First, total variance in life expectancy can be decomposed explicitly into more than two components to measure multiple macro and micro level stochastic elements in health inequalities. Secondly, the complete specification of the distribution of life expectancy by individual covariates, such as gender or social group, is provided in the fixed effects part of the model like any regression model, which determines means of gender or social groups. The variance components of random effects for each subgroup can also be estimated in the same model. The model may be of interest to both social scientists whose research is focused on the mean differences of social position on individuals and to those who are interested in the variation of life expectancy within subgroups of a population. Thirdly, area discrepancy in life expectancy can be studied based on random effects at the area level with statistical uncertainty around each area mean. Differences in the characteristics of residences across areas can be adjusted for in the same model in order to achieve a better comparison between areas, which is more efficient than the traditional method of calculating life expectancy for each area separately to study area differences. Finally, the effects of interactions between covariates at different levels, such as municipality income and gender in our case, can be studied directly in the same model without resulting in biased estimates of standard errors. Similarly, joint effects of individual and area income on life expectancy of individuals can be studied readily in the analysis.

The proposed method can also be applied to data with a shorter period of follow-up. The main difference in data between the longer and shorter time of follow-up is that for the latter design there would be more individuals who were still alive and be considered as censored data by the time for data analysis, and the mortality rate could be too low. An alternative approach is to fit a multilevel proportional hazard model, or multilevel Cox model to fit hazard rate based on death or live indicator at each time point [23].

For many cross-sectional studies, health inequalities in mortality measure would be the main interest of researchers. Given impacts of individual and area socioeconomic factors on mortalities, multilevel logistic models have been widely used to study health inequalities [19, 22, 24]. The same analytic strategies presented in the present study can be applied.

The calculation of summary measures of absolute and relative health inequalities across individuals and across areas, an important area of interest to health policy makers, can be studied based on multilevel models. The study showed that, for subgroups of individuals, the model-based change pattern in coefficient variance of the specific relative health inequalities was consistent with that of the raw CV but also sensitive to the variables controlled in the model, and that the magnitude of the measure depends on the estimates of the mean and variance for each subgroup. Simple comparison of such model-based coefficients between countries of different datasets with different covariates in the analysis for each country may not be sensible. This problem can be overcome, however, by pooling data from all countries together in one model for the analysis. In the absolute measure of inter-individual differences, our model based measure II[1,1] showed a similar pattern to the Gini coefficient (based on ecological analysis for subgroups of individuals of different household SES), a similar pattern between gender, and less inequalities by the model based estimates than a conventional simple approach as expected. The main reason for the difference could be that the Gini coefficient was calculated using the unadjusted survival probability of the Kaplan-Meier (KM) estimator within each subgroup, while the model based II[1,1] was calculated from the survival probability of the fitted Model E after adjusting for other variables, as well as taking into account variation between households and between parishes in the model. The model based approach enables us to understand what macro-level or micro-level factors might be associated with, or accountable for, how much variation in the health measure of interest (life expectancy or mortality) and estimate comparable inequality measures at both individual and area levels after controlling for known factors or confounding.

The model-based Gini coefficient and the coefficient of variance across parishes as summary measures of inequalities in the LE among parishes appeared comparable to results from an ecological analysis. The between-parish differences in the LE could be compared graphically and numerically based on random effects across parishes. For small samples at the area level, simulation based methods such as MCMC are available for estimating the distribution and confidence intervals of random effects [16, 18] for each parish.

Being focused on illustrating the methodology, this study used limited data which did not enable us to examine association between income at different levels or socio-economic position and inequalities in life expectancy in substantive depth. However, several important substantive findings have emerged in the study. Firstly, for this cohort most of the variation in remaining life expectancy after 1970 was among households and individuals. This suggests that in Sweden, where the socio-economic and health services tended to be distributed rather evenly among geographic areas, the health attributes and behaviours of household and individuals could be much more important in determining health outcome. However, this finding should not be generalised to other countries where wider regional social inequalities exist and cannot ignored. Multilevel modelling would be the adequate approach to analyse such data. Secondly, there were marked health inequalities among individual groups at different household socio-economic levels, especially within the poorest group, although household socio-economic differences explained only a small amount (1.8%) of variability in life expectancy at this level. Other factors such as ethnicity, lifestyle and morbidity rate of households could be further examined for their contribution to household health inequalities. Thirdly, comparable to other studies in the US and UK [15], males tended to have a larger coefficient of variation of life expectancy than females in this age group.

Multilevel survival models for longitudinal observation of survival time with censored data can handle complex data with geographic structure or many levels of hierarchy. They have the flexibility and efficiency of estimating mean differences between subgroups of individuals with different characteristics, adjusting for confounding, explicitly quantifying variation of the outcome at different levels, and measuring impacts of covariates on variability of the outcome at various levels simultaneously in one single model. The model can be extended to survival data collected in small areas by directly pooling data together in a 2-level model weighted by the variance of the outcome variable of the individual to produce variance estimates of small areas for the comparison of any summary measures of health inequalities between areas [15, 24]. The model illustrated can also be extended to an international model incorporating country as a separate level if many countries have small area health outcomes. If the outcome of interest is the risk of death or hazard rate instead of survival time or the LE, multilevel hazard models or multilevel discrete time models for death are readily applicable [25].

Much variability in life expectancy of individuals occurred at household and individual levels in this elderly Swedish cohort. Future research on determinants of health inequalities in the LE of the particular population should be focused on the household factors such as ethnicity, family life, contextual risks, health history of the family and individual factors such as genetics, socio-economic status, lifestyle, use of health services and health history etc. Further analysis on family and individual factors could reveal some important causes over the marked inequalities between individuals by their household poverty level reported in this study. Similar analysis on other age cohorts of the Swedish population could be interesting in establishing the magnitude of inequalities between age cohorts to identify possible intervention strategies to reduce such inequalities.