Spatio-Temporal Smoothing, Interpolation, and Prediction of Income Distributions Based on Grouped Data

Privacy concerns typically restrict access to information regarding peoples’ incomes. Rather, income data are typically available as grouped data, wherein individual incomes are grouped into a set of known income classes (e.g. annual incomes between JPY $3$ and $4$ million), followed by counting the number of individuals per class. As the amount of information in grouped data is considerably limited because of the grouping mechanism, the underlying structure cannot be expressly revealed. However, numerous studies in the literature on income distribution estimation explored various statistical methods and income models to estimate income distributions (see Kleiber and Kotz (2003), Chotikapanich (2008), and Jorda et al. (2021) for comprehensive reviews).

In Japan, the Housing and Land Survey (HLS) includes the municipality-level grouped income data of varying numbers of income classes over time. The HLS data can provide detailed information for effective policymaking with fine granularity. However, as HLS does not cover every municipality, the non-sampled areas must be spatially interpolated to create a complete map of an income or inequality measure. Additionally, HLS is conducted every five years; the public can access the results after a year of the survey. Therefore, the temporal prediction of the latest income status is necessary for an arbitrary time point to allow policymaking based on the latest information on income status. The proposed model can effectively overcome these challenges through the spatio-temporal heterogeneity embedded in the mixing proportions, as well as a set of covariates in the component distributions.

The extant studies on income distribution predominantly focused on estimating the income distribution of a single area at a single time point, such as the national-level income distribution estimated from a single round of an income survey. Examples of the existing methods (and countries to which they have been applied) include the maximum-likelihood estimation for the flexible five parameter beta distribution (McDonald and Xu, 1995) (the United States), generalised method of moments estimation for grouped income data (Hajargasht et al., 2012; Griffiths and Hajargasht, 2015) (China, India, Russia, Pakistan Poland, and Brazil) and its generalisation (Chen, 2018) (the United States and China), Bayesian estimation of parametric distributions for grouped income data (Kakamu and Nishino, 2019; Eckernkemper and Gribisch, 2021) (Japan, India, Peru, Ethiopia, and Iraq), maximum-likelihood and Bayesian estimations of a parametric Lorenz curve based on the Dirichlet likelihood (Chotikapanich and Griffiths, 2002, 2005) (Sweden and Brazil) and finite mixture of log-normal distributions for individual incomes (Flachaire and Nuñez, 2007; Lubrano and Ndoye, 2016) (the United Kingdom), grouped-data-based estimation of the parametric Lorenz curves via approximate Bayesian computation (Kobayashi and Kakamu, 2019) (Japan), and finite mixture of the beta distributions of the second kind (Chotikapanich et al., 2007) (Hong Kong, Japan, Malaysia, Philippines, Singapore, Korea, Taiwan, and Thailand).

The national-level or single-area income distributions and their associated inequalities may be estimated with some degree of accuracy using grouped data from a survey based on many households or individuals. However, they do not provide useful information for local-level policymaking, such as state- or city-level policymaking, which is very crucial for effective policymaking with fine granularity. Further, even with local-level data (as is the case for HLS), separately estimating the income distribution of each area is challenging because of the loss of information during grouping. Therefore, a joint statistical model that borrows the strength of each dataset is desired to increase the statistical efficiency of the estimation.

The joint estimation of income distributions over multiple areas and/or periods represents a new stream in the literature on the statistical analysis of income distributions. It is advantageous because of the recent availability of abundant data for statistical analyses. However, only a few studies have considered this approach despite its statistical and economic significance. For example, Nishino et al. (2012) and Nishino and Kakamu (2015) considered the log-normal income model comprising time-varying scale parameters. Kobayashi et al. (2022) considered the grouped-data-based state-space model for the Lorenz curve to stably estimate the Lorenz curves and Gini indices. These models could also be used for temporal predictions. Sugasawa et al. (2020) considered spatially varying income distributions based on a parametric family of the income distribution. This approach produces spatially smooth estimates of income distributions and enabled the interpolation of non-sampled areas. Although the cited studies demonstrated how joint estimation improves efficiency, they only considered spatial or temporal dependence, partly because of data availability. Additionally, the use of a parametric family can be restrictive because a suitable family may differ with the area or period. The four studies analysed grouped Japanese income data.

