# Random Forests in Count Data Modelling: An Analysis of the Influence of Data Features and Overdispersion on Regression Performance

Authors: Ciza Arsène Mushagalusa; Adandé Belarmain Fandohan; Romain Glèlè Kakaï
Journal: Journal of Probability and Statistics
 (2022)

---

## Abstract

Machine learning algorithms, especially random forests (RFs), have become an integrated part of the modern scientific methodology and represent an efficient alternative to conventional parametric algorithms. This study aimed to assess the influence of data features and overdispersion on RF regression performance. We assessed the effect of types of predictors (100, 75, 50, and 20% continuous, and 100% categorical), the number of predictors (p = 816 and 24), and the sample size (N = 50, 250, and 1250) on RF parameter settings. We also compared RF performance to that of classical generalized linear models (Poisson, negative binomial, and zero-inflated Poisson) and the linear model applied to log-transformed data. Two real datasets were analysed to demonstrate the usefulness of RF for overdispersed data modelling. Goodness-of-fit statistics such as root mean square error (RMSE) and biases were used to determine RF accuracy and validity. Results revealed that the number of variables to be randomly selected for each split, the proportion of samples to train the model, the minimal number of samples within each terminal node, and RF regression performance are not influenced by the sample size, number, and type of predictors. However, the ratio of observations to the number of predictors affects the stability of the best RF parameters. RF performs well for all types of covariates and different levels of dispersion. The magnitude of dispersion does not significantly influence RF predictive validity. In contrast, its predictive accuracy is significantly influenced by the magnitude of dispersion in the response variable, conditional on the explanatory variables. RF has performed almost as well as the models of the classical Poisson family in the presence of overdispersion. Given RF’s advantages, it is an appropriate statistical alternative for counting data.


---

## Body

## 1. Introduction

Interest in modelling count data has increased significantly over the past two decades. The Poisson distribution is still the most widely used distribution for modelling count data in many research areas, despite the violation of its well-known property that the mean and variance are equal, conditional on explanatory variables. However, the Poisson distribution is inadequate to model overdispersed count data [1, 2]. As a result, a multitude of alternative count data models have been developed to address the shortcomings of the Poisson regression in the presence of overdispersion, including quasi-Poisson [3], generalized Poisson [4], negative binomial [5–7], zero-inflated Poisson [8], and zero-truncated Poisson models [9].Nevertheless, traditional count data models remain limited as they may not be capable of detecting the presence of complex nonlinear interactions between explanatory variables and outcomes [10]. In addition, traditional statistical models used to predict overdispersed data from selected predictor variables are not suitable for “small-n-large-p” problems [11].The random forest (RF) is a good machine learning approach to improving prediction accuracy and model interpretation [12, 13]. RF is a “data-driven statistical method.” It is an ensemble learning approach developed to increase classification accuracy and regression tree prediction by combining many decision trees [14]. The algorithm takes advantage of two powerful techniques: random subspace selection at each split (“classification and regression trees (CART)-split criterion”) [15] and bagging (an abbreviation of “bootstrap aggregating”) of unpruned decision tree learners [16]. The nonlinear nature of RF gives it an edge over linear algorithms [14, 17]. RF algorithms provide variable importance measures for variable selection purposes. Moreover, they involve complex high-order interaction effects and are user-friendly because they have few parameters to set and are less demanding in preprocessing [18]. Over the past 25 years, RF has been successfully applied to a wide range of prediction problems, mainly for classification tasks [19–21]. RF has become an effective data analysis tool that performs well compared to many standard methods [22].However, in a regression problem, the range of predictions that RF can make is limited by the highest and lowest values in the training dataset. When training and testing data ranges differ, such as when predicting overdispersed response variables, this behaviour becomes troublesome [23]. Overdispersion is one of the important constraints in statistical analysis and may lead a variable to appear as a significant predictor when it is not, by deflating or underestimating the standard errors of the estimates [24]. As a result, it is essential to fully understand the data’s nature. Furthermore, most of the existing machine learning methods for regression tasks are designed for continuous data, and their performance for counted data is not well known [13].The primary motivation for this study is thus to evaluate the effectiveness of computer-based statistical approaches known as machine learning, such as RF, in the presence of overdispersed response variables. Since the performance of an algorithm depends mainly on the characteristics of datasets to be analysed (type of predictors, number of attributes in the dataset, sample size, noise, and dispersion), high performance requires the identification of the most appropriate algorithm for a given problem and dataset [25]. RFs are considered data type-specific, insensitive to unimportant and noisy predictors, and dependent on the number of informative variables [18, 19, 26].The RF method has shown comparable or even better prediction performance than other learning methods, such as neural networks, partial least squares regression, and support vector machines [14, 19, 22]. Unfortunately, the impact of overdispersion on the RF performance for different data features has rarely been assessed [27–29]. Thus, determining the impact of data features on RF parameters’ setting and accuracy is an avenue to explore as efforts are needed to improve the accuracy of RF estimates [12]. Therefore, this study aimed to assess the influence of data features on RF regression performance and compare RF regression to count data classic regressions. In addition, datasets with varying magnitudes of overdispersion from two real case studies are also analysed to examine RF performance when the response variable is overdispersed.


## 2. Materials and Methods

### 2.1. Statistical Models

#### 2.1.1. Random Forest Technique (RF)

RF methods deal with both supervised classification and regression tasks [14]. This research was focused on regression analysis.(i) Basic Principles. Assuming that we are given a training sample Dn=X1,Y1;…;Xn,Yn of i.i.d. 0,1d×R-valued random variables d≥2 with the same distribution as the independent prototype pair D defined by X,Y of size n satisfying the EY2<∞ condition. The space 0,1d is equipped with the standard Euclidean metric. For a given x∈0,1d, the regression function estimate mx can be obtained by(1)mx=EYX=x.Regarding this, the regression function estimatemn for the dataset Dn is assumed consistent if(2)EmnX−mX2⟶0,asn⟶∞.In RF, the final predictor is an average ofM independently randomized regression trees. For a given tree j randomly generated, the prediction value at the query point x is denoted mnX;Θj,Dn, where Θ1,Θ2,…,ΘM are the outputs of a randomized variable Θ (also known as independent random variables).In practice, the randomized variableΘ is used to resample the training set used in successive individual trees building [18]. The finite forest estimate (combination of the M trees grown) is obtained by(3)mM,nX;Θ1,…,ΘM,Dn=1M∑j=1MmnX;Θj,Dn.In the randomForest package, the default value ofM is 500. Since M may be chosen arbitrarily large (“limited only by available computing resources”), this estimate can be generalized from a modelling point of view, letting M tend to infinity by an infinite forest estimate [15], obtained by(4)m∞,nX;Dn=EΘmnX;Θ,Dn.It is noteworthy in practice that the above expectation is evaluated by Monte Carlo, i.e., by generatingM (usually large) random trees and taking the average of individual outcomes. This procedure is justified by the law of large numbers [14, 30] which asserts that the almost surely conditional on Dn is expressed as follows:(5)limM⟶∞mM,nX;Θ1,…,ΘM,Dn=m∞,nX;Dn.(ii) Algorithm. The whole process of the RF algorithm can be summarized as follows:LetM be the number of trees to grow in the forest and mtry be the number of variables to select at each node. The following steps are considered:(i)
We resample the training set to create new bootstrap samples.(ii)
For each bootstrap sample, we create an unpruned tree that is grown until all nodes contain observations no more than the maximal terminal node size defined by the user (a prespecified parameter).(iii)
At each split node, a randomly selected subset of predictors is uniformly chosen. In regression tasks forp predictors, a value equal to p/3 is utilized as the default value in the R package randomForest. Then, the best split using only these predictors is determined or returned. In addition, all observations are not used in tree construction, and a subset an is picked from the whole sample Dn. If an=n, the algorithm runs in a bootstrap mode, whereas an<n corresponds to subsampling (with replacement).(iv)
The RF overall prediction value is the average (in regression) of the results obtained from allM given by individual trees [15].(iii) Parameter Tuning. To date, there are limited studies focusing on RF parameter tuning [31]. According to Probst and Boulesteix [32], the most important parameters in RF methods are as follows:(i)
Number of variables to possibly split in each node(ii)
Minimal terminal node size(iii)
Fraction of observations to sample in each tree(iv)
Number of trees


#### 2.1.2. Poisson Regression Extensions Used to Predict Count Data

The prediction performance of RF was compared with that of the Poisson model and its extensions. It is suggested that they be used for assessing count data. The following models are utilized in this investigation and are briefly described as follows:(i) Linear Model. It is a linear function of the form:(6)Y=Xβ+ε;ε∼N0,σε2,where Y is the vectorized form of the response variable, X is a covariate matrix, β is a vector of coefficients, and ε is a vector of normally distributed errors with 0 mean and a variance of σε2 [33].(ii) Poisson Model (Poisson). It is a special case of the generalized linear model (GLM) with a log link function. It is the standard GLM for count data [34]. We assume that the response variable yi follows the Poisson distribution with the mean μi. Its probability mass function is written as(7)PY=yi;μi=e−μiμiyiyi!=expyilogμi−μi−logyi!,where the conditional mean is obtained by μi=Eyi/xi=expx′β with x′ being the vector of covariates and β being the vector of unknown parameters [35–37].(iii) Negative Binomial Model (Neg.Bin). In the negative binomial, a parameter of dispersion K is introduced to accommodate unobserved heterogeneity in the count data. Its probability density function is written as(8)fyixi=τyi+Kτyi+1τKKK+μiKμiK+μiyi,where K is a constant, τ is a function of yi, μi=expx′β, x′ is the vector of covariates, and β is the vector of unknown parameters [38].(iv) Quasi-Poisson. The generalized linear model is used to define the quasi-Poisson model. Let Y be a random variable such that(9)EY=μ,VarY=ϕ×μ,where Var (Y) is Y’s variance, μ>0, ϕ is the dispersion parameter and ϕ>0, and EY is Y’s expectation.(v) Zero-Inflated Poisson model (ZIP). This model is appropriate when there are zeros in excess. Its probability mass function is written as [34, 36](10)PY=yixi,zi=wizi+1−wiziPoissonμi;0xi,ifyi=0,1−wiPoissonμi;yixi,ifyi>0,where zi is the vector of covariates, wi is the probability that the response value is zero, Poissonμi;0xi=exp−μi, and Poissonμi;yixi=exp−μiμiyi/yi!.


### 2.2. Simulation Studies

#### 2.2.1. Effect of Dataset Characteristics on RF Algorithm Parameters’ Setting and Accuracy

Three sample sizes (N=50,250, and 1250), three numbers of predictors (p=8, 16, and 24), and five types of predictors (100, 75, 50, and 25% continuous, and 100% categorical) were considered. Schönbrodt and Perugini [39] stated that n should approach 250 for stable estimates. Biau [18] showed that n>500 is sufficiently large, while 50 represents a small population size. Altogether, 45 data types were set. For each case, 1000 datasets were created. Three RF parameters were considered to determine the best combination, while the number of trees was maintained constant (500):(i)
Number of variables to be randomly selected at each split (mtry): mtry= 2, 3, …, p, where p= the number of predictors(ii)
Proportion of samples for model training (sample size) with four levels (55,63.2,70, and 80%)(iii)
Minimal number of samples within each terminal node (node size= 2, 3, 4, …, 9)Depending on the number of predictors, the combination of different parameters under study resulted in 224, 480, and 736 different scenarios for each data type. Each scenario was run 1000 times to find optimal parameter combinations. For all datasets, the generated explanatory variables were not correlated and were not informative. The dependent variable was generated randomly from the Poisson distribution with lambda equal to one, which is denoted asPoissonλ=1. Categorical predictors were randomly generated by varying the number of levels from 2 to 8. Uncorrelated quantitative covariates were generated from a multivariate normal distribution (MVN) y∼Nn0,1, whose density function in the case of n dimensions is written as(11)fyiμ,Σ=fyi1,yi2,…,yin′μ=μ1,…,μn′,Σ=σ12σ1n⋱σn1σn2=12πnΣexp−12yi−μ′Σ−1yi−μ.The symbolΣ refers to the determinant of the matrix Σ. The matrix Σ must be positive semidefinite to assure that the most likely point is μ=μ1,μ2,…,μn and that, as yi moves away from μ in any direction, then the probability of observing yi declines. The denominator in the formula for the MVN is a normalizing constant, which assures that the distribution integrates to 1.0 [40].


#### 2.2.2. Effect of Overdispersion on RF Parameters’ Setting and Accuracy

In order to assess the effect of overdispersion on the setting and performance of RF parameters, negative binomial-distributedYnb and quasi-Poisson-approximated Yqp variables were randomly generated. The variance-to-mean relationship of Ynb is quadratic, while the variance-to-mean relationship of Yqp is linear such as(12)EYnb=μ,EYqp=μ,VarYnb=μ+ϕ×μ2,VarYqp=μ+ϕ×μ,where μ>0, ϕ is the dispersion parameter, and ϕ>0. For this objective, μ=3 with varying ϕ=1, 3, and 5. Three different samples (N=50,250, and 1250) and five types of predictors as in the first simulation design were considered. Figure 1 shows the shape and the distribution of randomly generated overdispersed outcomes.Figure 1
Histograms and density plots of some randomly generated overdispersed outcomes.


#### 2.2.3. RF Performance Compared to Poisson Model Extensions

Negative binomial-distributed variablesYi were generated with the linear and quadratic variance-to-mean relationships. We let the mean μ vary as a function of two informative covariates m=2 and different levels of the dispersion parameter ϕ. The mean of the outcome Yi of 1000 observations was obtained according to the model logEYim=y|Xim=logμi=α+∑m−12βmXim, where β is a collection of parameters β1=1,β2=0.5 and α=0.145. In addition, three noninformative variables were added from a multivariate normal distribution (MVN) ∼Nn0,1.All datasets were divided into two subsets: training and validation. 70% of the data were randomly selected for model calibration and the remaining 30% for model validation.Poisson regression, negative binomial, zero-inflated Poisson, and zero-inflated negative binomial were compared to the RF regression technique and linear regression applied to the log-transformed response variablelogY+1, as it is considered to be one of the alternatives when dealing with overdispersion or count data [41, 42]. The analysis of variance (ANOVA) test and Tukey’s honest significant difference test (Tukey HSD) were conducted to compare the performance of methods.


#### 2.2.4. Criteria for Methods’ Accuracy Assessment

