2.1. Mortality data and relative risk

In this study, we model the relative mortality risk in an individual LSOA $i \in \{1, \ldots, N\}$ compared to the average mortality in England. To define our measure of relative risk, we first define a baseline death rate $m^b_{ta}$ for year t and age a for the whole of England in the usual way:

(2.1) \begin{equation} m^b_{ta} = \frac{\sum_{i=1}^N D_{ita}}{\sum_{i=1}^N E_{ita}}.\end{equation}

The model will be fitted using data from years $\mathcal{T}$ and age range $\mathcal{A}$ . Without any additional information, the expected^{Footnote 1} total number of deaths $\hat{D}^0_i$ across all ages $a \in \mathcal{A}$ and years $t \in \mathcal{T}$ in LSOA i is given by

\[\hat{D}^0_i = \sum_{t \in \mathcal{T}, \, a \in \mathcal{A}} m^b_{ta} E_{ita} \mbox{ for all }i=1,\ldots,N,\]

and we define the observed relative risk of death $R^0_i$ for an individual living in LSOA i as the ratio of the actual number of deaths to the expected number of deaths in that LSOA, that is,

(2.2) \begin{eqnarray}R^0_i &=& \frac{\sum_{t \in \mathcal{T}, \, a \in \mathcal{A}} D_{ita}}{\hat{D}^0_i} \mbox{ for all } i=1,\ldots, N .\end{eqnarray}

With our definition, the realised relative risk ${R^0}$ in any neighbourhood is a random variable since the realised number of deaths is random. In the following, we are interested in modelling the conditional expectation of ${R^0}$ given a vector of socio-economic characteristics.

Note that the relative risk $R^0_i$ is not age and year specific. However, as mentioned above, we will calculate and model the relative risk using mortality data for different age ranges $\mathcal{A}$ , see Section 5 for details.

2.2. Socio-economic characteristics

In Section 3, we will construct an index that explains differences in the mortality rates in different LSOAs based on differences in their socio-economic characteristics. In Wen (Reference Wen2022), a large universe of predictive variables for LSOA-specific mortality rates were considered. Based on findings there, we restrict our attention in this paper to 12 variables. They are listed in Table 1. Further details about those 12 variables, including data sources, can be found in the supplementary material published online.

Table 1.

Predictive variables used in our study to model the relative mortality risk, ${R^0}$ . Variables $x_1, \ldots, x_9$ are standardised using a N(0,1) distribution function, $x_{11}$ and $x_{12}$ take values in [0, 1] and $x_{10}$ is a categorical variable taking one of five values explained in Table 2.

The possible values of the first nine numerical variables $x_1, \ldots, x_9$ in Table 1 are on very different scales. For the purpose of visualisation, we standardise them to have mean zero and variance one. Details of the standardisation procedure can be found in the supplementary material.

Variable $x_{10}$ is a categorical variable representing the urban–rural class of an LSOA and taking one of five values listed in Table 2.

Table 2.

Five categories for the urban–rural class (predictive variable $X_{10}$ ).

In summary, the socio-economic characteristics of any neighbourhood are given as a vector taking values in the $K=12$ dimensional space

(2.3) \begin{equation} L_{0} = \mathbb{R}^9 \times \{1, \ldots , 5\} \times [0,1]^2 .\end{equation}

Note, that our urban–rural class indicator $x_{10}$ distinguishes between urban conurbation in London and outside London. We have introduced that distinction as we found in previous research that mortality rates in London are rather different from mortality rates in other parts of England, see Cairns et al. (Reference Cairns, Kleinow and Wen2021) and Wen (Reference Wen2022).

It is to be expected that the covariates in Table 1 are correlated. We report the empirical correlations in Table 3. We observe in Table 3 that there are some strong correlations, but we argue that none of the observed correlations is so strong that a variable should be removed. The few strong correlations we observe might be seen as problematic for parametric models as the parameters for individual covariates might not be identifiable when one variable can act as a proxy for another variable. However, since we have 32,844 observations, the inclusion of highly correlated variables is still meaningful and significantly improves the fit. Additionally, in the random forest (RF) algorithm each decision only depends on one variable, and therefore this algorithm still works effectively when correlations are high – that is one of the reasons for choosing the RF estimator.

Table 3.

Correlations between the covariates, see Table 1 for details about the covariates. Empirical correlations have been calculated using all LSOAs in England and Wales regardless of their urban–rural classification. Note that $x_{10}$ (urban–rural classification) is not included in the table.