Some recent studies on the small area estimation (SAE) of some income measures are also worth mentioning. For instance, Marhuenda et al. (2013) explored the Fay–Herriot model exhibiting a spatio-temporal correlation structure to analyse the Italian income data. As the Fay–Herriot model is a model for a direct estimator, such as the median income, it cannot be used to treat grouped data. Kawakubo and Kobayashi (2023) considered a version of the area-level linear mixed model, and Walter et al. (2021) considered a unit-level small area model for grouped data. These methods were applied to treat Japanese and Mexican data, respectively. Although their method could also predict general small area parameters in non-sampled areas, it had no spatial consideration. Additionally, as their method was designed for a single period, it is not appropriate for temporal predictions. Furthermore, as the method of Walter et al. (2021) requires individual-level grouped income data, it cannot be applied to the analyses of HLS data. Gardini et al. (2022) considered a finite mixture of log-normal distributions for the unit-level small area model to stably estimate Italy’s income and poverty measures.

Based on the foregoing, we present a novel grouped-data-based framework for the spatio-temporal smoothing, interpolation, and prediction of areal income distributions by considering spatio-temporal changes. To overcome the limitations of the parametric approach, we considered a flexible finite mixture model. A unique feature of our modelling framework is that all the areas share common component distributions that do not change with time or space, apart from the covariate values. The area-wise heterogeneity of income distributions is characterised by area-wise mixing proportions that are hierarchically modelled. This approach is advantageous because of the possibility of suppressing the number of parameters and latent variables subject to estimation via the common component distributions and the possibility of stably estimating the model, even from grouped data. Further, it also offers an interesting interpretation in which the income distribution of each area and period is represented by a set of template income distributions that are mixed by weights change over space and time. Sugasawa et al. (2019) and Sugasawa (2021) adopted a similar mixture-modelling strategy, though both studies did not consider the spatio-temporal setting. Through the spatio-temporal mixing proportions, the proposed model can provide spatial interpolation and temporal prediction for an arbitrary spatial unit and at an arbitrary time point, which are the pre-requisites for analysing HLS data. As aforementioned, Flachaire and Nuñez (2007) and Lubrano and Ndoye (2016) also considered the finite mixture approach for modelling income distributions. However, their models were designed for national-level individual income data and not for grouped income data or municipality-level analyses. Furthermore, although they considered the income data of the United Kingdom over multiple periods, they estimated the mixture model independently for each period as their models did not include any temporal structure.

Some relevant extant studies on mixture models under spatio-temporal settings are as follows. Fernández and Green (2002) considered a spatial mixture model for areal data with a variable number of components and spatially varying mixing proportions. Neelon et al. (2014) considered a multivariate normal mixture model for areal data. The model included spatial random effects in the component distributions and mixing proportions. Considering the Poisson mixture, Hossain et al. (2014) proposed the algorithms for variable and model selections, as well as for relabeling posterior simulation. Viroli (2011) proposed a general finite mixture model for three-way arrays based on matrix normal distributions. A model for analysing spatio-temporal multivariate data presents a special case in which the spatial effects are incorporated in the mixing proportions. Paci and Finazzi (2018) and Vanhatalo et al. (2021) considered the spatio-temporal mixture models for point-referenced data. Several studies employed the Bayesian nonparametric mixture approaches to analyse spatio-temporal data. These approaches are built around or connected to the Dirichlet process mixture models (Kottas et al., 2008; Zhang et al., 2016; Youngmin and Kim, 2020; Wang et al., 2022).

The remainder of this paper is organised as follows. Section 2 describes the challenges of analysing the HLS data of Japan. Section 3 describes the spatio-temporal mixture modelling for the grouped income data. Section 4 presents the simulation studies to demonstrate the performance of the proposed method. Section 5 presents the application of the proposed method to the analyses of HLS data. Finally, Section 6 provides the summarised conclusions.

HLS^†††The HLS data are available at https://www.stat.go.jp/english/data/jyutaku/index.html is conducted by Japan’s Statistics Bureau to obtain primary data for formulating housing-related policies and to obtain information on the actual status and ownership of houses and buildings, as well as the status of households living in the houses and utilising the land. The HLS data also contain information on household incomes in the form of municipality-level grouped data, where the income classes vary from one round to another. Notably, HLS does not sample every municipality. The set of sampled municipalities varies with each round.