The root mean square errorRMSE and biases were computed for each model. These evaluation criteria are calculated as follows:(13)RMSE=1n∑i=1nPi−Oi2,Bias=1n∑i=1nEPi−Oi,where Pi and Oi are the predicted and observed values in the subunit i, respectively, n is the number of samples (subunits) [43], and EPi is the average of the predicted values [44]. The RMSE represents the overall quality of the prediction. Predictions become increasingly optimal as the RMSE approaches zero [43, 45]. All analyses were computed using R 4.0.3 software [46].


### 2.3. Motivating Real Datasets

We used two datasets from real case studies to examine RF performance at two magnitudes of overdispersion caused by outliers or zero inflation. The first example is tick burdens observed on two local goat breeds. Ticks were counted monthly between birth and the age of eight months on ten billy goats belonging to two local breeds [47]. The results are shown in Table 1. The dispersion statistic for this dataset is 3.050, caused by outliers.Table 1
Average tick counts on two local goat breeds.
Age (months)Breed 1Breed 2Billy 1Billy 2Billy 3Billy 4Billy 5Billy 6Billy 7Billy 8Billy 9Billy 10113414344222336187686634371131541014845716947128618516896191319391217613431169161614271241181810324379817312941112361018The second example is the dataset “NMES1988”: Demand for Medical Care in the National Medical Expenditure Survey (NMES) conducted in 1987 and 1988 [48]. This dataset contains 4406 observations on 19 variables. In this example, we examined the relationship between the number of physician office visits (response variable) and health (factor indicating self-perceived health status), age in years (divided by 10), gender, married (factor indicating that the individual is married), income (family income in USD 10000), and insurance (factor indicating that the individual is covered by private insurance) as predictors. The dispersion statistic for this subset is 7.345, caused by outliers and zero inflation.


## 2.1. Statistical Models

### 2.1.1. Random Forest Technique (RF)

RF methods deal with both supervised classification and regression tasks [14]. This research was focused on regression analysis.(i) Basic Principles. Assuming that we are given a training sample Dn=X1,Y1;…;Xn,Yn of i.i.d. 0,1d×R-valued random variables d≥2 with the same distribution as the independent prototype pair D defined by X,Y of size n satisfying the EY2<∞ condition. The space 0,1d is equipped with the standard Euclidean metric. For a given x∈0,1d, the regression function estimate mx can be obtained by(1)mx=EYX=x.Regarding this, the regression function estimatemn for the dataset Dn is assumed consistent if(2)EmnX−mX2⟶0,asn⟶∞.In RF, the final predictor is an average ofM independently randomized regression trees. For a given tree j randomly generated, the prediction value at the query point x is denoted mnX;Θj,Dn, where Θ1,Θ2,…,ΘM are the outputs of a randomized variable Θ (also known as independent random variables).In practice, the randomized variableΘ is used to resample the training set used in successive individual trees building [18]. The finite forest estimate (combination of the M trees grown) is obtained by(3)mM,nX;Θ1,…,ΘM,Dn=1M∑j=1MmnX;Θj,Dn.In the randomForest package, the default value ofM is 500. Since M may be chosen arbitrarily large (“limited only by available computing resources”), this estimate can be generalized from a modelling point of view, letting M tend to infinity by an infinite forest estimate [15], obtained by(4)m∞,nX;Dn=EΘmnX;Θ,Dn.It is noteworthy in practice that the above expectation is evaluated by Monte Carlo, i.e., by generatingM (usually large) random trees and taking the average of individual outcomes. This procedure is justified by the law of large numbers [14, 30] which asserts that the almost surely conditional on Dn is expressed as follows:(5)limM⟶∞mM,nX;Θ1,…,ΘM,Dn=m∞,nX;Dn.(ii) Algorithm. The whole process of the RF algorithm can be summarized as follows:LetM be the number of trees to grow in the forest and mtry be the number of variables to select at each node. The following steps are considered:(i)
We resample the training set to create new bootstrap samples.(ii)
For each bootstrap sample, we create an unpruned tree that is grown until all nodes contain observations no more than the maximal terminal node size defined by the user (a prespecified parameter).(iii)
At each split node, a randomly selected subset of predictors is uniformly chosen. In regression tasks forp predictors, a value equal to p/3 is utilized as the default value in the R package randomForest. Then, the best split using only these predictors is determined or returned. In addition, all observations are not used in tree construction, and a subset an is picked from the whole sample Dn. If an=n, the algorithm runs in a bootstrap mode, whereas an<n corresponds to subsampling (with replacement).(iv)
The RF overall prediction value is the average (in regression) of the results obtained from allM given by individual trees [15].(iii) Parameter Tuning. To date, there are limited studies focusing on RF parameter tuning [31]. According to Probst and Boulesteix [32], the most important parameters in RF methods are as follows:(i)
Number of variables to possibly split in each node(ii)
Minimal terminal node size(iii)
Fraction of observations to sample in each tree(iv)
Number of trees


### 2.1.2. Poisson Regression Extensions Used to Predict Count Data

The prediction performance of RF was compared with that of the Poisson model and its extensions. It is suggested that they be used for assessing count data. The following models are utilized in this investigation and are briefly described as follows:(i) Linear Model. It is a linear function of the form:(6)Y=Xβ+ε;ε∼N0,σε2,where Y is the vectorized form of the response variable, X is a covariate matrix, β is a vector of coefficients, and ε is a vector of normally distributed errors with 0 mean and a variance of σε2 [33].(ii) Poisson Model (Poisson). It is a special case of the generalized linear model (GLM) with a log link function. It is the standard GLM for count data [34]. We assume that the response variable yi follows the Poisson distribution with the mean μi. Its probability mass function is written as(7)PY=yi;μi=e−μiμiyiyi!=expyilogμi−μi−logyi!,where the conditional mean is obtained by μi=Eyi/xi=expx′β with x′ being the vector of covariates and β being the vector of unknown parameters [35–37].(iii) Negative Binomial Model (Neg.Bin). In the negative binomial, a parameter of dispersion K is introduced to accommodate unobserved heterogeneity in the count data. Its probability density function is written as(8)fyixi=τyi+Kτyi+1τKKK+μiKμiK+μiyi,where K is a constant, τ is a function of yi, μi=expx′β, x′ is the vector of covariates, and β is the vector of unknown parameters [38].(iv) Quasi-Poisson. The generalized linear model is used to define the quasi-Poisson model. Let Y be a random variable such that(9)EY=μ,VarY=ϕ×μ,where Var (Y) is Y’s variance, μ>0, ϕ is the dispersion parameter and ϕ>0, and EY is Y’s expectation.(v) Zero-Inflated Poisson model (ZIP). This model is appropriate when there are zeros in excess. Its probability mass function is written as [34, 36](10)PY=yixi,zi=wizi+1−wiziPoissonμi;0xi,ifyi=0,1−wiPoissonμi;yixi,ifyi>0,where zi is the vector of covariates, wi is the probability that the response value is zero, Poissonμi;0xi=exp−μi, and Poissonμi;yixi=exp−μiμiyi/yi!.


## 2.1.1. Random Forest Technique (RF)

RF methods deal with both supervised classification and regression tasks [14]. This research was focused on regression analysis.(i) Basic Principles. Assuming that we are given a training sample Dn=X1,Y1;…;Xn,Yn of i.i.d. 0,1d×R-valued random variables d≥2 with the same distribution as the independent prototype pair D defined by X,Y of size n satisfying the EY2<∞ condition. The space 0,1d is equipped with the standard Euclidean metric. For a given x∈0,1d, the regression function estimate mx can be obtained by(1)mx=EYX=x.Regarding this, the regression function estimatemn for the dataset Dn is assumed consistent if(2)EmnX−mX2⟶0,asn⟶∞.In RF, the final predictor is an average ofM independently randomized regression trees. For a given tree j randomly generated, the prediction value at the query point x is denoted mnX;Θj,Dn, where Θ1,Θ2,…,ΘM are the outputs of a randomized variable Θ (also known as independent random variables).In practice, the randomized variableΘ is used to resample the training set used in successive individual trees building [18]. The finite forest estimate (combination of the M trees grown) is obtained by(3)mM,nX;Θ1,…,ΘM,Dn=1M∑j=1MmnX;Θj,Dn.In the randomForest package, the default value ofM is 500. Since M may be chosen arbitrarily large (“limited only by available computing resources”), this estimate can be generalized from a modelling point of view, letting M tend to infinity by an infinite forest estimate [15], obtained by(4)m∞,nX;Dn=EΘmnX;Θ,Dn.It is noteworthy in practice that the above expectation is evaluated by Monte Carlo, i.e., by generatingM (usually large) random trees and taking the average of individual outcomes. This procedure is justified by the law of large numbers [14, 30] which asserts that the almost surely conditional on Dn is expressed as follows:(5)limM⟶∞mM,nX;Θ1,…,ΘM,Dn=m∞,nX;Dn.(ii) Algorithm. The whole process of the RF algorithm can be summarized as follows:LetM be the number of trees to grow in the forest and mtry be the number of variables to select at each node. The following steps are considered:(i)
We resample the training set to create new bootstrap samples.(ii)
For each bootstrap sample, we create an unpruned tree that is grown until all nodes contain observations no more than the maximal terminal node size defined by the user (a prespecified parameter).(iii)
At each split node, a randomly selected subset of predictors is uniformly chosen. In regression tasks forp predictors, a value equal to p/3 is utilized as the default value in the R package randomForest. Then, the best split using only these predictors is determined or returned. In addition, all observations are not used in tree construction, and a subset an is picked from the whole sample Dn. If an=n, the algorithm runs in a bootstrap mode, whereas an<n corresponds to subsampling (with replacement).(iv)
The RF overall prediction value is the average (in regression) of the results obtained from allM given by individual trees [15].(iii) Parameter Tuning. To date, there are limited studies focusing on RF parameter tuning [31]. According to Probst and Boulesteix [32], the most important parameters in RF methods are as follows:(i)
Number of variables to possibly split in each node(ii)
Minimal terminal node size(iii)
Fraction of observations to sample in each tree(iv)
Number of trees


## 2.1.2. Poisson Regression Extensions Used to Predict Count Data

The prediction performance of RF was compared with that of the Poisson model and its extensions. It is suggested that they be used for assessing count data. The following models are utilized in this investigation and are briefly described as follows:(i) Linear Model. It is a linear function of the form:(6)Y=Xβ+ε;ε∼N0,σε2,where Y is the vectorized form of the response variable, X is a covariate matrix, β is a vector of coefficients, and ε is a vector of normally distributed errors with 0 mean and a variance of σε2 [33].(ii) Poisson Model (Poisson). It is a special case of the generalized linear model (GLM) with a log link function. It is the standard GLM for count data [34]. We assume that the response variable yi follows the Poisson distribution with the mean μi. Its probability mass function is written as(7)PY=yi;μi=e−μiμiyiyi!=expyilogμi−μi−logyi!,where the conditional mean is obtained by μi=Eyi/xi=expx′β with x′ being the vector of covariates and β being the vector of unknown parameters [35–37].(iii) Negative Binomial Model (Neg.Bin). In the negative binomial, a parameter of dispersion K is introduced to accommodate unobserved heterogeneity in the count data. Its probability density function is written as(8)fyixi=τyi+Kτyi+1τKKK+μiKμiK+μiyi,where K is a constant, τ is a function of yi, μi=expx′β, x′ is the vector of covariates, and β is the vector of unknown parameters [38].(iv) Quasi-Poisson. The generalized linear model is used to define the quasi-Poisson model. Let Y be a random variable such that(9)EY=μ,VarY=ϕ×μ,where Var (Y) is Y’s variance, μ>0, ϕ is the dispersion parameter and ϕ>0, and EY is Y’s expectation.(v) Zero-Inflated Poisson model (ZIP). This model is appropriate when there are zeros in excess. Its probability mass function is written as [34, 36](10)PY=yixi,zi=wizi+1−wiziPoissonμi;0xi,ifyi=0,1−wiPoissonμi;yixi,ifyi>0,where zi is the vector of covariates, wi is the probability that the response value is zero, Poissonμi;0xi=exp−μi, and Poissonμi;yixi=exp−μiμiyi/yi!.


## 2.2. Simulation Studies

### 2.2.1. Effect of Dataset Characteristics on RF Algorithm Parameters’ Setting and Accuracy

Three sample sizes (N=50,250, and 1250), three numbers of predictors (p=8, 16, and 24), and five types of predictors (100, 75, 50, and 25% continuous, and 100% categorical) were considered. Schönbrodt and Perugini [39] stated that n should approach 250 for stable estimates. Biau [18] showed that n>500 is sufficiently large, while 50 represents a small population size. Altogether, 45 data types were set. For each case, 1000 datasets were created. Three RF parameters were considered to determine the best combination, while the number of trees was maintained constant (500):(i)
Number of variables to be randomly selected at each split (mtry): mtry= 2, 3, …, p, where p= the number of predictors(ii)
Proportion of samples for model training (sample size) with four levels (55,63.2,70, and 80%)(iii)
Minimal number of samples within each terminal node (node size= 2, 3, 4, …, 9)Depending on the number of predictors, the combination of different parameters under study resulted in 224, 480, and 736 different scenarios for each data type. Each scenario was run 1000 times to find optimal parameter combinations. For all datasets, the generated explanatory variables were not correlated and were not informative. The dependent variable was generated randomly from the Poisson distribution with lambda equal to one, which is denoted asPoissonλ=1. Categorical predictors were randomly generated by varying the number of levels from 2 to 8. Uncorrelated quantitative covariates were generated from a multivariate normal distribution (MVN) y∼Nn0,1, whose density function in the case of n dimensions is written as(11)fyiμ,Σ=fyi1,yi2,…,yin′μ=μ1,…,μn′,Σ=σ12σ1n⋱σn1σn2=12πnΣexp−12yi−μ′Σ−1yi−μ.The symbolΣ refers to the determinant of the matrix Σ. The matrix Σ must be positive semidefinite to assure that the most likely point is μ=μ1,μ2,…,μn and that, as yi moves away from μ in any direction, then the probability of observing yi declines. The denominator in the formula for the MVN is a normalizing constant, which assures that the distribution integrates to 1.0 [40].


### 2.2.2. Effect of Overdispersion on RF Parameters’ Setting and Accuracy

In order to assess the effect of overdispersion on the setting and performance of RF parameters, negative binomial-distributedYnb and quasi-Poisson-approximated Yqp variables were randomly generated. The variance-to-mean relationship of Ynb is quadratic, while the variance-to-mean relationship of Yqp is linear such as(12)EYnb=μ,EYqp=μ,VarYnb=μ+ϕ×μ2,VarYqp=μ+ϕ×μ,where μ>0, ϕ is the dispersion parameter, and ϕ>0. For this objective, μ=3 with varying ϕ=1, 3, and 5. Three different samples (N=50,250, and 1250) and five types of predictors as in the first simulation design were considered. Figure 1 shows the shape and the distribution of randomly generated overdispersed outcomes.Figure 1
Histograms and density plots of some randomly generated overdispersed outcomes.


