Mixed-effect models with trees

Tree-based regression models are a class of statistical models for predicting continuous response variables when the shape of the regression function is unknown. They naturally take into account both non-linearities and interactions. However, they struggle with linear and quasi-linear effects and assume iid data. This article proposes two new algorithms for jointly estimating an interpretable predictive mixed-effect model with two components: a linear part, capturing the main effects, and a non-parametric component consisting of three trees for capturing non-linearities and interactions among individual-level predictors, among cluster-level predictors or cross-level. The first proposed algorithm focuses on prediction. The second one is an extension which implements a post-selection inference strategy to provide valid inference. The performance of the two algorithms is validated via Monte Carlo studies. An application on INVALSI data illustrates the potentiality of the proposed approach.


Introduction
Mixed-effect or multilevel models (Snijders and Bosker 2012;Pinheiro and Bates 2006) are a valuable class of models able to deal with hierarchical/clustered data. Typical hierarchical data consist of statistical units (level 1 units) nested into clus-ters (level 2 units). Classic examples are students clustered within schools (individual cross-sectional data) or children's growth evaluated at several time points (repeated measures). Mixed-effect models consider unit clustering, including both fixed and random effects in the model. The standard linearity assumption for the fixed effects is too stringent in many situations, and a more flexible model specification, including non-linear effects and interactions, might be required. A worthwhile approach exploits regression trees (CART) (Breiman et al. 1984) to capture non-linear fixed effects and interactions via a piece-wise constant regression function. Several tree-based algorithms have been proposed in the literature to improve CART performance in different ways. For instance, see Loh (2002), Hothorn et al. (2006), Dusseldorp et al. (2010) and the subsequent literature to ensemble learning.
Tree-based models typically assume iid data and therefore ignore a clustered data structure when present. As it is well known that ignoring the clustering structure can lead to biased inference (see, e.g. Bryk and Raudenbush 2001, for linear models), some proposals have been suggested in the literature to adapt tree-based models to clustered, multilevel and longitudinal data. In the framework of longitudinal data, Segal (1992) was, up to our knowledge, the first to deal with this topic, generalizing the CART algorithm and its loss function to the case of correlated multiple binary responses. Zhang (1998) discussed the case of multiple binary response variables using as impurity measure the generalized entropy criterion, linked to the log-likelihood of a specific exponential family distribution. Other contributions are due to Abdolell et al. (2002), Loh and Zheng (2013), Eo and Cho (2014), among many. Being specific for longitudinal data, these solutions cannot be adopted when the clustered data do not have a regular structure.
In the general framework of multilevel data, regression trees have been extended to clustered data by Hajjem et al. (2011) and Sela and Simonoff (2012) modelling fixed effects with a decision tree while automatically accounting for random effects with a linear mixed model in a separate step. In particular, Hajjem et al. (2011) proposed mixed-effect regression trees (MERT) where CART selects the tree that models the fixed effect part, and a node-invariant linear structure is used to model the random effects. The algorithm is implemented within the framework of the expectationmaximization algorithm. Sela and Simonoff (2012) independently proposed a similar solution, called Random Effects/EM trees (RE-EM tree). It is shown that random effect regression trees are less sensitive to parametric assumptions and provide substantially improved predictive power compared to linear models with random effects and regression trees without random effects. The literature has grown with variants and extensions (e.g. Hajjem et al. 2014;Fu and Simonoff 2015;Hajjem et al. 2017;Miller et al. 2017;Pellagatti et al. 2021). Other solutions, including trees and clustered data, have been proposed in the framework of subgroups analysis. See for instance Fokkema et al. (2018) and Seibold et al. (2019).
This paper proposes a semi-parametric mixed-effect model where trees are added to a linear component for capturing interactions and non-linear effects. The idea of adding a linear component to a tree is not new in the literature. While regression trees can easily capture complex dependencies by reconstructing the entire regression function, they need many fortuitous splits to recreate a linear or quasi-linear function (Friedman et al. 2001). For this reason, tree-based models are said to struggle in fitting linear or additive dependencies. To overcome this issue, Dusseldorp and Meulman (2004) were the first to propose to combine a linear part and a regression tree in a single model, called regression trunk model, that can be efficiently estimated by a Simultaneous Threshold Interaction Modelling Algorithm (STIMA) proposed by Dusseldorp et al. (2010). Here we extend regression trunk models to the case of the mixed-effect models. Moreover, we include multiple trees in the tree component instead of one. For the particular application we have in mind, we consider here the case of three trees. The first tree captures non-linear relationships and interactions among level 1 predictors. The second tree captures non-linear relationships and interactions among level 2 predictors, while the third is devoted to detecting cross-level interactions. This last kind of interaction is of particular interest in multilevel analysis (see, for instance, Bauer and Curran 2005).
Following the reasoning of Efron (2020), the proposed approach lies in between a pure predictive algorithm and a more traditional regression method, providing an easy to interpret mixed-effect model with good predictive properties. The linear component helps to maintain the three trees as short as possible, while the three trees act as weak learners from the point of view of prediction. The resulting model is more easily interpretable than a single tree or a random forest. Moreover, it has a better predictive performance than a linear mixed-effect model. We call this class Three-tree mixedeffect (3Trees) models.
For the estimation of 3Trees models, we propose two different algorithms. The first algorithm relies on a backfitting-like procedure for selecting the three trees component, and then it jointly estimates the linear and tree-based parts. It provides good predictions of the response Y for new units in already observed or new clusters. Therefore, this algorithm is useful when a study has predictive purposes only. The second algorithm modifies the first one by introducing a post-selection inference procedure, which provides a pruning method based on hypothesis tests and valid confidence intervals for the predicted values. Specifically, we apply a split sample procedure (Cox 1975) on clusters. We evaluate the proposed algorithms through simulations. It turns out that both algorithms are computationally efficient and have satisfactory performance.
The remainder of the paper is organized as follows. Section 2 describes the proposed Three-tree mixed-effect model and its interpretation. Section 3 presents the two estimating algorithms. Section 4 reports a simulation study to evaluate the performance of the two algorithms in the case of a random intercept model. Section 5 illustrates an application of the 3Trees model to analyze the Maths scores achieved by Italian children at a national standardized test. The last section offers concluding remarks and outlines directions for future work.

