Identifying Linear Models in Multi-Resolution Population Data using Minimum Description Length Principle to Predict Household Income

TL;DR

This paper introduces a framework using regression and MDL to identify large regions with common income predictors, aiding targeted policy-making in resource-limited settings.

Contribution

It formalizes and solves the problem of finding largest regions with shared income models using MDL, demonstrated on real-world Thai data.

Findings

Framework outperforms baseline methods

Identifies key income indicators for Thai regions

Supports policy-making for income improvement

Abstract

One shirt size cannot fit everybody, while we cannot make a unique shirt that fits perfectly for everyone because of resource limitation. This analogy is true for the policy making. Policy makers cannot establish a single policy to solve all problems for all regions because each region has its own unique issue. In the other extreme, policy makers also cannot create a policy for each small village due to the resource limitation. Would it be better if we can find a set of largest regions such that the population of each region within this set has common issues and we can establish a single policy for them? In this work, we propose a framework using regression analysis and minimum description length (MDL) to find a set of largest areas that have common indicators, which can be used to predict household incomes efficiently. Given a set of household features, and a multi-resolution partition that represents administrative divisions, our framework reports a set C* of largest subdivisions that have a common model for population-income prediction. We formalize a problem of finding C* and propose the algorithm as a solution. We use both simulation datasets as well as a real-world dataset of Thailand's population household information to demonstrate our framework performance and application. The results show that our framework performance is better than the baseline methods. We show the results of our method can be used to find indicators of income prediction for many areas in Thailand. By increasing these indicator values, we expect people in these areas to gain more incomes. Hence, the policy makers can plan to establish the policies by using these indicators in our results as a guideline to solve low-income issues. Our framework can be used to support policy makers to establish policies regarding any other dependent variable beyond incomes in order to combat poverty and other issues.

Figures

Tables

Table 1. Table 1. Comparison between Mixture of Regression and our proposed framework.

Properties\Frameworks	Mixture of Regression	Proposed framework
Assign individuals into partitions	Yes	Yes
Input: multi-resolution partitions	Not required	Required
Input: number of clusters	Required	Not required
Homogeneity of model of clusters	Not report	Report
Time complexity of exact solution	No upper bound	Have upper bound

Table 2. Table 2. Notations and symbols

Term and notation	Description
$n$	Number of individuals
$d$	Number of independent variables
$l$	Number of parameters
$D=\{1,\dots ,n\}$	Set of individual indices
$C$	Cluster where $C\subseteq D$
$𝒞=\left\{C_{j,k}\right\}$	Set of Multi-resolution clusters where $C_{j,k}$ is a $j$th cluster at $k$th layer.
${𝒞}^{\ast }=\left\{C\right\}$	Maximal homogeneous partition
$X,x$	Independent variable and its realization
$Y,y$	Dependent variable and its realization
$S=\{(x_1,y_1),\dots ,(x_n,y_n)\}$	Dataset contains realizations of $X$ and $Y$
$ℒ\left(h\right)$	Number of bits of representation of a model $h$
$ℒ\left(S|h\right)$	Number of bits of representation of dataset $S$ given a model $h$
${ℒ}_{ℝ}\left(x\right)$	Number of bits of representation of real-number vector $x$
$\gamma$	Homogeneity threshold
$H,\mathscr{H},ℍ$	Function class, set of models, and hypothesis space
$ℐ({𝒞}^{\prime },H_1,H_2)$	Model Information Reduction Ratio where $𝒞$ is a partition of $D$
$ℐ({𝒞}_1^{\prime },{𝒞}_2^{\prime },H)$	Cluster Information Reduction Ratio where ${𝒞}_1^{\prime },{𝒞}_2^{\prime }$ are partitions of $D$

Table 3. Table 3. The official dimensions of MPI that policy makers of Thailand currently use to design policies that are related to poverty issues.

Main dimensions	Subdimensions
	Birth weight records
	Hygiene & healthy diet
	Accessing to necessary medicines
Health	Working out habits
	Living in a reliable house
	Accessing to clean water
	Getting enough water for consumption
Living conditions	Living in a tidy house
	Children as a pre-school age are prepared for a school
	Children as a school age can attend to mandatory education
	Everyone in household can attend at least high-school education
Education	Everyone in household can read
	Adults (age 15-59) have reliable jobs
	Seniors (age 60+) have incomes
Financial status	Average income of household members
	Seniors can access public services in need
Access to public services	People with disabilities can access public services in need

Table 4. Table 4. The average RMSE values of approaches on each dataset type. Each element in the table is the average of RMSE from 100 datasets per dataset type. Each row represents results from a specific type of datasets. Each column represents a different method.

	OPT	Greedy	Linear Regression	Mixture of Regression	$y_i-\overline{y}$
Type-1 Dataset					$6.30$
Type-2 Dataset			$5.30$	$0.56$	$6.12$
Type-3 Dataset			$5.94$	$1.90$	$6.35$
Type-4 Dataset			$5.58$	$12.67$	$6.45$

Table 5. Table 5. The average precision, recall, and F1 score values of approaches on each dataset type. Each element in the table is the average value from 100 datasets per dataset type.

	OPT			Greedy			Mixture of Regression
	Precision	Recall	F1	Precision	Recall	F1	Precision	Recall	F1
Type-1 Dataset	1	1	1	0.01	0.01	0.01	1	1	1
Type-2 Dataset	1	1	1	0.22	0.22	0.22	0.99	1	0.99
Type-3 Dataset	1	1	1	1	1	1	0.96	1	0.97
Type-4 Dataset	1	1	1	0.33	0.45	0.38	0.76	0.97	0.82

Table 6. Table 6. The average RMSE values of approaches on nonlinear datasets: polynomial and exponential datasets. Each element in the table is the average of RMSE from 100 datasets per dataset type. Each row represents results from a specific type of datasets. Each column represents a different method. The (LIN) in rows implies that we used the linear function for all methods except y i − y ¯ subscript 𝑦 𝑖 ¯ 𝑦 y_{i}-\bar{y} to fit the data while the (EXP) in rows implies we used exponential function to fit the data except y i − y ¯ subscript 𝑦 𝑖 ¯ 𝑦 y_{i}-\bar{y} .

	OPT	Greedy	Linear Regression	Mixture of Regression	$y_i-\overline{y}$
EXP Dataset (LIN)	8.41	7.42	12.91	5.00	13.99
Poly Dataset (LIN)	16.64	13.46	22.74	12.13	24.73
EXP Dataset (EXP)			12.26	5.16	14.20
Poly Dataset (EXP)	14.59	13.45	22.43	12.40	24.68

Table 7. Table 7. The average precision, recall, and F1 score values of approaches on nonlinear datasets: polynomial and exponential datasets. Each element in the table is the average value from 100 datasets per dataset type. The (LIN) in the rows implies that we used the linear function for all methods to fit the data while the (EXP) in rows implies we used exponential function to fit the data.

	OPT			Greedy			Mixture of Regression
	Precision	Recall	F1	Precision	Recall	F1	Precision	Recall	F1
EXP Dataset (LIN)	0.83	0.87	0.85	0.30	0.43	0.35	0.19	0.16	0.17
Poly Dataset (LIN)	0.53	0.57	0.55	0.30	0.43	0.35	0.16	0.16	0.16
EXP Dataset (EXP)	1.00	1.00	1.00	0.33	0.44	0.38	0.19	0.18	0.19
Poly Dataset (EXP)	0.92	0.94	0.93	0.30	0.43	0.35	0.16	0.15	0.15

Table 8. Table 8. The RMSE values of approaches on each province population. Each element in the table is the RMSE value. Each row represents results from a specific province. Each column represents a different method.

	OPT	Greedy	Linear Regression	$y_i-\overline{y}$
Khon Kaen	69,620.3	67,210	82,571	84,994.6
Chiang Mai	89,067.1	87,727.5	104,938	106,511

Table 9. Table 9. The number of clusters that are the members of output partition in each method.

	OPT				Greedy
	1st	2nd	3rd	4th	1st	2nd	3rd	4th
Khon Kaen	0	2	55	1,785	0	0	1	2,643
Chiang Mai	0	1	41	1,595	0	0	0	2,223

Equations

$$S=\{(x_1,y_1),\dots ,(x_n,y_n)\}.$$

$$Y_i=h^{\ast }(X_i)+ϵ.$$

$$\mathcal{L}(S{)}_{\mathbb{H}}=\underset{h\in \mathbb{H}}{\min }\mathcal{L}(h)+\mathcal{L}(S|h)$$

$$\mathcal{L}(h)=\sum _{j=1}^l{\mathcal{L}}_{\mathbb{R}}(p_j)$$

$$\mathcal{L}(S|h)=\sum _{i=1}^n{\mathcal{L}}_{\mathbb{R}}(x_i)+\sum _{i=1}^n{\mathcal{L}}_{\mathbb{R}}(y_i-h(x_i))$$

$$\mathcal{L}(\mathcal{C}^{\prime},H)=\sum_{C\in\mathcal{C}^{\prime}}\big{(}\mathcal{L}(H)_{C}+\mathcal{L}(C|H)\big{)}$$

$$\mathcal{L}(H{)}_C={\mathcal{L}}_{\mathbb{R}}({\hat{\beta }}_1^C)+{\mathcal{L}}_{\mathbb{R}}({\hat{\beta }}^C)$$

$$\mathcal{L}(C|H)=\sum _{i\in C}{\mathcal{L}}_{\mathbb{R}}(x_i)+\sum _{i\in C}{\mathcal{L}}_{\mathbb{R}}(y_i-{\hat{\beta }}_1^C-x_i^{\text{T}}{\hat{\beta }}^C)$$

$$\sum _{j=1}^k\sum _{i\in C_j}{\mathcal{L}}_{\mathbb{R}}(y_i-g_j(x_i))=\sum _{j=1}^k\sum _{i\in C_j}{\mathcal{L}}_{\mathbb{R}}(y_i-{\hat{\beta }}_{1,j}-x_i^{\text{T}}{\hat{\beta }}_j))\le |C^{\prime }|{\mathcal{L}}_{\mathbb{R}}({ϵ}^{\ast }).$$

$$\sum _{j=1}^k\sum _{i\in C_j}{\mathcal{L}}_{\mathbb{R}}(y_i-g_j(x_i))\le |C^{\prime }|{\mathcal{L}}_{\mathbb{R}}({ϵ}^{\ast })<\sum _{i\in C^{\prime }}{\mathcal{L}}_{\mathbb{R}}(y_i-g^{\prime }(x_i)).$$

$$\displaystyle\mathcal{L}(\mathcal{C}^{\prime},H)-\mathcal{L}(\mathcal{C}^{*},H)=\sum_{C^{\prime}_{i}\in\mathcal{C}^{\prime}}\big{(}\mathcal{L}(S_{C^{\prime}_{i}}|H)+\mathcal{L}(H)_{C^{\prime}_{i}}\big{)}-\mathcal{L}(S_{C}|H)-\mathcal{L}(H)_{C}$$

$$\mathcal{L}({\mathcal{C}}^{\prime },H)-\mathcal{L}({\mathcal{C}}^{\ast },H)=(|{\mathcal{C}}^{\prime }|-1)\mathcal{L}(H{)}_C>0$$

$$\mathcal{L}({\mathcal{C}}^{\prime },H)-\mathcal{L}({\mathcal{C}}^{\ast },H)=\mathcal{L}(S_{C^{\prime }}|H)+\mathcal{L}(H{)}_{C^{\prime }}-\sum _{C_i\in \mathcal{C}}(\mathcal{L}(S_{C_i}|H)+\mathcal{L}(H{)}_{C_i})$$