Our analysis features all the available HLS rounds of 1998, 2003, 2008, 2013, and 2018 ($T=5$) HLS. Table 1 presents the income classes for each round. During the 20 years between the first and last rounds, the municipal units were subjected to numerous changes, partly to ensure efficient local taxation. Some municipalities were absorbed into larger ones. Municipality groups were merged to create new municipalities, and others were simply upgraded from villages to cities. We tracked all these changes and compiled the grouped income data, thereby aggregating the observations in the merged and absorbed municipalities with the list of 1741 municipalities for 2018 as the standard. The resulting numbers of sampled municipalities were 814, 1040, 1117, 1110, and 1086 for 1998, 2003, 2008, 2013, and 2018, respectively.

Figure 1 presents the proportions of the nine income classes in the 2018 HLS data. The figure reveals that many inland municipalities were not sampled and that their data were not available. Based on the grouped data, the average incomes, as well as other income or poverty measures, could be computed crudely, at least for the sampled municipalities. However, as the top-income class is open-ended, the construction of a crude estimate is quite arbitrary, and its uncertainty quantification may be impossible.

Most municipalities have the largest proportions of residents in the third (between JPY 2 million and 3 million JPY) and fourth (between JPY 3 million and 4 million JPY) income classes. The municipalities, which are far from urban areas, such as those in the southwestern and northern parts of Japan, tend to have large proportions of the first (below 1 million JPY) and second (between 1 million and 2 million JPY) income classes. The plots for the other rounds of HLS are presented in the Supplementary Material.

HLS is conducted every five years, and the results are not made publicly available after at least one year. For flexible and rapid policymaking with fine granularity in time and space, it is desirable to make the information regarding the income states of all the municipalities available and updated periodically, e.g. annually.

In summary, although HLS collects some detailed data on household incomes across Japan, an analyst would face the following challenges:

1. 1.

   the data would be in the form of grouped data from which income and poverty measures cannot be readily obtained;

2. 2.

   the data would include many non-sampled municipalities;

3. 3.

   the survey would be conducted every five years even though one would desire to know the latest income status, for example, annually.

These challenges can be overcome by the spatio-temporal mixture model proposed in the following section.

The proposed method is first demonstrated using simulated data. The $M=200$ areal units on the two-dimensional space are generated by sampling each coordinate from $(d_{i1},d_{i2})\sim U(-1,1)^{2}$. The neighbours of an area are defined as those with coordinates within a $0.2$ Euclidean distance, and the average number of neighbours is $5.77$. In each area, $N_{it}$ is generated from $U(100,500)$ for $21$ periods. The groups are represented by $[0,1)$, $[1,2)$, $[2,3)$, $[3,4)$, $[4,5)$, $[5,7)$, $[7,10)$, $[10,15)$, and $[15,\infty)$ with $G=9$ for all $t$. These corresponded to those in the latter three HLS rounds (Table 1). In this simulation study, the grouped data are generated from the mixture model, (1), with $K=3$. The parameters of the component distribution are set, as follows: ${\text{\boldmath$\beta$}}_{1}=(0.5,0,1)^{t}$, ${\text{\boldmath$\beta$}}_{2}=(-0.5,1,0)^{t}$, ${\text{\boldmath$\beta$}}_{3}=(2.0,-1,-1)^{t}$, $(\sigma^{2}_{1},\sigma_{2}^{2},\sigma_{3}^{2})=(0.5,0.5,0.5)$. Each covariate, $x_{itj}$, is generated from the uniform distribution.

Using the same component distributions, the following two settings for the mixing proportions are considered. In the first setting, the mixing proportions are generated from (3) with $(\mu_{2},\mu_{3})=(-0.2,0.1)$, $(\rho_{2},\rho_{3})=(0.8,0.8)$, $(\tau_{2},\tau_{3})=(0.1,0.1)$, $(\alpha_{2},\alpha_{3})=(0.2,0.2)$, $(\eta_{02},\eta_{03})=(0,0)$. In the second setting, the spatial and temporal effects are determined by the following deterministic sequences: $u_{i2}=3(d_{i1}-d_{i2})$, $u_{i3}=3(d_{i2}-d_{i1})$, $\eta_{t2}=t/3-(T+1)/6$, $\eta_{t3}=t/6-T/12$. Furthermore, the two-dimensional space, $(-1,1)^{2}$, is partitioned into $5\times 5$ square blocks. The areas in the same block exert the same effect on the mixing proportions: $a_{h_{i}k}\sim N(0,0.2^{2}),\ h_{i}\in\{1,\dots,25\},\ i=1,\dots,200,\ k=2,3$. In this setting, the area, $i$, is a sub-area of the $h_{i}$th block. Thus, the mixing proportions for the $i$th area corresponded to the $h_{i}$th block are proportional to $\exp(\mu_{k}+u_{ik}+\eta_{tk}+a_{h_{i}k})$ for $k=2,3$. The average number of sub-areas in each block is $8$.