The Three-tree mixed-effect model
We assume that the true data generating process has the form where x i j is the vector of p 1 individual-level predictors, β the associated main fixed effect coefficients, including β 0 as intercept, while z j is the vector of p 2 cluster-level predictors and γ the associated main fixed effect coefficients. Let p = p 1 + p 2 . Here g(·) is an unknown function R p → R, ruling the non-linear fixed-effects, assumed to not vary across clusters. The level 1 errors ε i j are assumed to be iid N (0, σ 2 ε ). The level 2 errors (random effects) u j are supposed to be iid N (0, σ 2 u ). Moreover, the level 1 and level 2 errors are assumed to be independent. Thus, the model assumes that responses of units belonging to the same cluster are positively correlated, with correlation equal to although responses of units belonging to different clusters are uncorrelated. In model (1), for each cluster j, the random quantity b j = β 0 + u j is the so-called random intercept, varying across clusters, with b j ⊥ ⊥ b j for all j = j ∈ {1, . . . , J }.
As the non-linear function g(·) is unknown, we propose a nonparametric approximation based on regression trees. In particular, considering a two-level structure, we propose to use three trees. The first tree accounts for non-linearities at level 1, the second tree accounts for non-linearities at level 2, and the third one accounts for nonlinearities among all predictors. The primary purpose of the third tree is to account for cross-level interactions, namely interactions among level 1 and level 2 predictors. We call this model the Three-tree mixed-effect (3Trees) model In the 3Trees model, the function g(·) of Eq. (1) is approximated by the sum of three regression trees, T t (t = 1, 2, 3), each constructed over a subset of predictors. Specifically, h 1i j ⊆ x i j (level 1 predictors), h 2i j ⊆ z j (level 2 predictors) and h 3i j ⊆ (x i j ∪ z i j ) (level 1 and level 2 predictors). The model is additive in its components, and the tree component stands as a regionspecific intercept. Model (3) can be written equivalently as where R t1 , . . . , R t M t is the partition of the predictor space corresponding to the tree T t . From this point of view, each tree acts as a factor with M t categories. Each category corresponds to a tree leaf, representing a region of the selected partition of the predictor space. The region R tm can be identified by multiplying all the dummy variables defined by the binary splits along the path from the root node to each leaf in the tree T t . Therefore, some classical identifiability constraints, such as sum-to-zero or corner constraints, must be included for each of these three factors. When the unknown regression function can be assumed to be quasi-linear (Wermuth and Cox 1998), the number of leaf nodes M t (t = 1, 2, 3) can be kept small. In this situation, the proposed class of models is most convenient, providing a good balance between performance and interpretability. The presence of three trees instead of one, e.g. as used for iid data in STIMA (Dusseldorp et al. 2010), allows us to achieve the same goodness of fit with trees having smaller depth. In addition, devoting the trees to specific sets of predictors improves the ability of the algorithm to detect different kinds of non-linear effects.
The choice of the subsets of variables h 1i j h 2i j and h 3i j to be included in each tree is driven by subject-matter considerations. As a general practice, we recommend including all the potentially relevant variables and leaving it to the algorithm to select the most pertinent. Notice that the presence of the linear component attenuates the typical issue of the tree-based algorithms with correlated predictors.
It is worth noting that if the regression function has strong non-linearities, with abrupt inversions in the sign of dependence, model (3) may require deep trees to achieve good performance. One can include flexible functions in the linear part to overcome this issue, such as polynomials or splines. In this case, the so-called linear component captures such strong non-linearities, while the tree component mostly picks up the interaction effects.
The main difference of our procedure with respect to previous proposals (Hajjem et al. 2011;Sela and Simonoff 2012), is the inclusion of the linear component x i j β + z j γ in the mixed-effect model (3). This inclusion allows to avoid overfitting and helps interpretation.