$$\displaystyle=\sum_{i\in C^{\prime}}\mathcal{L}_{\mathbb{R}}(y_{i}-g^{\prime}(x_{i}))-\sum_{j=1}^{k}\sum_{i\in C_{j}}\Big{(}\mathcal{L}_{\mathbb{R}}(y_{i}-g(x_{i}))\Big{)}-(|\mathcal{C}^{\prime}|-1)\mathcal{L}(H)_{C^{\prime}}$$

$$\displaystyle\mathcal{L}(\mathcal{C}^{\prime},h)-\mathcal{L}(\mathcal{C}^{*},h)=\sum_{i\in C^{\prime}}\mathcal{L}_{\mathbb{R}}(y_{i}-g^{\prime}(x_{i}))-\sum_{j=1}^{k}\sum_{i\in C_{j}}\Big{(}\mathcal{L}_{\mathbb{R}}(y_{i}-g_{j}(x_{i}))\Big{)}\geq 0$$

$$h^{\ast }=\underset{h\in H}{\mathrm{argmin}}\sum _{i=1}^n(y_i-h(x_i){)}^2$$

$$h(x_i)={\vec{x}}_i^{\text{T}}\vec{\beta }+ϵ.$$

$$\hat{\beta }=\underset{\vec{\beta }\in {\mathbb{R}}^{d+1}}{\mathrm{argmin}}\sum _{i=1}^n(y_i-{\vec{x}}_i^{\text{T}}\vec{\beta }{)}^2$$

$$\hat{\beta }=({\text{X}}_1^{\text{T}}{\text{X}}_1{)}^{-1}{\text{X}}_1^{\text{T}}\vec{y}.$$