First, we fit the proposed spatio-temporal mixture model described in Section 3 with different component numbers: $K=2,\dots,6$. The proposed Gibbs sampler is run for 20,000 iterations, with an initial burn-in period of 10,000 iterations. Figure 2 shows the matching fractions for the proposed mixture models under both settings. The matching fractions are one as long as $K\leq 3$ in both settings. However, the figure shows that the matching fractions become less than one for the over-specified models with $K=4$, $5$, and $6$, indicating the overfitting of the mixture models. The figure indicates that the monitoring of the matching fraction from a small value of $K$ and the selection of an appropriate $K$ at the kink are reasonable approaches.

Next, we compare the proposed spatio-temporal mixture model with four alternative models. The first alternative model is the mixture model based on (1) and (3). The mixing proportions are modelled with the two-way independent random effects that represent the heterogeneity in space and time: $u_{ik}\sim N(0,\tau_{k})$, $\eta_{tk}\sim N(0,\alpha_{k})$, respectively. As this model exerts spatial and temporal effects, it can be used for spatial interpolation and temporal prediction. The second alternative model is the mixture of the log-normal distributions with the mixing proportions that include only $\mu_{k}$ and $u_{ik}$. It is a purely spatial model, which is independently fitted for each $t$ and used for spatial interpolation only. This model may be considered as a spatial mixture extension of the small area model for grouped data proposed by Kawakubo and Kobayashi (2023). The third alternative is the simple log-normal model for grouped data. It is independently estimated using the maximum likelihood for each area and period. The fourth alternative is an SAE model. More specifically, we consider the linear mixed model proposed by Torabi and Rao (2014), as recommended by one of the reviewers. It is given by the following:

$$\hat{y}_{bj}={\text{\boldmath$x$}}^{t}_{bj}{\text{\boldmath$\beta$}}+v_{b}+u_{bj}+e_{bj},\quad b=1,\dots,25,\quad j=1,\dots,M_{b},$$

where $M_{b}$ is the number of sub-areas in the area, $b$, such that $\sum_{b=1}^{25}M_{b}=200$, $\hat{y}_{bj}$ is the crude average income described below, $v_{b}\sim N(0,\sigma^{2}_{v})$ is the area effect, $u_{bj}\sim N(0,\sigma^{2}_{u})$ is the sub-area effect, and $e_{bj}\sim N(0,\sigma^{2}_{ej})$ is the sampling error. The average incomes are computed based on $N_{it}^{-1}\sum_{g=1}^{G}N_{itg}{\rm mid}_{tg}$, where ${\rm mid}_{tg}$ is the mid-point of the income classes. For the top-income class ($g=G$) with an open-ended interval, ${\rm mid}_{tG}$ is set to be $20$. This model is fitted in Setting 2, whereas the model without $v_{b}$ is fitted in Setting 1 as the true model does not have a sub-area structure. The small area models are estimated independently per period.

The data for the $m=150$ areas and $T=20$ periods are used for model estimation. The data for the remaining $m^{*}=50$ areas and the last period are subjected to spatial interpolation and temporal prediction.

The models are first compared based on the root mean squared errors (RMSE) and coverage of the 95% credible or prediction intervals for the average incomes in the in-sample estimation, spatial interpolation, and temporal prediction. For example, in the in-sample estimation case, RMSE is defined by $\sqrt{\frac{1}{mT}\sum_{i=1}^{m}\sum_{t=1}^{T}(\widehat{\rm AI}_{it}-{\rm AI}_{it})^{2}}$, and the coverage is defined by $\frac{1}{mT}\sum_{i=1}^{m}\sum_{t=1}^{T}I\left\{L_{it}\leq\text{AI}_{it}\leq U_{it}\right\}$, where ${\rm AI}_{it}$ is the truth, $\widehat{\rm AI}_{it}$ is the posterior mean or maximum likelihood estimate for the sampled areas and $L_{it}$ and $U_{it}$ are the endpoints of the 95% credible interval under the Bayesian models. Regarding the spatial interpolation and temporal prediction cases, $\widehat{\rm AI}_{it}$ corresponds to the posterior predictive mean for the non-sampled areas with $i=m+1,\dots,M$ and $t=1,\dots,T$; it also corresponds to the posterior predictive mean for all the areas with $t=T+1$. Similarly, $L_{it}$ and $U_{it}$ correspond to the endpoints of the 95% prediction intervals.