Two estimating algorithms
We propose two iterative procedures to select the trees and estimate the model parameters. While previous proposals for trees with clustered data were based on a kind of EM algorithm (see, for instance, Hajjem et al. 2011), we propose a procedure similar to backfitting for the selection of the trees.
Recalling that we denoted the total number of sample units as n tot = J j=1 n j , let us call X the n tot × p 1 matrix of the individual-level predictors, with rows x i j . Let Z be the n tot × p 2 matrix of the cluster-level predictors, with rows z i j = z j for all the units of the same cluster, i = 1, . . . , n j . Denote H L = (X, Z) the n tot × p matrix, with p = p 1 + p 2 of all the predictors included into the linear component, H 1 the matrix n tot × h 1 , h 1 ≤ p 1 , whose columns are selected from X to be considered for the tree T 1 , H 2 the matrix n tot × h 2 , h 2 ≤ p 2 , whose columns are selected from Z to be considered for the tree T 2 and H 3 the matrix n tot × h 3 , h 3 ≤ p, whose columns are selected from H L for T 3 . Finally, call y the n tot × 1 vector of the responses.
The first algorithm we are proposing, called 3Trees-Alg and whose pseudo-code can be found in Algorithm 1, has predictive purposes. It consists of two main steps: L from the fitted model using both fixed and random effects sub-step Tree1: 8 Compute the vector of partial residuals for T 1 as r Fit a regression tree with the CART algorithm for r (l) 1 on H 1 ;

10
Predict y (l) T 1 from the fitted tree sub-step Tree2: 11 Compute the vector of partial residuals for T 2 as r Fit a regression tree with the CART algorithm of r (l) 2 on H 2 ;

13
Predict y (l) T 2 from the fitted tree sub-step Tree3: 14 Compute the vector of partial residuals for T 3 as r Fit a regression tree with the CART algorithm of r the Selection step and the Estimation step. The Selection step aims to find the best trees conditionally on the linear component. The algorithm initialises the sum of the three trees equal to the sample mean and the linear component at the standard least squared estimates. This choice provides a neutral starting point that, in our experience, facilitates algorithm convergence. Then the procedure iteratively runs the CART algorithm (Breiman et al. 1984) on the partial residuals with respect to the other trees and the mixed effect linear component. At the end of each iteration, the algorithm computes the Mean Squared Error (MSE) based on the predictions for each tree and the mixed-effect model. The iterative procedure stops when the convergence criterion is met, namely, the difference between MSE in two successive iterations is below a given threshold, or the maximum number of iterations is achieved. Notice that as the regression function is not smooth, the algorithm can fall into a periodic loop. Hiabu et al. (2021), in a different framework, proved the validity of backfitting in the presence of non-smoothness. Therefore, the Selection step of 3Trees-Alg returns the trees corresponding to the minimum MSE. For the Selection step, one has also to set the tuning parameters of the CART algorithm, such as the maximum depth of the tree, the minimum number of units in a leaf node and the complexity parameter CART. If this latter parameter is set to a value different from zero, then the CART algorithm within 3Trees-Alg performs tree pruning at each iteration via cross-validation.
In the second step of the 3Trees-Alg algorithm, the Estimation step, each tree selected in the Selection step is specified as a factor, and a mixed-effect model with the linear component and the three new factors as in (4) is fitted via maximum likelihood. The step returns the estimates of all the parameters μ tm , with m = 1, . . . , M t , t = 1, 2, 3, for the tree component, the vectors of coefficients for the linear components β and γ, and the variances σ 2 u , σ 2 ε . It is worth pointing out that the tree component is selected non-parametrically. The CART procedure to choose the trees is greedy, and therefore it addresses the splittings of the predictors providing the local largest decrease in MSE. This implies a lack of order invariance in the tree component. Namely, if we permute the order of the trees in the Selection step, the resulting selected trees can be different. The order invariance can be guaranteed when the subsets of variables included in the trees are disjoint, i.e. each variable can participate in the construction of only one tree, and the variables are independent between trees. However, in multilevel settings, we are often interested in detecting cross-level interactions, which requires defining a third tree that includes all predictors. Thus the trees cannot be kept disjoint. In principle, one could try all six possible orderings of trees and choose by predictive performance comparison. However, we suggest the order dictated by the standard strategy of model building in multilevel analysis (Snijders and Bosker 2012), namely first select level 1 predictors, then level 2 predictors, and finally cross-level interactions.
The 3Trees-Alg algorithm, described in Algorithm 1, is implemented in a user written R code (avalaible from the authors), using the package lme4 (Bates et al. 2015a) for the estimation of the mixed models and on rpart (Therneau and Atkinson 2019) for selecting the trees via the CART algorithm.