$$\mathcal{L}_{\mathbb{R}}(y)=\left\{\begin{array}[]{@{}ll@{}}\lceil\text{log}_{2}|y|\rceil+1,&|y|\geq 1\\ 1,&\text{otherwise}\end{array}\right.$$

$$\mathcal{L}(\mathcal{C}^{\prime},H_{1})-\mathcal{L}(\mathcal{C}^{\prime},H_{2})\propto\sum_{C\in\mathcal{C}^{\prime}}\big{(}\sum_{i\in C}\mathcal{L}_{\mathbb{R}}(y_{i}-h^{C}_{1}(x_{i}))-\sum_{i\in C}\mathcal{L}_{\mathbb{R}}(y_{i}-h^{C}_{2}(x_{i}))+(\mathcal{L}(H_{1})_{C}-\mathcal{L}(H_{2})_{C})\big{)}.$$

$$\mathcal{I}(\mathcal{C}^{\prime},H_{1},H_{2})=\frac{\sum_{C\in\mathcal{C}^{\prime}}\big{(}\sum_{i\in C}\mathcal{L}_{\mathbb{R}}(y_{i}-h^{C}_{1}(x_{i}))-\sum_{i\in C}\mathcal{L}_{\mathbb{R}}(y_{i}-h^{C}_{2}(x_{i}))+(\mathcal{L}(H_{1})_{C}-\mathcal{L}(H_{2})_{C})\big{)}}{\sum_{C\in\mathcal{C}^{\prime}}\big{(}\sum_{i\in C}\mathcal{L}_{\mathbb{R}}(y_{i}-h^{C}_{1}(x_{i}))+\mathcal{L}(H_{1})_{C}\big{)}}.$$

$$\mathcal{I}(\mathcal{C}^{\prime}_{1},\mathcal{C}^{\prime}_{2},H)=\frac{\sum_{C_{1}\in\mathcal{C}^{\prime}_{1}}\big{(}\sum_{i\in C_{1}}\mathcal{L}_{\mathbb{R}}(y_{i}-h^{C_{1}}(x_{i}))\big{)}-\sum_{C_{2}\in\mathcal{C}^{\prime}_{2}}\big{(}\sum_{i\in C_{2}}\mathcal{L}_{\mathbb{R}}(y_{i}-h^{C_{2}}(x_{i}))\big{)}+c_{H}(|\mathcal{C}^{\prime}_{1}|-|\mathcal{C}^{\prime}_{2}|)}{\sum_{C_{1}\in\mathcal{C}^{\prime}_{1}}\big{(}\sum_{i\in C_{1}}\mathcal{L}_{\mathbb{R}}(y_{i}-h^{C_{1}}(x_{i}))\big{)}+c_{H}|\mathcal{C}^{\prime}_{1}|}.$$

$$\eta(C)=\frac{\sum_{C^{\prime}\subset C}\big{(}corr(y_{C^{\prime}},h_{C-C^{\prime}}(x_{C^{\prime}}))\big{)}^{2}}{|\{C^{\prime}:C^{\prime}\subset C\}|}$$

$${\gamma }^{\prime }=\underset{C^{\prime }\in {\mathcal{C}}^{\prime }}{\min }\eta (C^{\prime }).$$

$$Y_j=c_1\cdot X_j+c_2$$

$$Y_j=c_1\cdot e^{X_j}+c_2$$

$$Y_j=c_1\cdot {X_j}^d+c_2$$

$$M_0=q_0\times a_0.$$

$${\beta }^{\ast }(\lambda )={\text{argmin}}_{\beta \in {\mathbb{R}}^l}\frac{1}{n}\sum _{i=1}^{n}loss(y_i,f_{\beta }(x_i))+\lambda {\Vert \beta \Vert }_0$$

$$\mathcal{L}(h,\lambda )=\sum _{i=1}^l{\mathcal{L}}_{\mathbb{R}}(p_i)+\lambda {\mathcal{L}}_{\mathbb{R}}({\Vert \beta \Vert }_0)$$

Full text

multi-resolution data, regression analysis, minimum description length, population data, model selection

††journalyear: 2021††copyright: rightsretained††journal: TKDD††journalvolume: 15††journalnumber: 2††article: 15††publicationmonth: 1††doi: 10.1145/3424670††isbn: 978-1-4503-9999-9/18/06††ccs: Computing methodologies Machine learning approaches††ccs: Information systems Data mining††ccs: Applied computing Sociology

1. Introduction

Ending poverty in all its form everywhere has been recognized as the greatest global challenge in the 2030 Agenda for Sustainable Development (Assembly, 2015). Annually, there are at least 18 millions human deaths that caused by poverty (Pogge, 2005). A common goal is to increase income beyond poverty line (e.g. 1.90 USD per day). However, poverty alleviation often requires comprehensive measures. In practice, implementing poverty alleviation programs would depend on the realities, capabilities and level of development of each nation. In the past, there was a tendency to use a one-size-fit-all policy as a solution in the international level (Lahn and Stevens, 2017). Even though one-size-fit-all policies enable uncomplicated deployment of government resources, in reality, they are not suitable in alleviating poverty as each region possesses different problems and socioeconomic characteristics (Berdegue et al., 2002). The solution in one region might not work for another region even if they share some properties (Lahn and Stevens, 2017; Rioja and Valev, 2004). If one were to drill down to the household-level poverty, one would often find that each household faces different problems in respect to demographic, health, education, living standards, and accessibility to available resources. Subsequently, lifting those household members above the poverty line in sustainable way would often require targeted poverty alleviation with unique solutions. However, in reality, many developing countries often do not possess enough resources. Hence, it is challenging for policy makers to be able to optimize between cheap one-size-fit-all and expensive targeted individual-based approaches.

2. Related work

To combat poverty, one should find root causes of issues and all related factors that contribute to poverty of populations. There are many factors that can contribute to poverty, such as health issues, lacking of education, inaccessibility of public service, severe living condition, lacking of job opportunities, etc  (Strotmann and Volkert, 2018; Alkire et al., 2019). To measure the degree of poverty, several works proposed to use the poverty indices, such as Human Poverty Index (HPI) (Anand and Sen, 1997), Multidimensional Poverty Index (MPI) (Alkire and Santos, 2010; Alkire et al., 2019), etc. However, a specific poverty index covers only some aspects of population, which might not include the actual root cause predictors (Alkire and Santos, 2010; Strotmann and Volkert, 2018). Additionally, different areas typically have different issues and root causes of poverty (Commins, 2004; Pringle et al., 2000).

To infer which factors are the main causes of poverty in the area from population data, the first step is to find the associations between a dependent variable (e.g. poverty index), and independent variables (potential predictors). This is because the causal inference requires association relations among variables to be known before we can find causal directions (Pearl, 2009; Peters et al., 2017). One of the approaches that is widely used in social science to find association relations is Regression Analysis (Nathans et al., 2012). It has been used as a part of framework in poverty indices analysis by the works in (Njuguna and McSharry, 2017; Jean et al., 2016). Linear regression has been extended and used in various types of data beyond linear-static data, such as time series data (Granger, 1969; Udny Yule, 1927; Atukeren et al., 2010), non-linear data (An et al., 2007; Kostopoulos et al., 2018), etc.

In this work, we are interested to find variables that associate with population household incomes, which is one of the main components in several poverty indices (Anand and Sen, 1997; Alkire and Santos, 2010; Alkire et al., 2019).

However, different areas resolutions (e.g. provinces vs. country levels) might have different predictors that are associated with household incomes. In the perspective of policy makers, placing the wrong policies to the areas might not be able to solve the issues. For example, urban and rural areas have different issues of poverty (Commins, 2004). Hence, it is important to find the way to select the model from the right resolutions that truly represents issues of areas. According to the Occam’s razor principle, we prefer a policy that is simple and covers issues of larger population as an efficient policy. One of the model selection approaches that is widely used is Minimum Description Length (MDL) (Hansen and Yu, 2001), which was introduced in the field of data analytics by (Rissanen, 1978). It is used in linear regression setting as a feature selection criteria to find the informative features (Schmidt and Makalic, 2012). Nevertheless, to the best of our knowledge, there is no framework that compares the regression models from different areas resolution using MDL and asks which one is a better model as a predictor for household-income prediction.

2.1. Comparison with the clustering problem

Given a set of multi-resolution partitions as an input of our framework along with the population data (realizations of dependent and independent variables). We propose the framework that is able to infer a subset of multi-resolution partitions (Fig. 1), which is a set of largest partition s.t. each partition has a model that represents population well in the term of predictive performance (see Section 4.3). In contrast, in the clustering problem, the main goal is to find a set of clusters or partitions that is fitting the population data well w.r.t. some similarity function. Typically, the inputs of clustering are the number of clusters and/or some threshold(s) that define clusters, as well as population data.

Hence, in a dataset that has only population data without a set of multi-resolution partitions, we can use some data clustering approach to infer the hierarchical clusters (multi-resolution partitions). Then, our proposed framework can be used to find the optimal partition in the hierarchical clusters that covers the entire population s.t. the partition models represent the population well in the aspect of prediction of the target dependent variable.

Policy makers might consider each cluster in the optimal partition from our framework as a group of individuals that shares the same indicators that can predict the target variable well. If policy makers want to change the value of target variable in the population, then they might establish a policy for each cluster in our optimal partition instead of either making a single one-size-fit-all policy or having a single unique policy for each individual.

2.2. Comparison with mixture models

The concept of mixture models has been widely used for modeling complex distributions by using several mixture components for representing “subdistribution” or “subpopulations” (McLachlan and Peel, 2000).

Generally, an individual subdistribution is modeled using a “simpler” distribution, e.g., Gaussian distribution. Given a specific number of mixture components and a distribution type of mixture components, a finite mixture model will estimate degree of membership of an individual datapoint to all mixture components as well as parameters of the components. Data points belong to the same partition if they have the highest membership degree to the same mixture component. Mixture models can be generalized for regression problems (de Veaux, 1989; Wedel and DeSarbo, 1995). Particularly, a regression function is used, instead of a distribution function, for representing a mixture component. Datapoints belong to the same component would share similar “behavior” according to the regression function and its variables.

In contrary to traditional data clustering techniques, a mixture model can form a data partition without any similarity function defined. The model, however, requires an assumption of subpopulations that are distributed or characterized according to a certain class of distribution or regression function. Ideally, the mixture model framework would naturally fit to our problem if a set of policies is intended to apply to subpopulations, which can be clearly defined by a certain class of distribution or a regression function. Since it is not straightforward to impose structural constraints such as administrative regions (e.g., a pair of administrative regions at the same level cannot belong to the same subpopulation) into the inference algorithm, one can simply let the mixture model to flexibly “search” for subpopulation partitions and subsequently impose the administrative constraints if some partition violates such constraints. Our approach, however, explicitly takes the structural constraints into account by “searching” and evaluating the data partition along the administrative hierarchy. Consequently, the data partition from our approach always complies to the constraints and requires no additional step – the constraint checking – as in the mixture model case.

Table 1 briefly illustrates the difference between Mixture of Regression and our proposed framework. Both approaches assign individuals into a partition. While Mixture of Regression requires the number of clusters (components) as a parameter, our approach uses the multi-resolution partitions as an input. Additionally, instead of only assigning individuals into a partition, our approach can also evaluate how strong the individuals in the same cluster share the same model (see Section 4.3). Moreover, the time complexity of our approach is bounded for the exact solution, while the bound of Mixture approach is unknown; but there is an algorithm (Li and Liang, 2018) that can provide the approximate solution with the bound $\mathcal{O}(Nd)$ where $N$ is a number of individuals and $d$ is a number of dimensions.

2.3. Our contributions

To fill the gaps mentioned above, in this paper, we aim to formalize the multi-resolution model selection problem in regression analysis using MDL principle and propose the framework (Fig. 1) as a solution of our model selection problem. The proposed framework is capable of:

• $\bullet$

Inferring the optimal-resolution partition: inferring the best fit partition that effectively represents population in the aspect of model prediction and MDL; and

• $\bullet$

Inferring the optimal model for each area: inferring the best fit model among candidate functions given an area population data as well as choices of functions in a linear function class.

After formalizing the new problem, we provide the algorithm that guarantees soundness and completeness of finding the optimal partition, which is the output of the new problem. We evaluate our approach performance and demonstrate its application using both simulation and real-world datasets of Thailand population households. Our framework can be generalized beyond the context of poverty analysis.

3. Problem formalization

Given a vector $\vec{y}=(y_{1},\dots,y_{n})$, an $n\times d$ matrix $\textbf{X}=(x_{1},\dots,x_{n})^{\text{T}}$, and a set of Multi-resolution clusters $\mathcal{C}=\{C_{j,k}\}$ (Def. 3.1). $y_{i}\in\mathbb{R}$ is a realization of dependent random variable $Y$ of individual $i$, $\textbf{X}(i,j)\in\mathbb{R}$ is a realization of independent random variable $X_{j}$ of individual $i$, and $C_{j,k}$ is a $j$ cluster in $k$th layer of multi-resolution cluster set where individual $i$ is a member of $C_{j,k}$. In this work, the main purpose is to solve Problem 2 where a hypothesis space $H$ is a linear function class. The details of definitions and the problem formalization can be found the this section.

Specifically, in Section 3.1, we explore the concept of Minimum Description Length (MDL). Then, in Section 3.2, we provide definitions as well as the formalization of Problem 2, which is the main problem we aim to solve in this paper. In Section 3.3, we provide the properties of linear model in the MDL setting. We also provide the Table 2 as a reference for the notations and symbols used in this paper.

3.1. Minimum description length on regression models

Given a set of observation data below:

[TABLE]

Where $(x_{i},y_{i})$ are realizations of random variables $(X_{i},Y_{i})$ s.t. $x_{i}\in\mathbb{R}^{d}$ and $y_{i}\in\mathbb{R}$. The random variables $X_{1},\dots,X_{n}$ are i.i.d. with an unknown distribution $P_{X}$ or $X_{i}\sim P_{X}$. The random variable $Y_{i}$ has a relation with $X_{i}$ below:

[TABLE]

Where $h^{*}(X_{i})$ is an unknown function that is a member of some hypothesis space $\mathbb{H}$ and $\epsilon$ is a noisy constant s.t. $X_{i}$ and $\epsilon$ are statistically independent or X_{i}\rotatebox[origin={c}]{90.0}{\models}\epsilon.

In MDL, we would like to find the number of bits of the shortest representation that we can have for any given data based on some hypothesis space.

[TABLE]

Where $\mathcal{L}(h)$ is a number of bits of representation for a predictive model $h$ and $\mathcal{L}(S|h)$ is a number of bits of representation for data $S$ given the predictive model $h$. In our case, $\mathcal{L}(h)$ is the number of bits we encode parameters and related information of function $h$. The $\mathcal{L}(S|h)$ is the number of bits we encode $S$ using $h$.

[TABLE]

[TABLE]

Where $l$ is a number of parameters of $h$, $p_{j}$ is $j$th parameter of $h$, $\mathcal{L}_{\mathbb{R}}(x)$ is a function that returns a number of bits we need to encode a real number $x$, and $n$ is a number of data points in $S$.

The model complexity in Eq. 3 can be controlled by $\mathcal{L}(h)$. If two models perform equally, then the Eq. 3 chooses the less complex model that has lower $\mathcal{L}(h)$ as a solution. See Section 7.2 for more discussion regarding how model complexity can be implemented and used to compared functions from difference classes.

In the sparsity regularization aspect, we can modify $\mathcal{L}(h)$ to add the higher bits for a model representation (penalty) when the model $h$ has a higher number of parameters (nonzero coefficients). See Section 7.3 for more discussion.

In this paper, we the assumption below:

Assumption 1.

For any $x_{1},x_{2}\in\mathbb{R}$ s.t. $|x_{1}|\geq|x_{2}|$, we have $\mathcal{L}_{\mathbb{R}}(x_{1})\geq\mathcal{L}_{\mathbb{R}}(x_{2})$.

We have Assumption 1 to represent the fact that the larger number contains similar amount of or more information that we need to encode. In the function $\mathcal{L}(S|h)$, there are two parts: the part that we encode $x$ or $\sum_{i=1}^{n}\mathcal{L}_{\mathbb{R}}(x_{i})$ and the part that we encode $y$ by $h\in\mathbb{H}$ or $\sum_{i=1}^{n}\mathcal{L}_{\mathbb{R}}(y_{i}-h(x_{i}))$. If $h$ predicts $y_{i}$ using $x_{i}$ perfectly, then we expect to encode no information of $y_{i}$, which makes $\mathcal{L}_{\mathbb{R}}(y_{i}-h(x_{i}))=0$ bits. However, if $h^{\prime}\in\mathbb{H}$ makes larger error of predicting $y_{i}$ from $x_{i}$ than $h$, we expect to use more bits of encoding for $y_{i}$ using $h^{\prime}$ and $x_{i}$ than $h$. Hence, Assumption 1 naturally represents this intuitive notion of having more bits of encoding for a function $h^{\prime}$ that makes more error of predicting $y_{i}$ from $x_{i}$.

Now, we are ready to formalize MDL-function inference problem.

3.2. Model selection on multi-resolution clusters

Given a set of individual indices $D=\{1,\dots,n\}$ where $i\in D$ refers to $i$th individual in a population, we can define a set of multi-resolution clusters below.

Definition 3.1 (Multi-resolution cluster set (MRC set)).

A set of multi-resolution clusters $\mathcal{C}=\{C_{j,k}\}$ of $n_{c}$ layers, where $k\in[1,n_{c}]$, and $C_{j,k}\subseteq D$ is a $j$th cluster at $k$th layer, has a following properties.

• $\bullet$

If $k>1$, then $\forall C_{j,k},\exists C_{j^{\prime},k-1},C_{j,k}\subseteq C_{j^{\prime},k-1}$.

• $\bullet$

All individuals belong to some cluster in each layer or $\forall k,\bigcup_{j}C_{j,k}=D$.

• $\bullet$

Clusters in the same layer are disjoint or $\forall k,\bigcap_{j}C_{j,k}=\emptyset$.

An example of multi-resolution cluster set is a national administrative division of regions in a country where the first layer is a national level, the second layer is a province level, and the last layer is a village level. Each cluster in each layer is a specific group of individuals that are governed by the same local government in a particular level. For instance, in the province level, each cluster $C$ represents citizens or households within a province $C$. Next, we have a concept of clusters that all members share the same joint distribution $P_{X,Y}$.

Definition 3.2 (Homogeneous cluster).

Given a cluster $C\subseteq D$, and a set of population data $S_{\mathcal{C}}=\{(x_{1},y_{1}),\dots,(x_{n},y_{n})\}$ where $x_{i}\in\mathbb{R}^{d}$ and $y_{i}\in\mathbb{R}$. A cluster $C$ is a homogeneous cluster of $S_{C}=\{(x_{i},y_{i})|i\in C\}$ where $x_{i},y_{i}$ are realizations of $X_{i},Y_{i}$ respectively, if the following conditions hold.

$\exists P_{X},\forall i\in C,X_{i}\sim P_{X}$.

• 2.

$Y_{i}=h^{*}(X_{i})+\epsilon$ for some unknown function $h^{*}$.

Where $\epsilon\in\mathbb{R}$ is a noise random variable from some unknown truncated distribution with the bound $[a,b]$ s.t. $\epsilon^{*}=max(|a|,|b|)$, $\mathbb{E}[\epsilon]=0$, and $|\epsilon|\leq\epsilon^{*}$. All subsets of $C$ are also homogeneous clusters.

The concept of homogeneous cluster represents a specific area (e.g. province, village) that share a common relationship between a dependent variable and independent variables. For instance, a village $A$ has the health issue that affects household incomes of all villagers. Suppose $Y$ is a dependent variable of income and $X$ is an independent variable of degree of health-issue severeness. If a village $A$ is a homogeneous cluster, then $|y_{i}|\propto|x_{i}|$ for any household $i$ in this village. Because all individual households share the similar properties, policy makers should make a single set of policies to support the entire population within a homogeneous cluster.

Next, we need the concept of a cluster set that covers all individuals.

Definition 3.3 (MRC Partition).

Given a multi-resolution cluster set $\mathcal{C}$. A set of clusters $\mathcal{C}_{p}$ is an MRC partition if $\mathcal{C}_{p}\subseteq\mathcal{C}$, all clusters in $\mathcal{C}_{p}$ are disjoint, and $\bigcup_{C_{i}\in\mathcal{C}_{p}}C_{i}=D$.

An MRC Partition represents a set of administrative divisions that covers all individuals within a country. A set of all provinces is an example of the MRC Partition.

However, we need the MRC Partition concept that can include the sub-areas from mixed layers of MRC set but still covers the entire population. This is because some provinces might be homogeneous clusters that share the same properties, while, in other provinces, we need to break down to the village layer to get the homogeneous clusters. Our goal is to use the concept of MRC Partition to define a set of maximal homogeneous clusters (largest possible homogeneous cluster for a specific area) that covers the entire population.

Definition 3.4 (Maximal homogeneous cluster).

Given a multi-resolution cluster set $\mathcal{C}$, and a set of population data $S_{\mathcal{C}}=\{(x_{1},y_{1}),\dots,(x_{n},y_{n})\}$ where $x_{i}\in\mathbb{R}^{d}$ and $y_{i}\in\mathbb{R}$. A cluster $C^{*}\in\mathcal{C}$ is a maximal homogeneous cluster if $C$ is a homogeneous cluster, and there is no other homogeneous cluster $C^{\prime}$ s.t. $C^{\prime}\in\mathcal{C}$ and $C^{*}\subset C^{\prime}$.

The numbers of elements within MRC Partition might represent a number of sets of policies that policy makers should consider to make in order to support each area common issues. Next, we need to define a concept to call datasets that have the maximal homogeneous partition, which is the MRC Partition that contains only maximal homogeneous clusters as members.

Definition 3.5 (Multi-resolution-bivariate set (MRB)).

Let $\mathcal{C}$ be a multi-resolution cluster set, and $S_{\mathcal{C}}$ be a set of $\mathcal{C}$’s population data. $S_{\mathcal{C}}$ is a Multi-resolution-bivariate set if $S_{\mathcal{C}}$ satisfies the following properties.

• $\bullet$

There exists an MRC partition of homogeneous clusters, $\mathcal{C}^{*}=\{C\}$, s.t. for each $C\in\mathcal{C}^{*}$, $C$ is a maximal homogeneous cluster. We call $\mathcal{C}^{*}$ as the “Maximal Homogeneous Partition”.

• $\bullet$

Each homogeneous $C_{i}\in\mathcal{C}^{*}$ has a separate function $h^{*}_{i}$ that might be different from any other function $h^{*}_{j}$ of $C_{j}\in\mathcal{C}^{*}$ where $C_{i}\neq C_{j}$.

Suppose we have an MRB set $S_{\mathcal{C}}=\{(x_{1},y_{1}),\dots,(x_{n},y_{n})\}$, and a function class $H$ (e.g. linear class, polynomial class, convex class). We can use the data in each cluster $C_{j,k}$ to estimate the parameters of function in $H$ to build a set of models $H_{j,k}\subset H$. Each model $h\in H_{j,k}$ has parameter values estimated from the data points $\{(x_{i},y_{i})\}$ in the cluster $C_{j,k}$. Let $\mathcal{H}=\bigcup H_{j,k}$ be a set of models of $\mathcal{C}$, we can formalize our new computational problem: MDL-function inference on multi-resolution data problem, in Problem 2.

3.3. Linear function properties

Suppose we have $H$ as a linear function class, and $\mathcal{C}^{\prime}$ as an MRC partition of $D$. We can have a measure of how well the models derived from $\mathcal{C}^{\prime}$ represent the set of population data $S$ in the aspect of MDL below.

[TABLE]

[TABLE]

[TABLE]

Where $\hat{\beta}_{1}^{C}\in\mathbb{R},\hat{\beta}^{C}\in\mathbb{R}^{d}$ are the estimators of linear function parameters based on the data points in $C$. We have the next proposition to show that the sum of encoding bits of residuals from homogeneous clusters are less than the sum of encoding bits of residuals of function that its parameters are approximated by the union of these homogeneous clusters.

Proposition 3.6.

Let $C=\{C_{1},\dots,C_{k}\}$ be a set of homogeneous clusters of $S_{\mathcal{C}}$ but the cluster $C^{\prime}=\bigcup_{C_{i}\in C}C_{i}$ is not a homogeneous cluster, and $H$ is a linear function class. Assuming that the functions that generate $Y$ from $X$ in $S_{\mathcal{C}}$ for each cluster $C_{i}\in C$ are in $H$.

If $g^{\prime}$ is a function that its parameters are approximated by $C^{\prime}$ and $g_{1}\dots,g_{k}$ are functions that their parameters are approximated by $C_{1},\dots,C_{k}$ respectively, where $g^{\prime}$ and $g_{1}\dots,g_{k}$ are in $H$. We have $\sum_{i\in C^{\prime}}\mathcal{L}_{\mathbb{R}}(y_{i}-g^{\prime}(x_{i}))\geq\sum_{j=1}^{k}\sum_{i\in C_{j}}\mathcal{L}_{\mathbb{R}}(y_{i}-g_{j}(x_{i}))$.

Proof.

Because each cluster in $C$ is homogeneous and the relation between $Y$ and $X$ is $Y_{i}\approx\beta_{1}+X_{i}^{\text{T}}\beta+\epsilon$. We can approximate $\beta_{1},\beta$ using some estimator approach. Suppose $\hat{\beta}_{1},\hat{\beta}$ are estimators for $\beta_{1},\beta$, we can have

[TABLE]

In contrast, if two clusters $C_{1},C_{2}\in C$ have a different functions $g_{1},g_{2}$, for example, $g_{1}(X)=-X^{\text{T}}\beta_{c}\ +\epsilon$ while $g_{2}(X)=-g_{1}(X)$, then by approximating $\beta_{1},\beta$ for $g^{\prime}(X)$ using data points from both $C_{1},C_{2}$, the sum of residuals of data points from both clusters must be higher than $\epsilon^{*}$. Precisely, according to Def. 3.2, $|\epsilon|\leq\epsilon^{*}$ for any homogeneous cluster. Suppose $g_{1}\neq g_{2}$ and $C_{1},C_{2}$ are homogeneous clusters of $g_{1}$ and $g_{2}$ respectively. If we approximate the parameters from data points of both $C_{1},C_{2}$, then we would have $g^{\prime}$ that is not the same as both $g_{1},g_{2}$ but $g^{\prime}$ is something between $g_{1},g_{2}$. For $g^{\prime}$, the sum of encoding bits of residuals of $C_{1}$ is $\sum_{i\in C_{1}}\mathcal{L}_{\mathbb{R}}(y_{i}-g^{\prime}(x_{i}))\leq|C_{1}|\mathcal{L}_{\mathbb{R}}(\epsilon^{*})+\sum_{i\in C_{1}}\mathcal{L}_{\mathbb{R}}(g^{\prime}(x_{i})-g_{1}(x_{i}))$ where the term $\sum_{i\in C_{1}}\mathcal{L}_{\mathbb{R}}(g^{\prime}(x_{i})-g_{1}(x_{i}))>0$ since $g^{\prime}\neq g_{1}$. The same applies for $C_{2}$.

Hence, the sum of errors using $g^{\prime}$ to predict $y$ must be greater than $|C^{\prime}|\mathcal{L}_{\mathbb{R}}(\epsilon^{*})$.

Therefore, we have

[TABLE]

Hence, $\sum_{i\in C^{\prime}}\mathcal{L}_{\mathbb{R}}(y_{i}-g^{\prime}(x_{i}))\geq\sum_{j=1}^{k}\sum_{i\in C_{j}}\mathcal{L}_{\mathbb{R}}(y_{i}-g_{j}(x_{i}))$. ∎

We have the next theorem to state that the maximal homogeneous partition has the lowest encoding bits compared to any other MRC partitions.

Theorem 3.7.

Let $D$ be a set of individual indices, $\mathcal{C}$ be a multi-resolution cluster set , $S_{\mathcal{C}}$ be an MRB set, $H$ be a linear function class, and $\mathcal{C}^{*}$ be the maximal homogeneous partition of $S_{\mathcal{C}}$. For any $\mathcal{C}^{\prime}$ that is an MRC partition, $\mathcal{L}(\mathcal{C}^{*},H)\leq\mathcal{L}(\mathcal{C}^{\prime},H)$.

Proof.

Suppose $\mathcal{L}(\mathcal{C}^{*},H)\neq\mathcal{L}(\mathcal{C}^{\prime},H)$, because both $\mathcal{C}^{*}$ and $\mathcal{C}^{\prime}$ are MRC partitions, all cluster members in both partitions are also members of MRC set $\mathcal{C}$. It is obvious that $\mathcal{L}(\mathcal{C}^{*},H)\neq\mathcal{L}(\mathcal{C}^{\prime},H)$ because some clusters between $\mathcal{C}^{\prime}$ and $\mathcal{C}^{*}$ are different.

Let $A=\mathcal{C}^{*}\cap\mathcal{C}^{\prime}\subset\mathcal{C}$. There are some regions $\mathcal{C}^{*}-A$ that make both partitions different. There are two cases that we need to consider.

Case 1: suppose $C\in\mathcal{C}^{*}-A$ and $C^{\prime}=\{C^{\prime}_{1},\dots,C^{\prime}_{k}\}\subseteq\mathcal{C}^{\prime}-A$ where $\bigcup_{C^{\prime}_{i}\in C^{\prime}}C^{\prime}_{i}=C$. Since $C$ is a homogeneous cluster, the subsets of $C$ which are $\{C^{\prime}_{1},\dots,C^{\prime}_{k}\}$ must be the homogeneous clusters from the same $X\sim P_{X}$ and $Y=h(X)+\epsilon$. Because $\bigcup_{C_{i}\in C^{\prime}}C_{i}=C$ implies that we estimate the parameters of $h$ from the same region, hence, $\sum_{C^{\prime}_{i}\in\mathcal{C}^{\prime}}\mathcal{L}(S_{C^{\prime}_{i}}|H)\approx\mathcal{L}(S_{C}|H)$.

[TABLE]

Since $H$ is a linear class and all homogeneous clusters in this case share the same $h(X)$, then $\mathcal{L}(H)_{C}\approx\mathcal{L}(H)_{C^{\prime}_{i}}$.

[TABLE]

Case 2: suppose $C=\{C_{1},\dots,C_{k}\}\subseteq\mathcal{C}^{*}-A$ and $C^{\prime}\in\mathcal{C}^{\prime}-A$ where $\bigcup_{C_{i}\in C}C_{i}=C^{\prime}$. Since $\mathcal{C}^{*}$ is a maximal partition, $C^{\prime}$ cannot be a homogeneous cluster. In fact, different homogeneous cluster might have different function $h^{*}$ where $Y_{i}=h^{*}(X_{i})+\epsilon$. Let $g^{\prime}$ be a model that parameters are approximated by $C^{\prime}$ and $g_{1}\dots,g_{k}$ be models that parameters are approximated by $C_{1},\dots,C_{k}$ respectively. Again, we assume that $\mathcal{L}(H)_{C^{\prime}}\approx\mathcal{L}(H)_{C_{i}}$.

[TABLE]

Typically, the number of parameters of linear functions is less than the number of individuals. Hence, it is safe to assume that $(|\mathcal{C}^{\prime}|-1)\mathcal{L}(H)_{C^{\prime}}\ll\sum_{i\in C^{\prime}}\mathcal{L}_{\mathbb{R}}(y_{i}-g^{\prime}(x_{i}))$. According to Proposition 3.6, we have $\sum_{i\in C^{\prime}}\mathcal{L}_{\mathbb{R}}(y_{i}-g^{\prime}(x_{i}))\geq\sum_{j=1}^{k}\sum_{i\in C_{j}}\mathcal{L}_{\mathbb{R}}(y_{i}-g_{j}(x_{i}))$. Hence,

[TABLE]

In both cases, $\mathcal{L}(\mathcal{C}^{\prime},h)-\mathcal{L}(\mathcal{C}^{*},h)\geq 0$. Therefore, $\mathcal{L}(\mathcal{C}^{*},h)\leq\mathcal{L}(\mathcal{C}^{\prime},h)$.

∎

4. Methods

In this section, we provide the details of method that is used to solve Problem 2. Fig. 1 illustrates the overview of our proposed framework. In the first step, the framework performs the parameter estimation (Section 4.1). Then, we use MDL (Section 4.2) and cross-validation scheme (Section 4.3) to infer the maximal homogeneous partition in Section 4.4. We provide the time complexity of our method in Section 4.5.

4.1. Parameter estimation

In linear models, to minimize the magnitude of residuals of prediction, the optimal $h^{*}$, which minimizes the residuals, can be inferred using the least squares technique. Given $H$ is a linear function class, we have an optimization problem using least square loss below.

[TABLE]

The linear function is in the form:

[TABLE]

Where $\vec{x}_{i}$ is a vector of $x_{i}$ s.t. $\vec{x}_{i}(1)=1$ and $\vec{x}_{i}(j)=x_{i}(j-1)$, and $\epsilon\sim\mathcal{N}(0,\tau^{*})$ is a noisy value where $\tau^{*}$ is an unknown variance. Hence,

[TABLE]

We call $\hat{\beta}$ an ordinary-least-squares (OLS) estimator of $\vec{\beta}$. Suppose all columns in X are linear independent, the optimal solution for Eq. 11 can be solved by the closed-form equation as follows:

[TABLE]

Where $\textbf{X}_{1}$ is an $n\times d+1$ matrix s.t. the 1st column of $\textbf{X}_{1}$ contains only ones, and for any other column $j>1$, $\textbf{X}_{1}(i,j)=\textbf{X}(i,j-1)$.

In this step, for each cluster $C\in\mathcal{C}$, we estimate $\hat{\beta}_{C}$ using data points in cluster $C$ ($(x_{i},y_{i})$ s.t. $i\in C$ ).

4.2. Minimum-description-length measure of models

To estimate $\mathcal{L}(\mathcal{C}^{\prime},H)$ in Eq. 6, we need to compute $\mathcal{L}(H)_{C}$ in Eq. 7 and $\mathcal{L}(C|H)$ in Eq. 8 that both functions require the estimation of $\mathcal{L}_{\mathbb{R}}(x)$.

Suppose we have two integer numbers $y_{1}$ and $y_{2}$ where $|y_{1}|\ll|y_{2}|$. When we compress $y_{1},y_{2}$ into a computer memory space, because the magnitude of $y_{1}$ is a lot smaller than $y_{2}$, if we need $q$ bits to represent $y_{2}$, then we can approximately represent $y_{1}$ without using all $q$ bits.

For integer numbers, we need at least $\mathcal{L}_{\mathbb{I}}(y)=\lceil\text{log}_{2}|y|\rceil+1$ bits to represent $y$ in a binary representation where $+1$ bit is for the sign representation of $y$. For a real number representation, the work by  (Lindstrom et al., 2018) provided the framework to implement the universal coding for real numbers s.t. for any $y\in\mathbb{R}$, $|\mathcal{E}(y)|\propto\lceil\text{log}_{2}|y|\rceil$ where $|\mathcal{E}(y)|$ is a number of bits of real-number representation in (Lindstrom et al., 2018)’s work. Hence, for simplicity, for any real number $y$, we can approximately calculate $\mathcal{L}_{\mathbb{R}}(y)$ below.

[TABLE]

If Y is a vector or matrix, then $\mathcal{L}_{\mathbb{R}}(\textbf{Y})=\sum_{i}\mathcal{L}_{\mathbb{R}}(y_{i})$ where $y_{i}$ is a real number element in Y.

4.2.1. Comparing two linear models in the same data points

Given $H$ is a linear function class and $H_{1},H_{2}\subset H$ where $H_{1}\neq H_{2}$. To compare $H_{1}$ and $H_{2}$ using $\mathcal{L}(\mathcal{C}^{\prime},H)$ in Eq. 6 that their parameters are estimated from the same set of data points and the same set of clusters $\mathcal{C}^{\prime}$, as well as $\sum_{i\in C}\mathcal{L}_{\mathbb{R}}(x_{i})$ in both $\mathcal{L}(C|H_{1})$ and $\mathcal{L}(C|H_{2})$ are the same, we have a following relation.

[TABLE]

Where $h^{C}_{1}(x_{i})\in H_{1}$ is a linear function that its $\hat{\beta}$ is estimated using data points in $C$, the same applies for $h^{C}_{2}(x_{i})\in H_{2}$. The normalization of the equation above is as follows:

[TABLE]

Where $\mathcal{I}(\mathcal{C}^{\prime},H_{1},H_{2})\in[-\infty,1]$, and $\mathcal{L}(H_{1})_{C},\mathcal{L}(H_{2})_{C}$ are the numbers of bits needed for the compression of parameters of a function that is approximated by data points in $C$ using $H_{1}$ and $H_{2}$ respectively. We call $\mathcal{I}(\mathcal{C}^{\prime},H_{1},H_{2})$ “Model Information Reduction Ratio”. The function $\mathcal{I}(\mathcal{C}^{\prime},H_{1},H_{2})$ represents how well we can reduce the size of our data using $H_{2}$ instead of $H_{1}$. For example, if $\mathcal{I}(\mathcal{C}^{\prime},H_{1},H_{2})=0.5$, it implies that $H_{2}$ reduces the size of space we need to encode the data at $50\%$ compared to the size we encode the same data using $H_{1}$. In contrast, if $\mathcal{I}(\mathcal{C}^{\prime},H_{1},H_{2})=-0.4$, it means using $H_{1}$ is still a better option since $H_{2}$ increases the size of space we need to encode the same data around $40\%$.

4.2.2. Comparing two sets of clusters in the same data points

Given $H$ is a linear function class, and $\mathcal{C}^{\prime}_{1},\mathcal{C}^{\prime}_{2}$ are two sets of clusters where they cover the same data points or $\bigcup_{C\in\mathcal{C}^{\prime}_{1}}=\bigcup_{C\in\mathcal{C}^{\prime}_{2}}$. Assuming that $\mathcal{L}(H)_{C}$ for each cluster have the same size $c_{H}$, we can have “Cluster Information Reduction Ratio ” $\mathcal{I}(\mathcal{C}^{\prime}_{1},\mathcal{C}^{\prime}_{2},H)$ below:

[TABLE]

Where $h^{C_{1}},h^{C_{2}}$ are functions that their parameters are approximated using data points in cluster $C_{1}$ and $C_{2}$ respectively. The function $\mathcal{I}(\mathcal{C}^{\prime}_{1},\mathcal{C}^{\prime}_{2},H)\in[-\infty,1]$. $\mathcal{I}(\mathcal{C}^{\prime}_{1},\mathcal{C}^{\prime}_{2},H)$ represents how well the different sets of clusters can encode the same information. For example, if $\mathcal{I}(\mathcal{C}^{\prime}_{1},\mathcal{C}^{\prime}_{2},H)=0.5$, it implies that $\mathcal{C}_{2}$ reduces the size of space we need to encode the data at $50\%$ compared to the size we can encode the same data using $\mathcal{C}_{1}$. In contrast, if $\mathcal{I}(\mathcal{C}^{\prime}_{1},\mathcal{C}^{\prime}_{2},H)=-0.4$, it means using $\mathcal{C}_{1}$ is still a better option since $\mathcal{C}_{2}$ increases the size of space we need to encode for the same data around $40\%$.

4.3. Homogeneity measure of clusters

In order to infer the degree of homogeneity of clusters, we deploy the cross-validation concept. It is typically used as a methodology to measure the generalization of model. If the model is generalized well, then it is capable of performing well in any given dataset outside of a training data that is used to approximate the model parameters (Abu-Mostafa et al., 2012). Given a cluster $C$ and its subsets $\{C_{1},\dots,C_{k}\}$, we use the squared correlation between predicted and real $y$ of cross validation among $C$ subsets as a measure of homogeneity of cluster $C$.

[TABLE]

Where $h_{C-C^{\prime}}$ is a function that its parameters are approximated using data points in $C-C^{\prime}$, and $corr(y_{C^{\prime}},h_{C-C^{\prime}}(x_{C^{\prime}}))$ is a correlation between true values of $y$ in $C^{\prime}$ and the predicted values of $y$ using $x_{C}^{\prime}$ in $C^{\prime}$ by $h_{C-C^{\prime}}$.

We estimate $h_{C-C^{\prime}}$’s parameters using data points from the rest of subset clusters in $C$ except $C^{\prime}$, then using $h_{C-C^{\prime}}$ to predict $y$ using data points of $x$ in $C^{\prime}$. If $C$ has no subsets in a multi-resolution cluster set $\mathcal{C}$, then we can use $10$-fold cross validation scheme to generate subsets of $C$ by uniformly separating members of $C$ into 10 subsets that each has equal size.

Note that the Eq. 10 in (Browne, 2000), which is used as a cross validation evaluation, is similar to \big{(}corr(y_{C^{\prime}},h_{C-C^{\prime}}(x_{C^{\prime}}))\big{)}^{2} in Eq. 16. See (Browne, 2000) for more details regarding other indices that can be used to evaluate a linear model in a cross validation scheme.

Next, we explore a theoretical property between $\eta(C)$ and homogeneity of clusters.

Definition 4.1 ($\gamma$-correlation-linear-learnable dependency relation).

Given a function a set of data points $S=\{(x_{1},y_{1}),\dots,(x_{n},y_{n})\}$ where $x_{i}\in\mathbb{R}^{l}$ is a realization of $X\sim P_{X}$, $y_{i}$ is a realization of $Y$ where $Y=h(X)$ s.t. $h(X)$ is a linear function, and a threshold $\gamma\in[0,1]$. Let $S_{1}$ and $S_{2}$ be disjoint subsets of $S$. Assuming that all dimensions in $x_{i}$ are linearly independent.

We say that $S_{1}$ and $S_{2}$ are $\gamma$-correlation-linear-learnable dependent if the following conditions are satisfied

1. There exists the ordinary least square estimator $\vec{\beta}$ of the linear function $g(x)=\vec{x}^{\text{T}}\vec{\beta}$ s.t. $\vec{\beta}$ is estimated using all data points from $S_{1}$ and $corr(y_{S_{2}},g(x_{S_{2}}))\geq\gamma$ where $corr(y_{S_{2}},g(x_{S_{2}}))$ is a correlation between the true values of $y_{i}$ in $S_{2}$ and predicted values of $y_{i}$ using $g(x)$ and $x_{i}$ in $S_{2}$.

2. The same must be true if we use $S_{2}$ for training $g^{\prime}$ function. We have $corr(y_{S_{1}},g^{\prime}(x_{S_{1}}))\geq\gamma$.

We denote $S_{1}\approxeq_{\gamma}S_{2}$.

Lemma 4.2.

A relation in Def. 4.1 is a dependency relation (Aalbersberg and Rozenberg, 1988), which possesses the reflexive and symmetric properties.

Proof.

First, we show that $S_{1}\approxeq_{\gamma}S_{1}$. The Theorem 4 in (Turney, 1994) states that if all dimensions in $x$ are linearly independent, then, using OLS and data points in $S_{1}$ to train a function $g$, $y_{i}=g(x_{i})$ for all $(x_{i},y_{i})$ in $S_{1}$. This implies $corr(y_{S_{1}},g(x_{S_{1}}))=1$. Hence, the relation in Def. 4.1 is reflexive. For the symmetric property, for any $S_{1}\approxeq_{\gamma}S_{2}$, it implies $S_{2}\approxeq_{\gamma}S_{1}$ by definition. ∎

By using $\eta(C)$ in Eq. 16, we can define the version of homogeneous cluster with the homogeneous degree property.

Definition 4.3 ($\gamma$-Homogeneous cluster).

Given a set of individual indices $D=\{1,\dots,n\}$, and a set of population data $S=\{(x_{1},y_{1}),\dots,(x_{n},y_{n})\}$ where $x_{i}\in\mathbb{R}^{d}$ and $y_{i}\in\mathbb{R}$. A cluster $C\subseteq D$ is a homogeneous cluster of $S$ if $\eta(C)\geq\gamma$.

Proposition 4.4.

Given $C$ as a homogeneous cluster and a set of disjoint subsets $\mathcal{C}=\{C_{1},\dots,C_{m}\}$ where $C_{i}\subset C$ and $\bigcup_{C_{i}\in\mathcal{C}}C_{i}=C$. For any $C_{i}\in\mathcal{C}$, suppose $S_{C_{i}},S_{C-C_{i}}$ are sets of data points in $C_{i}$ and $C-C_{i}$ respectively. If $S_{C_{i}}\approxeq_{\gamma}S^{\prime}_{C-C_{i}}$ for any $C_{i}\in\mathcal{C}$, then $\eta(C)\geq\gamma$, which implies $C$ is a $\gamma$-homogeneous cluster.

Proof.

Since $C$ is a homogeneous cluster, all data points $x\sim X$ in $C$ are generated from the same distribution $X\sim P_{X}$ and the relation $Y=h(X)$ is the same for $Y$ and $X$ random variables in $C$. For any $C_{i}\in\mathcal{C}$, because $S_{C_{i}}\approxeq_{\gamma}S^{\prime}_{C-C_{i}}$, it implies that $corr(y_{C_{i}},g_{C-C_{i}}(x_{C_{i}}))\geq\gamma$, where $g_{C-C_{i}}$ is a function trained by data points in $C-C_{i}$, and $y_{C_{i}},x_{C_{i}}$ are data points in $C_{i}$. Therefore, by averaging all $corr(y_{C_{i}},g_{C-C_{i}}(x_{C_{i}}))$ from all $C_{i}$, $\eta(C)\geq\gamma$. ∎

The Proposition 4.4 suggests that if we know that all pairs of subsets $C_{i}$ and $C-C_{i}$ of $C$ have a $\gamma$-correlation-linear-learnable dependency relation, we can imply that $C$ is a $\gamma$-homogeneous cluster. The dependency relation in Def. 4.1 can be checked by running the real data.

4.4. Inferring maximal homogeneous partition

In this section, we propose the algorithm to solve Problem 2. To determine whether linear models provide informative solution, we compare a linear model of each cluster $C$ with the null model, which is the model that we predict $y$ using the average of $y$, $\bar{y}=\frac{\sum_{i\in C}y_{i}}{|C|}$. We use Eq. 14 (Model Information Reduction Ratio $\mathcal{I}(\mathcal{C}^{\prime},H_{1},H_{2})$) to compare linear and null models. We set $\mathcal{C}^{\prime}=\{C\}$, $H_{1}=\{h(x)=\bar{y}\}$, $H_{2}=\{h(x)=\vec{x}_{i}^{\text{T}}\vec{\beta}_{C}\}$ where $\vec{\beta}_{C}$ is estimated using data points in $C$. Hence, $\mathcal{L}(H_{1})_{C}=\mathcal{L}_{\mathbb{R}}(\bar{y})$ and $\mathcal{L}(H_{2})_{C}=\mathcal{L}_{\mathbb{R}}(\vec{\beta}_{C})$.

Next, we show that our proposed algorithm (Algorithm 3) provides the solution for Problem 2.

Proposition 4.5.

Given a multi-resolution cluster set $\mathcal{C}=\{C_{j,k}\}$, MRB set $\mathcal{S}_{\mathcal{C}}$, and the threshold $\gamma=0$. Algorithm 3 always returns the maximal homogeneous partition.

Proof.

For soundness, given $\mathcal{S}_{\mathcal{C}}$ and $\mathcal{C}=\{C_{j,k}\}$, we prove that Algorithm 3 always provides a maximal homogeneous partition as an output. According to Theorem 3.7, the maximal homogeneous partition $\mathcal{C}^{*}$ always has the $\mathcal{L}(\mathcal{C}^{*},H)$ lower or equal any MRC partition. First, we show that Algorithm 3 provides MRC partition.

For line 4 to 5, the algorithm seeks the clusters from a top layer to a bottom one. It implies that if there exists a cluster $C\in\mathcal{C}^{*}$ in the above layer, then it is included to $\mathcal{C}^{*}$ before its subsets. The condition in line 7 prevents the algorithm to add any subsets of the cluster members of $\mathcal{C}^{*}$, hence, all clusters in $\mathcal{C}^{*}$ are disjoint. For some cluster $C$ that is not the last layer member of $\mathcal{C}$, it is either included to $\mathcal{C}^{*}$ by the line 8-11 or its sub clusters are included in $\mathcal{C}^{*}$ for some later iteration of the loop. If there is no sub clusters of $C$ are included until the last layer, then all sub clusters of $C$ at the last layer are included into $\mathcal{C}^{*}$ by default. Hence, $\mathcal{C}^{*}$ covers all individuals because the union of the first layer clusters must cover all individuals and later layer clusters are subsets of some first layer cluster. Hence, the $\mathcal{C}^{*}$ is MRC partition.

Second, we show that $\mathcal{C}^{*}$ is the maximal homogeneous partition. In line 11, we include the cluster $C$ into $\mathcal{C}^{*}$ only if $\mathcal{I}(\{C_{1},\dots,C_{k}\},\{C\},H_{\text{lin}})>0$. For the homogeneous cluster $C$ and its subsets $C_{1},\dots,C_{k}$, by Case 1 in Theorem 3.7, $\mathcal{I}(\{C_{1},\dots,C_{k}\},\{C\},H_{\text{lin}})>0$. Since our algorithm performs top-down searching, it always gets a homogeneous cluster $C$ in the highest possible layers before its subsets by the line 11. Hence, Algorithm 3 provides the maximal homogeneous partition.

For completeness, given $\mathcal{C}^{*}$ as any maximal homogeneous partition, we prove that there is only one possible unique $\mathcal{C}^{*}$ that Algorithm 3 provides and no other $\mathcal{C}^{\prime}$ exists. Suppose $\mathcal{C}^{*}_{1},\mathcal{C}^{*}_{2}$ are maximal homogeneous partitions of $\mathcal{S}_{\mathcal{C}}$ and $\mathcal{C}=\{C_{j,k}\}$. Let $A=\mathcal{C}^{*}_{1}\cap\mathcal{C}^{*}_{2}$, there are some clusters in $\mathcal{C}^{*}_{1}-A$ that are different from clusters in $\mathcal{C}^{*}_{2}-A$. Let assume that $C\in\mathcal{C}^{*}_{1}-A$ and $\{C_{1},\dots,C_{k}\}\subseteq\mathcal{C}^{*}_{2}-A$ where $\bigcup C_{i}=C$. According to Case 1 in Theorem 3.7, because $C$ and its subsets are homogeneous, the length of encoding by $C$ must smaller than using $\{C_{1},\dots,C_{k}\}$. Hence, $\mathcal{C}^{*}_{2}$ is not maximal homogeneous partition, which is a contradiction!

Therefore, Algorithm 3 always returns the maximal homogeneous partition which is unique. ∎

The $\eta(C_{j,k})_{\text{cv}}\geq\gamma$ condition in line guarantees that all homogeneous clusters are $\gamma$-homogeneous clusters. By setting $\gamma$, Algorithm 3 provides the maximal homogeneous partition s.t. all homogeneous clusters are $\gamma$-homogeneous clusters. Otherwise, it provides the set of MRC partition that contains the $\gamma$-homogeneous clusters from the highest possible layer that can be found, and the clusters from the last layers for the population that has no $\gamma$-homogeneous cluster in any layer.

Let $\mathcal{C}^{\prime}$ be the homogeneous partition with some none-$\gamma$-homogeneous clusters generated by the algorithm at the first time. We can have

[TABLE]

After we set the threshold $\gamma=\gamma^{\prime}$ and run Algorithm 3 for the second time, then the result of the algorithm is the maximal homogeneous partition with all $\gamma^{\prime}$-homogeneous clusters. This is true since we know that all clusters in $\mathcal{C}^{\prime}$ are $\gamma^{\prime}$-homogeneous clusters. By running the algorithm again at the second time, if the result $\mathcal{C}^{\prime}_{2}$ is not the same as the first running result $\mathcal{C}^{\prime}$, by the restriction as line 11, the different parts must be $\gamma^{\prime}$-homogeneous clusters. However, the experts in the field should determine whether $\gamma^{\prime}$ is appropriate for their problems.

4.5. Time complexity

The least square approach has the time complexity as $\mathcal{O}(n^{2}d)$ where $n$ is a number of individuals and $d$ is a number of $X$ dimensions. Given $n_{max}$ is a number of individuals in the largest cluster and $|\mathcal{C}|$ is a total number of clusters from all layers. Algorithm 3 has a time complexity as $\mathcal{O}(|\mathcal{C}|n_{max}^{2}d)$ where $n_{max}\leq n$. The lower bound is $\Omega(n^{2}d)$.

5. Experimental setup

We use both simulation and real-world datasets to evaluate our method performance.

5.1. Simulation data

Given a multi-resolution cluster set $\mathcal{C}=\{C_{j,k}\}$, we assign some subsets of $\mathcal{C}$ to be a maximal homogeneous partition $\mathcal{C}^{*}=\{C\}$. We generate realizations of independent random variables $X_{1},\dots,X_{K}$ from normal distribution where $X_{k}\sim\mathcal{N}(0,1)$. For all members in each cluster $C_{j}\in\mathcal{C}^{*}$, their random variables $Y_{j},X_{j}$ have the same joint distribution. We have a following relation between $X_{j}$ and $Y_{j}$.

[TABLE]

Where $c_{1},c_{2}\in\mathbb{R}$, and \forall k\neq j,Y_{j}\rotatebox[origin={c}]{90.0}{\models}X_{k}. This guarantees that each $C_{j}\in\mathcal{C}^{*}$ has a different function $Y=h(X)$.

In our simulation setting, we have $X_{1},\dots X_{20}$. Each individual $i$ has a pair of values $(y_{i},\vec{x}_{i})$ where $\vec{x}_{i}=(x_{1,i},\dots,x_{20,i})$. We have four types of simulation datasets for linear models. Each type has a different $\mathcal{C}=\{C_{j,k}\}$ and $\mathcal{C}^{*}=\{C\}$ (see Fig. 2). Each cluster in the last layer of any dataset have 10,000 individuals as members. Hence, a type-1/2/3 dataset has $\vec{y}=(y_{1},\dots,y_{80000})$, and $80000\times 20$ matrix X, while a type-4 dataset has $\vec{y}=(y_{1},\dots,y_{240000})$, and $240000\times 20$ matrix X.

Moreover, we also generated two types of datasets that has $Y=h(X)$ as a nonlinear function. The first one is the dataset type of exponential function.

[TABLE]

The second type is the polynomial function.

[TABLE]

Where $c_{1},c_{2}\in\mathbb{R}$, and \forall k\neq j,Y_{j}\rotatebox[origin={c}]{90.0}{\models}X_{k}. We set $d=3$ as the polynomial degree. We generated these nonlinear datasets base on a multi-resolution cluster set of type-4 dataset in Fig. 2. The parameter settings of both types are the same as the type-4 datasets except there are 100 individuals per clusters for the last layer.

We use these simulation datasets to evaluate whether our framework can infer the correct $\mathcal{C}^{*}$. We generated 100 datasets for each type and used them to report the averages of framework performance for each dataset type. We define true positive cases (TP) of prediction as a number of individuals of clusters within the ground-truth maximal homogeneous partition s.t. these clusters are also within the predicted maximal homogeneous partition. The false negative cases (FN) of prediction is a number of individuals in clusters within the ground-truth maximal homogeneous partition that are not in the predicted maximal homogeneous partition. The false positive cases (FP) of prediction is a number of individuals that belong to clusters outside the maximal homogeneous partition but the predicted results claim that these clusters are in the maximal homogeneous partition. We use TP, FN, and FP to compute precision, recall, and F1-score values in the result section.

Lastly, the distributions of all types of simulation data are shown in Figure 3.

5.1.1. Motivation of the simulation design

The Multidimensional Poverty Index (MPI) (Alkire and Santos, 2010; Alkire et al., 2019) is used by United Nations Development Programme (UNDP) to measure acute poverty around the world. By design, MPI considers several factors that can contribute to poverty of people. In Thailand, there are five main dimensions that government policy makers consider in MPI: health, living condition, education, financial status, and access to the public services (see Table 3 for more details). Note that MPI typically analyzes in the household level.

Given $Q$ as a matrix derived from surveys where $Q(i,j)$ is a binary value of whether household $i$ fails the dimension $j$ in MPI: one for a fail status and zero for a pass status. In other words, a row $i$ of $Q$ represents a vector of poverty statuses of household $i$. For example, suppose $j$ is the dimension of household average income, if household $i$ has the average income below the government threshold, then $Q(i,j)=1$. In the analysis, suppose $i$ fails $d_{i}$ dimensions. $i$ is considered to be poor if $\bar{d}_{i}=d_{i}/n_{0}$ is below a government threshold where $n_{0}$ is a number of dimensions. Now, we are ready to define the MPI index.

[TABLE]

Where $q_{0}\in[0,1]$ is a ratio of a number of poor people divided by a number of total population, $a_{0}$ is the average of $\bar{d}_{i}$ among poor people, and $M_{0}\in[0,1]$ is the MPI index. The higher $M_{0}$ implies higher issues of poverty.

Based on the MPI concept, however, MPI cannot be used to find resolutions of common problems and dependency among dimensions (e.g. income and health issues). Hence, our aim of study is to identify whether other dimensions (e.g. health, education) affect incomes of households. We are interested in identifying the resolutions of common issues from specific dimensions that affect people’s incomes. Given that there are different linear dependencies between synthetic poverty dimensions ($X$ variable) and income ($Y$ variable) in a simulation data in various resolutions, we want to firstly test whether our proposed framework can identify the correct resolution of models that share common dependencies between income and independent dimensions. Hence, the main purpose of our simulation is to be used as a sanity check whether our framework can correctly identify resolution w.r.t. our ground truth before deploying the proposed framework to the real world datasets.

5.2. Real-world data: Thailand’s population household information

We obtained the dataset of Thailand household-population surveys from Thai government. The surveys were collected in 2019. The surveys were used for estimating the Multidimensional Poverty Index (MPI) (Alkire and Santos, 2010; Alkire et al., 2019), which is the main poverty index that UN is currently using. We used the surveys from two provinces, Khon Kaen and Chiang Mai, to perform our analysis. There are 353,910 households for Khon Kaen, and 378,465 households for Chiang Mai. For each individual household, there are 30 dimensions that describe the characteristics of each household. Each dimension has three possible values: 1 means a good condition, 0 means no data, and -1 means a bad condition. These 30 dimensions can be categorized into five aspects: health, financial status, education, access of public services, and the living conditions. We remove the average household income dimension from the set of independent variable and make it as a dependent variable for prediction. There are five layers of Thailand administrative divisions: 1) the nation (Thailand), 2) provinces (e.g. Khon Kaen and Chiang Mai), 3) amphoes, 4) tambons, and 5) villages. We created a Multi-resolution cluster set from these administrative division. Our task is to infer which layer of these administrative divisions are predictive for income prediction given the household information.