### 2.2.3. RF Performance Compared to Poisson Model Extensions

Negative binomial-distributed variablesYi were generated with the linear and quadratic variance-to-mean relationships. We let the mean μ vary as a function of two informative covariates m=2 and different levels of the dispersion parameter ϕ. The mean of the outcome Yi of 1000 observations was obtained according to the model logEYim=y|Xim=logμi=α+∑m−12βmXim, where β is a collection of parameters β1=1,β2=0.5 and α=0.145. In addition, three noninformative variables were added from a multivariate normal distribution (MVN) ∼Nn0,1.All datasets were divided into two subsets: training and validation. 70% of the data were randomly selected for model calibration and the remaining 30% for model validation.Poisson regression, negative binomial, zero-inflated Poisson, and zero-inflated negative binomial were compared to the RF regression technique and linear regression applied to the log-transformed response variablelogY+1, as it is considered to be one of the alternatives when dealing with overdispersion or count data [41, 42]. The analysis of variance (ANOVA) test and Tukey’s honest significant difference test (Tukey HSD) were conducted to compare the performance of methods.


### 2.2.4. Criteria for Methods’ Accuracy Assessment

The root mean square errorRMSE and biases were computed for each model. These evaluation criteria are calculated as follows:(13)RMSE=1n∑i=1nPi−Oi2,Bias=1n∑i=1nEPi−Oi,where Pi and Oi are the predicted and observed values in the subunit i, respectively, n is the number of samples (subunits) [43], and EPi is the average of the predicted values [44]. The RMSE represents the overall quality of the prediction. Predictions become increasingly optimal as the RMSE approaches zero [43, 45]. All analyses were computed using R 4.0.3 software [46].


## 2.2.1. Effect of Dataset Characteristics on RF Algorithm Parameters’ Setting and Accuracy

Three sample sizes (N=50,250, and 1250), three numbers of predictors (p=8, 16, and 24), and five types of predictors (100, 75, 50, and 25% continuous, and 100% categorical) were considered. Schönbrodt and Perugini [39] stated that n should approach 250 for stable estimates. Biau [18] showed that n>500 is sufficiently large, while 50 represents a small population size. Altogether, 45 data types were set. For each case, 1000 datasets were created. Three RF parameters were considered to determine the best combination, while the number of trees was maintained constant (500):(i)
Number of variables to be randomly selected at each split (mtry): mtry= 2, 3, …, p, where p= the number of predictors(ii)
Proportion of samples for model training (sample size) with four levels (55,63.2,70, and 80%)(iii)
Minimal number of samples within each terminal node (node size= 2, 3, 4, …, 9)Depending on the number of predictors, the combination of different parameters under study resulted in 224, 480, and 736 different scenarios for each data type. Each scenario was run 1000 times to find optimal parameter combinations. For all datasets, the generated explanatory variables were not correlated and were not informative. The dependent variable was generated randomly from the Poisson distribution with lambda equal to one, which is denoted asPoissonλ=1. Categorical predictors were randomly generated by varying the number of levels from 2 to 8. Uncorrelated quantitative covariates were generated from a multivariate normal distribution (MVN) y∼Nn0,1, whose density function in the case of n dimensions is written as(11)fyiμ,Σ=fyi1,yi2,…,yin′μ=μ1,…,μn′,Σ=σ12σ1n⋱σn1σn2=12πnΣexp−12yi−μ′Σ−1yi−μ.The symbolΣ refers to the determinant of the matrix Σ. The matrix Σ must be positive semidefinite to assure that the most likely point is μ=μ1,μ2,…,μn and that, as yi moves away from μ in any direction, then the probability of observing yi declines. The denominator in the formula for the MVN is a normalizing constant, which assures that the distribution integrates to 1.0 [40].


## 2.2.2. Effect of Overdispersion on RF Parameters’ Setting and Accuracy

In order to assess the effect of overdispersion on the setting and performance of RF parameters, negative binomial-distributedYnb and quasi-Poisson-approximated Yqp variables were randomly generated. The variance-to-mean relationship of Ynb is quadratic, while the variance-to-mean relationship of Yqp is linear such as(12)EYnb=μ,EYqp=μ,VarYnb=μ+ϕ×μ2,VarYqp=μ+ϕ×μ,where μ>0, ϕ is the dispersion parameter, and ϕ>0. For this objective, μ=3 with varying ϕ=1, 3, and 5. Three different samples (N=50,250, and 1250) and five types of predictors as in the first simulation design were considered. Figure 1 shows the shape and the distribution of randomly generated overdispersed outcomes.Figure 1
Histograms and density plots of some randomly generated overdispersed outcomes.


## 2.2.3. RF Performance Compared to Poisson Model Extensions

Negative binomial-distributed variablesYi were generated with the linear and quadratic variance-to-mean relationships. We let the mean μ vary as a function of two informative covariates m=2 and different levels of the dispersion parameter ϕ. The mean of the outcome Yi of 1000 observations was obtained according to the model logEYim=y|Xim=logμi=α+∑m−12βmXim, where β is a collection of parameters β1=1,β2=0.5 and α=0.145. In addition, three noninformative variables were added from a multivariate normal distribution (MVN) ∼Nn0,1.All datasets were divided into two subsets: training and validation. 70% of the data were randomly selected for model calibration and the remaining 30% for model validation.Poisson regression, negative binomial, zero-inflated Poisson, and zero-inflated negative binomial were compared to the RF regression technique and linear regression applied to the log-transformed response variablelogY+1, as it is considered to be one of the alternatives when dealing with overdispersion or count data [41, 42]. The analysis of variance (ANOVA) test and Tukey’s honest significant difference test (Tukey HSD) were conducted to compare the performance of methods.


## 2.2.4. Criteria for Methods’ Accuracy Assessment

The root mean square errorRMSE and biases were computed for each model. These evaluation criteria are calculated as follows:(13)RMSE=1n∑i=1nPi−Oi2,Bias=1n∑i=1nEPi−Oi,where Pi and Oi are the predicted and observed values in the subunit i, respectively, n is the number of samples (subunits) [43], and EPi is the average of the predicted values [44]. The RMSE represents the overall quality of the prediction. Predictions become increasingly optimal as the RMSE approaches zero [43, 45]. All analyses were computed using R 4.0.3 software [46].


## 2.3. Motivating Real Datasets

We used two datasets from real case studies to examine RF performance at two magnitudes of overdispersion caused by outliers or zero inflation. The first example is tick burdens observed on two local goat breeds. Ticks were counted monthly between birth and the age of eight months on ten billy goats belonging to two local breeds [47]. The results are shown in Table 1. The dispersion statistic for this dataset is 3.050, caused by outliers.Table 1
Average tick counts on two local goat breeds.
Age (months)Breed 1Breed 2Billy 1Billy 2Billy 3Billy 4Billy 5Billy 6Billy 7Billy 8Billy 9Billy 10113414344222336187686634371131541014845716947128618516896191319391217613431169161614271241181810324379817312941112361018The second example is the dataset “NMES1988”: Demand for Medical Care in the National Medical Expenditure Survey (NMES) conducted in 1987 and 1988 [48]. This dataset contains 4406 observations on 19 variables. In this example, we examined the relationship between the number of physician office visits (response variable) and health (factor indicating self-perceived health status), age in years (divided by 10), gender, married (factor indicating that the individual is married), income (family income in USD 10000), and insurance (factor indicating that the individual is covered by private insurance) as predictors. The dispersion statistic for this subset is 7.345, caused by outliers and zero inflation.


## 3. Results and Discussion

### 3.1. Results

#### 3.1.1. Effect of Data Features on RF Regression’s Parameters’ Tuning

The effect of predictor types on RF parameter settings for a response variable following a Poisson distribution with the parameter isλ=1.In most scenarios, as shown in Figure2, the RF algorithm performs better with the number of variables randomly selected at each split (mtry), varying between two and four when there are many observations and a few predictors. Figure 2 shows that the number of variables to be randomly selected at each split is independent of the type of data. Nonetheless, the combination of various types of predictors causes some variability in the mtry values that yield the optimal predictive performance. When the sample size is small n=50 and the number of predictors is large p=24, there is a considerable fluctuation in mtry values (Figure 2(g)).Figure 2
Best RF mtry as a function of the population size, number, and type of predictors (Numeric=100% are quantitative predictors, Mixt1_3=25% are quantitative predictors, Mixt1_2=50% are quantitative predictors, Mixt3_1=75% are quantitative predictors, and Categorical=100% are categorical predictors). (a) N=50,P=8, (b) N=250,P=8, (c) N=1250,P=8, (d) N=50,P=16, (e) N=250,P=16, (f) N=1250,P=16, (g) N=50,P=24, (h) N=250,P=24, and (i) N=1250,P=24.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure3 indicates that when the sample size is large (n = 1250) and the number of predictors is small (8), the RF algorithm performs better in the majority of cases with the highest minimum number of observations in each terminal node (between seven and nine). Furthermore, when the fraction of qualitative predictors was high, even for large sample sizes and a small number of predictors, the proportion of scenarios that achieved the optimal predictive performance with high values of the minimum number of observations in each terminal node (best node size ≥) also increased. Small sample sizes n=50 showed a considerable variation in the number of observations in terminal nodes.Figure 3
Best RF node size as a function of the population size, number, and type of predictors (Numeric=100% are quantitative predictors, Mixt1_3=25% are quantitative predictors, Mixt1_2=50% are quantitative predictors, Mixt3_1=75% are quantitative predictors, and Categorical=100% are categorical predictors). (a) N=50,P=8, (b) N=250,P=8, (c) N=1250,P=8, (d) N=50,P=16, (e) N=250,P=16, (f) N=1250,P=16, (g) N=50,P=24, (h) N=250,P=24, and (i) N=1250,P=24.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure4 shows that, in most circumstances, the RF method performs better when 55 percent of the whole sample size is used, followed by 63.2 percent. However, it turns out that the variation in this parameter depends on the ratio of the sample size to the number of predictors. The greater the ratio, the better the algorithm performs with fewer observations, whereas the lower the ratio, the more observations are used.Figure 4
Best RF sample size as a function of the population size, number, and type of predictors (Numeric=100% are quantitative predictors, Mixt1_3=25% are quantitative predictors, Mixt1_2=50% are quantitative predictors, Mixt3_1=75% are quantitative predictors, and Categorical=100% are categorical predictors). (a) N=50,P=8, (b) N=250,P=8, (c) N=1250,P=8, (d) N=50,P=16, (e) N=250,P=16, (f) N=1250,P=16, (g) N=50,P=24, (h) N=250,P=24, and (i) N=1250,P=24.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


#### 3.1.2. Effect of Overdispersion on the RF Parameter Settings for Various Predictor Types

As shown in Table2, the RF algorithm works similarly for all types of predictors and magnitudes of dispersion. In most cases, regardless of the number of predictors or the extent of overdispersion, the optimal performance of RF is obtained by setting the number of variables to randomly sample at each split equal to two (mainly a small number of variables between two and four). Overall, for all types of predictors and observation sizes, it is observed that dispersion does not influence RF parameters’ tuning.Table 2
Effect of predictor types and dispersion amplitude on the number of variables randomly selected at each split.
Data typesVariance-to-mean relationshipBest mtryN = 50 (%)N = 250 (%)N = 1250 (%)p=8p=24p=8p=24p=8p=24ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5CategoricalLinear2818990565865989799808486100100100919294>2≤p/3936272730131171514000986>p/3≤2∗p/38749115000000000000>2∗p/3210840100310000000Quadratic2898888627474989999818592100100100979998>2≤p/378428161820118158000312>p/3≤2∗p/3231476010100000000>2∗p/321763200000000000025% of predictors are quantitativeLinear20010001001001001001001000100100100100100100100>2≤p/300000000001000000000>p/3≤2∗p/31000010000000000000000>2∗p/301000000000000000000Quadratic210010010010010001001001001000100100100100100100100>2≤p/300000000001000000000>2∗p/30000010000000000000050% of predictors are quantitativeLinear210010010000010010000100100100100100100100100>2≤p/3000001000010010000000000>p/3≤2∗p/300010000000000000000>2∗p/300001000000000000000Quadratic210010010010010010010010010010010010010010010010010010075% of predictors are quantitativeLinear210010010001001001001001001000100100100100100100100>2≤p/30001000000001000000000Quadratic210001000100100100100100100100100100100100100100100>2≤p/30100010000000000000000QuantitativeLinear2100010010010010001001000100100100100100100100100>2≤p/3010000001000010000000000Quadratic21001000100100100100100100100100100100100100100100100>2≤p/300100000000000000000According to Table3, when predictors are categorical, for the same number of observations, the RF algorithm tends to use smaller samples to train the model as the number of predictors decreases. Thus, when eight predictors were considered, regardless of their type or level of overdispersion, the best performance was obtained in most cases using 55% of the total number of observations. In most scenarios, this proportion is followed by 63.2%. Furthermore, as the number of predictors increases, the 55% proportion loses dominance. As a result, the likelihood of obtaining the best performance with a large proportion of samples to train the model increases.Table 3
Effect of overdispersion on the RF optimal size of the sample to draw.
Data typesVariance-to-mean relationshipsample sizeN = 50 (%)N = 250 (%)N = 1250 (%)p=8p=24p=8p=24p=8p=24ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5CategoricalLinear0.555646563049316371592737398989743439430.632183018241724191226322527109222725250.7131114211623111262216211242321160.8131312251822759192213000161516Quadratic0.555762574036375868624138377970864636430.6322113212125302819232718251926101932260.712121019221398121722201402315170.810131220172055315221810412171425% of predictors are quantitativeLinear0.5510000000100010010001001000100010000.632001001000000000001000100000.70000000100001000000001000.8010000100100000000000000Quadratic0.551000100000001001000100100001001001000.632000100001000001000001000000.7010000100100010000000000000.80000000000000100000050% of predictors are quantitativeLinear0.551001001001000100010010000010010000000.632000000000010000010010001000.70000001000010001000000000.80000100000000000001000Quadratic0.550100000010010000001001000001000.63210001000100000010000001000000.700000000100000000010000.8000100010000001001000001000075% of predictors are quantitativeLinear0.550000100001001000001001000100000.6321001000100000001000100001000000.700100001000000000001001000.80000001000001000000000Quadratic0.550010000010000001001000100001000.632010001000100010001000000010010000.7100000100000001000010000000.800000000100000000000QuantitativeLinear0.551000100100001001001001000100010010001001000.63201000000000010000000000.700001000000000000100000.80000010000000010000000Quadratic0.5500100000010001000010010010001001000.63201000000100000000000000.70001000100001000000000000.81000001000000010010000010000However, when all predictors are quantitative continuous or mixed with qualitative predictors, no trend was observed for the best size of the sample to draw and the minimum size of terminal nodes in RF, regardless of the complexity of the response variable.In most cases, the best performance of RF is obtained with a minimum terminal node size greater than five, the default value used for regression. Table4 shows that, regardless of the degree of overdispersion of the dependent variable, the RF algorithm performs best with a minimum terminal node size greater than five. For qualitative predictors, the preponderance of the high number of observations per node decreases as the number of predictors increases. However, as the number of predictors increases, the fraction of cases, where the best performance is achieved with a large number of observations per node, decreases. Thus, it turns out that globally, irrespective of the type of predictors, for the same number of observations, the higher the number of predictors, the smaller the minimum size of the terminal nodes. However, it should be noted that the number of predictors is not always inversely proportional to the minimum size of terminal nodes. In some scenarios, the minimum size of terminal nodes does not vary with the number of predictors. In contrast, in others, the best performance is obtained with small minimum sizes of terminal nodes, even when there are few predictors. Thus, when there are a few observations and many predictors, no trend is predictable, and the parameter will depend only on the complexity of the relationship between the data F.Table 4
The impact of overdispersion and data features on the minimal terminal node size tuning.
Data typesVariance-to-mean relationshipNode sizeN = 50 (%)N = 250 (%)N = 1250 (%)p=8p=24p=8p=24p=8p=24ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5CategoricalLinear≤52526183231345112272625000181517>5≤9121016293126312128474839221628435652>9636466393840646870262636788472392931Quadratic≤5221619312731465301818100151010>5≤9151711262214221625335245141610495751>963677043515574787037303785849036333925% of predictors are quantitativeLinear≤5100001001000000010010000010000>5≤90010000100100001000001001000100100>901000000010010000010000000Quadratic≤50000010000000000010000>5≤9010000100010010000010001000000>91000100100000010010010001000100010010050% of predictors are quantitativeLinear≤5000001000000010000010000>5≤90010010000100001000010010000100100>910010000100001001000100000100000Quadratic≤50100000100000010010000001000>5≤91000010000100000000010010000>9001000100001001001000010010000010075% of predictors are quantitativeLinear≤510001000010010010001000000010000>5≤9000000000010010010010000100100>90100010010000010000000100000Quadratic≤501000100001000100010000010001000>5≤90000000001000100100100010000>9100010001001000100000000000100QuantitativeLinear≤50100100001000100100100010010000100100100>5≤9000100100000001000000000>910000000100000000100100000Quadratic≤51000100001000010001001000001000100>5≤9000100100000010000010010001000>901000000100100000010000000