Algorithm based on post-selection inference
The proposed tree-based mixed-effect model inherits the theoretical properties of parametric models if the model is well specified. If the selected regression function is not assumed to be the true one but rather a good approximation of it, inference is referred to the projection parameters (see, for instance Buja et al. 2019), where the regression function is optimally projected into the space of the semilinear functions in (4). The tree component identification utilizes the CART algorithm that is proved to provide a locally optimal solution (Breiman et al. 1984).
In a context where the model is selected using the data, as in our procedure, classical tools for inference, such as p-values and confidence intervals, are invalid (see, for instance Benjamini 2010; Berk et al. 2013). In particular, as confidence intervals are computed after model selection, the classical procedures provide narrower confidence intervals than due at the nominal level, with an actual coverage below the nominal one. Interesting literature on post-selection inference proposes methods mainly in the context of ordinary least square and Lasso-type estimators for iid data. A simple postselection inference procedure, proper when the interest is to discover dependencies and influential predictors in a data-driven model, is the split sampling procedure (Cox 1975). It involves the random division of the sample into two parts, one used to select the model and one used for inference. This procedure allows deriving conditional (post-selection) tests or confidence intervals. In addition, when required, inference based on sample splitting followed by bootstrap provides assumption-lean, robust confidence intervals for the projection parameters (Rinaldo et al. 2019).
To obtain adequate inference for the parameters of a 3Trees model, we propose adopting a split sampling procedure. At this aim, we need to assume that the maximum size of the selected model is under control. This corresponds to fixing the maximum depth of each tree.
For multilevel data, the standard assumption of iid data applies at the cluster-level. Therefore, sample splitting is obtained by randomly partitioning the J clusters into two subsets, approximately of equal size. This partition induces a partition of the original data set D into two subsets, say D S and D E . The first subset, D S , is used to select the three trees using the Selection step of Algorithm 1. The second subset, D E , is used for inference on model parameters using the standard linear mixed effect procedures and software. The proposed algorithm, called 3Trees-Alg-PSI, is summarized in Algorithm 2.
Algorithm 2: 3Trees-Alg-PSI Backfitting algorithm for the Three-tree linear mixed models using sample splitting procedure

Estimation step:
Compute parameter estimates for model M using the standard maximum likelihood estimator for mixed effect models on the subset of data Applying 3Trees-Alg-PSI, we obtain an interpretable predictive model laying in between exploratory and confirmatory data analysis. The confidence intervals computed using 3Trees-Alg-PSI can also be used to prune a tree when the complexity parameter in the CART sub-step is set to zero or kept very small. For this purpose, consider the two regions R tm and R tm corresponding to the left and right terminal leaves, produced by the same variable splitting of the tree T t . Let μ tm and μ tm be the corresponding coefficients in the 3Trees model as in (4). Pruning the tree corresponds to joining R tm and R tm . Therefore, the tree can be safely pruned whenever the null hypothesis H 0 : μ tm = μ tm vs H A : μ tm = μ tm cannot be rejected or, equivalently, when the confidence intervals for the two parameters overlap. To perform a deeper pruning, one has to check the equality of all the coefficients involved in the pruning. The pruned 3Trees model is nested in the larger 3Trees model by imposing equality constraints on the parameters. Therefore, standard techniques for nested model parameters, such as the Likelihood Ratio test, can be safely utilized.

Simulation studies
In this section, we propose two simulation studies to evaluate the algorithms presented in Sect. 3 for fitting 3Trees models.
The first simulation study aims at evaluating the predictive performance of Algorithm 1 (3Trees-Alg). In contrast, the second one is designed to assess the inferential accuracy in terms of confidence intervals of Algorithm 2 (3Trees-Alg-PSI).
We consider three scenarios with different data generating processes. Each scenario includes two individual-level predictors, X 1 and X 2 , having bivariate Normal distribution, N 2 (0, Σ), and two cluster-level predictors, Z 1 and Z 2 , with distribution N 2 (0, Γ ). We set the following values for the covariance matrices The random components are u j ∼ N (0, σ 2 u ), with σ 2 u = 3, and ε i j ∼ N (0, σ 2 ε ), with σ 2 ε = 1. Those values imply a high Intraclass Correlation Coefficient (I CC = 0.75), which allows us to clearly show the relevance of using the random effects in predicting the response. The data generating processes of the three scenarios differ in the shape of the regression function, as described in the following.
Scenario 3: we use the non-linear model of Hajjem et al. (2014), except that we replace the third and fourth level 1 predictors (X 3 and X 4 ) with level 2 predictors, denoted as Z 1 and Z 2 .
In each scenario, some parameters of the linear part are set to zero to provide a glimpse of the ability of the proposed methods to correctly identify predictors with no linear effects in this small dimensional setting.

