Generalized Estimating Equations Boosting (GEEB) machine for correlated data

Rapid development in data science enables machine learning and artificial intelligence to be the most popular research tools across various disciplines. While numerous articles have shown decent predictive ability, little research has examined the impact of complex correlated data. We aim to develop a more accurate model under repeated measures or hierarchical data structures. Therefore, this study proposes a novel algo‑ rithm, the Generalized Estimating Equations Boosting (GEEB) machine, to integrate the gradient boosting technique into the benchmark statistical approach that deals with the correlated data, the generalized Estimating Equations (GEE). Unlike the pre‑ vious gradient boosting utilizing all input features, we randomly select some input features when building the model to reduce predictive errors. The simulation study evaluates the predictive performance of the GEEB, GEE, eXtreme Gradient Boosting (XGBoost), and Support Vector Machine (SVM) across several hierarchical structures with different sample sizes. Results suggest that the new strategy GEEB outperforms the GEE and demonstrates superior predictive accuracy than the SVM and XGBoost in most situations. An application to a real‑world dataset, the Forest Fire Data, also revealed that the GEEB reduced mean squared errors by 4.5% to 25% compared to GEE, XGBoost, and SVM. This research also provides a freely available R function that could implement the GEEB machine effortlessly for longitudinal or hierarchical data.


Introduction
Correlated data measure the dependent variable repeatedly across multiple dimensions, such as longitudinal, clustered, spatial, or multilevel data [17].Correlated data frequently occur in medicine, public health, and other research fields, requiring specialized statistical approaches to handle the complex correlation structure, avoid potential estimation biases, and ensure the accuracy of estimations.
Generalized Estimating Equations, also known as GEE [6,12,13], is a statistical method initially proposed by Liang and Zeger [11].It extends the framework of Generalized Linear Models (GLMs) and overcomes the assumption of independence among observations, making it particularly useful for handling correlated data.One of the strengths of GEE is that it only assumes a "working correlation matrix" to describe the correlation structure among observations.This characteristic reduces the need for such restricted distribution assumptions and keeps parameter estimation consistent even when the working correlation matrix is misspecified under mild regularity conditions [7].The Mixed-effects model is another standard statistical method for correlated data [10,15].The difference between GEE and the mixed-effects model is that GEE estimates the population average effects, and the mixed-effects model estimates individual random effects.
The widespread adoption of modern technology and digitization has created vast amounts of data in recent decades.Coupled with the general availability of digital services and advancements in storage technologies, a massive accumulation of data has occurred.This phenomenon has driven the necessity for big data analysis, which has fueled the flourishing development of machine learning (ML).ML aims to construct computer models that could automatically optimize algorithms based on past experiences to predict future outcomes.Conceptually, one can think of ML as having numerous settings of candidate models and using a large amount of past experiential data to guide the computer in finding the model setting that optimizes performance indicators [9].
Currently, some popular supervised ML models include eXtreme Gradient Boosting (XGBoost) [3], Random Forest [8], and Support Vector Machine (SVM) [4].As one of the fastest-growing technologies, data scientists applied ML in various fields such as finance, marketing, computer vision, aerospace, biomedicine, etc.Every discipline is increasingly utilizing ML for prediction and decision support.Over the past two decades, ML has made significant progress and achievements, from academic research to the most popular commercial applications [16].
Thus, in the era of big data, in addition to traditional statistical methods, ML has provided advanced choices for data analysis.Numerous studies have compared statistical methods to ML, with some articles indicating that ML outperforms statistical methods [2,14,19].However, limited research discusses more complex data structures, such as correlated or hierarchical ones.Therefore, we propose a novel algorithm, the Generalized Estimating Equations Boosting (GEEB) machine, to integrate the gradient boosting technique from ML algorithms into the GEE.Under such hybrid algorithms, we aim to create a new ML model to deal with correlated data, avoid biased estimates, and provide a more accurate prediction.