#### 3.1.3. Impact of Overdispersion on the RF Predictive Performance for Different Types of Predictors

There was no difference among the types of covariates considered to be part of this simulation study. Using the best combination of parameters, RF performed similarly for all types of covariates regardless of the level of dispersion in the response variable. However, the performance of the RF for the sample sizen=50 varied enormously for different levels of dispersion, whereas the values were stable for large sample sizes (in our case n=1250). RF performance was impacted by the magnitude of overdispersion, but relatively more accurate results were obtained when the population size was large in comparison to a small population size, where the RMSE varied enormously for different levels of dispersion (Figures 5(a) and 5(b)).Figure 5
Impact of the response variable overdispersion on RF accuracy (relation is the outcome’s variance-to-mean relationship). (a)RMSE and (b) BIAS.
(a)(b)


#### 3.1.4. Assessment of the RF Performance Compared to Traditional Approaches for Count Data Prediction

In this section, we analyse the performance of RF regression compared to classical log-linear, Poisson, quasi-Poisson, negative binomial, and zero-inflated Poisson models for different sample sizes and overdispersion degrees. Table5 and Figure 5 show that linear regression applied to the log-transformed response variable has the lowest performance. However, in most scenarios, the biases and RMSE obtained using the RF approach to predict the results of the overdispersed count are not statistically different from the biases and RMSE obtained using traditional Poisson family methods. Furthermore, no statistical difference was observed when the variance-mean relationship was quadratic, with some exceptions for low overdispersion. However, for the linear variance-to-mean relationship, the mean RMSE of RF is higher than that of GLMs, i.e., less accurate with higher variability. As a result, while predicting overdispersed outcomes, the most accurate and precise method will have to be identified case by case, depending on the available data. In some situations, RF will be the best performing method, while in others, it will not. Figures 6(a) and 6(b) present RF compared to the biases of GLM family models and RMSE regression lines as a function of dispersion and indicate the 95% confidence interval (range), represented by the grey zone in which the biases and RMSE would be if the experiment were repeated. Globally, the mean of the RMSE and the biases of RF tend to increase as dispersion increases. Consequently, the slope of the regression of the RMSE in the function of dispersion and the range of RMSE values are lower for outcomes with linear variance-to-mean relationships than for those with quadratic variance-to-mean relationships.Table 5
Description of methods’ predictive performance and Tukey’s HSD pairwise comparisons for different levels of overdispersion.
Variance-to-mean relationshipMethodsϕ=1ϕ=2ϕ=3ϕ=4ϕ=5ϕ=6ϕ=7ϕ=8ϕ=9ϕ=10BIASLinearLog_linear1.53±0.18b1.62±0.23c1.65±0.24c1.58±0.25b1.58±0.26d1.58±0.26d1.56±0.29c1.56±0.29c1.70±0.34c1.67±0.32cNeg.Bin0.12±0.10a0.17±0.15ab0.18±0.13a0.23±0.17a0.26±0.22a0.26±0.22a0.28±0.22a0.53±0.44a0.34±0.26a0.32±0.25aPoisson0.12±0.09a0.15±0.11a0.18±0.13a0.20±0.16a0.22±0.17b0.22±0.17b0.22±0.17b0.27±0.20b0.27±0.21b0.27±0.20bQuasi-Poisson0.12±0.09a0.15±0.11a0.18±0.13a0.20±0.16a0.22±0.17b0.22±0.17b0.22±0.17b0.27±0.20b0.27±0.21b0.27±0.20bRF0.14±0.11a0.18±0.13b0.23±0.16b0.23±0.18a0.26±0.20ac0.26±0.20ac0.25±0.19ab0.29±0.22b0.30±0.23ab0.30±0.23abZIP0.12±0.09a0.16±0.12ab0.19±0.14a0.21±0.16a0.23±0.17bc0.23±0.17bc0.24±0.18b0.27±0.21b0.27±0.22b0.28±0.22bQuadraticLog_linear1.65±0.28c1.67±0.57c1.67±0.32b1.99±0.66b1.90±0.58b1.90±0.58b1.93±0.66b1.46±0.50b2.06±0.72b2.02±0.91bNeg.Bin0.23±0.17a0.41±0.43ab0.34±0.27a0.61±0.50a0.59±0.46a0.59±0.46a0.59±0.49a0.40±0.33a0.71±0.53a0.75±0.69aPoisson0.24±0.18a0.41±0.42ab0.34±0.26a0.63±0.54a0.60±0.47a0.60±0.47a0.57±0.49a0.38±0.32a0.73±0.57a0.66±0.66aQuasi-Poisson0.24±0.18a0.41±0.42ab0.34±0.26a0.63±0.54a0.60±0.47a0.60±0.47a0.57±0.49a0.38±0.32a0.73±0.57a0.66±0.66aRF0.30±0.22b0.49±0.51a0.38±0.28a0.66±0.50a0.65±0.49a0.65±0.49a0.63±0.51a0.43±0.35a0.80±0.61a0.73±0.68aZIP0.24±0.18a0.40±0.44b0.33±0.25a0.60±0.50a0.58±0.46a0.58±0.46a0.56±0.48a0.39±0.32a0.71±0.53a0.65±0.63aRMSELinearLog_linear3.98±0.56c4.68±0.89d5.19±0.93d4.99±0.87c4.97±0.75d5.23±0.82d5.23±1d5.54±0.95c5.95±1.14d6.03±1.16cNeg.Bin2.15±0.2a2.94±0.67a3.14±0.3a3.71±0.51a3.81±0.47ab4.18±0.64a4.26±0.56a6.02±2.27a5.29±0.93a5.19±0.91aPoisson2.13±0.18a2.76±0.27b3.14±0.29a3.61±0.49a3.75±0.45a4.03±0.54b4.09±0.58b4.68±0.73b5.01±0.87b4.98±0.87bQuasi-Poisson2.13±0.18a2.76±0.27b3.14±0.29a3.61±0.49a3.75±0.45a4.03±0.54b4.09±0.58b4.68±0.73b5.01±0.87b4.98±0.87bRF2.42±0.31b3.24±0.55c3.88±0.61b3.86±0.66b4.03±0.49c4.33±0.61c4.43±0.78c4.86±0.76b5.19±0.93ac5.28±0.93aZIP2.14±0.21a2.9±0.38a3.39±0.51c3.69±0.61a3.87±0.49b4.11±0.58ab4.25±0.74a4.75±0.77b5.04±0.93bc5.11±0.96abQuadraticLog_linear5.93±1.45c8.13±6.74a6.54±1.29c8.31±2.71a10.45±4.29a10.46±4.41a10.4±5.3a6.96±2.51a11.96±6.01a11.88±8.47aNeg.Bin4.44±1a7.19±5.85a5.93±1.03ab8.05±2.38a10.42±4.09a10.23±3.91a10.32±4.81a6.85±2.16a11.74±5.25a12.02±7.39aPoisson4.52±0.98a7.12±5.37a5.88±1.07ab8.1±2.4a10.59±4.14a10.26±3.95a10.15±4.67a6.82±2.27a11.89±5.29a11.79±7.56aQuasi-Poisson4.52±0.98a7.12±5.37a5.88±1.07ab8.1±2.4a10.59±4.14a10.26±3.95a10.15±4.67a6.82±2.27a11.89±5.29a11.79±7.56aRF5.1±1.11b7.68±6.11a6.06±1.03a8.15±2.35a10.62±3.88a10.49±4.01a10.46±4.62a7.03±2.4a12.08±5.41a12.04±8aZIP4.57±1.08a7.11±5.71a5.81±1.12b7.99±2.4a10.38±4.01a10.14±3.97a10.08±4.71a6.81±2.24a11.74±5.34a11.73±7.49aFigure 6
Effect of overdispersion on RF predictive performance compared to GLM models (relation is the outcome’s variance-to-mean relationship). (a)Bias and (b) RMSE.
(a)(b)


#### 3.1.5. Real Dataset Results

RF yielded the lowest mean of the absolute values of biases (1.51± 1.16) compared to linear regression applied to the logarithmic transformation of tick burden (8.24 ± 1.43), negative binomial (1.56 ± 1.19), quasi-Poisson (1.52 ± 1.17), and Poisson (1.52 ± 1.17) regressions for the dataset with overdispersion caused by outliers (Figure 7(a)). However, the difference between the mean absolute values of biases is not significant. Quasi-Poisson and Poisson (7.72 ± 1.32) have the lowest RMSE values, followed by the RF (7.75 ± 1.34) and negative binomial (7.82 ± 1.29) but are not significantly different according to the Tukey HSD test (Figures 7(a) and 7(b)).Figure 7
Predictive performance of RF compared to GLMs for tick burden and NMES1988 datasets. (a) Tick burden dataset:BIAS, (b) tick burden dataset: RMSE, (c) NMES1988 dataset: BIAS, and (d) NMES1988 dataset: RMSE.
(a)(b)(c)(d)For the NMES1988 dataset with overdispersion caused by outliers and zero inflation, quasi-Poisson, Poisson, zero-inflated Poisson, and negative binomial produced the lowest absolute value of biases (0.18± 0.14), followed by RF (0.19 ± 0.14). Furthermore, the zero-inflated Poisson model (6.56 ± 0.39) yielded the lowest RMSE, followed by quasi-Poisson (6.57 ± 0.39), negative binomial (6.57 ± 0.39), and RF (6.59 ± 0.4). As for the tick burden dataset, no statistically significant difference was observed between RF and traditional Poisson family models (Figures 7(c) and 7(d)).


### 3.2. Discussion

#### 3.2.1. Impact of Data Features on RF Regression