Predictive performance
For each scenario described in the previous section, we generate 500 data sets with J = 500 clusters of n j = 100 units, for a total sample size of 50 000 units. We randomly split each cluster into two equal parts to form the train and the test sets, each composed of J = 500 clusters of n j = 50 units. We compare the predictive performance of Algorithm 1 (3Trees-Alg) for fitting the three-tree models with the following commonly used methods: RE-EM trees (Sela and Simonoff 2012), regression trees (CART) (Breiman et al. 1984), and a linear mixedeffect model with main effects only (Linear-ME) fitted with maximum likelihood. In addition, we compare the performance of mixed effect random forest as proposed by Hajjem et al. (2014). As a benchmark, we also report the results for a mixed-effect model specified as the data generating model (True-ME). Notice that in Scenario 1, the True-ME model and the Linear-ME model coincide.
We implement the 3Trees-Alg in R (R Core Team 2020), using the packages rpart (Therneau and Atkinson 2019) for the tree component and lme4 (Bates et al. 2015b) for the mixed-effect linear component. As tuning parameters for the tree component we set rpart maxdepth tuning parameters at 2 for the first two scenarions, and at 3 for the third one. The rpart complexity parameter cp is set at 0.0001. The rpart package has also been used for CART, while lme4 has also been used to fit linear mixed-effect models Linear-ME and True-ME. The RE-EM algorithm has been implemented using the R package REEMtree (Sela and Simonoff 2021) with the complexity parameter for pruning at 0.01 and the number of standard errors used in pruning set at 1. For the mixed-effect random forest, we used two different R functions. The first function, MixRF, that we used with the default setting, is implemented in the package MixRF (Wang et al. 2016). The second function, MERF, is implemented in the package LongituRF (Capitaine et al. 2021). We used the default setting for this function, but without the stochastic longitudinal component, not adequate for the scenarios considered in this section.
In clustered data, we can distinguish two types of predicted values for the response variable: (i) the value of the response for a statistical unit of a hypothetical (new) cluster and (ii) the value of the response for a statistical unit belonging to a cluster already included in the sample. The two types of prediction differ in the value assigned to the random effect. In the first type, the random effect is set at its expectation, namely zero, thus predicting the response with the fixed part only. On the other hand, in the second type, the random effect is set at the BLUP prediction for the observed cluster, thus predicting the response with the fixed part plus the BLUP prediction (Robinson 1991;Skrondal and Rabe-Hesketh 2009).
To evaluate the predictive performance of the methods under comparison, we compute the Mean Squared Error on the train data and the Predictive Mean Squared Error on the test data (new data). We denote with MSE and PMSE, respectively, the measures referred to the predictions with the fixed part only. In addition, we denote with clusMSE and clusPMSE the measures when the prediction refers to a new unit in an observed cluster, i.e. including the random effect. Table 1 reports the Monte Carlo averages and standard deviations of the mean squared errors for the methods under comparison. The predictive performance of the proposed 3Trees-Alg algorithm, as measured by clusPMSE and PMSE, is very close to the benchmark of the correctly specified mixed-effect model (True-ME) in Scenario 1 and Scenario 2. As expected, in the highly non-linear model of Scenario 3, the performance of 3Trees-Alg does not reach the benchmark. However, it is still better than the other algorithms, apart from the two mixed-effect random forest algorithms. It is worth noting that the performance of the 3Trees-Alg improves when the depth of the trees increases. As an example, for Scenario 3, we report the performance of 3Trees-Alg with trees of depth 6. This performance is substantially improved but at the cost of a less easy interpretation. An interesting aspect is the comparison of MSE and PMSE, which are pretty similar for the 3Trees algorithm, suggesting that it is never overfitting.
As expected, the error measures for predicting a unit in a new cluster (MSE and PMSE) are higher than the corresponding measures relative to the prediction for a unit of an observed cluster (clusMSE and clusPMSE). The magnitude is much larger (about three times) due to the high proportion of the between-cluster variation in the simulated data (ICC=0.75). It is worth noting that clusMSE and clusPMSE are obtained using BLUP predictions of level 2 errors, which are crucial for obtaining an accurate prediction for a unit of an observed cluster.
The standard deviation of the random effects σ u is interesting per se and, in addition, it plays a key role in the prediction for an observed cluster since it affects the BLUP of the random effect. We compare the performance of the competing methods in estimating the parameter σ u by plotting the Monte Carlo distributions (except for CART, which does not provide an estimate of σ u ). In Scenario 1 (Fig. 1), there are no relevant differences, except for MERF, which seems to underestimate the random effect variance.
On the other hand, Scenarios 2 and 3 are noteworthy, and the results are reported in Fig. 2. In all cases, the proposed method (3Trees-Alg) yields estimates of σ u similar to the benchmark (True-ME), whereas the performance of RE-EM is not satisfactory. Indeed, RE-EM and MixRF tend to overestimate σ u , with a large bias in Scenario 3,   while MERF underestimates it in Scenario 2. The 3Trees algorithm performs better than RE-EM and Linear-ME also in estimating the level 1 standard deviation σ ε , as shown by the plots of the Monte Carlo distributions reported in Appendix A (Figs. 6 and 7).
Another aspect to consider is the estimate of the parameters of the linear component, β 1 , β 2 , γ 1 and γ 2 . Notice that the proposed algorithm is not meant to recover the true data generating process but to predict the response by selecting a quasi-linear model closest to the true data generating process. When the tree component adequately approximates interactions and non-linearities, the parameters of the linear part should be close to the true ones. Table 2 reports the Monte Carlo averages and standard deviations of such estimates for the linear component. In both Scenario 1 and 2, the algorithm 3Trees-Alg provides estimates as good as the benchmark algorithm for null and non-null parameters. In Scenario 3, the parameters for the main effects of the level 1 predictors are correctly captured by the 3Trees-Alg algorithm, together with the null coefficient of the cluster-level predictor Z 2 . Conversely, the coefficient of Z 1 turns out to be underestimated. This behaviour is due to the particular shape of the non-linear component, i.e. the interaction term Z 1 j ln |X 1i j |, which is probably inaccurately accounted for by the tree component. This behaviour confirms that this scenario would probably require a deeper tree to capture such a particular non-linear shape.
The simulation results are similar when we consider a smaller number of clusters, and smaller clusters (J = 50 and n j = 15), with a test set for prediction evaluation of equal size. For these simulations, we also consider the predictive performance of a single regression trunk model for iid data (Dusseldorp et al., 2010) estimated via the R package stima with maxsplit set to 2. Table 7 of Appendix A reports the results for this case for the three scenarios. Simulations indicate that the behaviour of the proposed algorithm with respect to the benchmark is confirmed in terms of average clusPMSE, with a slight increase in Monte Carlo variability. Further settings are considered in Table 8 of Appendix A. Specifically, we report simulation results for five variations of Scenario 2. In Scenario 2A, we set the predictors' covariance at 0.75 instead of 0.4 to check the algorithm's behaviour in the case of highly correlated predictors. Scenario 2B considers the case of independent predictors instead. The 3Trees-Alg algorithm seems not influenced by the correlation between predictors, reporting similar performance in the two scenarios. Scenario 2C considers the case of a lower value of the intraclass correlation coefficient. In this case, the performance remains unchanged in terms of clusPMSE, but it is improved in terms of PMSE due to the reduced impact of the random component. Finally, in Scenarios 2D and 2E, the coefficients of the linear part increase to 5 or decrease to 0.5. Overall, the performance of the 3Trees algorithm is essentially unchanged. To check the behaviour of the procedures with more sample sizes, Table 9 in Appendix A presents the results of the same study as Table 2, where other eight non-significant explanatory variables were added to Scenario 2. The procedure is still able to capture which variables have no linear effect in this quasi-linear scenario.