Lastly, the distributions of real datasets are shown in Figure 4.

5.3. Baseline method: greedy algorithm

To the best of our knowledge, since there is no method that we can compare against our approach directly, we compare our Algorithm 3 performance with the greedy algorithm in Algorithm 4 as a main baseline. The greedy algorithm infer the partition by greedily selecting the clusters that have the lowest root mean square error (RMSE) into the output partition. Hence, no algorithm can provide the partition that has the total RMSE as low as Algorithm 4. On the line 3 of Algorithm 4, the RMSE of the cluster $C_{j,k}$ is calculated by fitting all data points in $C_{j,k}$ to make a model using linear regression. Then, we use the predicted values of the dependent variable from the inferred model and the true values of the dependent variable to compute RMSE. In the result section, we will show that our approach in Algorithm 3 provides the maximal homogeneous partition that has the total RMSE as low as Algorithm 4 partition’s.

The least square approach has the time complexity as $\mathcal{O}(n^{2}d)$ where $n$ is a number of individuals and $d$ is a number of $X$ dimensions. Given $n_{max}$ is a number of individuals in the largest cluster and $|\mathcal{C}|$ is a total number of clusters from all layers. Algorithm 4 has a time complexity as $\mathcal{O}(|\mathcal{C}|n_{max}^{2}d)$ where $n_{max}\leq n$. The time complexity is dominated by the code line 1-4 in Algorithm 4 that infers the model from data using the least square approach, which also occurs in Algorithm 3. This makes both Algorithm 3 and Algorithm 4 have the same time complexity. The lower bound is $\Omega(n^{2}d)$.