One of the characteristics of many databases that practitioners consider is the complex relationship among data, which typically yields enormous noise. In some situations, this complexity may directly or indirectly result from variables at the hand and affect many statistical procedures. Computer-based statistical methods are more appropriate than traditional statistical models for effectively solving this complex problem and dealing with large-scale issues. However, despite an increasing interest and practical use, computer-based methods have rarely been evaluated [18]. There is a growing trend of replacing traditional statistical models with machine learning in several research areas to improve predictive accuracy. In addition, these techniques are now used to model data that do not fit into the “small-n-large-p” category.The RF method works well for all covariates (categorical, continuous, or their mixture) in regression with or without overdispersion in the response variable. No original predictor transformation is required for the predictive task when encoding categorical covariates. Furthermore, RF accuracy is not significantly affected by covariate types. Singh et al. [25] found similar results in classification, suggesting that RF performs equally well with numeric-only or text-only datasets. However, fluctuations in RF parameter settings were observed for exclusively categorical predictors when the number of observations or the relation (ratio) between the number of observations and the number of predictors was low.Frequently, empirical studies on classification find that the default RF’s parameter values yield good prediction accuracy [19, 49].However, among all possible combinations, in most scenarios, the best RF performance was not observed with the default parameters in Breiman and Cutler’s random forests for the classification and regression package (randomForest), widely used in R software [44]. For regression tasks, the number of variables to select in each split should equal the number of predictors divided by 3. Our results support the findings of Han and Kim [50], which imply that the default mtry of p/3 cannot guarantee the highest level of accuracy in regression. Additionally, Strobl et al. [51] argued that various mtry values should be taken into account.For all types of covariates, RF performed better when a small number of variables were selected at each splitmtry=2, regardless of the response variable overdispersion, while for some scenarios, mtry=p/3 was too small. RF prediction accuracy is expected to be higher for smaller values of mtry in general, according to Strobl et al. [51]. Smaller mtry values are preferred if the dataset’s features have similar relevance. If not, greater mtry values are preferred [50]. Furthermore, a small number of predictors in each split can better utilize the available information in many informative predictors with different strengths [52].However, for a variable to have a chance to appear in a sufficient number of trees, a large number of preselected predictors are required when there are many predictor variables. Therefore, the average variable importance measure will be accurate and not just a random fluctuation based on enough trials [51]. However, Svetnik et al. [49] claimed that RF performance changes very little over a wide range of mtry values, except for extremes (mtry=1 and mtry=p).RF was relatively insensitive to changes in the minimum number of observations in each terminal node. The minimum node size was adjusted to speed up the computation for large datasets. Our findings suggest that the minimum node size behaviour is influenced by the relationship between the total sample size and the number of predictors. The best RF performance was obtained with many observations in each end node for largen and small p (between seven and nine). When dealing with situations with many observations and a few variables, increasing the terminal node size with the sample size improved the convergence rate of the prediction accuracy [53]. There was no discernible pattern for small n and large p. Although the minimum node size appears to have a minor impact on predictive accuracy, it is necessary to consider this parameter when tuning RF. When determining the ideal combination of the best parameters of RF, the prevalence of specific values can indicate the stability of estimates, avoiding the effect of random fluctuations. According to Lin and Jeon [53], growing large trees, that is, using small observations in each terminal node, generally results in the highest performance for “small-n-large-p” situations.We observed that sample size behaviour is also influenced by the relationship between the number of observations and predictors. When the ratio of the number of observations to the number of predictors is high, RF produces the best results with fewer observations (55 percent or 63.2 percent of the total sample size). Buja and Stuetzle [54] argued that half-sampling might be more suited because bootstrap sampling shares several theoretical traits. Wösthoff [55] claimed that the subsample size is not a critical parameter as the results of sizes 0.5 and the default value (0.632) are nearly identical. For permutation importance estimation in classification, significant predictor variables may not be recognized if the sample size is too small. The variability of the importance measure tends to increase if the subsample size is so large that the out-of-bag sample becomes very small [55].Considering overdispersion, we realized that, for small sizes (less than 250), the best mtry value was randomly selected among possible values. In contrast, for large population sizes, this value was maintained constant (mtry = 2) for all types of covariates. According to Svetnik et al. [49], this parameter is the main RF tuning parameter. Furthermore, the sample size influenced the RF accuracy level. For small sample sizes (less than 250), the effect of overdispersion resulted in an erratic trend, challenging to interpret. In contrast, the trend was linear for large population sizes, illustrating the overdispersion effect.Although RF parameters seemed insensitive to the degree of overdispersion of the response variable, RF accuracy was affected. Thus, despite being an essential step in RF analyses, parameter tuning did not affect the best parameter selection when the response variable was overdispersed. Some authors have recommended data transformation for RF analyses, especially population estimation. Transforming the response variable allows the RF algorithm to find good splits in data as the relationship between distance-based covariates and the population is, in most cases, more uniformly distributed [56]. Transformation methods may be applied to observed response variables before RF model fitting for higher prediction accuracy. However, O’Hara and Kotze [57] stated that we should not transform count data because, in addition to performing poorly, they can lead to impossible predictions.


#### 3.2.2. RF Performance Compared to Poisson Model Extensions

According to Anderson et al. [33], one of the most critical decisions in the population modelling process is understanding and selecting the appropriate model structure based on the available data. RF regression is an efficient and credible alternative in prediction tasks, as reported by previous studies [12]. The results obtained in this study suggest that RF is as accurate as Poisson, zero-inflated Poisson, and negative binomial regressions when dealing with overdispersed count data. The tick burden dataset confirmed this statement. Novel methods such as RF are more efficient for predicting spatial patterns in species distribution models than more traditional linear models [58]. Furthermore, they consistently outperform more established methods. RF’s consistency and predictive performance could be an additional point to explore and reorient the debate on modelling count data. Some researchers suggest using the Poisson model and its extensions [57]. Others recommend using the linear model when the objective is to test the significance of estimates [41]. On the contrary, others think there is no best solution because it depends on the database in the presence [42]. However, these models have shown similar performance for different levels of overdispersion. Considering GLMs, RF is robust to misspecification of the relationship between the mean and variance [57]. Anderson et al. [33] claimed that traditional regression-based models may not be the most effective model, depending on the available census data (for example, presence of sampling bias, missing data, outliers, and imperfect detection). Given the predictive accuracy demonstrated by RF, its flexibility in dealing with covariates of various natures, and other recognised advantages such as the ability to deal with complex interaction structures, highly correlated variables and the measurement of variable significance, RF is an effective alternative method for predicting count outcomes [52]. RF may be a tool of capital importance for decision and policymakers working in limited data areas or with unreliable census data, especially in developing countries. Furthermore, producing accurate and reliable estimates is a problem of interest for resource allocation, public health, outbreak early detection, food security, conservation policy, environmental planning, and others [33, 56, 59].


## 3.1. Results

### 3.1.1. Effect of Data Features on RF Regression’s Parameters’ Tuning

The effect of predictor types on RF parameter settings for a response variable following a Poisson distribution with the parameter isλ=1.In most scenarios, as shown in Figure2, the RF algorithm performs better with the number of variables randomly selected at each split (mtry), varying between two and four when there are many observations and a few predictors. Figure 2 shows that the number of variables to be randomly selected at each split is independent of the type of data. Nonetheless, the combination of various types of predictors causes some variability in the mtry values that yield the optimal predictive performance. When the sample size is small n=50 and the number of predictors is large p=24, there is a considerable fluctuation in mtry values (Figure 2(g)).Figure 2
Best RF mtry as a function of the population size, number, and type of predictors (Numeric=100% are quantitative predictors, Mixt1_3=25% are quantitative predictors, Mixt1_2=50% are quantitative predictors, Mixt3_1=75% are quantitative predictors, and Categorical=100% are categorical predictors). (a) N=50,P=8, (b) N=250,P=8, (c) N=1250,P=8, (d) N=50,P=16, (e) N=250,P=16, (f) N=1250,P=16, (g) N=50,P=24, (h) N=250,P=24, and (i) N=1250,P=24.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure3 indicates that when the sample size is large (n = 1250) and the number of predictors is small (8), the RF algorithm performs better in the majority of cases with the highest minimum number of observations in each terminal node (between seven and nine). Furthermore, when the fraction of qualitative predictors was high, even for large sample sizes and a small number of predictors, the proportion of scenarios that achieved the optimal predictive performance with high values of the minimum number of observations in each terminal node (best node size ≥) also increased. Small sample sizes n=50 showed a considerable variation in the number of observations in terminal nodes.Figure 3
Best RF node size as a function of the population size, number, and type of predictors (Numeric=100% are quantitative predictors, Mixt1_3=25% are quantitative predictors, Mixt1_2=50% are quantitative predictors, Mixt3_1=75% are quantitative predictors, and Categorical=100% are categorical predictors). (a) N=50,P=8, (b) N=250,P=8, (c) N=1250,P=8, (d) N=50,P=16, (e) N=250,P=16, (f) N=1250,P=16, (g) N=50,P=24, (h) N=250,P=24, and (i) N=1250,P=24.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure4 shows that, in most circumstances, the RF method performs better when 55 percent of the whole sample size is used, followed by 63.2 percent. However, it turns out that the variation in this parameter depends on the ratio of the sample size to the number of predictors. The greater the ratio, the better the algorithm performs with fewer observations, whereas the lower the ratio, the more observations are used.Figure 4
Best RF sample size as a function of the population size, number, and type of predictors (Numeric=100% are quantitative predictors, Mixt1_3=25% are quantitative predictors, Mixt1_2=50% are quantitative predictors, Mixt3_1=75% are quantitative predictors, and Categorical=100% are categorical predictors). (a) N=50,P=8, (b) N=250,P=8, (c) N=1250,P=8, (d) N=50,P=16, (e) N=250,P=16, (f) N=1250,P=16, (g) N=50,P=24, (h) N=250,P=24, and (i) N=1250,P=24.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 3.1.2. Effect of Overdispersion on the RF Parameter Settings for Various Predictor Types

As shown in Table2, the RF algorithm works similarly for all types of predictors and magnitudes of dispersion. In most cases, regardless of the number of predictors or the extent of overdispersion, the optimal performance of RF is obtained by setting the number of variables to randomly sample at each split equal to two (mainly a small number of variables between two and four). Overall, for all types of predictors and observation sizes, it is observed that dispersion does not influence RF parameters’ tuning.Table 2
Effect of predictor types and dispersion amplitude on the number of variables randomly selected at each split.
Data typesVariance-to-mean relationshipBest mtryN = 50 (%)N = 250 (%)N = 1250 (%)p=8p=24p=8p=24p=8p=24ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5CategoricalLinear2818990565865989799808486100100100919294>2≤p/3936272730131171514000986>p/3≤2∗p/38749115000000000000>2∗p/3210840100310000000Quadratic2898888627474989999818592100100100979998>2≤p/378428161820118158000312>p/3≤2∗p/3231476010100000000>2∗p/321763200000000000025% of predictors are quantitativeLinear20010001001001001001001000100100100100100100100>2≤p/300000000001000000000>p/3≤2∗p/31000010000000000000000>2∗p/301000000000000000000Quadratic210010010010010001001001001000100100100100100100100>2≤p/300000000001000000000>2∗p/30000010000000000000050% of predictors are quantitativeLinear210010010000010010000100100100100100100100100>2≤p/3000001000010010000000000>p/3≤2∗p/300010000000000000000>2∗p/300001000000000000000Quadratic210010010010010010010010010010010010010010010010010010075% of predictors are quantitativeLinear210010010001001001001001001000100100100100100100100>2≤p/30001000000001000000000Quadratic210001000100100100100100100100100100100100100100100>2≤p/30100010000000000000000QuantitativeLinear2100010010010010001001000100100100100100100100100>2≤p/3010000001000010000000000Quadratic21001000100100100100100100100100100100100100100100100>2≤p/300100000000000000000According to Table3, when predictors are categorical, for the same number of observations, the RF algorithm tends to use smaller samples to train the model as the number of predictors decreases. Thus, when eight predictors were considered, regardless of their type or level of overdispersion, the best performance was obtained in most cases using 55% of the total number of observations. In most scenarios, this proportion is followed by 63.2%. Furthermore, as the number of predictors increases, the 55% proportion loses dominance. As a result, the likelihood of obtaining the best performance with a large proportion of samples to train the model increases.Table 3
Effect of overdispersion on the RF optimal size of the sample to draw.
Data typesVariance-to-mean relationshipsample sizeN = 50 (%)N = 250 (%)N = 1250 (%)p=8p=24p=8p=24p=8p=24ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5CategoricalLinear0.555646563049316371592737398989743439430.632183018241724191226322527109222725250.7131114211623111262216211242321160.8131312251822759192213000161516Quadratic0.555762574036375868624138377970864636430.6322113212125302819232718251926101932260.712121019221398121722201402315170.810131220172055315221810412171425% of predictors are quantitativeLinear0.5510000000100010010001001000100010000.632001001000000000001000100000.70000000100001000000001000.8010000100100000000000000Quadratic0.551000100000001001000100100001001001000.632000100001000001000001000000.7010000100100010000000000000.80000000000000100000050% of predictors are quantitativeLinear0.551001001001000100010010000010010000000.632000000000010000010010001000.70000001000010001000000000.80000100000000000001000Quadratic0.550100000010010000001001000001000.63210001000100000010000001000000.700000000100000000010000.8000100010000001001000001000075% of predictors are quantitativeLinear0.550000100001001000001001000100000.6321001000100000001000100001000000.700100001000000000001001000.80000001000001000000000Quadratic0.550010000010000001001000100001000.632010001000100010001000000010010000.7100000100000001000010000000.800000000100000000000QuantitativeLinear0.551000100100001001001001000100010010001001000.63201000000000010000000000.700001000000000000100000.80000010000000010000000Quadratic0.5500100000010001000010010010001001000.63201000000100000000000000.70001000100001000000000000.81000001000000010010000010000However, when all predictors are quantitative continuous or mixed with qualitative predictors, no trend was observed for the best size of the sample to draw and the minimum size of terminal nodes in RF, regardless of the complexity of the response variable.In most cases, the best performance of RF is obtained with a minimum terminal node size greater than five, the default value used for regression. Table4 shows that, regardless of the degree of overdispersion of the dependent variable, the RF algorithm performs best with a minimum terminal node size greater than five. For qualitative predictors, the preponderance of the high number of observations per node decreases as the number of predictors increases. However, as the number of predictors increases, the fraction of cases, where the best performance is achieved with a large number of observations per node, decreases. Thus, it turns out that globally, irrespective of the type of predictors, for the same number of observations, the higher the number of predictors, the smaller the minimum size of the terminal nodes. However, it should be noted that the number of predictors is not always inversely proportional to the minimum size of terminal nodes. In some scenarios, the minimum size of terminal nodes does not vary with the number of predictors. In contrast, in others, the best performance is obtained with small minimum sizes of terminal nodes, even when there are few predictors. Thus, when there are a few observations and many predictors, no trend is predictable, and the parameter will depend only on the complexity of the relationship between the data F.Table 4
The impact of overdispersion and data features on the minimal terminal node size tuning.
Data typesVariance-to-mean relationshipNode sizeN = 50 (%)N = 250 (%)N = 1250 (%)p=8p=24p=8p=24p=8p=24ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5CategoricalLinear≤52526183231345112272625000181517>5≤9121016293126312128474839221628435652>9636466393840646870262636788472392931Quadratic≤5221619312731465301818100151010>5≤9151711262214221625335245141610495751>963677043515574787037303785849036333925% of predictors are quantitativeLinear≤5100001001000000010010000010000>5≤90010000100100001000001001000100100>901000000010010000010000000Quadratic≤50000010000000000010000>5≤9010000100010010000010001000000>91000100100000010010010001000100010010050% of predictors are quantitativeLinear≤5000001000000010000010000>5≤90010010000100001000010010000100100>910010000100001001000100000100000Quadratic≤50100000100000010010000001000>5≤91000010000100000000010010000>9001000100001001001000010010000010075% of predictors are quantitativeLinear≤510001000010010010001000000010000>5≤9000000000010010010010000100100>90100010010000010000000100000Quadratic≤501000100001000100010000010001000>5≤90000000001000100100100010000>9100010001001000100000000000100QuantitativeLinear≤50100100001000100100100010010000100100100>5≤9000100100000001000000000>910000000100000000100100000Quadratic≤51000100001000010001001000001000100>5≤9000100100000010000010010001000>901000000100100000010000000


### 3.1.3. Impact of Overdispersion on the RF Predictive Performance for Different Types of Predictors