Inferential accuracy
The Algorithm 3Trees-Alg-PSI (Algorithm 2 in Sect. 3.1) is devised to provide valid confidence intervals. We run a further simulation study applying the 3Trees-Alg-PSI Algorithm to evaluate its performance, using 500 data sets of J = 500 clusters of size n j = 100. For comparison, we compute the confidence intervals for the parameters of interest using a correctly-specified mixed effect model (true model) fitted either on the entire data set (True-ME) and on the same split sample used in the Estimation step of Algorithm 3Trees-Alg-PSI, labelled Split-True-ME. We compare the Monte Carlo averages of the interval length and the Monte Carlo actual coverage. We are reporting here the results for Scenario 2 (Table 3), while those for Scenario 1 and 3 are in Appendix A (Table 10).
Note that the 3Trees-Alg-PSI is directly comparable with Split-True-ME, as they both are fitted on half of the data set. In summary, the length of the confidence intervals obtained with the 3Trees-Alg-PSI is in line with a mixed effect model fitted on the halved sample (Split-True-ME). Instead, the confidence intervals obtained by True-ME are shorter since they are obtained using the whole data set. Therefore, the loss of efficiency is attributable to the reduced sample size and not to the tree selection procedure. Moreover, 3Trees-Alg-PSI provides confidence intervals of the expected coverage, in line with correctly specified mixed effect models, as if the regression function was known. The results for Scenario 1 (Table 10) are in line with those of Scenario 2. As expected, for Scenario 3 the performances are worse. Deeper trees are required for this data generating process complexity, as highlighted by simulations in Table 1.