GEE
The core of the new machine is the GEE.Here, we briefly introduce the fundamentals of the GEE.Assume that y ij , i = 1tok and j = 1ton i represent the jth response of the ith subject, which has a vector of covariates x ij .There are n i measurements on subject i, and the maximum number of measurements per subject is T. Let the responses of the ith subject be y i = [y i1 , . . ., y ini ]′ with corresponding meansµ i = [µ i1 , . . ., µ ini ]′.
The marginal mean µ ij of the response y ij is related to a linear predictor through a link function g(µ ij ) = x ij ′β , and the variance of y ij depends on the mean through a variance function ν(µ ij ) for generalized linear models (GLM).
Solving the generalized estimating equations, we could obtain the estimate of the parameter: where V i is the working covariance matrix of Y i We only require the mean and the covariance of Y i in the GEE method, we do not need the full specification of the joint distribution of the correlated responses.This feature of the GEE is desirable and leads to a convenient way of analysis since the joint distribution for noncontinuous outcome variables involves high-order associations and is complicated to specify.In addition, the regression parameter estimates are consistent even when the working covariance is incorrectly specified.However, the GEE approach can lead to biased estimates when missing responses depend on previous responses.The "Weighted Generalized Estimating Equations under the MAR Assumption" can provide an unbiased estimate.
Working correlation matrix Suppose R i (α) is an n i × n i "working" correlation matrix specified by the vector of parameters.The covariance matrix of Y i is modeled as: where A i is a diagonal matrix ( n i × n i ) whose jth diagonal element is ν(µ ij ) and W i is a diagonal matrix ( n i × n i ) whose jth diagonal is w ij , where w ij is a variable indicating the weight.If not weighted, w ij = 1 for all i and j.If R i (α) is the true correlation matrix of Y i , then V i is the true covariance matrix of Y i .
In practice, the working correlation matrix is usually unknown, which must be estimated in the iterative fitting process by using the current value of the parameter vector β to compute appropriate functions of the Pearson residual: If the working correlation matrix is the identity matrix (I), the GEE reduces to the independence estimating equations.The table from SAS [20] demonstrates the working correlation structure [21].

GEEB
GEEB has a hybrid design with GEE and gradient boosting.The strength of the GEEB algorithm lies in its ability to handle complex relationships in correlated data while optimizing the algorithm further through the gradient-boosting technique.
Gradient Boosting is a prevailing ML algorithm that applies the idea of gradient descent to ensemble learners.In the gradient boosting framework, each iteration builds a new i learner based on the prediction errors calculated by the learner in the previous iteration.When the iterations meet the stopping rule, all the predictions from the iterations are weighted and summed to obtain the final prediction.More specifically, because the loss function is the difference between the predicted and actual values, the goal of gradient boosting is to progressively move towards the minimum value of the loss function by minimizing the prediction errors in each iteration.Negative gradients of the loss function computed at the previous iteration and learning rate determine the direction and range of progressive movement in the current iteration.In other words, the algorithm minimizes the loss function and improves the overall prediction accuracy by updating the model based on the negative gradients.
The GEEB algorithm has four components: an initial setting and three computational steps.The initial stage defines the input dataset and the loss function.The input dataset contains n samples with some input features ( x i ) and a continuous output fea- ture ( y i ), represented as . When the dependent variable ( y ) is continu- ous, the algorithm defines the loss function as a modified version of the mean squared error: Here, L(, ) represents the loss function, y i denotes the actual outcome of the i th data point, and F (x i ) represents the model's predicted outcome for the i th data point.Modifying the loss function is crucial as it facilitates more straightfor- ward computation of gradients in subsequent steps.
After defining the input dataset and the loss function, the first step is to compute the initial prediction values of the model.The initial prediction value is a constant number chosen to start the iterations at the most efficient point.In this case, the initial prediction value, denoted as F 0 (x) , is defined as F 0 (x) = argmin F (x i ) n i=1 L(y i , F (x i )) .Thus, when the loss func- tion is defined as L y i , F (x i ) = 1 2 (y i − F (x i )) 2 , the initial prediction value is set to the mean of the feature values F 0 (x) = n i=1 y i n .
The second step of the algorithm involves the iterative model updates for M times.Each iteration consists of five parts: (A), (B), (C), (D), and (E).The iterations continue until the M th iteration converges.Unlike the conventional gradient boosting machines that incorporate all input features, the GEEB randomly selects some input features when building the model to reduce predictive errors.In part (A), a subset of the data is created from the dataset by randomly selecting some features to generate the model.This subset is denoted as , where x i ′ represents the selected features and y i represents the target feature.Part (B) involves calculating the residuals between the true and predicted outcomes of the subset (A).The residuals are computed as In part (C), the residuals calculated in (B) are used to fit the generalized estimating equations, obtaining the coefficients for this iteration.Part (D) uses the coefficients obtained in (C) to predict the residuals for the entire dataset in this iteration.Finally, in part (E), the progress of this iteration's prediction for the model is updated.The predicted values of the residuals for this iteration, denoted asp i,m , are multiplied by the learning rate ( ν ) and added to the previous overall predicted valueF m−1 (x i ) .The resulting computation represents the pre- diction for this iteration,F m (x i ).
After M iterations, the third step is to output the overall prediction results of the model, denoted as F M (x i ).
The following presents the GEEB algorithm.