There was no difference among the types of covariates considered to be part of this simulation study. Using the best combination of parameters, RF performed similarly for all types of covariates regardless of the level of dispersion in the response variable. However, the performance of the RF for the sample sizen=50 varied enormously for different levels of dispersion, whereas the values were stable for large sample sizes (in our case n=1250). RF performance was impacted by the magnitude of overdispersion, but relatively more accurate results were obtained when the population size was large in comparison to a small population size, where the RMSE varied enormously for different levels of dispersion (Figures 5(a) and 5(b)).Figure 5
Impact of the response variable overdispersion on RF accuracy (relation is the outcome’s variance-to-mean relationship). (a)RMSE and (b) BIAS.
(a)(b)


### 3.1.4. Assessment of the RF Performance Compared to Traditional Approaches for Count Data Prediction

In this section, we analyse the performance of RF regression compared to classical log-linear, Poisson, quasi-Poisson, negative binomial, and zero-inflated Poisson models for different sample sizes and overdispersion degrees. Table5 and Figure 5 show that linear regression applied to the log-transformed response variable has the lowest performance. However, in most scenarios, the biases and RMSE obtained using the RF approach to predict the results of the overdispersed count are not statistically different from the biases and RMSE obtained using traditional Poisson family methods. Furthermore, no statistical difference was observed when the variance-mean relationship was quadratic, with some exceptions for low overdispersion. However, for the linear variance-to-mean relationship, the mean RMSE of RF is higher than that of GLMs, i.e., less accurate with higher variability. As a result, while predicting overdispersed outcomes, the most accurate and precise method will have to be identified case by case, depending on the available data. In some situations, RF will be the best performing method, while in others, it will not. Figures 6(a) and 6(b) present RF compared to the biases of GLM family models and RMSE regression lines as a function of dispersion and indicate the 95% confidence interval (range), represented by the grey zone in which the biases and RMSE would be if the experiment were repeated. Globally, the mean of the RMSE and the biases of RF tend to increase as dispersion increases. Consequently, the slope of the regression of the RMSE in the function of dispersion and the range of RMSE values are lower for outcomes with linear variance-to-mean relationships than for those with quadratic variance-to-mean relationships.Table 5
Description of methods’ predictive performance and Tukey’s HSD pairwise comparisons for different levels of overdispersion.
Variance-to-mean relationshipMethodsϕ=1ϕ=2ϕ=3ϕ=4ϕ=5ϕ=6ϕ=7ϕ=8ϕ=9ϕ=10BIASLinearLog_linear1.53±0.18b1.62±0.23c1.65±0.24c1.58±0.25b1.58±0.26d1.58±0.26d1.56±0.29c1.56±0.29c1.70±0.34c1.67±0.32cNeg.Bin0.12±0.10a0.17±0.15ab0.18±0.13a0.23±0.17a0.26±0.22a0.26±0.22a0.28±0.22a0.53±0.44a0.34±0.26a0.32±0.25aPoisson0.12±0.09a0.15±0.11a0.18±0.13a0.20±0.16a0.22±0.17b0.22±0.17b0.22±0.17b0.27±0.20b0.27±0.21b0.27±0.20bQuasi-Poisson0.12±0.09a0.15±0.11a0.18±0.13a0.20±0.16a0.22±0.17b0.22±0.17b0.22±0.17b0.27±0.20b0.27±0.21b0.27±0.20bRF0.14±0.11a0.18±0.13b0.23±0.16b0.23±0.18a0.26±0.20ac0.26±0.20ac0.25±0.19ab0.29±0.22b0.30±0.23ab0.30±0.23abZIP0.12±0.09a0.16±0.12ab0.19±0.14a0.21±0.16a0.23±0.17bc0.23±0.17bc0.24±0.18b0.27±0.21b0.27±0.22b0.28±0.22bQuadraticLog_linear1.65±0.28c1.67±0.57c1.67±0.32b1.99±0.66b1.90±0.58b1.90±0.58b1.93±0.66b1.46±0.50b2.06±0.72b2.02±0.91bNeg.Bin0.23±0.17a0.41±0.43ab0.34±0.27a0.61±0.50a0.59±0.46a0.59±0.46a0.59±0.49a0.40±0.33a0.71±0.53a0.75±0.69aPoisson0.24±0.18a0.41±0.42ab0.34±0.26a0.63±0.54a0.60±0.47a0.60±0.47a0.57±0.49a0.38±0.32a0.73±0.57a0.66±0.66aQuasi-Poisson0.24±0.18a0.41±0.42ab0.34±0.26a0.63±0.54a0.60±0.47a0.60±0.47a0.57±0.49a0.38±0.32a0.73±0.57a0.66±0.66aRF0.30±0.22b0.49±0.51a0.38±0.28a0.66±0.50a0.65±0.49a0.65±0.49a0.63±0.51a0.43±0.35a0.80±0.61a0.73±0.68aZIP0.24±0.18a0.40±0.44b0.33±0.25a0.60±0.50a0.58±0.46a0.58±0.46a0.56±0.48a0.39±0.32a0.71±0.53a0.65±0.63aRMSELinearLog_linear3.98±0.56c4.68±0.89d5.19±0.93d4.99±0.87c4.97±0.75d5.23±0.82d5.23±1d5.54±0.95c5.95±1.14d6.03±1.16cNeg.Bin2.15±0.2a2.94±0.67a3.14±0.3a3.71±0.51a3.81±0.47ab4.18±0.64a4.26±0.56a6.02±2.27a5.29±0.93a5.19±0.91aPoisson2.13±0.18a2.76±0.27b3.14±0.29a3.61±0.49a3.75±0.45a4.03±0.54b4.09±0.58b4.68±0.73b5.01±0.87b4.98±0.87bQuasi-Poisson2.13±0.18a2.76±0.27b3.14±0.29a3.61±0.49a3.75±0.45a4.03±0.54b4.09±0.58b4.68±0.73b5.01±0.87b4.98±0.87bRF2.42±0.31b3.24±0.55c3.88±0.61b3.86±0.66b4.03±0.49c4.33±0.61c4.43±0.78c4.86±0.76b5.19±0.93ac5.28±0.93aZIP2.14±0.21a2.9±0.38a3.39±0.51c3.69±0.61a3.87±0.49b4.11±0.58ab4.25±0.74a4.75±0.77b5.04±0.93bc5.11±0.96abQuadraticLog_linear5.93±1.45c8.13±6.74a6.54±1.29c8.31±2.71a10.45±4.29a10.46±4.41a10.4±5.3a6.96±2.51a11.96±6.01a11.88±8.47aNeg.Bin4.44±1a7.19±5.85a5.93±1.03ab8.05±2.38a10.42±4.09a10.23±3.91a10.32±4.81a6.85±2.16a11.74±5.25a12.02±7.39aPoisson4.52±0.98a7.12±5.37a5.88±1.07ab8.1±2.4a10.59±4.14a10.26±3.95a10.15±4.67a6.82±2.27a11.89±5.29a11.79±7.56aQuasi-Poisson4.52±0.98a7.12±5.37a5.88±1.07ab8.1±2.4a10.59±4.14a10.26±3.95a10.15±4.67a6.82±2.27a11.89±5.29a11.79±7.56aRF5.1±1.11b7.68±6.11a6.06±1.03a8.15±2.35a10.62±3.88a10.49±4.01a10.46±4.62a7.03±2.4a12.08±5.41a12.04±8aZIP4.57±1.08a7.11±5.71a5.81±1.12b7.99±2.4a10.38±4.01a10.14±3.97a10.08±4.71a6.81±2.24a11.74±5.34a11.73±7.49aFigure 6
Effect of overdispersion on RF predictive performance compared to GLM models (relation is the outcome’s variance-to-mean relationship). (a)Bias and (b) RMSE.
(a)(b)


### 3.1.5. Real Dataset Results

RF yielded the lowest mean of the absolute values of biases (1.51± 1.16) compared to linear regression applied to the logarithmic transformation of tick burden (8.24 ± 1.43), negative binomial (1.56 ± 1.19), quasi-Poisson (1.52 ± 1.17), and Poisson (1.52 ± 1.17) regressions for the dataset with overdispersion caused by outliers (Figure 7(a)). However, the difference between the mean absolute values of biases is not significant. Quasi-Poisson and Poisson (7.72 ± 1.32) have the lowest RMSE values, followed by the RF (7.75 ± 1.34) and negative binomial (7.82 ± 1.29) but are not significantly different according to the Tukey HSD test (Figures 7(a) and 7(b)).Figure 7
Predictive performance of RF compared to GLMs for tick burden and NMES1988 datasets. (a) Tick burden dataset:BIAS, (b) tick burden dataset: RMSE, (c) NMES1988 dataset: BIAS, and (d) NMES1988 dataset: RMSE.
(a)(b)(c)(d)For the NMES1988 dataset with overdispersion caused by outliers and zero inflation, quasi-Poisson, Poisson, zero-inflated Poisson, and negative binomial produced the lowest absolute value of biases (0.18± 0.14), followed by RF (0.19 ± 0.14). Furthermore, the zero-inflated Poisson model (6.56 ± 0.39) yielded the lowest RMSE, followed by quasi-Poisson (6.57 ± 0.39), negative binomial (6.57 ± 0.39), and RF (6.59 ± 0.4). As for the tick burden dataset, no statistically significant difference was observed between RF and traditional Poisson family models (Figures 7(c) and 7(d)).


## 3.1.1. Effect of Data Features on RF Regression’s Parameters’ Tuning

The effect of predictor types on RF parameter settings for a response variable following a Poisson distribution with the parameter isλ=1.In most scenarios, as shown in Figure2, the RF algorithm performs better with the number of variables randomly selected at each split (mtry), varying between two and four when there are many observations and a few predictors. Figure 2 shows that the number of variables to be randomly selected at each split is independent of the type of data. Nonetheless, the combination of various types of predictors causes some variability in the mtry values that yield the optimal predictive performance. When the sample size is small n=50 and the number of predictors is large p=24, there is a considerable fluctuation in mtry values (Figure 2(g)).Figure 2
Best RF mtry as a function of the population size, number, and type of predictors (Numeric=100% are quantitative predictors, Mixt1_3=25% are quantitative predictors, Mixt1_2=50% are quantitative predictors, Mixt3_1=75% are quantitative predictors, and Categorical=100% are categorical predictors). (a) N=50,P=8, (b) N=250,P=8, (c) N=1250,P=8, (d) N=50,P=16, (e) N=250,P=16, (f) N=1250,P=16, (g) N=50,P=24, (h) N=250,P=24, and (i) N=1250,P=24.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure3 indicates that when the sample size is large (n = 1250) and the number of predictors is small (8), the RF algorithm performs better in the majority of cases with the highest minimum number of observations in each terminal node (between seven and nine). Furthermore, when the fraction of qualitative predictors was high, even for large sample sizes and a small number of predictors, the proportion of scenarios that achieved the optimal predictive performance with high values of the minimum number of observations in each terminal node (best node size ≥) also increased. Small sample sizes n=50 showed a considerable variation in the number of observations in terminal nodes.Figure 3
Best RF node size as a function of the population size, number, and type of predictors (Numeric=100% are quantitative predictors, Mixt1_3=25% are quantitative predictors, Mixt1_2=50% are quantitative predictors, Mixt3_1=75% are quantitative predictors, and Categorical=100% are categorical predictors). (a) N=50,P=8, (b) N=250,P=8, (c) N=1250,P=8, (d) N=50,P=16, (e) N=250,P=16, (f) N=1250,P=16, (g) N=50,P=24, (h) N=250,P=24, and (i) N=1250,P=24.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure4 shows that, in most circumstances, the RF method performs better when 55 percent of the whole sample size is used, followed by 63.2 percent. However, it turns out that the variation in this parameter depends on the ratio of the sample size to the number of predictors. The greater the ratio, the better the algorithm performs with fewer observations, whereas the lower the ratio, the more observations are used.Figure 4
Best RF sample size as a function of the population size, number, and type of predictors (Numeric=100% are quantitative predictors, Mixt1_3=25% are quantitative predictors, Mixt1_2=50% are quantitative predictors, Mixt3_1=75% are quantitative predictors, and Categorical=100% are categorical predictors). (a) N=50,P=8, (b) N=250,P=8, (c) N=1250,P=8, (d) N=50,P=16, (e) N=250,P=16, (f) N=1250,P=16, (g) N=50,P=24, (h) N=250,P=24, and (i) N=1250,P=24.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 3.1.2. Effect of Overdispersion on the RF Parameter Settings for Various Predictor Types