We would also argue that there is no surprise in the correlation table: for example $x_1$ is positively correlated with $x_2$ , but negatively correlated with $x_5$ and $x_8$ meaning that the higher the level of income deprivation, the higher the level of employment deprivation and the smaller are the houses, and the fewer are working in management positions. For completeness, we also report correlation tables for individual urban–rural classes in the supplementary material published online. The general conclusions from those tables are similar to those obtained from correlations across all LSOAs. However, correlations with the more-minor variables $x_3$ to $x_9$ do vary more between different urban–rural classes.

3.1. Modelling relative mortality risk

As a starting point, we model the conditional expectation of the relative mortality risk ${R^0}$ given characteristics x:

(3.1) \begin{equation} f(x) \,{:\!=}\, \mathbb{E}[{R^0} | x] \mbox{ for any } x \in L_0\end{equation}

where $x = (x_1, \ldots, x_K)$ is the vector of predictive socio-economic variables taking values in a K-dimensional space $L_0$ of possible realisations of $X_i$ . For the 12 variables in our empirical study, (see Table 1) $L_0$ is the 12-dimensional space defined in (2.3).

To estimate the regression function f, we will use a supervised machine learning algorithm called a RF, and we will denote this estimator of f by $${\hat f^{{\text{RF}}}}$$ . Details of the estimation procedure are given in Section 4. Let us mention that other estimators for f could be used. For example, Wen (Reference Wen2022) compares the RF estimator with a local linear regression estimator.

3.2. Care homes

When we construct the LIFE index as an estimator for the relative risk in any LSOA, we need to take into account that in our sample of LSOAs there are some with care homes and some without care homes. Clearly, if a significant proportion of individuals in any LSOA are living in a care home, then this will increase the mortality rate in that LSOA and, therefore, increase the relative risk. However, this does not then properly reflect the main socio-economic characteristics of the LSOA ( $x_1,\ldots, x_{10}$ ).

To offset this effect when constructing the LIFE index in the next section, we will make assumptions about the proportion of people living in a care home for any LSOA rather than using the actual proportion of people living in care homes in that LSOA. In other words, we are trying to answer the question: What would be the relative risk of dying in LSOA i if we kept all socio-economic variables to the values observed in that LSOA, but changed the proportion of people living in care homes to the average for the whole of England, or to some other chosen value.

3.3. The LIFE index

Based on the discussion so far, we now define our Longevity Index for England as the value of f for specific neighbourhoods using the socio-economic characteristics of this neighbourhood but replacing the proportion of people living in care homes with the average for the whole of England. More precisely, we define the LIFE index for LSOA i as

(3.2) \begin{equation} R_i = f(\tilde X_i) \mbox{ with } \tilde X_i = \left(X_{i,1}, \ldots, X_{i,9},X_{i,10}, \bar X_{11}, \bar X_{12} \right)\end{equation}

where $\bar X_{11}$ and $\bar X_{12}$ denote the average values of the proportion of an LSOA’s population living in care homes with nursing and care homes without nursing, respectively. In replacing the true with the mean values of $\bar X_{11}$ and $\bar X_{12}$ , it is helpful to note from Table 3 that $x_{11}$ and $x_{12}$ have a very low correlation with other socio-economic variables. First, this implies that care homes are not concentrated in neighbourhoods with particular socio-economic characteristics. Second, the lack of correlation means that when we replace $X_{i,11}$ and $X_{i,12}$ with their mean values, we do not need to alter the values of other predictive variables to compensate: that is, the presence of a care home does not artificially inflate or deflate an LSOA’s other socio-economic, predictive variables.

Note that the index could be constructed with other adjustments to the care home variables. For example, we could choose to calculate the index based on setting $X_{i,11} = X_{i,12} = 0$ . That would also be a good choice to model the relative mortality risk of the popoulation not living in care homes. However, we prefer the setting in (3.2) as we will use the index to calculate life expectancies as a function of socio-economic characteristics. Setting $X_{i,11} = X_{i,12}=0$ would implicitly assume that no individuals will ever be in a care home, which is, of course, not reasonable. Additionally, base mortality $m^b_{ta}$ incorporates excess care home deaths, and so we prefer to reflect this also in the index values.

Let us mention that our index is based on the conditional expectation in (3.1) and therefore, we can calculate the relative risk of dying for any values of the socio-economic variables, they do not need to be those observed in a specific LSOA. In other words, we can calculate the LIFE score for fictional neighbourhoods with specified socio-economic characteristics. This allows us to investigate how sensitive the relative mortality risk is with respect to changes in certain socio-economic variables. We will return to that point in Section 8.

