How to deal with hierarchical / nested data in machine learning

Question

I'll explain my problem with an example. Suppose you want to predict the income of an individual given some attributes: {Age, Gender, Country, Region, City}. You have a training dataset like so

train <- data.frame(CountryID=c(1,1,1,1, 2,2,2,2, 3,3,3,3), 
             RegionID=c(1,1,1,2, 3,3,4,4, 5,5,5,5), 
             CityID=c(1,1,2,3, 4,5,6,6, 7,7,7,8), 
             Age=c(23,48,62,63, 25,41,45,19, 37,41,31,50), 
             Gender=factor(c("M","F","M","F", "M","F","M","F", "F","F","F","M")),
             Income=c(31,42,71,65, 50,51,101,38, 47,50,55,23))
train
   CountryID RegionID CityID Age Gender Income
1          1        1      1  23      M     31
2          1        1      1  48      F     42
3          1        1      2  62      M     71
4          1        2      3  63      F     65
5          2        3      4  25      M     50
6          2        3      5  41      F     51
7          2        4      6  45      M    101
8          2        4      6  19      F     38
9          3        5      7  37      F     47
10         3        5      7  41      F     50
11         3        5      7  31      F     55
12         3        5      8  50      M     23

Now suppose I want to predict the income of a new person who lives in City 7. My training set has a whopping 3 samples with people in City 7 (assume this is a lot) so I can probably use the average income in City 7 to predict the income of this new individual.

Now suppose I want to predict the income of a new person who lives in City 2. My training set only has 1 sample with City 2 so the average income in City 2 probably isn't a reliable predictor. But I can probably use the average income in Region 1.

Extrapolating this idea a bit, I can transform my training dataset as

    Age Gender CountrySamples CountryIncome RegionSamples RegionIncome CitySamples CityIncome
 1:  23      M              4         52.25             3        48.00           2    36.5000
 2:  48      F              4         52.25             3        48.00           2    36.5000
 3:  62      M              4         52.25             3        48.00           1    71.0000
 4:  63      F              4         52.25             1        65.00           1    65.0000
 5:  25      M              4         60.00             2        50.50           1    50.0000
 6:  41      F              4         60.00             2        50.50           1    51.0000
 7:  45      M              4         60.00             2        69.50           2    69.5000
 8:  19      F              4         60.00             2        69.50           2    69.5000
 9:  37      F              4         43.75             4        43.75           3    50.6667
10:  41      F              4         43.75             4        43.75           3    50.6667
11:  31      F              4         43.75             4        43.75           3    50.6667
12:  50      M              4         43.75             4        43.75           1    23.0000

So, the goal is to somehow combine the average CityIncome, RegionIncome, and CountryIncome while using the number of training samples for each to give a weight/credibility to each value. (Ideally, still including information from Age and Gender.)

What are tips for solving this type of problem? I prefer to use tree based models like random forest or gradient boosting, but I'm having trouble getting these to perform well.

UPDATE

For anyone willing to take a stab at this problem, I've generated sample data to test your proposed solution here.

Answer 1

I have been thinking about this problem for a while, with inspirations from the following questions on this site.

• How can I include random effects into a randomForest?

• Random forest on grouped data

• Random Forests / adaboost in panel regression setting

• Random forest for binary panel data

• Modelling clustered data using boosted regression trees

Let me first introduce the mixed-effects models for hierarchical/nested data and start from a simple two-level model (samples nested within cities). For the $j$-th sample in the $i$-th city, we write the outcome $y_{ij}$ as a function of covariates $\boldsymbol x_{ij}$ (a list of variables including gender and age), $$ y_{ij}=f(\boldsymbol x_{ij})+{u_i}+\epsilon_{ij},$$ where ${u_i}$ is the random intercept for each city, $j=1,\ldots,n_i$. If we assume $u_i$ and $\epsilon_{ij}$ follow normal distributions with mean 0 and variances $\sigma^2_u$ and $\sigma^2$, the empirical Bayesian (EB) estimate of $u_i$ is $$\hat{u}_i=\frac{\sigma^2_u}{\sigma^2_u+\sigma^2/n_i}(\bar{\mathbf{y}}_{i.}-f(\bar{\boldsymbol x}_{i.})),$$ where $\bar{\mathbf{y}}_{i.}=\frac{1}{n_i}\sum_i^{n_i}y_{ij}$, $f(\bar{\boldsymbol x}_{i.})=\frac{1}{n_i}\sum_i^{n_i}f(\boldsymbol x_{ij}).$ If we treat $(\bar{\mathbf{y}}_{i.}-f(\bar{\boldsymbol x}_{i.}))$ as the OLS (ordinary least square) estimate of $u_i$, then the EB estimate is a weighted sum of 0 and the OLS estimate, and the weight is an increasing function of the sample size $n_i$. The final prediction is $$\hat{f}(\boldsymbol x_{ij})+\hat{u}_{i},$$ where $\hat{f}(\boldsymbol x_{ij})$ is the estimate of the fixed effect from linear regression or machine learning method like random forest. This can be easily extended to any level of data, say samples nested in cities and then regions and then countries. Other than the tree-based methods, there is a method based on SVM.

For random-forest-based method, you can try MixRF() in our R package MixRF on CRAN.

Answer 2

Given that you only have two variables and straightforward nesting, I would echo the comments of others mentioning a hierarchical Bayes model. You mention a preference for tree-based methods, but is there a particular reason for this? With a minimal number of predictors, I find that linearity is often a valid assumption that works well, and any model mis-specification could easily be checked via residual plots.