As shown in Table2, the RF algorithm works similarly for all types of predictors and magnitudes of dispersion. In most cases, regardless of the number of predictors or the extent of overdispersion, the optimal performance of RF is obtained by setting the number of variables to randomly sample at each split equal to two (mainly a small number of variables between two and four). Overall, for all types of predictors and observation sizes, it is observed that dispersion does not influence RF parameters’ tuning.Table 2
Effect of predictor types and dispersion amplitude on the number of variables randomly selected at each split.
Data typesVariance-to-mean relationshipBest mtryN = 50 (%)N = 250 (%)N = 1250 (%)p=8p=24p=8p=24p=8p=24ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5CategoricalLinear2818990565865989799808486100100100919294>2≤p/3936272730131171514000986>p/3≤2∗p/38749115000000000000>2∗p/3210840100310000000Quadratic2898888627474989999818592100100100979998>2≤p/378428161820118158000312>p/3≤2∗p/3231476010100000000>2∗p/321763200000000000025% of predictors are quantitativeLinear20010001001001001001001000100100100100100100100>2≤p/300000000001000000000>p/3≤2∗p/31000010000000000000000>2∗p/301000000000000000000Quadratic210010010010010001001001001000100100100100100100100>2≤p/300000000001000000000>2∗p/30000010000000000000050% of predictors are quantitativeLinear210010010000010010000100100100100100100100100>2≤p/3000001000010010000000000>p/3≤2∗p/300010000000000000000>2∗p/300001000000000000000Quadratic210010010010010010010010010010010010010010010010010010075% of predictors are quantitativeLinear210010010001001001001001001000100100100100100100100>2≤p/30001000000001000000000Quadratic210001000100100100100100100100100100100100100100100>2≤p/30100010000000000000000QuantitativeLinear2100010010010010001001000100100100100100100100100>2≤p/3010000001000010000000000Quadratic21001000100100100100100100100100100100100100100100100>2≤p/300100000000000000000According to Table3, when predictors are categorical, for the same number of observations, the RF algorithm tends to use smaller samples to train the model as the number of predictors decreases. Thus, when eight predictors were considered, regardless of their type or level of overdispersion, the best performance was obtained in most cases using 55% of the total number of observations. In most scenarios, this proportion is followed by 63.2%. Furthermore, as the number of predictors increases, the 55% proportion loses dominance. As a result, the likelihood of obtaining the best performance with a large proportion of samples to train the model increases.Table 3
Effect of overdispersion on the RF optimal size of the sample to draw.
Data typesVariance-to-mean relationshipsample sizeN = 50 (%)N = 250 (%)N = 1250 (%)p=8p=24p=8p=24p=8p=24ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5CategoricalLinear0.555646563049316371592737398989743439430.632183018241724191226322527109222725250.7131114211623111262216211242321160.8131312251822759192213000161516Quadratic0.555762574036375868624138377970864636430.6322113212125302819232718251926101932260.712121019221398121722201402315170.810131220172055315221810412171425% of predictors are quantitativeLinear0.5510000000100010010001001000100010000.632001001000000000001000100000.70000000100001000000001000.8010000100100000000000000Quadratic0.551000100000001001000100100001001001000.632000100001000001000001000000.7010000100100010000000000000.80000000000000100000050% of predictors are quantitativeLinear0.551001001001000100010010000010010000000.632000000000010000010010001000.70000001000010001000000000.80000100000000000001000Quadratic0.550100000010010000001001000001000.63210001000100000010000001000000.700000000100000000010000.8000100010000001001000001000075% of predictors are quantitativeLinear0.550000100001001000001001000100000.6321001000100000001000100001000000.700100001000000000001001000.80000001000001000000000Quadratic0.550010000010000001001000100001000.632010001000100010001000000010010000.7100000100000001000010000000.800000000100000000000QuantitativeLinear0.551000100100001001001001000100010010001001000.63201000000000010000000000.700001000000000000100000.80000010000000010000000Quadratic0.5500100000010001000010010010001001000.63201000000100000000000000.70001000100001000000000000.81000001000000010010000010000However, when all predictors are quantitative continuous or mixed with qualitative predictors, no trend was observed for the best size of the sample to draw and the minimum size of terminal nodes in RF, regardless of the complexity of the response variable.In most cases, the best performance of RF is obtained with a minimum terminal node size greater than five, the default value used for regression. Table4 shows that, regardless of the degree of overdispersion of the dependent variable, the RF algorithm performs best with a minimum terminal node size greater than five. For qualitative predictors, the preponderance of the high number of observations per node decreases as the number of predictors increases. However, as the number of predictors increases, the fraction of cases, where the best performance is achieved with a large number of observations per node, decreases. Thus, it turns out that globally, irrespective of the type of predictors, for the same number of observations, the higher the number of predictors, the smaller the minimum size of the terminal nodes. However, it should be noted that the number of predictors is not always inversely proportional to the minimum size of terminal nodes. In some scenarios, the minimum size of terminal nodes does not vary with the number of predictors. In contrast, in others, the best performance is obtained with small minimum sizes of terminal nodes, even when there are few predictors. Thus, when there are a few observations and many predictors, no trend is predictable, and the parameter will depend only on the complexity of the relationship between the data F.Table 4
The impact of overdispersion and data features on the minimal terminal node size tuning.
Data typesVariance-to-mean relationshipNode sizeN = 50 (%)N = 250 (%)N = 1250 (%)p=8p=24p=8p=24p=8p=24ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5ϕ=1ϕ=3ϕ=5CategoricalLinear≤52526183231345112272625000181517>5≤9121016293126312128474839221628435652>9636466393840646870262636788472392931Quadratic≤5221619312731465301818100151010>5≤9151711262214221625335245141610495751>963677043515574787037303785849036333925% of predictors are quantitativeLinear≤5100001001000000010010000010000>5≤90010000100100001000001001000100100>901000000010010000010000000Quadratic≤50000010000000000010000>5≤9010000100010010000010001000000>91000100100000010010010001000100010010050% of predictors are quantitativeLinear≤5000001000000010000010000>5≤90010010000100001000010010000100100>910010000100001001000100000100000Quadratic≤50100000100000010010000001000>5≤91000010000100000000010010000>9001000100001001001000010010000010075% of predictors are quantitativeLinear≤510001000010010010001000000010000>5≤9000000000010010010010000100100>90100010010000010000000100000Quadratic≤501000100001000100010000010001000>5≤90000000001000100100100010000>9100010001001000100000000000100QuantitativeLinear≤50100100001000100100100010010000100100100>5≤9000100100000001000000000>910000000100000000100100000Quadratic≤51000100001000010001001000001000100>5≤9000100100000010000010010001000>901000000100100000010000000


## 3.1.3. Impact of Overdispersion on the RF Predictive Performance for Different Types of Predictors

There was no difference among the types of covariates considered to be part of this simulation study. Using the best combination of parameters, RF performed similarly for all types of covariates regardless of the level of dispersion in the response variable. However, the performance of the RF for the sample sizen=50 varied enormously for different levels of dispersion, whereas the values were stable for large sample sizes (in our case n=1250). RF performance was impacted by the magnitude of overdispersion, but relatively more accurate results were obtained when the population size was large in comparison to a small population size, where the RMSE varied enormously for different levels of dispersion (Figures 5(a) and 5(b)).Figure 5
Impact of the response variable overdispersion on RF accuracy (relation is the outcome’s variance-to-mean relationship). (a)RMSE and (b) BIAS.
(a)(b)


## 3.1.4. Assessment of the RF Performance Compared to Traditional Approaches for Count Data Prediction

In this section, we analyse the performance of RF regression compared to classical log-linear, Poisson, quasi-Poisson, negative binomial, and zero-inflated Poisson models for different sample sizes and overdispersion degrees. Table5 and Figure 5 show that linear regression applied to the log-transformed response variable has the lowest performance. However, in most scenarios, the biases and RMSE obtained using the RF approach to predict the results of the overdispersed count are not statistically different from the biases and RMSE obtained using traditional Poisson family methods. Furthermore, no statistical difference was observed when the variance-mean relationship was quadratic, with some exceptions for low overdispersion. However, for the linear variance-to-mean relationship, the mean RMSE of RF is higher than that of GLMs, i.e., less accurate with higher variability. As a result, while predicting overdispersed outcomes, the most accurate and precise method will have to be identified case by case, depending on the available data. In some situations, RF will be the best performing method, while in others, it will not. Figures 6(a) and 6(b) present RF compared to the biases of GLM family models and RMSE regression lines as a function of dispersion and indicate the 95% confidence interval (range), represented by the grey zone in which the biases and RMSE would be if the experiment were repeated. Globally, the mean of the RMSE and the biases of RF tend to increase as dispersion increases. Consequently, the slope of the regression of the RMSE in the function of dispersion and the range of RMSE values are lower for outcomes with linear variance-to-mean relationships than for those with quadratic variance-to-mean relationships.Table 5
Description of methods’ predictive performance and Tukey’s HSD pairwise comparisons for different levels of overdispersion.
Variance-to-mean relationshipMethodsϕ=1ϕ=2ϕ=3ϕ=4ϕ=5ϕ=6ϕ=7ϕ=8ϕ=9ϕ=10BIASLinearLog_linear1.53±0.18b1.62±0.23c1.65±0.24c1.58±0.25b1.58±0.26d1.58±0.26d1.56±0.29c1.56±0.29c1.70±0.34c1.67±0.32cNeg.Bin0.12±0.10a0.17±0.15ab0.18±0.13a0.23±0.17a0.26±0.22a0.26±0.22a0.28±0.22a0.53±0.44a0.34±0.26a0.32±0.25aPoisson0.12±0.09a0.15±0.11a0.18±0.13a0.20±0.16a0.22±0.17b0.22±0.17b0.22±0.17b0.27±0.20b0.27±0.21b0.27±0.20bQuasi-Poisson0.12±0.09a0.15±0.11a0.18±0.13a0.20±0.16a0.22±0.17b0.22±0.17b0.22±0.17b0.27±0.20b0.27±0.21b0.27±0.20bRF0.14±0.11a0.18±0.13b0.23±0.16b0.23±0.18a0.26±0.20ac0.26±0.20ac0.25±0.19ab0.29±0.22b0.30±0.23ab0.30±0.23abZIP0.12±0.09a0.16±0.12ab0.19±0.14a0.21±0.16a0.23±0.17bc0.23±0.17bc0.24±0.18b0.27±0.21b0.27±0.22b0.28±0.22bQuadraticLog_linear1.65±0.28c1.67±0.57c1.67±0.32b1.99±0.66b1.90±0.58b1.90±0.58b1.93±0.66b1.46±0.50b2.06±0.72b2.02±0.91bNeg.Bin0.23±0.17a0.41±0.43ab0.34±0.27a0.61±0.50a0.59±0.46a0.59±0.46a0.59±0.49a0.40±0.33a0.71±0.53a0.75±0.69aPoisson0.24±0.18a0.41±0.42ab0.34±0.26a0.63±0.54a0.60±0.47a0.60±0.47a0.57±0.49a0.38±0.32a0.73±0.57a0.66±0.66aQuasi-Poisson0.24±0.18a0.41±0.42ab0.34±0.26a0.63±0.54a0.60±0.47a0.60±0.47a0.57±0.49a0.38±0.32a0.73±0.57a0.66±0.66aRF0.30±0.22b0.49±0.51a0.38±0.28a0.66±0.50a0.65±0.49a0.65±0.49a0.63±0.51a0.43±0.35a0.80±0.61a0.73±0.68aZIP0.24±0.18a0.40±0.44b0.33±0.25a0.60±0.50a0.58±0.46a0.58±0.46a0.56±0.48a0.39±0.32a0.71±0.53a0.65±0.63aRMSELinearLog_linear3.98±0.56c4.68±0.89d5.19±0.93d4.99±0.87c4.97±0.75d5.23±0.82d5.23±1d5.54±0.95c5.95±1.14d6.03±1.16cNeg.Bin2.15±0.2a2.94±0.67a3.14±0.3a3.71±0.51a3.81±0.47ab4.18±0.64a4.26±0.56a6.02±2.27a5.29±0.93a5.19±0.91aPoisson2.13±0.18a2.76±0.27b3.14±0.29a3.61±0.49a3.75±0.45a4.03±0.54b4.09±0.58b4.68±0.73b5.01±0.87b4.98±0.87bQuasi-Poisson2.13±0.18a2.76±0.27b3.14±0.29a3.61±0.49a3.75±0.45a4.03±0.54b4.09±0.58b4.68±0.73b5.01±0.87b4.98±0.87bRF2.42±0.31b3.24±0.55c3.88±0.61b3.86±0.66b4.03±0.49c4.33±0.61c4.43±0.78c4.86±0.76b5.19±0.93ac5.28±0.93aZIP2.14±0.21a2.9±0.38a3.39±0.51c3.69±0.61a3.87±0.49b4.11±0.58ab4.25±0.74a4.75±0.77b5.04±0.93bc5.11±0.96abQuadraticLog_linear5.93±1.45c8.13±6.74a6.54±1.29c8.31±2.71a10.45±4.29a10.46±4.41a10.4±5.3a6.96±2.51a11.96±6.01a11.88±8.47aNeg.Bin4.44±1a7.19±5.85a5.93±1.03ab8.05±2.38a10.42±4.09a10.23±3.91a10.32±4.81a6.85±2.16a11.74±5.25a12.02±7.39aPoisson4.52±0.98a7.12±5.37a5.88±1.07ab8.1±2.4a10.59±4.14a10.26±3.95a10.15±4.67a6.82±2.27a11.89±5.29a11.79±7.56aQuasi-Poisson4.52±0.98a7.12±5.37a5.88±1.07ab8.1±2.4a10.59±4.14a10.26±3.95a10.15±4.67a6.82±2.27a11.89±5.29a11.79±7.56aRF5.1±1.11b7.68±6.11a6.06±1.03a8.15±2.35a10.62±3.88a10.49±4.01a10.46±4.62a7.03±2.4a12.08±5.41a12.04±8aZIP4.57±1.08a7.11±5.71a5.81±1.12b7.99±2.4a10.38±4.01a10.14±3.97a10.08±4.71a6.81±2.24a11.74±5.34a11.73±7.49aFigure 6
Effect of overdispersion on RF predictive performance compared to GLM models (relation is the outcome’s variance-to-mean relationship). (a)Bias and (b) RMSE.
(a)(b)


## 3.1.5. Real Dataset Results

RF yielded the lowest mean of the absolute values of biases (1.51± 1.16) compared to linear regression applied to the logarithmic transformation of tick burden (8.24 ± 1.43), negative binomial (1.56 ± 1.19), quasi-Poisson (1.52 ± 1.17), and Poisson (1.52 ± 1.17) regressions for the dataset with overdispersion caused by outliers (Figure 7(a)). However, the difference between the mean absolute values of biases is not significant. Quasi-Poisson and Poisson (7.72 ± 1.32) have the lowest RMSE values, followed by the RF (7.75 ± 1.34) and negative binomial (7.82 ± 1.29) but are not significantly different according to the Tukey HSD test (Figures 7(a) and 7(b)).Figure 7
Predictive performance of RF compared to GLMs for tick burden and NMES1988 datasets. (a) Tick burden dataset:BIAS, (b) tick burden dataset: RMSE, (c) NMES1988 dataset: BIAS, and (d) NMES1988 dataset: RMSE.
(a)(b)(c)(d)For the NMES1988 dataset with overdispersion caused by outliers and zero inflation, quasi-Poisson, Poisson, zero-inflated Poisson, and negative binomial produced the lowest absolute value of biases (0.18± 0.14), followed by RF (0.19 ± 0.14). Furthermore, the zero-inflated Poisson model (6.56 ± 0.39) yielded the lowest RMSE, followed by quasi-Poisson (6.57 ± 0.39), negative binomial (6.57 ± 0.39), and RF (6.59 ± 0.4). As for the tick burden dataset, no statistically significant difference was observed between RF and traditional Poisson family models (Figures 7(c) and 7(d)).


## 3.2. Discussion

### 3.2.1. Impact of Data Features on RF Regression