Materials
This chapter describes the data source, the detailed background of simulation studies, and the application to a real-world dataset, the Forest Fire Data.

Data source
The Institute of Digital Research and Education (IDRE) at the University of California, Los Angeles, published the HDP simulated data in July 2012 [1].The HDP data is based on a large-scale lung cancer-related study.The correlations exist in its hierarchical structure, which consists of three nested levels: Doctors are nested within hospitals, and patients are nested within doctors.Researchers could adjust the number of hospitals, doctors, and patients according to their research requirements.
The simulated data includes nine different outcomes.For this research, we select tumor size, which follows a Gaussian distribution, as the target output feature.The patient-related features include age (Age), marital status (Married), family history (Fami-lyHx), smoking history (SmokingHX), sex (Sex), cancer stage (CancerStage), length of stay in hospital (LengthofStay), white blood cell count (WBC), red blood cell count (RBC), body mass index (BMI), interleukin-6 (IL6), and C-reactive protein (CRP).At the doctor level, there is doctor ID (DID), the experience of the doctor (Experience), the quality of the school doctors trained (School), and the number of lawsuits (Lawsuits).
Note that the variable "School" is divided into two categories (top vs. average).Due to the highly imbalanced distribution, the "school" variable may have only one group that introduces errors in estimating the GEE function with the R package.Therefore, we did not include the "school" variable in simulation studies.The hospital-related features include hospital ID (HID) and Medicaid at the given hospital (Medicaid).Consequently, there are 17 predictors in the simulation study.Note that not all 17 features are related to the target response, and these features are noises to the predictive models.

Real-world data
Cortez and Morais [5] published the Forest Fire Data.This dataset covers the period from January 2000 to December 2003 and includes records of forest fires in the Montesinho Natural Park in northeastern Portugal.Multiple institutions collected the data and encompassed numerous variables, such as the Fire Weather Index (FWI) [22], spatial, temporal, and weather-related information.
We generated a new feature for "season" to construct the third-level hierarchical structure.In this way, the day is nested within the month, and the month is nested within season.Note that "season" was derived from the "month" variable.The four seasons are (1) Spring, from March to May; (2) Summer, from June to August; (3) Autumn, from September to November; and (4) Winter, from December to February.As a result, there are 14 variables (refer to Table 1), and the sample size is 517.

Experiment
We examine the consistency and accuracy of (1) the new machine GEEB, (2) the statistical method GEE, (3) the SVM, and (4) XGBoost under different hierarchical structures and sample sizes.Regarding the hyperparameter of SVM, the kernel is a radial basis function (RBF).Hyperparameters of XGB are objective = 'reg:squarederror' , nrounds = 50, and verbose = 0. Other settings of hyperparameters yielded similar results.When developing simulation studies, we tuned the SVM and XGB with different parameter settings, such as the max_depth and learning_rate for the XGB.The results could be better or worse.In each scenario, there are 1000 repetitions.Each repetition could find its best parameter setting, but the comparisons are similar.Therefore, we used the most common settings for SVM and XGB.