An illustrative example on invalsi tests in italian schools
We apply a 3Trees model to data from the Italian Invalsi test (see, e.g. Cardone et al. 2019). These data concern students who participated in the Invalsi Math tests while attending 5 th grade in 2013-2014 and then attending the 8 th grade in 2016-2017. The Invalsi data have a multilevel structure where the individual-level units are represented

MATH5
Test score at 5th grade ( by pupils and the cluster-level units by the schools, grouping the students. The data set contains 409 528 pupils in 5 777 Italian schools, with an average number of tested students per school of 70.94. Here the aim is to predict the Math achievement at 8 th grade, having some clues on understanding why results differ. The student and school-level predictors are described in Table 4. We exploit the 3Trees model (4) with pupils as level 1 units and schools as level 2 units, where the random effect part is aimed to capture unobserved factors at the school-level. We first use the 3Trees-Alg (Algorithm 1) to predict the Math score at the 8 th grade. Table 5 compares the performance of the 3Trees model using Algorithm 1 (with tree depths set at 3) and the same algorithms considered in the simulation study. The 3Trees-Alg has the best performance for predicting the math score of a student in an observed school (clusPMSE) and predicting the response of a new student in a new school (PMSE). The RE-EM provides the second-best prediction when cpmin= 0.0001, providing a tree with 55 terminal nodes. The performance when cpmin= 0.001 is worse, but the resulting tree, with 15 terminal nodes, is more interpretable. The CART algorithm, ignoring the two-level structure of the data, performs worse than 3Trees and Linear-ME but avoids overfitting, thanks to the pruning procedure based on cross-validation.
The three methods accounting for clustering show a reduction of the PMSE when using the random effects to predict the response for a new unit in an observed cluster (clusPMSE). For the 3Trees model, this reduction is about 12%, which is much smaller than what was observed in the simulations (Table 1) due to the lower value of the ICC (0.13 versus 0.75).
In order to make inference on the coefficients of the predictors, we apply the 3Trees-Alg-PSI (Algorithm 2), which randomly splits the schools into two independent sub-samples for post-selection inference. Estimating the contextual effects of socio-economic status (SES) and mathematical background (MATH5) of pupils belonging to the same school is relevant for inferential purposes. To this end, we add to the predictors the school averages of these two variables, named gr.M(SES) and gr.M(MATH5). To facilitate the interpretation of the results and keep control of the maximum size of the selected model, we set the maximum tree depth at three. We first describe the tree component in the 3Trees model selected by 3Trees-Alg-PSI. The algorithm represents the predictors' partition of each tree as a factor with M t labels corresponding to the partition regions. The factor is then parameterized as (M t − 1) dummy variables for the usual identifiability constraints. The first region plays as the reference. For simplicity, in the factor description, T stands for Tree and R for Region. The regions selected on the student-level predictors by the first tree are depicted in Fig. 3 and summarized as follows. The second tree, on the school-level predictors, is depicted in Fig. 4 and its selected regions are listed below. The parameter estimates and their confidence intervals are reported in Table 6 in the column labelled 3Trees model. As expected, the Math score at 5 th grade (MATH5) is a key predictor, with a significant main effect in the linear component, in addition to non-linear and interaction effects captured by the first and third trees. Concerning other student-level predictors, it is worth noting that one of the most relevant predictors is the socio-economic status, confirming an enduring issue and the necessity of additional support for children with low socio-economic status in order to overcome social inequalities. Regarding the school-level predictors, the coefficient of type of school (SCTYPE) is negative. Thus, conditionally on the other predictors, students of public schools have a lower performance.
Looking at the confidence intervals, we note that two school-level predictors, namely the number of classes in the school (CLSIZE) and the school's location (TOWN), have confidence intervals crossing zero. In addition, the tree for the schoollevel predictors has two regions (R2.12 and R2.13) whose confidence intervals for the coefficients largely overlap. Similarly, in the third tree, the regions R3.10 and R3.11. Consequently, these two trees can be pruned. Table 6 also reports the parameter estimates obtained from a reduced model where the trees are pruned and the nonsignificant predictors SCSIZE and TOWN are omitted. This reduced model has four parameters less than the full model and is nested in it. The pruned 3Trees model is preferred based on the Likelihood Ratio test comparing the two models (statistic = 4.67, df = 4, p-value = 0.322).

Conclusions
Tree-based regression models are a class of predictive algorithms that are conceptually simple and able to take into account both interactions and non-linearities. However, they struggle with linear dependencies and assume iid data. This paper proposes a multilevel extension of regression trunk models where the tree component consists of multiple trees. In particular, we call 3Trees model a mixed-effect model with a linear part and three trees. The trees are aimed to take into account nonlinearities and interactions among individual-level predictors, cluster-level predictors and cross-levels. The model can be easily extended to a larger number of trees.
The key idea behind our proposal is to develop a flexible class of models capable of accurately approximating the true regression function while preserving interpretability. Indeed, the limitation of a standard regression tree is that, even if it is able to reconstruct a regression function in the presence of non-linear effects and interactions, it does not directly reveal which effect is linear and which is not, which predictor is involved in an Table 6 Estimates and 95% confidence intervals for the 3Trees model using Algorithm 2 and the corre- Level 2 (school) std deviation 6.004 ( 5.830 , 6.186 ) 6.011 ( 5.836 , 6.192 ) Level 1 (pupil) std deviation 14.752 (14.706 , 14.797) 14.752 (14.706 , 14.797) interaction and which is not. As explained by Gottard et al. (2020), in most tree-based algorithms, including CART and random forests, the ordering of the splitting variables and their importance can be different from the ordering and the importance of the direct effects on the response. The interpretation of a tree is in terms of predictive power and not on attribution (Efron 2020) and direct effects as in traditional regression. From this point of view, 3Trees models constitute an interesting bridge between a proper statistical model and a machine learning algorithm.
To estimate the 3Trees model parameters, we introduce two iterative algorithms. The first algorithm is specific for predictive purposes and showed a good predictive performance in the simulation study and the applied study on Invalsi data. It produces transparent and interpretable predictions. However, this algorithm does not provide valid confidence intervals, as the shape of the regression function is learned during the estimation process. We propose a post-selection inference procedure that provides valid confidence intervals for the coefficients and predicted values to overcome this limitation. This procedure is based on the splitting sample procedure, which is easy to apply and shows good performance if the number of clusters in the data is not too small. Notice that, in likelihood-based inference for multilevel models, the number of clusters is crucial for estimating the level 2 variance σ 2 u , which in turn affects the standard errors of regression coefficients, especially those related to level 2 predictors (Elff et al. 2021). In the case of a few clusters, say less than 30 in each split sample, the level 2 variance is appreciably underestimated. Thus the confidence intervals are shorter than due. This issue can be addressed by REML estimation, though this comes at the cost of losing efficiency and it prevents using LR tests for pruning.
The proposed algorithms for 3Trees models can be easily extended to allow for predictors with random slopes (see Appendix B). Indeed, there is no need to modify the algorithm but to specify the random component properly. It is worth noting that the role of a random slope is to account for a large unexplained between-cluster variation in the regression coefficient of a level 1 predictor. This scenario is less likely in a 3Trees model due to the flexibility of the tree component, notably for what concerns crosslevel interactions. It can happen that the increase in complexity induced by random slopes may not be adequately rewarded in terms of predictive power. It will be up to the researcher to decide whether to adopt a random slope or leave cross-level interactions to be explained by the trees. In general, random slopes should be considered only for predictors playing a key role in the research aims.
Finally, further work is needed to apply the proposed 3Trees model and related algorithms to high dimensional settings. When the number of predictors increases, the proposed algorithms can be modified to deal with the high dimensionality issue. For instance, one could replace the likelihood-based estimation procedure with the lassotype estimator of Groll and Tutz (2014). In addition, recently Rügamer et al. (2022) proposed a procedure for post-selection inference for linear mixed-effect models and additive models, implementing a multi-stage selection procedure. Further research is needed to evaluate if their proposal can be adapted to the case of 3Trees models.
Funding Open access funding provided by Università degli Studi di Firenze within the CRUI-CARE Agreement.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.    Fig. 10 Third tree of the 3Trees model on Invalsi data using Algorithm 2 with random slope on SES algorithm is similar to the one for random intercept, but, both in the selection and in the estimation step, the specification of the random component in the mixed-effect model has been modified to include a random slope. In the following, we first describe the selected trees and then report the estimates table.
The regions selected on the student-level predictors by the first tree are depicted in Fig. 8 and listed below. Comparing this tree with the corresponding tree in the random intercept model (Fig. 3), we note that it involves the same predictors, namely MATH5 and SES. However, SES plays a less important role in the tree since it acts only when MATH ≥ 73. This is an expected consequence given that in the current model the effect of SES is modelled in a finer way by the random slope. The second tree, dealing with the school-level predictors, is depicted in Fig. 9 and its selected regions are summarized as follows. Estimates and confidence intervals are reported in the first column of Table 11 (3Trees model). The estimates for the pruned version are presented in the second column (Pruned 3Trees model). According to the LR test (statistic = 4.543, df = 3, p-value = 0.208) the pruned version is to be preferred.