From the construction of the LIFE index in (3.2), it is clear that the index is not explicitly age specific. Instead, it is an index summarising the socio-economic characteristics of all members of a small community regardless of their age. However, the LIFE index relies on an estimate of the relative risk function f in (3.2) that links socio-economic characteristics to death rates. The age range for death counts and exposures used to estimate the function f will of course have an impact on the obtained LIFE index values. We will investigate this further in Section 5.

4.1. Overview

Our RF algorithm consists of three stages. For each stage, we use R (R Core Team, 2021) and the R package randomForest (Liaw and Wiener, Reference Liaw and Wiener2002).

4.2. Fitting a single tree

A RF consists of $B>1$ regression trees also known as decision trees. We will here briefly discuss how each tree is constructed as a crude estimator of the regression function f in (3.1). In the next section, we will then turn to combining many trees into a RF.

For fitting an individual tree with index $b \in \{1, \ldots, B\},$ we only use a subset $\mathcal{S}^b$ of the LSOAs in the training set $\mathcal{S}^{\text{train}}$ in both stages. The procedure for growing a tree is the same for stages one and two. The choice of $\mathcal{S}^b$ is explained in Section 4.3.

Constructing an individual tree is an iterative procedure. We start with defining our initial estimator $\hat{f}^{(b)}_0$ as the average of all observed values of the relative risk ${R^0}$ of the LSOAs in the set $\mathcal{S}^b$ , that is,

\[ \hat{f}^{(b)}_0(x) = \frac{1}{|\mathcal{S}^b|} \sum_{i \in \mathcal{S}^b} R^0_i \qquad \mbox{ for all } x \in L_{0} \]

where $|\mathcal{S}^b|$ is the size of the data set $\mathcal{S}^b$ .

In the next step, we choose one explanatory variable, say $x_{k^*}$ , and a level $l^*$ , and split the initial node $L_{0}$ into the two disjoint subsets:

(4.1) \begin{align} L^{b}_{1,1} &= \{x \in L_0 \,{:}\, x_{k^*} < l^*\} \end{align}

(4.2) \begin{align} L^{b}_{1,2} &= \{x \in L_0\,{:}\, x_{k^*} \ge l^*\} \end{align}

This procedure is now repeated but, in addition to choosing an explanatory variable $x_{k^*}$ and a threshold level $l^*$ , we also choose one of the sets (nodes) $L^{b}_{1,1}$ and $L^{b}_{1,2}$ which we then split in the next step.

Starting with $s=1$ and the two nodes defined in (4.1) and (4.2), we now apply the following iterative procedure:

• Choose one subset $L^{b}_{s,j^*}$ ( $j^* \in 1, \ldots, s+1$ ) out of the $s+1$ subsets formed by the first s splits. Also, choose an explanatory variable $x_{k^*}$ and a threshold $l^*$

• Split $L^{b}_{s,j^*}$ into two subsets and leave all other subsets unchanged.

With this procedure, we have the following nodes available after $s+1$ splits:

(4.3) \begin{align}L^{b}_{s+1,j} &= L^{b}_{s,j} \mbox{ for } j=1, \ldots, s+1, \; j \neq j^*\end{align}

(4.4) \begin{align}L^{b}_{s+1,j^*} &= \{x \in L^{b}_{s,j^*} \,{:}\, x_{k^*} < l^*\}\end{align}

(4.5) \begin{align}L^{b}_{s+1,s+2} &= \{x \in L^{b}_{s,j^*} \,{:}\, x_{k^*} \geq l^*\}\end{align}

Equation (4.3) states that split $s+1$ does not affect any nodes other than $L^{b}_{s,j^*}$ . Equations (4.4) and (4.5) mean that all LSOAs with characteristics x in node $L^{b}_{s,j^*}$ for which $x_{k^*} \geq l^*$ are put into a new node $L^{b}_{s+1,s+2}$ so that only those LSOAs with characteristics in $L^{b}_{s,j^*}$ and $x_{k^*} < l^*$ remain in that node and are then contained in $L^{b}_{s+1,j^*}$ after $s+1$ splits.

We now define the estimator $\hat{f}^{(b)}_s(x)$ of the regression function f in (3.1) obtained from one tree b after splitting $L_0$ into $s+1$ nodes as

\[ \hat{f}^{(b)}_s(x) = \sum_{j = 1}^{s+1} r_j \mathbb{I}_{(x \in L^{b}_{s,j})} \]

where $r_j$ is the average observed relative risk ${R^0}$ for LSOAs $X_i$ in node j, that is,