5.4. Baseline method: finite mixtures of regression models

We compare our approach with the finite mixtures of regression models (MR) developed by  (Grün et al., 2007; Leisch, 2004; Grün and Leisch, 2006). The mixture-model implementation is in the R package “flexmix”. We set a number of components of mixture models (parameter $k$) corresponding to a number of homogeneous clusters within each dataset. We use a mixture model to show that even if mixtures of linear regression knows the number of homogeneous clusters, in complicated scenarios, its performance is quite unstable. In contrast, our framework that deploys the simple linear regression with multi-resolution partitions can perform better than mixture models in the same complicated datasets.

Additionally, given $n$ is a number of individuals and $k=|\mathcal{C}|$ is a number of clusters. Algorithm 3 has a time complexity as $\mathcal{O}(k\times n^{2}\times d)$ while the time complexity of mixtures models parameter estimating using Expectation-maximization (EM) algorithm (Dempster et al., 1977) (deployed by flexmix) is $\mathcal{O}(k\times t_{0}\times n\times d)$ where $t_{0}\in[1,\infty)$ is a number of time steps that EM algorithm converges w.r.t. some predefined threshold. This makes mixtures of regression models running time is unknown compared to our approach.

5.5. Software and hardware used in the analysis

The computer specification that we used in this experiment is the Lenovo Thinkpad T480s, with CPU Intel Core i7-8650U 1.9 GHz, and RAM 16 GB. The software we used to conduct the experiment is the R studio version 1.2.5033 based on the R version 3.6.2. The R packages we used in the analysis are igraph (Csardi and Nepusz, 2006), caret (Kuhn, 2020), and ggplot2 (Wickham, 2016). All experiments were conducted on Microsoft Window 10. The implementation of our framework is in the form of R package with documentation that can be found at (Amornbunchornvej, 2020).