Simulation studies
The core concept of the GEEB machine involves the random selection of features.Note that validation sets in the training data could find the optimal proportion.However, we randomly select eight proportions (30%, 40%, 50%, 60%, 70%, 80%, 90%, and 100%) of features in GEEB denoted as Model 1 to Model 8.The new approach is more promising if the GEEB outperforms other methods without an optimal proportion obtained by tenfold cross-validation.Table 2 presents detailed model settings.Next, we examine the impact of sample size on the predictions.Therefore, three sample sizes were defined: (A) small, (B) medium, and (C) large.Due to the random variation ±1 set for the number of doctors and patients, the sample size is approximately estimated as the mean value.The three estimated sizes with the minimum and maximum in the brackets are (A) small sample: 200 [72, 405], (B) medium sample: 500 [200, 720], (C) large sample: 1000 [650, 1500].Because the results are consistent from 72 to 1,500 patients, we did not increase the sample size after 1,500.Since a small sample size introduces more statistical issues than big data, the simulation study suggests that a minimum of 72 subjects is sufficient to implement the GEEB.
Additionally, we explore the impact of different hierarchical data structures and consider five scenarios in three different sample sizes: (1) a structure with a small number of hospitals, followed by some doctors, and then more patients in a ratio of 1:3:5.(2) A structure with a small number of hospitals, followed by more doctors and then a more significant number of patients, in a more disparate ratio of 1:5:9.(3) An equal number of hospitals, doctors, and patients in a balanced ratio 1:1:1.(4) A structure with many hospitals, followed by some doctors, and then a small number of patients in a ratio of 5:3:1.
(5) A structure with even more hospitals, followed by some doctors and a few patients, with a more extreme ratio of 9:5:1.Tables 3, 4, 5 show a detailed summary of the hierarchical structures.
Finally, the research framework diagram in Fig. 1 indicates the study flow.In the beginning, the input datasets undergo data preprocessing, and then the dataset is split into an 80% training set to build the models.The remaining 20% of the data is the testing set that evaluates the predictive performance.Lastly, we record the predictive performance in every scenario.
Regarding other hyperparameters, the learning rate for the GEEB model is set to 0.1.The number of iterations for the GEEB model is set to 100 since the convergence takes many iterations.Lastly, the simulation of the HDP dataset repeats 1000 times for each parameter setting.

Application to a real-world data
In the Forest Fire Data, the preprocessing step standardized all input features.We split the data into 80% training and 20% testing to evaluate the performance.Due to the stochastic nature of data splitting and feature selection, we will repeat the analysis one or 100 times.Tables 13, 14 and Fig. 2 reveal the analysis results and Feature Importance.

Evaluation metric
The Mean Square Error (MSE) measures the performance since the output feature is Gaussian.The MSE is defined as: The formula represents the expected value of the squared difference between the true values ( y i ) and the predicted values ( F (x i ) ) for each subject in a dataset of size n.There- fore, a smaller MSE indicates that the model's overall predicted results are closer to the actual values, meaning better performance, and vice versa.

Results
Computer simulations and applications to the Forest Fire Data are implemented by R version 4.2.1 R Core Team [18].R: A language and environment for statistical computing.R Foundation for Statistical Computing).A computer with 11th Gen Intel(R) Core(TM) i7-11700 @ 2.50 GHz and 16 GB of RAM in a 64-bit platform is used to implement all the experiments.

Simulation results
The new model GEEB performed well according to simulation results.The GEEB, with the selection of all features, is consistently superior to the benchmark GEE in Tables 6, 7, 8.With a suitable random feature selection proportion, the GEEB has a further improvement.
We discovered that three factors, the proportion of random feature selection, sample size, and hierarchical structure, impact model performance.The GEEB and GEE  perform better in A1-A2, B1-B2, and C1-C2 in Tables 6, 7, 8.These scenarios have fewer hospitals, a moderate number of doctors, and more patients.Models within the same hierarchical structure perform better as the dataset increases (B1 vs. A1, C1 vs. B1, …, etc.).Although the optimal feature selection proportion varies with different sample sizes and hierarchical structures, we suggest that Model 3 (with a 50% random feature selection) demonstrates consistent and satisfying predictive results (Tables 6, 7, 8).Therefore, without tenfold cross-validation searching for the optimal ratio, we adopted Model 3 in Tables 9, 10, 11, 12 for the GEEB.The GEEB with Model 3 shows a superior MSE than the SVM and XGBoost (Tables 9, 10, 11, 12).
In addition to the main results mentioned above, the subsequent sections have details in three aspects.First, we explore the impact of different random feature selection proportions in the GEEB.Secondly, we examine the influence of sample size.Lastly, we see how the predictive ability differs among the GEEB, GEE, XGBoost, and SVM.

