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))
   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.


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

  • 7
    $\begingroup$ Hierarchical bayesian models are very natural to exploit structure as described in your data. Check out the classical example on Radon contamination modeling mc-stan.org/documentation/case-studies/radon.html $\endgroup$ – Vladislavs Dovgalecs Jun 30 '16 at 6:09
  • $\begingroup$ Take a look at this Kaggle Getting Started tutorial challenge: kaggle.com/c/titanic. It deals with a similar problem, that is to predict if a person has survived the Titanic disaster given various attributes about the person such as Gender, ticket type etc. The best submitted solutions for this use advanced methods such as Gradient Boosting and Hierarchical Bayesian Models etc. $\endgroup$ – Vihari Piratla Jul 2 '16 at 5:05
  • 7
    $\begingroup$ @VihariPiratla Thanks for the input, but I'm familiar with the Titanic dataset and challenge, and I don't see how it relates to the nested data issue I asked about. $\endgroup$ – Ben Jul 2 '16 at 5:40
  • $\begingroup$ Using L2/L1 regularisation for those models eg logistic regression/SVM that perform it should help (poor mans hierarchical bayes). effectively you penalize coefficients, so unless a coefficient significantly ( ie on lots of data) improves error it will be set to close to zero. and you use crossvalidation to decide level of penalisation $\endgroup$ – seanv507 Jul 6 '16 at 7:52
  • $\begingroup$ Do you also want to be able to predict the income of a person who lives in city 9? $\endgroup$ – jan-glx Dec 8 '17 at 14:02

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

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.

  • $\begingroup$ Could you explain why the intercept is allowed to vary with $i$, but $f$ has a single set of parameters for all $i$? Is it a simplifying assumption to avoid overfitting or make the problem tractable? $\endgroup$ – user20160 Jul 5 '16 at 7:34
  • $\begingroup$ @user20160 Yes, we call $f$ part as fixed-effects part, and $u_i$ as the random effect. We can also have some random slopes like $\mathbf{x}_{ij}^{'} \mathbf{u}_{i},$ but should limit the number of random slopes because this would introduce a large number of variance components in the covariance matrix of the random-effects if it is unstructured. $\endgroup$ – Randel Jul 5 '16 at 15:24

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.

  • 1
    $\begingroup$ +1 for introducing the mboost package. $\endgroup$ – Randel Jul 6 '16 at 20:40
  • $\begingroup$ The real data I'm working with has a lot more than two variables which add a lot of real-world mess that isn't picked up by my simple example (e.g. non linearity, codependence, missing values, categorical values, etc.). In my experience, tree based learners do the best job at handling all the real-world mess, which is why I lean towards using them. (There's a reason XGBoost wins almost all the structured data competitions on Kaggle.) mboost sounds interesting - I'll give it a look. Thanks $\endgroup$ – Ben Jul 6 '16 at 22:38
  • 1
    $\begingroup$ I see, agreed that trees can certainly be useful in that situation. In that case sticking with @Randel's suggestion would be a good option. Mboost also has a tree base learner as well which might prove useful in combination with the random effects base learners. $\endgroup$ – dmartin Jul 6 '16 at 23:27

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 recent 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:


Also, you can find his dissertation here:


  • 3
    $\begingroup$ Thanks for the reference, Erik! One thing to keep in mind with my dissertation is that it only evaluated how forest methods perform "out of the box" when exposed to multi-level data structures. Methods that @Randel mentioned are discussed in Chapter 3 as a literature review, but the simulation study investigated mainly Breiman's original algorithm and conditional inference trees/forests. $\endgroup$ – dmartin Jul 6 '16 at 19:25

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:

 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 .

  • $\begingroup$ To anyone interested in this, RFcluster appears to be just a random forest cross-validated such that the id variable is taken into account when constructing the folds. It does not in any way take the id information into account in the fitting procedure. $\endgroup$ – Frans Rodenburg Apr 12 at 9:31

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


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


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.

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

# 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))
# Or the same thing using the lme4 package
mod <- lmer(Income ~ Age + Gender + (1|CountryID) + (1|RegionID) + (1|CityID), data=train, REML=FALSE)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.