6. Results

6.1. Simulation

In this section, we report the results of our analysis from simulation datasets (Section 5.1). We compared our approach (Algorithm 3) with the baseline methods from Section 5.3 and Section  5.4.

The result of output partition from our approach is in Fig. 5. The red nodes are clusters selected as members of maximal homogeneous partition. Comparing to the ground truth in Fig. 2, our approach can infer the maximal homogeneous partition $100\%$ correctly whether the maximal homogeneous partition consists of clusters at 1st layer (Type-1 datasets), 2nd layer (Type-2 datasets), 3rd layer (Type-3 datasets), or multiple layers (type-4 datasets). In contrast, the result of output partition from the greedy algorithm in Fig. 6 shows that it mostly selected the clusters in the last layers as the members of the partition. The greedy algorithm even included the lower clusters that are the subset of homogeneous cluster from the upper layer as the members of output partition. This is because the average RMSE of clusters from lower layer are typically lower than the upper layer clusters. However, the difference between average RMSE of lower and upper layer are not significant (see T1-T3 Datasets in Table 4). In fact, the lower-layer clusters trend to have a lower RMSE compared to their super-set cluster because the OLS seems to fit noise easier with less number of individuals. This indicates that the greedy algorithm is sensitive to noise, which implies that it has an over-fitting issue.