Table 2 presents the RMSEs and coverages for the average incomes under the five models. Generally, the proposed model outperform the other models regarding the RMSE. The coverages are around the targeted value of $0.95$. Regarding the two-way model, although the RMSEs in the in-sample estimation are the second smallest in both settings, following the proposed model, it produced high RMSEs in the spatial interpolation and temporal prediction. In addition, the coverages are lower than the target value of $0.95$. The RMSEs of the spatial model regarding the spatial interpolation are the second smallest in both settings, following the proposed model, with the coverages of approximately $0.95$. The results indicate the efficacy of the spatio-temporal effects, and the proposed model comprising both effects demonstrated superior performance. Estimating the grouped-data model for each area and period separately via maximum likelihood resulted in high RMSEs in the in-sample estimation because of the lack of borrowing information across space and time. The employed SAE method is not designed specifically for grouped data or a spatio-temporal setting. Thus, it produced the largest RMSEs for the in-sample estimation and spatial interpolation.

Finally, the proposed, two-way, and spatial models are compared based on the RMSE and coverage for the number of households. Table 3 presents the results for both settings. The proposed model with the spatial and temporal effects performs well, exhibiting the smallest RMSEs, except for the in-sample estimation in Setting 2 and the coverages around the targeted $0.95$. The two-way model yielded the largest RMSEs in all cases. The spatial model produced smaller RMSEs than the two-way model regarding the spatial interpolation. However, the proposed model that borrows information over time produced smaller RMSEs than those of the spatial model.

In this study, we proposed a novel spatio-temporal mixture model for analysing income distribution based on grouped data. Through the covariate information and spatio-temporal heterogeneous effects, the proposed model allows for the inference and prediction of the income distributions of sampled and non-sampled areas at any time point. When an appropriate covariate is available at a higher frequency than the grouped income data, a great insight into the latest income status can be obtained through predictions. By applying the proposed method to the HLS data, we obtained the complete maps of the average incomes, median incomes, and Gini indices. We observed that the overall average and median incomes decreased over the past 20 years and that the variation in the Gini indices decreased over the same period. We also demonstrated that adopting an appropriate covariate was essential to revealing the income statuses of the areas. This is particularly true for the non-sampled areas, where the covariate information is the only available information that can be used to infer income statuses.

The proposed method can reveal the income distributions and associated poverty measures of any area and at any period as long as their covariate information is available. While it was not shown in the present paper, it is also possible to estimate the proportions of the population in each area whose incomes are below any given threshold income levels. Exploring the results of the proposed method, policymakers can identify areas that require political intervention to mitigate poverty and inequality across the nation through efficient policy implementation and resource allocation. Such political interventions can be implemented based on the latest information on income distributions, as elucidated by the spatial interpolation and temporal prediction of the proposed model.

This study adopts the full MCMC approach to the posterior inferences. The advantage of the proposed Gibbs sampler is its non-requirement of algorithm tuning. As the inclusion of many latent variables may cause some inefficiencies, we may argue for the need for an alternative computational strategy. In fact, we also attempted to implement an MH algorithm and its adaptive version to sample ${\text{\boldmath$\beta$}}_{k}$ and $\sigma^{2}_{k}$ without generating individual household incomes. Put differently, (2) was directly used to contribute to the likelihood function. However, the MH algorithms failed to explore the parameter space and mostly stuck around the initial values. Some other faster posterior-computing methods, such as the integrated nested Laplace approximation (Rue et al., 2009) and the variational Bayes approximation (Blei et al., 2017), are appealing alternatives. Regarding the application to the HLS data, the computing time, e.g. for $K=4$ was approximately 40 hours based on non-optimised R code on the Apple M1 Ultra chip without any parallelisation. As we prefer to ensure the quality of the posterior inference and because the computing time was not the main issue in our application, we considered the present computational approach acceptable. However, developing a more computationally efficient algorithm would also be a relevant research direction that must be explored.

This work was supported by JSPS KAKENHI (#22K01421, #21K01421, #21H00699, #20H00080 and #18H03628) and Japan Center for Economic Research.