\[r_j = \frac{1}{\sum_{i=1}^N \mathbb{I}_{(X_i \in L^{b}_{s,j})} } \sum_{i=1}^N \mathbb{I}_{(X_i \in L^{b}_{s,j})} R^0_i \qquad \mbox{ for all } j = 1, \ldots, s+1.\]

As explained above, for each new split, we need to choose an existing node $L^{b}_{s,j^*}$ , an explanatory variable $x_{k^*}$ and a threshold $l^*$ . Those are chosen such that the fit of the new estimator $\hat{f}^{(b)}_{s+1}(x) = \hat{f}^{(b)}_{s+1}\left(x;\ j^*, x_{k^*}, l^*\right)$ to the observed values of the relative risk, $R^0_i$ for $i \in \mathcal{S}^b$ , is optimised. More specifically, our choice minimises the residual sum of square

(4.6) \begin{align} \mbox{RSS}^{b}_s \left(\,j, x_n, l\right) &= \sum_{i \in \mathcal{S}^b} \bigg(R^0_i - \hat{f}^{(b)}_{s+1}\left(X_i;\ j, x_n, l\right)\bigg)^2 \end{align}

(4.7) \begin{align} \big(\,j^*, n^*, l^*\big)_s^b &= \underset{j, n, l}{\mathrm{argmin}} \, \mbox{RSS}^b_s \left(\,j, x_n, l\right) \end{align}

Within the RF algorithm, for each split, rather than optimise over all K= 12 of the predictive variables, we optimise over a subset of m variables. This subset is chosen randomly for each new split (see, for example, James et al., Reference James, Witten, Hastie and Tibshirani2013). The purpose of this is to increase the variability and reduce the correlation between individual trees. Without restricting the set of variables we find that trees are very similar. The parameter m is a hyperparameter, and we explain its choice in Section 4.4.

Finally, we stop splitting nodes further as soon as any obtained node contains less than M observations where M is another hyperparameter. In our empirical study, we choose $M=200$ as that choice achieves a good balance between goodness-of-fit and overfitting. Díaz-Uriarte and Alvarez de Andrés (Reference Díaz-Uriarte and Alvarez de Andrés2006) provide a general discussion of the minimum node size M, principles of choosing it, and examples on its potential impact on the model performance and computation time. Wen (Reference Wen2022) has a more detailed discussion around M in the specific context of modelling relative mortality risk at neighbourhood level using the RF algorithm, including its impact on the complexity of underlying trees and the standard deviation of the outcomes produced by individual trees. Although the model’s out of sample performance does not appear to be very sensitive to the choice of M, choosing M to be 200 rather than, say, 5 or 50 significantly saves computation time without sacrificing the predictive power of the RF estimator in our application.

In Figure 1, we illustrate the construction of one tree using our data set of LSOAs. In this example, only two variables are considered for potentially splitting nodes: old-age income deprivation, $x_1$ and employment deprivation, $x_2$ . In total, $s=5$ splits have been performed leaving us with $s+1 = 6$ nodes, and our estimator for f is given by