Table 4 illustrates the RMSE comparison between different methods. In all datasets, our approach (OPT) results are the same as the greedy one. In contrast, linear regression perform poorly in all datasets except the Type-1 datasets. This is because the Type-1 datasets have only a single homogeneous cluster, which implies that there is only a single linear function $h(x)$ that generated $y$ from $x$. However, Type-2/3/4 datasets have more than one linear function $h(x)$. By fitting a single linear regression function to these datasets, we expect to have a huge RMSE for these datasets (see Proposition 3.6). For a mixture of regression (Section 5.4), it performs well for a dataset that possesses a low number of homogeneous clusters. Nevertheless, for the complicated datasets (Type-4) that has 13 homogeneous clusters in the different layers, the mixture of regression was unable to perform well. This indicates that the mixture of regression, despite of the generalization of linear regression, is quite unstable when the number of homogeneous clusters rises. The last column is the RMSE from the residual of difference between any $y_{i}$ with the mean of $y$ that we designate it as a null model. The results demonstrate that OPT and Greedy perform a lot better than the null model.

In the aspect of performance of inferring the maximal homogeneous partition, we reported the results of methods in Table 5. In the table, our method (OPT) can infer maximal homogeneous partition perfectly while the greedy performed mostly poorly except in Type-3 datasets. This is because the greedy algorithm reports mostly the clusters at the last layer to be members of output partition and Type-3 datasets have all clusters members of the maximal homogeneous partition at the last layer. For a mixture of regression, despite of no access to multi-resolution partitions, it can retrieve homogeneous clusters well for datasets that have few homogeneous clusters (Type-1/2/3). However, it performs fairly in the complicated Type-4 datasets.

We also extended our approach and other methods to be able to fit data using the exponential function. We reported the results of using nonlinear datasets in the analysis. Table 6 shows the average RMSE values of approaches on nonlinear datasets: polynomial and exponential datasets. Our method (OPT) is not the best compared to Greedy and Mixture of Regression methods in the aspect of reducing RMSE.

The reason that the Mixture of regression outperforms OPT is that it tries to find the fitting models without boundary constraints of partitions, while OPT fits models w.r.t. the given boundary constraints. It is possible that there are some partitions that are not the same as the given partitions from the input that have lower RMSEs when we fit models using them. If such partitions exist, then the Mixture of regression might use the unconstrained partitions to fit its models. Hence, the RMSEs of the Mixture regression can be lower than OPT’s in this case. Additionally, the assumption we implicitly made was that if all methods fit models under the partition constraints, then the Greedy has the lowest RMSE. Because the Mixture regression breaks this assumption, its RMSE can be lower than the RMSE of the greedy algorithm.

However, when we consider the average precision, recall, and F1 score values in Table 7 for the task of inferring the correct maximal homogeneous partitions, OPT performed better than all methods.

Comparing linear and exponential model fitting, in Table 7 (LIN) rows, linear model fitting of OPT performed well in exponential datasets but it did not perform well on the polynomial datasets. In contrast, the result of exponential model fitting in (EXP) rows shows that OPT performed well in both exponential and polynomial datasets. This indicates that our approach is able to be extended to infer maximal homogeneous partitions in nonlinear datasets well.

In summary, the results in this section indicate that our approach performance has the same RMSE as the greedy one (Table 4), nevertheless, its performance is a lot better in the aspect of inferring the maximal homogeneous partition (Table 5). For mixture of regression model, it performed well for Type-1/2/3 datasets, while its performance decreased for a Type-4 datasets, which have clusters in the maximal homogeneous partition from multiple layers. In fact, the homogeneous clusters that are inferred by a mixture of regression might not be consistent with a given multi-resolution partitions. However, we can deploy mixtures of regression to approximate homogeneous clusters when datasets come without multi-resolution partitions. But we still need our approach to estimate the degree of homogeneity of each cluster as well as using it to compare different ways of partitioning a population. Moreover, we can even use the mixture of regression as a kernel instead of using linear regression to predict our dependent variable.

6.2. Case study: inferring informative-administrative-level subdivision to predict population incomes

We use household population data (Section 5.2) to demonstrate the application of our framework to help the policy maker to make a policy. In reality, policy maker cannot create the unique policy for every village to solve the poverty issues because of the limited resource. On the other hand, making a single policy cannot solve all poverty issues since each region has their own unique problems. Hence, we propose to use our framework to find the largest level of administrative subdivision that have enough common issues (measuring by $\eta$ in Eq. 16) so that the policy makers do not need to establish a policy for each village. In our framework, a $\gamma$-Homogeneous cluster is considered to be an informative subdivision that the policy makers can establish a policy. In this case study, we set $\gamma=0.05$ which requires the correlation between dependent variable $Y$ and predicted $\hat{Y}$ around 0.22 that is closed to a moderate correlation (the correlation around 0.3 is considered as a moderate correlation (Cohen, 1988, 2013)).