If you did have a large number of predictors, the RF example based on the EM approach mentioned by @Randel would certainly be an option. One other option I haven't seen yet is to use model-based boosting (available via the mboost package in R). Essentially, this approach allows you to estimate the functional form of your fixed-effects using various base learners (linear and non-linear), and the random effects estimates are approximated using a ridge-based penalty for all levels in that particular factor. This paper is a pretty nice tutorial (random effects base learners are discussed on page 11).

I took a look at your sample data, but it looks like it only has the random effects variables of City, Region, and Country. In this case, it would only be useful to calculate the Empirical Bayes estimates for those factors, independent of any predictors. That might actually be a good exercise to start with in general, as maybe the higher levels (Country, for example), have minimal variance explained in the outcome, and so it probably wouldn't be worthwhile to add them in your model.

Answer 3

This is more of a comment or suggestion rather than an answer, but I think you ask an important question here. As someone who works exclusively with multilevel data, I can say that I have found very little about machine learning with multilevel data. However, Dan Martin, a PhD graduate in quantitative psychology at the University of Virginia, did his dissertation on the use of regression trees with multilevel data. Below is a link to an R package he wrote for some of these purposes:

https://github.com/dpmartin42/mleda/blob/master/README.md

Also, you can find his dissertation here:

http://dpmartin42.github.io/about.html

Answer 4

You can use a mixed effect model that models the ID variables as random effects. By doing so, you allow for information pooling: you use both data from the global average and from the group averages, and the less data you have per group, the more weight is given to the global average. If you want to use a machine learning model for the fixed effects part, you can, for instance, use tree-boosting.

The GPBoost library with Python and R packages builds on LightGBM and allows for combining tree-boosting and mixed effects models. Simply speaking it is an extension of linear mixed effects models where the fixed-effects are learned using tree-boosting. See this blog post and Sigrist (2020) for further information. Below is some sample R code on how to run the analysis in your case. Note that currently the R package is not yet on CRAN (hopefully soon) but the the Python package is available on PyPI.

Disclaimer: I am the author of the GPBoost library.

library(gpboost)
train <- data.frame(CountryID=c(1,1,1,1, 2,2,2,2, 3,3,3,3), 
                    RegionID=c(1,1,1,2, 3,3,4,4, 5,5,5,5), 
                    CityID=c(1,1,2,3, 4,5,6,6, 7,7,7,8), 
                    Age=c(23,48,62,63, 25,41,45,19, 37,41,31,50), 
                    Gender=factor(c("M","F","M","F", "M","F","M","F", "F","F","F","M")),
                    Income=c(31,42,71,65, 50,51,101,38, 47,50,55,23))
# Prepare data
X <- as.matrix(cbind(Gender=train[,"Gender"],Age=train[,"Age"]))# fixed effects data
group_data <- train[,c("CountryID","RegionID","CityID")]# grouping data for random effects
y <- train[,c("Income")]# response variable
# Define a random effects model
gp_model <- GPModel(group_data = group_data)
# Run boosting algorithm (this will not give meaningfull results as the data is too small)
bst <- gpboost(data = X, label = y, gp_model = gp_model, verbose = -1,
               objective = "regression_l2",  nrounds=10, learning_rate=0.1)
# Show estimated variance parameters
summary(gp_model)

# A linear mixed effects model also has problems with this small data and some of the variances are 0
gp_model <- fitGPModel(group_data=group_data, y=y, X=cbind(Intercept = rep(1,length(y)),X))
summary(gp_model)
# Or the same thing using the lme4 package
library(lme4)
mod <- lmer(Income ~ Age + Gender + (1|CountryID) + (1|RegionID) + (1|CityID), data=train, REML=FALSE)
summary(mod)

Answer 5

The function RFcluster() from the gamclass package for R "adapts random forests to work (albeit clumsily and inefficiently) with clustered categorical outcome data". The following example is from the help page for RFcluster:

 library(randomForest)
 library(gamclass)
 data(mlbench::Vowel)
 RFcluster(formula=Class ~., id = V1, data = Vowel, nfold = 15,
           tree=500, progress=TRUE, printit = TRUE, seed = 29)

This returns an OOB accuracy (where the "bags" are bags of speakers, not bags of individual speaker samples), that my machine gives as 0.57 .

Answer 6

R package glmertree allows for fitting decision trees to multilevel and longitudinal data (which would otherwise be modeled with a mixed-effects model).

It allows for specifying a random effects structure, and partitions the dataset into subgroups using level-I, or higher-level predictors.

For further reference, see the package vignette (tutorial): https://cran.r-project.org/web/packages/glmertree/vignettes/glmertree.pdf.

Fokkema, M., Smits, N., Zeileis, A., Hothorn, T., & Kelderman, H. (2018). Detecting treatment-subgroup interactions in clustered data with generalized linear mixed-effects model trees. Behavior research methods, 50(5), 2016-2034. https://doi.org/10.3758/s13428-017-0971-x

Fokkema, M., Edbrooke-Childs, J., & Wolpert, M. (2021). Generalized linear mixed-model (GLMM) trees: A flexible decision-tree method for multilevel and longitudinal data. Psychotherapy Research, 31(3), 329-341. https://doi.org/10.1080/10503307.2020.1785037

Answer 7

You may want to have a look at metboost: Miller PJ et al. metboost: Exploratory regression analysis with hierarchically clustered data.arXiv:1702.03994

Quote from the abstract: We propose an extension to boosted decision decision trees called metboost for hierarchically clustered data. It works by constraining the structure of each tree to be the same across groups, but allowing the terminal node means to differ. This allows predictors and split points to lead to different predictions within each group, and approximates nonlinear group specific effects. Importantly, metboost remains computationally feasible for thousands of observations and hundreds of predictors that may contain missing values.

It is implemented in the R package mvtboost