(4.8) \begin{equation} \hat{f}^{(b)}_5(x) = \left\{ \begin{array}{lcl}0.7117 & \mbox{ for } & x \in L^{b}_{5,1} =\{x\,{:}\, x_1 < -0.418\} \\ 0.9830 & \mbox{ for } & x \in L^{b}_{5,2} =\{x\,{:}\, -0.418 \leq x_1 < 0.351\} \\ 1.2070 & \mbox{ for } & x \in L^{b}_{5,3} =\{x\,{:}\, 0.351 \leq x_1 < 1.06 \mbox{ and } x_2 < 0.701\} \\ 1.4470 & \mbox{ for } & x \in L^{b}_{5,4} =\{x\,{:}\, 0.351 \leq x_1 < 1.06 \mbox{ and } 0.701 \leq x_2\} \\ 1.5660 & \mbox{ for } & x \in L^{b}_{5,5} =\{x\,{:}\, 1.06 \leq x_1 \mbox{ and } x_2 < 1.441 \} \\ 1.8750 & \mbox{ for } & x \in L^{b}_{5,6} =\{x\,{:}\, 1.06 \leq x_1 \mbox{ and } 1.441 \leq x_2 \}\end{array}\right.\end{equation}

Figure 1.

The values of the piecewise constant regression tree $\hat{f}^{(b)}_5(x)$ in (4.8) after five splits. Blue solid lines show the boundaries of the six nodes. Red numbers are the estimated relative risk for all LSOAs in each of the those nodes. Each of the gray dots represents the observed values of old-age income deprivation and employment deprivation for a single LSOA in the training set for this example.

In Figure 2, we show the order of the performed splits: the first three splits are all based on $x_1$ (old-age income deprivation). Only after three splits using $x_1$ , two of the obtained nodes are split using $x_2$ (employment deprivation). This clearly shows, that for this example $x_1$ has a higher explanatory power than $x_2$ since splitting the early nodes in our tree according to the old-age income deprivation score reduces the residual sum of squares $\mbox{RSS}$ more than early splits with respect to employment deprivation would achieve.

Figure 2.

A graphical representation of our example regression tree $\hat{f}^{(b)}_5(x)$ in (4.8), and the residual sum of squares, $\mbox{RSS}_s^b$ (4.6) as a function of the number of splits s for this tree.

To illustrate this specific example further, we also report the residual sum of squares, $\mbox{RSS}_s^b$ , in Figure 2. As expected, we find that the early splits result in the greatest reduction of $\mbox{RSS}_s^b$ .

4.3. Many trees form a RF

Having seen how an individual regression tree is fitted to observations in a set $\mathcal{S}^b$ , we now turn to describing how we choose the sets $\mathcal{S}^b$ and how we combine many trees to obtain our final estimator ${\hat{f}^{\text{RF}}}$ for the regression function f in (3.1).

For each tree $b \in \{1, \ldots, B\}$ , the set $\mathcal{S}^b$ is obtained by (see, for example, James et al., Reference James, Witten, Hastie and Tibshirani2013)

1. 1. sampling randomly with replacement from the training data set $\mathcal{S}^{\text{train}}$ to obtain a sample of the same size as $\mathcal{S}^{\text{train}}$ , and then

2. 2. removing all duplicates from that sample.^{Footnote 2}

Repeating this procedure B times, we obtain B subsets $\mathcal{S}^b \subseteq \mathcal{S}^{\text{train}} $ of the training data set.

To introduce more randomness in the construction of the RF estimator ${\hat{f}^{\text{RF}}}$ we also, as remarked before, restrict the predictive variables considered at each split in any individual tree. Rather than choosing a predictive variable $x_k$ out of all p variables when minimising the residual sum of squares $\mbox{RSS}_s^b$ in (4.6), we follow James et al. (Reference James, Witten, Hastie and Tibshirani2013) and choose $x_k$ from a subset of m predictive variables. As mentioned in Section 4.2, this subset of predictive covariates is randomly chosen for each split within each tree.

So, each tree $b = 1, \ldots B$ is fitted to a randomly chosen subset $\mathcal{S}^b$ of observations from the training data set, and $\mbox{RSS}_s^b$ is optimised with respect to m randomly chosen predictive variables. In this way, we obtain a total of B regression functions $\hat{f}^{(b)}$ . The number m of predictive variables considered for each split of nodes is a hyperparameter and will be chosen in stage one using cross-validation. As mentioned earlier, m is then fixed in stage two.

Our final RF estimator ${\hat{f}^{\text{RF}}}$ for the regression function f in (3.1) is obtained by taking the average over all individual regression trees $\hat{f}^{(b)}$ , that is,

(4.9) \begin{equation} {\hat{f}^{\text{RF}}}(x) = \frac{1}{B}\sum_{b=1}^B \hat{f}^{(b)}(x) \mbox{ for any } x \in L_0 \end{equation}

Note that ${\hat{f}^{\text{RF}}}$ is piecewise constant over the full range of values of $x\in L_0$ as it is an average over a finite number of piecewise constant regression tree functions $\hat{f}^{(b)}$ . However, ${\hat{f}^{\text{RF}}}$ can take many more values compared to any individual tree $\hat{f}^{(b)}$ .

4.4. Hyperparameter selection (stage 1)

With the minimum node size, $M=200$ , fixed, the regression function ${\hat{f}^{\text{RF}}}$ in (4.9) will depend on two further hyperparameters: the number m of predictive variables considered for each split, and the number of trees, B. Both of those parameters need to be chosen. One could argue that we can also choose the size $N^{\text{train}}$ of the training set $\mathcal{S}^{\text{train}}$ , but, for simplicity, we choose that set to include half of the available observations with the other half being included in the validation set $\mathcal{S}^{\text{val}}$ in stage one. For the $N=32,844$ LSOAs in our empirical study, we clearly have that both sets include $N^{\text{train}} = N^{\text{val}} = 16,422$ LSOAs.

The hyperparameters B and m are chosen in stage one in the following way: we fit ${\hat{f}^{\text{RF}}} = {\hat{f}^{\text{RF}}}_{B,m}$ to the data in the training set using different values for B and m. For each combination (B,m) considered, we then evaluate the fit of the obtained estimate ${\hat{f}^{\text{RF}}}_{B,m}$ to the data in the validation set using the mean squared error as criterion.

(4.10) \begin{equation} \mbox{MSE}(B,m) = \frac{1}{N^{\text{val}}}\sum_{i\in \mathcal{S}^{\text{val}}} \frac{\big(R^0_i-{\hat{f}^{\text{RF}}}_{B,m}(X_i)\big)^2}{{\hat{f}^{\text{RF}}}_{B,m}(X_i)/\hat{D}^0_i} = \frac{1}{N^{\text{val}}}\sum_{i\in \mathcal{S}^{\text{val}}} \frac{\big(D_i-\hat{D}^{RF}_i\big)^2}{\hat{D}^{RF}_i} \end{equation}

where $D_i = R^0_i \hat{D}^0_i=\sum_{t,a}D_{ita}$ is the observed total number of deaths across all ages a and years t in LSOA i, and $\hat{D}^{RF}_i = {\hat{f}^{\text{RF}}}_{B,m}(X_i) \hat{D}^0_i$ is the expected total number of deaths adjusted with the fitted relative risk ${\hat{f}^{\text{RF}}}_{B,m}(X_i)$ for LSOA i.

In Figure 3, we plot $\mbox{MSE}(B,m)$ in (4.10) for different values of B (with $m=4$ ) and different values of m (with $B=2500$ ). The figure shows that the out-of-sample performance of the RF is not worsening as more trees are grown, and we choose $B=2500$ as we think this will be a good compromise between computational effort and goodness-of-fit. However, we find that considering $m=4$ predicative variables in (4.6) leads to the smallest mean squared error.

Figure 3.

Validation MSE calculated following (4.10) and over LSOAs in the validation set $\mathcal{S}^{va}$ , of a random forest model trained using the relative risk of England males aged 70–79, and with different numbers of trees B (left plot, with m set as 4) and different numbers of variables considered per split m (right plot, with B set as 2500). The 12 predictive variables outlined in Table 1 are used.

Table 4 summarises our choice of hyper-parameters.

Table 4.

Settings of the final random forest model we use for creating the mortality index for England males.

4.5. Goodness-of-fit (stage 2)

In order to assess the out-of-sample performance of the proposed RF estimator applied to the mortality data for the 32,844 LSOAs, we move on to stage two as mentioned in Section 4.1.

To this end, we randomly split the set of all 32,844 LSOAs into two equally sized subsets: the training set $\mathcal{S}^{\text{train}}$ (different from before) and the test set $\mathcal{S}^{\text{test}}$ . The training set $\mathcal{S}^{\text{train}}$ chosen at this stage is a random sample of $\mathcal{S}$ and independent of the training set chosen in stage 1. The sets $\mathcal{S}^{\text{train}}$ and $\mathcal{S}^{\text{test}}$ are disjoint and $\mathcal{S}^{\text{train}}\cup\mathcal{S}^{\text{test}}=\mathcal{S}$ . We then construct an estimator ${\hat{f}^{\text{RF}}}$ using the data in $\mathcal{S}^{\text{train}}$ and the hyper-parameters in Table 4. To quantify the goodness-of-fit of the obtained estimator ${\hat{f}^{\text{RF}}}$ , we evaluate its out-of-sample fit to data in the test set by calculating the mean squared error as in (4.10) but now considering LSOAs in the test set rather than the validation set, that is,

(4.11) \begin{equation} \mbox{MSE}^{\text{test}} = \frac{1}{N^{\text{test}}}\sum_{i\in \mathcal{S}^{\text{test}}} \frac{\big(R^0_i-{\hat{f}^{\text{RF}}}(X_i)\big)^2}{{\hat{f}^{\text{RF}}}(X_i)/\hat{D}^0_i} \end{equation}

Clearly, the realised values of $\mbox{MSE}^{\text{test}}$ will depend on the randomly chosen LSOAs in $\mathcal{S}^{\text{train}}$ and $\mathcal{S}^{\text{test}}$ . To get an idea of how sensitive the results are to the randomised choice of $\mathcal{S}^{\text{train}}$ and $\mathcal{S}^{\text{test}}$ , we calculate $\mbox{MSE}^{\text{test}}$ for three different splits (rounds) of our data into training and test sets. The results are provided in Table 5 for data based on different age ranges. We find that there is not much variation in the values of $\mbox{MSE}^{\text{test}}$ between rounds. We also see that the goodness-of-fit of ${\hat{f}^{\text{RF}}}$ is much better when it is estimated from mortality data at younger ages.

Table 5.

Test set MSE of the proposed random forest model fitted to three randomly chosen training sets (rounds) for data from different age groups. The applied hyperparameters are listed in Table 4.

4.6. Robustness

The rather small variation of the test set MSEs over different randomly chosen training sets in Table 5 is an indication that the fitted relative risk, and therefore the LIFE index, is a robust estimator of the true underlying relative mortality risk. To investigate robustness further we now split the annual data for the observation period 2001–2018 into two subsets: data for even years 2002, 2004, … and data for odd years 2001, 2003, … This split leaves us with two subsets each consisting of nine years of observations. We chose to split the observation period in this way to avoid any impact of potential trends in the relative mortality risk over time.

We now apply the above methods to obtain estimates ${\hat{f}^{\text{RF}}}(x)$ with data from only one of the two observation subsets and then compare the results.

We present scatter plots of the estimated values of ${\hat{f}^{\text{RF}}}(x)$ based on odd years (horizontal axis) and even years (vertical axis) using mortality data for different age ranges in Figure 4. The plots clearly show that the estimated values of the relative risk for individual LSOAs are very similar when mortality data from different years are used, in particular, there seems to be no systematic differences – this is further evidence that the results of our RF estimator are robust. Any variation we see is most likely due to sampling variation in the deaths counts rather than systematic differences.

Figure 4.

Estimated relative risk over 16,422 LSOAs in $\mathcal{S}^{te}$ by the random forest model trained using LSOAs in $\mathcal{S}^{tr}$ and relative risks of the two year groups, $R^{0,odd}$ and $R^{0,even}$ . x-axis: model trained with odd years; y-axis: model trained with even years. Left to right: relative risk of age 60–69, 70–79 and 80–89.

4.7. Final index values (stage 3)

The hyperparameters have been chosen in stage 1 and goodness-of-fit and robustness assessed in stage 2, and it has been concluded that the RF algorithm has produced a good estimate of f with the chosen parameters. In the final stage, we simply rerun the RF algorithm but, instead of using only half of the data, we use the full set of 32,844 LSOAs.

9.1. Relative risk in different deciles

With the proposed LIFE index, we can also zoom in on specific populations identified by their mortality risk. For example, Figure 11 shows the LIFE index values where we only consider the 10% of LSOAs with the lowest relative risk (left plot) and the 10% of LSOAs with the highest relative risk (right plot).

Figure 11.

Heatmaps of observed relative risk values as a function of income deprivation at ages 65+ and employment deprivation. All other variables are unchanged. The LIFE index is fitted to mortality data for ages 60–69. In the left plot, only LSOAs with $g(i) = 1$ (highest relative risk) are shown, and for the right plot, only LSOAs with $g(i)=10$ (lowest) are used.

More specifically, we introduce a decile function g which identifies the decile k for any LSOA i with $g(i) = k$ for $k \in \{1, \ldots, 10\}$ . The 10% of LSOAs with the highest mortality risk are LSOAs i with $g(i) =1$ and the LSOAs with the lowest relative risk have $g(i) = 10$ .

The plots in Figure 11 clearly show that in both subpopulations the most important variable to explain mortality differences is old-age income deprivation while employment deprivation has less impact.

9.2. Explanatory power of the LIFE index

The LIFE index has been constructed with the aim to explain the mortality risk in an LSOA based on socio-economic variables. Therefore, the question arises whether the LIFE index can indeed differentiate between high and low mortality areas. Or, in other words, are the socio-economic variables and the constructed index able to predict the relative mortality risk in an LSOA. To investigate this, we calculate age standardised mortality rates (ASMR) and ADSMR.

The ASMR for any population in year t is calculated as follows:

(9.1) \begin{equation} ASMR_{gt} = \frac{\sum_{a\in\mathcal{X}}m_{gta}E^s_a}{\sum_{a\in\mathcal{X}}E^s_a} \end{equation}

where $\mathcal{X}$ refers to the age range used and $m_{gta} = D_{gta}/E_{gta}$ is the crude death rate for the underlying population g in year t for age a. The standard population $E^s_a$ used in our study is the European Standard Population^{Footnote 3} (ESP) in 2013.

To investigate how well the LIFE index can distinguish between different levels of mortality, we split the set of all LSOAs into 10 deciles as described above and plot the ASMRs for each decile over time. We then compare our results to ASMRs calculated for deciles obtained from the Index of Multiple Deprivation (IMD), see Wen et al. (Reference Wen, Cairns and Kleinow2021). The results are shown in Figure 12.

Figure 12.

ASMR by deprivation decile based on IMD scores (left) and LIFE scores (right) (log scale). The LIFE index and the ASMRs have been calculated for age groups 60–69 and 80–89. Similar figures for other age groups can be found in the supplementary material published online.

We can see in Figure 12 that the 10 deciles obtained from either index, LIFE or IMD, produce very different mortality rates. The figure also shows that the LIFE index leads to a wider spread of mortality rates meaning that it is better than the IMD in identifying low and high mortality on the basis of socio-economic variables. However, we must keep in mind that the IMD was not designed to predict mortality while the LIFE index was chosen to do so. We also find that both indices show a widening of the mortality gap, the difference between rates for the most deprived as compared to the least deprived. For the IMD, we discussed this issue in detail in Wen et al. (Reference Wen, Cairns and Kleinow2021).

9.3. Explaining regional differences

While the LIFE index only uses information about socio-economic covariates, including an urban–rural class, to predict the relative mortality risk in an LSOA, it might be that LSOAs in different regions in England have different mortality rates although they have the same socio-economic characteristics.

To investigate this question further, we group all LSOAs into nine geographical regions. Table 8 lists the regions and shows the percentage of a region’s LSOAs that belong to the risk groups defined in Section 7 for the LIFE index fitted to ages 60–69.

The table shows that 12.4% of the LSOAs in the North West of England belong to the group of LSOAs with the 5% highest mortality risk, closely followed by the North East. The lowest LIFE scores are observed in the South East with 11.4% of those LSOAs belonging to the 5% English LSOAs with the lowest mortality risk. The results for other age groups are very similar, and, therefore, not reported in this paper.

Table 8.

Distribution of LSOA into regions for different subpopulation classes. The LIFE index has been fitted to ages 60–69.

To obtain a more detailed picture of the regional differences, we consider age standardised mortality rates. The ASMRs by region for age 60–69 and 80–89 can be seen in the Figure 13. We find that there are apparently substantial inequalities between regions, but how much of that can be explained by differences in the socio-economic mix of the nine regions?

Figure 13.

ASMRs by region for mortality data for ages 60–69 and 80–89 (log scale). Similar figures for other age groups can be found in the supplementary material published online.

To address this question, we propose using what we call the Age and Deprivation Standardised Mortality Rate, ADSMR. In the same way that the basic ASMR removes differences between populations that have different age profiles, the ADSMR is designed to remove the impact of different deprivation profiles. Thus, the ADSMR in region r in year t is defined as

(9.2) \begin{equation} ADSMR_{rt} = \frac{1}{10}\sum_{k=1}^{10} ASMR_{rkt} \end{equation}

where $ASMR_{rkt}$ is the ASMR in year t of all LSOAs in region r and index decile k. Deciles are formed on the basis of the LIFE index values as described above, and, for comparison, on an IMD basis.

Assuming that all mortality differences are explained by the LIFE (or IMD) index based on socio-economic variables, the ADSMRs should be the same for different regions (since $ASMR_{rkt}$ would be independent of r). In Figure 14, we plot $ADSMR_{rt}$ for the nine regions where deciles have been obtained from using the IMD (left plots) and the LIFE index (right plots).

Figure 14.

ADSMRs on the basis of IMD deciles (left) and LIFE deciles (right) where mortality data for different age groups have been used (log scale). Similar figures for other age groups can be found in the supplementary material published online.

We clearly see that the purpose built LIFE index is better suited to explain mortality differences between regions, resulting in much smaller differences in the ADSMR between the regions than the IMD-based ADSMR. While this is true for both age ranges (and other ranges as shown in Figures published online in the supplementary material), it is for the oldest age group 80–89 that the LIFE index can explain most of the differences between regions – the ADSMRs in the bottom right plot are much more similar than in the bottom left plot. This can also be observed for age group 70–79, see the online supplementary material for the relevant graphs. The reason for the much better ability of the LIFE index to explain differences between regions for higher age groups as compared to the IMD might be the inclusion of old-age income deprivation in the LIFE index rather than general income deprivation (which is one of the domains of deprivation used for the IMD).

Another interesting feature we observe in Figure 14 and similar figures in the online supplementary material are the mortality improvement rates in London which seem to be much greater than in other regions. Considering the oldest age group 80–89 and measuring deprivation using the IMD suggests that mortality in London improved significantly more than in other areas even after accounting for deprivation and that this improvement continued after 2011 when other regions have experienced no or little improvements. However, measuring deprivation using the LIFE index changes this conclusion. The ADSMRs for the nine regions are very close together, and there is no London-effect visible for this age range.