Table 8 illustrates RMSE of each methods from two provinces: Khon Kaen and Chiang Mai. The result shows that our approach (OPT) has a bit higher RMSE compared to the greedy algorithm. In contrast, linear regression and the null model $y_{i}-\bar{y}$ have larger RMSE. Linear regression represents the approach of finding one policy for all regions within a province to predict the income of households, while our approach is trying to find a unique policy for each specific region. The result in Table 8 indicates that each region has its own problem. This is why fitting linear regression for the entire province population performed poorly.

Table 9 shows the number of clusters that are the members of output partition in each method. The first layer is the province level, the second layer is the amphoe level, the third layer is the tambon level, and the fourth layer is the village level. The result indicates that our approach can detect informative subdivisions beyond the last layer while the greedy approach can report only the clusters in the last layer.

After we got the maximal homogeneous partition, we need to validate whether the result is consistent with the ground truth of the government records.111 All reference data of government records used here are from https://www.tpmap.in.th/. We used the 2019 records of poor people as a ground truth in this result section. The following areas are examples of members of the maximal homogeneous partition.

First, in Khon Kaen province, the “Subsomboon” tambon, which is on the 3rd layer, has 11 villages in 4th layer. 8 out of 11 villages have the majority of poor people suffering from the lack of the education issue. Second, another tambon in Khon Kaen province is Khu Kham. It has 7 out of 8 villages faces the lack of education issue among poor people. Third, in Chiang Mai province, in San Sai amphoe (2nd layer) has 12 tambons. 10 out of 12 tambons has financial and education as leading issues.

We can see that, for each tambon in the examples above, policy makers can make a single policy for all villages because the majority of subareas in each homogeneous area have almost the same issues.

Next, the non-homogeneous cases are provided below.

In Khon Kaen province, Sila tambon (3rd layer) consists of 28 villages. There are 11 villages that have no poor people. There are 11 villages that have financial issues among the majority of poor people. There are 5 villages facing health issues among the majority of poor people in each village. One village has equal numbers of poor people facing financial and health issues. We can see that the policy makers need to group villages in this tambon w.r.t. the issue types before making policies that fit each type of issue.

In Chiang Mai province, Chiang Dao amphoe (2nd layer) has 7 tambons. Two tambons have financial issues among the majority of poor people. Two tambons have financial and education issues. Two tambons have financial and health issues. Lastly, one tambon has the education issue as a main problem. It is quite challenging for policy makers to place a single policy for the entire Chiang Dao amphoe.

The next one is the example of how policy makers can use our system to combat poverty. One of the tambons in the output maximal homogeneous partition of our approach from Khon Kaen province is “Subsomboon”. We used the t-test to find the variable importance of coefficients of regression that have values far from zero for this tambon. The null hypothesis is that the absolute of specific coefficient is zero, while the alternative hypothesis is that the absolute of coefficient is greater than zero. We reject the null hypothesis at $\alpha=0.01$.

Out of 30 coefficients, there are only four coefficients that we can reject the null hypothesis. However, the only positive coefficient is the indicator of whether a household has a saving account with a bank. This implies that having saving account is associated with income. The policy maker should consider why some households can have saving account to make a policy of combating poverty in Subsomboon tambon. For Chiang Mai province, San Sai amphoe is one of the informative clusters in the output partition. The variables that are important and have positive coefficients mostly are health issues in the household that might make the members of household cannot work efficiently. Hence, the policy maker should consider to make a policy that supports the health of people in the area.

In practice, the greedy algorithm performs well only when all homogeneous clusters are mostly in the last layer. Even in this case, its performance is not significantly better than OPT. Moreover, the time and space complexities of both algorithms are the same. Hence, OPT should always be used to find the maximal homogeneous partition.

7. Discussion

7.1. Inferring multiresolution partitions

In practice, the obvious case that a dataset of multiresolution partitions can be found is a dataset that is related to administrative subdivisions of country. In each sub-districts in a specific resolution, the data can be varied, such as household income, types of land use, distribution of labor force, etc. In the eyes of policy makers, knowing which area share similar models for target properties make them have an easier way to manage resource or to place policies.

However, in some case, there is no given multiresolution partitions. For example, suppose we divide a natural farming area into a grid and we have information about history of land use of each block in the grid (e.g. the remaining resource). Our goal is to infer partitions of multiple neighbor blocks that we can declare as either preserved areas to protect the land or non-preserved areas. In this situation, we can use hierarchical clustering (Sibson, 1973; Defays, 1977; Jr., 1963) (hclust function in R programming (R Core Team, 2019) ) to build a hierarchical tree of blocks. We can assign layers of partitions based on the distance of nodes from the nearest leaf. Another approach is to use the mixture regression method (Grün et al., 2007; Leisch, 2004; Grün and Leisch, 2006) in Section 5.4. We might run the mixture regression method to get the first-layer result, then, for each inferred partition, we might run the method again to find the sub-districts. Nevertheless, this approach requires users to set an appropriate number of clusters $k$.

7.2. Extension to non-linear models

It is possible to extend our work to non-linear models if we care only a number of bits and information we need to explain models. However, in machine learning, there are issues of underfitting and overfitting that we should consider.

To resolve underfitting issues, we try to find $h$ that minimizes $\mathcal{L}(S|h)$ in Eq. 5. Nevertheless, a complex model typically fits data well, which makes $\mathcal{L}(S|h)$ small. To resolve the overfitting issue, since we prefer simpler model that has its performance close to complicated model, we try to find $h$ that also minimizes $\mathcal{L}(h)$ in Eq. 4.

We can say that $\mathcal{L}(h)$ represents the model complexity. However, there are many challenges that we need to address in order to measure $\mathcal{L}(h)$ for any arbitrary $h$.

It is straightforward to compare $\mathcal{L}(h_{i})$ against $\mathcal{L}(h_{j})$ from the same function class (e.g. a linear function class, an exponential function class, etc) by measuring a number of bits we need to keep coefficients and terms. For example, the model complexity of one linear function that has less number of terms with small-size coefficients can be considered as a simpler model than a complicated linear function that has many terms with large coefficients. On the contrary, it is challenging to compare two models from different classes and declare that one model is less complicated than another (e.g. sine function vs. logarithmic function).

One of the possible ways to measure the complexity is to utilize the concept of the Vapnik-Chervonenkis (VC) dimension (Vapnik and Chervonenkis, 1971) in learning theory. The VC dimension of a function class $\mathbb{H}$ can be used to measure a sample complexity (Abu-Mostafa et al., 2012; Mohri et al., 2018); a model with a higher VC dimension requires more data to achieve the same generalization performance compared to a model with a lower VC dimension. In our context, we can give a higher number of bits to $\mathcal{L}(h)$ for a model from a function class with the higher VC dimension. Nevertheless, not all function classes have VC dimensions. Alternatively, in learning theory, we can use Rademacher complexity (Koltchinskii, 2001) that measures the richness of a function family $\mathbb{H}$ in term of how well $\mathbb{H}$ can fit random noise. The more complex function class can fit random noise better than the simpler function class (Mohri et al., 2018). Hence, we make $\mathcal{L}(h)$ to have more bits if a function class fits random noise better than another function class.

7.3. Sparsity regularization

In the context of sparsity regularization, we also want to get the model that has the minimum number of dimensions we used in the model along with the prediction performance $\mathcal{L}(S|h)$ (Eq. 5) and the model complexity $\mathcal{L}(h)$ (Eq. 4). Suppose $\beta=(p_{1},\dots,p_{l})$ is a vector of parameters for a model $h$, given a set of data $S=\{(x_{1},y_{1}),\dots,(x_{n},y_{n})\}$ and a penalty weight $\lambda$, we can have the following optimization problem to find the optimal $\beta$ from data that utilizes the sparsity regularization.

[TABLE]

Where $\left\|\beta\right\|_{0}$ is the number of parameters in $\beta$ that has nonzero values, $f_{\beta}(x_{i})$ is a function that returns predicted value of $y_{i}$ using $\beta$ and $x_{i}$, and $loss()$ is a loss function. Based on the optimization problem above, we can modify $\mathcal{L}(h)$ to utilize the sparsity regularization.

[TABLE]

Hence, we have

[TABLE]

The Eq. 22 above encourages the model $h$ that has a good performance (minimizing $\mathcal{L}(S|h)$) as well as using less parameters with the parameters that have small magnitudes (minimizing $\mathcal{L}(h,\lambda)$).

7.4. Encoding methods for real numbers

As a reminder, throughout this paper, we have Assumption 1 states as follows.

For any $x_{1},x_{2}\in\mathbb{R}$ s.t. $|x_{1}|\geq|x_{2}|$, we have

[TABLE]

This assumption is true for the Single-Precision Floating-point Format (the IEEE Standard 754 (Hough, 2019)) that is widely used in several well-known programming languages (e.g. Fortran, C, C++, C#, Java) to encode real numbers. However, the floating-point format makes every real number use an equal number of bits: 32 bits per number.

Nevertheless, the work by  (Lindstrom et al., 2018) provided the framework to implement the universal coding for real numbers. In this coding, a number of bits varies proportional to a real-number magnitude. Specifically, for any $y\in\mathbb{R}$,

[TABLE]

where $|\mathcal{E}(y)|$ is a number of bits of real-number representation in (Lindstrom et al., 2018)’s work. In practice, if a specific application has high variation of magnitudes of real numbers, then the encoding in  (Lindstrom et al., 2018) might be able to reduce spaces requires for keeping real numbers and our framework can obviously gain this benefit. For a specific application, one might check whether the encoding scheme for real numbers that complied with our assumption is suitable for being used under the application.

8. Conclusion

In this paper, we addressed the challenges encountered by many policy makers: to find the largest possible areas with common problems where similar policies can be implemented and limited resources are optimally utilized. Given a set of target variable and predictors, as well as a set of multi-resolution clusters that represent administrative divisions as inputs, we proposed a framework to infer a set of largest-informative clusters that have efficient predictors. Efficient predictors of each informative cluster cover the entire cluster population. We used both simulation datasets and real-world dataset of Thailand’s population household information to evaluate and to illustrate the application of our framework. The results showed that our framework performed better than all baseline methods. Moreover, in Thailand’s population household information, the framework can infer the predictors of people income for the specific areas, which can be used to guide policy makers to utilize this information to combat poverty. Particularly, simulated data shows that our approach can estimate homogeneity of each cluster and compare different ways of population partition, which would enhance the formulation of evidence-based policies as well as their assessments. This finding is confirmed with real-world data from the two major provinces of Thailand where it is possible to cluster regions with similar problems up to the 2nd (Amphoe) level. More importantly, our framework can be used beyond the context of income prediction but for any kind of regression analysis on multi-resolution clusters. Lastly, we also provide the R package, MRReg, which is the implementation of our framework in R language with the manual at (Amornbunchornvej, 2020). The official link for MRReg at the Comprehensive R Archive Network (CRAN) can be found at https://cran.r-project.org/package=MRReg.

Acknowledgment

The authors would like to thank the National Electronics and Computer Technology Center (NECTEC), Thailand, to provide our resources in order to successfully finish this work.

This paper was supported in part by the Thai People Map and Analytics Platform (TPMAP), a joint project between the office of National Economic and Social Development Council (NESDC) and the National Electronic and Computer Technology Center (NECTEC), National Science and Technology Development Agency (NSTDA), Thailand.