One of the characteristics of many databases that practitioners consider is the complex relationship among data, which typically yields enormous noise. In some situations, this complexity may directly or indirectly result from variables at the hand and affect many statistical procedures. Computer-based statistical methods are more appropriate than traditional statistical models for effectively solving this complex problem and dealing with large-scale issues. However, despite an increasing interest and practical use, computer-based methods have rarely been evaluated [18]. There is a growing trend of replacing traditional statistical models with machine learning in several research areas to improve predictive accuracy. In addition, these techniques are now used to model data that do not fit into the “small-n-large-p” category.The RF method works well for all covariates (categorical, continuous, or their mixture) in regression with or without overdispersion in the response variable. No original predictor transformation is required for the predictive task when encoding categorical covariates. Furthermore, RF accuracy is not significantly affected by covariate types. Singh et al. [25] found similar results in classification, suggesting that RF performs equally well with numeric-only or text-only datasets. However, fluctuations in RF parameter settings were observed for exclusively categorical predictors when the number of observations or the relation (ratio) between the number of observations and the number of predictors was low.Frequently, empirical studies on classification find that the default RF’s parameter values yield good prediction accuracy [19, 49].However, among all possible combinations, in most scenarios, the best RF performance was not observed with the default parameters in Breiman and Cutler’s random forests for the classification and regression package (randomForest), widely used in R software [44]. For regression tasks, the number of variables to select in each split should equal the number of predictors divided by 3. Our results support the findings of Han and Kim [50], which imply that the default mtry of p/3 cannot guarantee the highest level of accuracy in regression. Additionally, Strobl et al. [51] argued that various mtry values should be taken into account.For all types of covariates, RF performed better when a small number of variables were selected at each splitmtry=2, regardless of the response variable overdispersion, while for some scenarios, mtry=p/3 was too small. RF prediction accuracy is expected to be higher for smaller values of mtry in general, according to Strobl et al. [51]. Smaller mtry values are preferred if the dataset’s features have similar relevance. If not, greater mtry values are preferred [50]. Furthermore, a small number of predictors in each split can better utilize the available information in many informative predictors with different strengths [52].However, for a variable to have a chance to appear in a sufficient number of trees, a large number of preselected predictors are required when there are many predictor variables. Therefore, the average variable importance measure will be accurate and not just a random fluctuation based on enough trials [51]. However, Svetnik et al. [49] claimed that RF performance changes very little over a wide range of mtry values, except for extremes (mtry=1 and mtry=p).RF was relatively insensitive to changes in the minimum number of observations in each terminal node. The minimum node size was adjusted to speed up the computation for large datasets. Our findings suggest that the minimum node size behaviour is influenced by the relationship between the total sample size and the number of predictors. The best RF performance was obtained with many observations in each end node for largen and small p (between seven and nine). When dealing with situations with many observations and a few variables, increasing the terminal node size with the sample size improved the convergence rate of the prediction accuracy [53]. There was no discernible pattern for small n and large p. Although the minimum node size appears to have a minor impact on predictive accuracy, it is necessary to consider this parameter when tuning RF. When determining the ideal combination of the best parameters of RF, the prevalence of specific values can indicate the stability of estimates, avoiding the effect of random fluctuations. According to Lin and Jeon [53], growing large trees, that is, using small observations in each terminal node, generally results in the highest performance for “small-n-large-p” situations.We observed that sample size behaviour is also influenced by the relationship between the number of observations and predictors. When the ratio of the number of observations to the number of predictors is high, RF produces the best results with fewer observations (55 percent or 63.2 percent of the total sample size). Buja and Stuetzle [54] argued that half-sampling might be more suited because bootstrap sampling shares several theoretical traits. Wösthoff [55] claimed that the subsample size is not a critical parameter as the results of sizes 0.5 and the default value (0.632) are nearly identical. For permutation importance estimation in classification, significant predictor variables may not be recognized if the sample size is too small. The variability of the importance measure tends to increase if the subsample size is so large that the out-of-bag sample becomes very small [55].Considering overdispersion, we realized that, for small sizes (less than 250), the best mtry value was randomly selected among possible values. In contrast, for large population sizes, this value was maintained constant (mtry = 2) for all types of covariates. According to Svetnik et al. [49], this parameter is the main RF tuning parameter. Furthermore, the sample size influenced the RF accuracy level. For small sample sizes (less than 250), the effect of overdispersion resulted in an erratic trend, challenging to interpret. In contrast, the trend was linear for large population sizes, illustrating the overdispersion effect.Although RF parameters seemed insensitive to the degree of overdispersion of the response variable, RF accuracy was affected. Thus, despite being an essential step in RF analyses, parameter tuning did not affect the best parameter selection when the response variable was overdispersed. Some authors have recommended data transformation for RF analyses, especially population estimation. Transforming the response variable allows the RF algorithm to find good splits in data as the relationship between distance-based covariates and the population is, in most cases, more uniformly distributed [56]. Transformation methods may be applied to observed response variables before RF model fitting for higher prediction accuracy. However, O’Hara and Kotze [57] stated that we should not transform count data because, in addition to performing poorly, they can lead to impossible predictions.


### 3.2.2. RF Performance Compared to Poisson Model Extensions

According to Anderson et al. [33], one of the most critical decisions in the population modelling process is understanding and selecting the appropriate model structure based on the available data. RF regression is an efficient and credible alternative in prediction tasks, as reported by previous studies [12]. The results obtained in this study suggest that RF is as accurate as Poisson, zero-inflated Poisson, and negative binomial regressions when dealing with overdispersed count data. The tick burden dataset confirmed this statement. Novel methods such as RF are more efficient for predicting spatial patterns in species distribution models than more traditional linear models [58]. Furthermore, they consistently outperform more established methods. RF’s consistency and predictive performance could be an additional point to explore and reorient the debate on modelling count data. Some researchers suggest using the Poisson model and its extensions [57]. Others recommend using the linear model when the objective is to test the significance of estimates [41]. On the contrary, others think there is no best solution because it depends on the database in the presence [42]. However, these models have shown similar performance for different levels of overdispersion. Considering GLMs, RF is robust to misspecification of the relationship between the mean and variance [57]. Anderson et al. [33] claimed that traditional regression-based models may not be the most effective model, depending on the available census data (for example, presence of sampling bias, missing data, outliers, and imperfect detection). Given the predictive accuracy demonstrated by RF, its flexibility in dealing with covariates of various natures, and other recognised advantages such as the ability to deal with complex interaction structures, highly correlated variables and the measurement of variable significance, RF is an effective alternative method for predicting count outcomes [52]. RF may be a tool of capital importance for decision and policymakers working in limited data areas or with unreliable census data, especially in developing countries. Furthermore, producing accurate and reliable estimates is a problem of interest for resource allocation, public health, outbreak early detection, food security, conservation policy, environmental planning, and others [33, 56, 59].


## 3.2.1. Impact of Data Features on RF Regression

One of the characteristics of many databases that practitioners consider is the complex relationship among data, which typically yields enormous noise. In some situations, this complexity may directly or indirectly result from variables at the hand and affect many statistical procedures. Computer-based statistical methods are more appropriate than traditional statistical models for effectively solving this complex problem and dealing with large-scale issues. However, despite an increasing interest and practical use, computer-based methods have rarely been evaluated [18]. There is a growing trend of replacing traditional statistical models with machine learning in several research areas to improve predictive accuracy. In addition, these techniques are now used to model data that do not fit into the “small-n-large-p” category.The RF method works well for all covariates (categorical, continuous, or their mixture) in regression with or without overdispersion in the response variable. No original predictor transformation is required for the predictive task when encoding categorical covariates. Furthermore, RF accuracy is not significantly affected by covariate types. Singh et al. [25] found similar results in classification, suggesting that RF performs equally well with numeric-only or text-only datasets. However, fluctuations in RF parameter settings were observed for exclusively categorical predictors when the number of observations or the relation (ratio) between the number of observations and the number of predictors was low.Frequently, empirical studies on classification find that the default RF’s parameter values yield good prediction accuracy [19, 49].However, among all possible combinations, in most scenarios, the best RF performance was not observed with the default parameters in Breiman and Cutler’s random forests for the classification and regression package (randomForest), widely used in R software [44]. For regression tasks, the number of variables to select in each split should equal the number of predictors divided by 3. Our results support the findings of Han and Kim [50], which imply that the default mtry of p/3 cannot guarantee the highest level of accuracy in regression. Additionally, Strobl et al. [51] argued that various mtry values should be taken into account.For all types of covariates, RF performed better when a small number of variables were selected at each splitmtry=2, regardless of the response variable overdispersion, while for some scenarios, mtry=p/3 was too small. RF prediction accuracy is expected to be higher for smaller values of mtry in general, according to Strobl et al. [51]. Smaller mtry values are preferred if the dataset’s features have similar relevance. If not, greater mtry values are preferred [50]. Furthermore, a small number of predictors in each split can better utilize the available information in many informative predictors with different strengths [52].However, for a variable to have a chance to appear in a sufficient number of trees, a large number of preselected predictors are required when there are many predictor variables. Therefore, the average variable importance measure will be accurate and not just a random fluctuation based on enough trials [51]. However, Svetnik et al. [49] claimed that RF performance changes very little over a wide range of mtry values, except for extremes (mtry=1 and mtry=p).RF was relatively insensitive to changes in the minimum number of observations in each terminal node. The minimum node size was adjusted to speed up the computation for large datasets. Our findings suggest that the minimum node size behaviour is influenced by the relationship between the total sample size and the number of predictors. The best RF performance was obtained with many observations in each end node for largen and small p (between seven and nine). When dealing with situations with many observations and a few variables, increasing the terminal node size with the sample size improved the convergence rate of the prediction accuracy [53]. There was no discernible pattern for small n and large p. Although the minimum node size appears to have a minor impact on predictive accuracy, it is necessary to consider this parameter when tuning RF. When determining the ideal combination of the best parameters of RF, the prevalence of specific values can indicate the stability of estimates, avoiding the effect of random fluctuations. According to Lin and Jeon [53], growing large trees, that is, using small observations in each terminal node, generally results in the highest performance for “small-n-large-p” situations.We observed that sample size behaviour is also influenced by the relationship between the number of observations and predictors. When the ratio of the number of observations to the number of predictors is high, RF produces the best results with fewer observations (55 percent or 63.2 percent of the total sample size). Buja and Stuetzle [54] argued that half-sampling might be more suited because bootstrap sampling shares several theoretical traits. Wösthoff [55] claimed that the subsample size is not a critical parameter as the results of sizes 0.5 and the default value (0.632) are nearly identical. For permutation importance estimation in classification, significant predictor variables may not be recognized if the sample size is too small. The variability of the importance measure tends to increase if the subsample size is so large that the out-of-bag sample becomes very small [55].Considering overdispersion, we realized that, for small sizes (less than 250), the best mtry value was randomly selected among possible values. In contrast, for large population sizes, this value was maintained constant (mtry = 2) for all types of covariates. According to Svetnik et al. [49], this parameter is the main RF tuning parameter. Furthermore, the sample size influenced the RF accuracy level. For small sample sizes (less than 250), the effect of overdispersion resulted in an erratic trend, challenging to interpret. In contrast, the trend was linear for large population sizes, illustrating the overdispersion effect.Although RF parameters seemed insensitive to the degree of overdispersion of the response variable, RF accuracy was affected. Thus, despite being an essential step in RF analyses, parameter tuning did not affect the best parameter selection when the response variable was overdispersed. Some authors have recommended data transformation for RF analyses, especially population estimation. Transforming the response variable allows the RF algorithm to find good splits in data as the relationship between distance-based covariates and the population is, in most cases, more uniformly distributed [56]. Transformation methods may be applied to observed response variables before RF model fitting for higher prediction accuracy. However, O’Hara and Kotze [57] stated that we should not transform count data because, in addition to performing poorly, they can lead to impossible predictions.


## 3.2.2. RF Performance Compared to Poisson Model Extensions

According to Anderson et al. [33], one of the most critical decisions in the population modelling process is understanding and selecting the appropriate model structure based on the available data. RF regression is an efficient and credible alternative in prediction tasks, as reported by previous studies [12]. The results obtained in this study suggest that RF is as accurate as Poisson, zero-inflated Poisson, and negative binomial regressions when dealing with overdispersed count data. The tick burden dataset confirmed this statement. Novel methods such as RF are more efficient for predicting spatial patterns in species distribution models than more traditional linear models [58]. Furthermore, they consistently outperform more established methods. RF’s consistency and predictive performance could be an additional point to explore and reorient the debate on modelling count data. Some researchers suggest using the Poisson model and its extensions [57]. Others recommend using the linear model when the objective is to test the significance of estimates [41]. On the contrary, others think there is no best solution because it depends on the database in the presence [42]. However, these models have shown similar performance for different levels of overdispersion. Considering GLMs, RF is robust to misspecification of the relationship between the mean and variance [57]. Anderson et al. [33] claimed that traditional regression-based models may not be the most effective model, depending on the available census data (for example, presence of sampling bias, missing data, outliers, and imperfect detection). Given the predictive accuracy demonstrated by RF, its flexibility in dealing with covariates of various natures, and other recognised advantages such as the ability to deal with complex interaction structures, highly correlated variables and the measurement of variable significance, RF is an effective alternative method for predicting count outcomes [52]. RF may be a tool of capital importance for decision and policymakers working in limited data areas or with unreliable census data, especially in developing countries. Furthermore, producing accurate and reliable estimates is a problem of interest for resource allocation, public health, outbreak early detection, food security, conservation policy, environmental planning, and others [33, 56, 59].


## 4. Conclusion

RF has become a popular analysis tool in many application fields and will probably remain relevant in the future due to its high flexibility and many advantages. This study evaluated the impact of different characteristics of the dataset and the magnitude of overdispersion on the settings of RF parameters and predictive performance through simulated and real-life datasets. We found that the RF algorithm performs relatively well for all types of predictors. Data features do not influence the setting of RF parameters and, consequently, its predictive accuracy. However, the ratio of observations to the number of predictors influences the stability of the optimal RF parameters for a given dataset. Although the degree of overdispersion of the response variable did not significantly influence RF predictive validity, the magnitude of dispersion influenced RF predictive accuracy. The RF algorithm performed almost equally well as standard count data models such as quasi-Poisson, zero-inflated Poisson, and negative binomial models for the outcomes with linear and quadratic variance-to-mean relationships. However, for some scenarios, RF sometimes had a higher bias or RMSE, sometimes a lower bias or RMSE compared to Poisson family models. Furthermore, for the same scenario, depending on the selected subsample, a model may appear to perform better when it does not perform well overall. Thus, the best model choice should be based on data. RF is a reliable method when the assumptions for GLMs are not satisfied. Therefore, despite the sensitivity of RF to the overdispersion of the response variable, it is an efficient method for predicting count outcomes. Thereby, it is a tool of capital importance for decisions and policymakers working in limited data areas. Therefore, research should focus on decreasing the sensitivity of predictive accuracy to overdispersion.



---

*Source: 2833537-2022-12-01.xml*