The proportion of random feature selection
In the small sample scenarios (A1 to A5), Table 6 shows that the GEEB consistently outperforms the GEE even when all features are included (Model 8).Therefore, the boosting technique improved the accuracy compared to the conventional statistical approach.The optimal model is identified as Model = 95.51834122 ).These results demonstrate the impact of the proportion of random feature selection on the GEEB across various sample sizes.
In Tables 6, 7, 8, Model 8 (GEEB with all features) outperforms GEE even without random feature selection, showing better predictive results than the GEE model.Furthermore, in the small sample scenarios (Table 6), the MSE in Models 1 to 7 is smaller than the numbers in Model 8.The comparisons demonstrate the improved GEEB through random feature selection.Depending on the sample size, the optimal random feature selection proportion falls between 30 to 60%.Table 7 (B2-B4) and 8 (C1-C5) show a curved pattern when the percentage is too small.Models 1 or 2 may encounter a higher MSE compared to Model 8 of the GEEB.For example, in B2, the best scenario is Model 2 (40% features, MSE M2,B2 = 93.35856425), whereas a lower selection proportion Model 1 (30% features, MSE M1,B2 = 93.53285239 ) performs worse than Model 2 and even under- performs the GEE ( MSE GEE,B2 = 93.43940338).Similarly, in C1, the best situation is Model 4 (60% features, MSE M4,C1 = 93.71736127), whereas a lower selection proportion Model 1 (30% features, MSE M1,C1 = 94.02394485 ) performs worse than Model 4 and even underperforms the GEE ( MSE GEE,C1 = 93.72922513).We think that the informa- tion content of 30% or 40% of the features is insufficient to provide accurate predictions.
In conclusion, although each sample size has its optimal range of random feature selection proportions, we recommend the following: (1) random feature selection is a hyperparameter for the proposed GEEB machine.The optimal selection proportion falls between 30 to 60% through simulation studies.Hence, one can find the optimal hyperparameter through techniques such as a validation set or k-fold cross-validation.(2) A 50% selection proportion consistently demonstrates stable and excellent performance across all scenarios.In cases where it is impossible to employ any validation technique, this research suggests setting the random feature selection proportion to 50%.Thus, the GEEB function in the R language has a default random feature selection proportion set to 50%.
The following two sections (Tables 9, 10, 11, 12) will explore the impact of sample size and hierarchical structure using the GEEB with Model 3.

Sample size
In Table 9, all models, including the GEEB, GEE, SVM, and XGBoost, perform better as the sample size increases.The GEEB with Model 3 in the small sample A1 has MSE M3,A1 = 100.08083977 .In the medium sample sce- nario B1, it is MSE M3,B1 = 94.62586832 .In the large sample scenario C1, it is MSE M3,C1 = 93.71743381 .The pattern indicates an improvement in reduced errors as the sample size increases.For GEE, the MSE in A1 is MSE GEE,A1 = 100.53386987, in B1, it is MSE GEE,B1 = 94.71798990, and in C1, it is MSE GEE,C1 = 93.72922513 .For SVM, in A1, it has MSE SVM,A1 = 112.12767754, in B1, it has MSE SVM,B1 = 97.14783528, and in C1, it has MSE SVM,C1 = 90.40457748 .For XGBoost, in A1, it has MSE XGBoostA1 = 112.45269936, in B1, it has MSE XGBoost,B1 = 97.25146238, and in C1, it has MSE XGBoost,C1 = 89.93082567 .In summary, under the same hierarchical structure, increasing the sample size leads to a decreasing trend in MSE, indicating improved predictive ability.Besides, SVM and XGBoost are more sensitive to sample size variations than the GEEB and GEE.

Hierarchical structure
Five types of hierarchical structures show the impact on the MSE.In Tables 10, 11, 12, there is a significant contrast between the first (ratio = 1:3:5) and fourth (ratio = 5:3:1), as well as the second (ratio = 1:5:9) and fifth (ratio = 9:5:1) hierarchical structures across all scenarios with the same dataset size.All models perform better in the first and second hierarchical structures, where the setting presents fewer hospitals, a moderate number of doctors, and the highest number of patients.In contrast, they do not perform well in the fourth and fifth hierarchical structures, when the data have more hospitals and minimal patients.Regarding the third hierarchical structure, where the data involves an equal number of hospitals, doctors, and patients, the predictive capability lies between all models.
Note that the second and fifth hierarchical structures demonstrate more extreme ratios, implying a more significant disparity in the size of hospitals, doctors, and patients.These situations investigate whether the models would exhibit more extreme MSEs.However, only in the small sample scenario A could we observe relatively notable differences.We did not see significant differences in the medium and large sample scenarios.The XGBoost and SVM are more sensitive hierarchical structures than the GEEB and GEE.
In the small sample size under the first and fourth structure, the GEEB with Model 3 yields MSE M3,A1 = 100.).Thus, we observe signif- icant differences between the first and fourth scenarios when using the SVM and XGBoost.It may be because the two ML models are not well-suited to handle hierarchical data, leading to their inferior performance in the fourth hierarchical structure, which has relatively enormous clusters.
Furthermore, according to the increasing number of clusters, we observed that SVM and XGBoost show a clear inverse relationship in both medium and large samples.Besides, as the number of clusters increases, the predictive performance of SVM and XGBoost decreases, indicating worse performance with more clusters and more pronounced inter-cluster correlations.In contrast, the GEEB and GEE demonstrate consistent and satisfying predictions.
The SVM and XGBoost can outperform the GEEB in scenarios C1 and C2 for larger sample sizes (Table 12).The reason may be with fewer clusters and relatively large data sizes, SVM and XGBoost can overlook the inter-cluster correlation structure and treat the data as independent.

Results of the Forest Fire Data
According to the simulation study, the GEEB model with 50% random feature selection demonstrates consistent and improved predictive performance.Therefore, we adopt the GEEB with Model 3 as the default model in real-world data analysis.
The Forest Fire Data analyses for each method are shown in Table 13, indicating that the GEEB with Model 3 exhibits the minimum MSE compared to the GEE, SVM, and XGBoost.The MSE of GEE is approximately 4.5% higher than the GEEB.The XGBoost is about 25.2% higher than the GEEB.Therefore, the GEEB has a decent improvement compared to the most famous statistical model for correlated data and the most promising ML approaches, SVM and XGBoost.Feature Importance of the GEEB with Model 3 is in Table 14, and the visualization is in Fig. 2.

Discussions
In this study, we propose a new ML strategy named the Generalized Estimating Equations Boosting (GEEB) machine.This method integrates the gradient boosting technique with the gold standard model for correlated data, the GEE.Computer simulations confirmed that the GEEB outperforms the GEE.In most situations, GEEB performs better than the famous SVM and XGBoost.Besides, the GEEB demonstrates the best prediction for the Forest Fire Data.Therefore, our findings suggest: (1) the gradient boosting technique enables the GEEB to outperform the GEE model.(2) Although the XGBoost and SVM are known for their excellent predictive ability, they may not perform well This research also provides the code that computes all research results.The geebm() is an R function that implements the GEEB machine.This function has seven arguments: formula, id, iteration, feature_rate, lrate, standardize, and data.Note that formula must be specified in the format "response ~ predictors" to list the predictors (input features) and response variable (output feature) in the dataset.id is a vector that identifies the clusters and can support multiple levels arranged in the order of multilayer structure.iteration is an integer representing the number of iterations, set to default at 100 iterations.feature_rate represents the proportion of random feature selection.When set to 1, it uses all features; by default, it is set to 0.5, using half of the features.lrate is a hyperparameter for the learning rate, with a default value of 0.1.standardize determines whether features are standardized, and the default does not perform standardization.data is used to input the training dataset.For example, when training the model with the Forest Fire Data in this study, the function would be: geebm(area ~ X + Y + FFMC + DMC + DC + ISI + temp + RH + wind + rain + day, id = c("season","month"), iteration = 100, feature_ rate = 0.5, lrate = 0.1, standardize = T, data = Dataset).
The GEEB is also inspired by the Random Forest that incorporates Bootstrap while randomly selecting features.However, the results were not satisfying.Our research aims to compare the GEEB with other benchmark ML and statistical models in correlated data.When deciding which ML models to include, we primarily considered models that are widely discussed and used in academia and industry and frequently win in various data science competitions.Therefore, we included the XGBoost, SVM, and Random Forest.However, when conducting the simulation studies, we discovered that the random-Forest package in R cannot handle datasets with more than 53 categories.Since each doctor within each hospital is treated as a separate category and there are other categorical features such as gender and cancer stage, this number exceeds the limitation.Therefore, we must exclude Random Forest in the comparison as it does not apply to the hierarchical dataset.
Integrating the concept of gradient boosting and using statistical model GEE as the base learner, combined with a random feature selection, the proposed novel approach GEEB has several advantages.Compared to the ML model XGBoost, which also utilizes Gradient Boosting, GEEB performs better in most scenarios.GEEB can handle such data more effectively, resulting in improved predictive performance.Furthermore, compared to using the GEE model alone, after conducting 1000 simulations, we observed that GEEB achieves more accurate predictions.

Limitations
There are some limitations in this study.Firstly, the simulated data used in this study is based on the publicly available HDP dataset from UCLA, and the investigation of the impact of the level of variable correlations, such as weak to high correlations, has not been further explored.
Secondly, in this study, we investigate the predictions of tumor size.The underlining techniques of GEEB are GEE and gradient boosting, both of which support classification tasks.However, this research focused on regression tasks only.The performance of the GEEB under other types of output features is unknown.
Here, we only roughly categorized the structures into three types: (1) fewer hospitals, followed by some doctors and more patients; (2) more hospitals, followed by some doctors and fewer patients; and (3) an equal number of hospitals, doctors, and patients.The study also considered varying sample sizes, including more extreme cases.Consequently, we examined five hierarchical structures: 1:3:5, 1:5:9, 1:1:1, 5:3:1, and 9:5:1.While this design provides initial insights, we could explore more detailed hierarchical structures in future works.

Future research topics
The corresponding theoretical work and simulation studies are great topics for future research with dichotomous, ordinal, or categorical nominal correlated datasets.
learning rate, with a default value of 0.1.standardize determines whether features are standardized, and the default does not perform standardization.data is used to input the training dataset.For example, when training the model with the Forest Fire Data in this study, the function would be: geebm(area~X+Y+FFMC+DMC+DC+ISI+temp+RH+ wind+rain+day, id=c("season", "month"), iteration=100, feature_rate=0.5,lrate=0.1,standardize=T, data=Dataset).

Table 1
The preprocessed Forest Fire Data attributes

Table 2
Settings of Model 1-Model 8

Table 3
Parameter settings of a small sample size

Table 4
Parameter settings of a medium sample size

Table 5
Parameter settings of a large sample size

Table 6
Simulation results of GEEB and GEE with a small sample size

Table 7
Simulation results of GEEB and GEE with a medium sample size

Table 8
Simulation results of GEEB and GEE with a large sample size

Table 9
Simulation results of GEEB with Model 3, GEE, XGBoost, and SVM concerning the hierarchical structure

Table 10
Simulation results of GEEB with Model 3, GEE, SVM, and XGBoost with a small sample size

Table 11
Simulation results of GEEB with Model 3, GEE, SVM, and XGBoost with a medium sample size

Table 13
Applications to the Forest Fire Data with hierarchical data.Treating subjects as independent failed to capture the correlation structure.

Table 14
Feature importance of the Forest Fire Data (single run/averaged 100 runs)