Advancing Psychological Research With Random Forests: A Review of Methods,Tools,and Applications

Decision trees

Decision trees are nonparametric statistical-learning methods that automatically take care of complex nonlinear relations and interactions. They are supervised-learning techniques for predicting categorical outcomes (classification trees; e.g., depression status) or continuous outcomes (regression trees; e.g., work burnout score). The tree is constructed by recursive partitioning, a step-by-step procedure in which the predictor space (i.e., feature space) X is repeatedly partitioned into subspaces to form more homogeneous subgroups of observations for the purpose of prediction. Its structure resembles an upside-down tree, featuring a root at the top and branches with leaves growing downward, as illustrated in Figure 1a. This structural property is why they are referred to as “decision trees.” It first starts with the root node (i.e., the entire data set), ℜ 0 = X = { x ∈ ℝ d } ; at each subsequent step, a splitting variable X_j is selected along with a splitting criterion θ _j ; on this basis, the parent node is divided into nonoverlapping regions accordingly: ℜ child 1 = { x ∈ ℜ parent | X j ≤ θ j } and ℜ child 2 = { x ∈ ℜ parent | X j > θ j } , as illustrated in Figure 1b. These subregions, also called “child nodes,” are expected to be more homogeneous in terms of the outcome Y. The most commonly used partitioning algorithm is binary recursive partitioning, in which at each split, the parent node is being divided strictly into two child nodes. The splitting stops when it hits a predetermined stopping criterion (e.g., each subregion should contain no fewer than 10 observations). The regions or nodes that are not further split are called “terminal nodes” ( ℜ t ). All terminal nodes altogether define a partition of the data (see the dashed lines in Fig. 1b).

OPEN IN VIEWER

Fig. 1.

A hypothetical conceptual example of classification tree. This tree has a depth of 3: From the root node to ℜ 1 or ℜ 2 , the path goes through three splits. (a) A classification tree for predicting depression status. (b) The corresponding partition of the predictor space.

Following the partition, each observation in the data should fall into one and only one terminal node. With a single tree, each terminal node is assigned a particular predicted outcome value, Y ^ t . For categorical outcomes, the predicted value is usually computed as the majority vote in the terminal node; for continuous outcomes, the predicted value is usually computed as the average in each terminal node. Therefore, with a single tree, all observations that fall into the same terminal node will have the same predicted value. On the other hand, it is also possible for different terminal nodes to have identical predicted value, and thus, observations ending up in different terminal nodes can be predicted to have the same outcome.

In the following sections, we mainly focus on one of the most widely adopted decision-tree algorithms—classification and regression trees (CARTs; Breiman et al., 1984)—as the base learner of RFs (methodological limitations are briefly discussed in the Discussion section). Note that other alternative decision-tree algorithms also exist, such as the C4.5 classifier (Quinlan, 2014) and conditional inference trees (CTrees; Hothorn et al., 2006). These algorithms share the fundamental decision-tree framework with CART but differ in certain technical aspects, including the criterion for selecting a splitting variable and the internal procedure for handling missing data. Some of these alternatives were reported to better address the inherent biases of CART (Hothorn et al., 2006; Strobl et al., 2008, 2009). Nonetheless, CART still remains the most commonly used base learner in current software implementations of RFs. For a more thorough treatment of the tree-based methods and different algorithms, readers can refer to Mienye and Jere (2024) and Strobl et al. (2009).

Figure 1 provides a simple hypothetical example of classification trees built with CART. Here, the outcome of predictive interest is the depression status (depressed vs. not depressed). The two covariates—sleep quality and stress level—are all measured on a 5-point Likert scale; high values indicate higher levels in the corresponding measure. The tree starts growing from the entire sample space. At each step, the CART algorithm finds the best variable to serve as the splitting variable by considering every possible covariate and every possible split for each covariate. All these possible splits are evaluated and compared simultaneously. The one single split that yields the best fit is then selected for further growing the tree. For classification tree, a possible way to evaluate the fit is to compute the node-impurity measure, such as entropy

E t = − p log p − ( 1 − p ) log ( 1 − p )

(1)

or Gini index

Gini t = 2 p ( 1 − p ) ,

(2)

where p = P(Y = depressed) is the observed proportion of depressed observations in node t. The goodness of a possible split can therefore be defined as the reduction in the chosen node-impurity measure after splitting the parent node. The split that leads to the largest reduction in node impurity (i.e., the largest improvement in node purity) is selected. ¹ In this example, for instance, sleep quality is chosen as the first splitting variable, and the split occurs at a value of 3 ( { x ∈ ℜ 0 | X SleepQuality < 3 } and { x ∈ ℜ 0 | X SleepQuality ≥ 3 } ), which means this particular split at the root node yields the best result overall compared with other possible splits across all variables and all parameter spaces.

Now suppose the tree algorithm goes through the recursive-partitioning process and renders the tree model shown in Figure 1 as the best predictive model. After partitioning, the feature space is divided into a total of five terminal nodes ( ℜ 1 − ℜ 5 ; each contains 36%, 8%, 14%, 7%, and 35% of the observations, respectively). Recall that the goal of recursive partitioning is to create more homogeneous subgroups after splitting. In this example, 42% of the observations are diagnosed as “depressed,” and 58% of the observations are “not depressed” at the root node, which translates into a Gini index of 0.487. After the initial split at X SleepQuality = 3 , 62% of the observations in the left child node are reported depressed, yielding a reduced Gini index of 0.471; in the right child node, only 14% of the observations are reported depressed, resulting in a greatly reduced Gini index of 0.241. Following further recursive partitioning, the node-purity measure is greatly improved, with 100%, 0%, 0%, 84%, and 0% of observations diagnosed as “depressed” in each terminal node, respectively. This partitioning results in zero variation in the observed outcome across most terminal nodes. ²

In this classification tree, the predicted outcome is the majority depression status observed among the training cases in each terminal node. To predict the depression status for a new individual, we simply follow this decision tree to send them down through the branches, guided by the covariates’ values accordingly. In this example, to predict an individual’s depression status, the first question to ask is whether this person has a sleep quality score below 3. If yes, the person moves down the left branch; if not, the person moves to the right. Once in the left child node, the next question is if the stress level is above 3. An answer of yes sends this individual further down to the left branch, ending up in terminal node ℜ 1 ; in contrast, an answer of no sends the person to the right branch, reaching terminal node ℜ 2 instead. If this individual is ultimately sent to terminal node ℜ 2 , for instance, given all training observations in that node are not depressed (100% not depressed), our prediction for this new individual should be “not depressed.” Consider another individual with a sleep-quality score of 1 and a stress level of 4. Following the same logic, this person will be sent down to terminal node ℜ 1 , where all training observations are depressed. Consequently, the second individual is predicted to be depressed.

How does one evaluate the prediction accuracy of a decision tree? With decision trees, the prediction accuracy is usually assessed by examining the model’s performance on a separate test set or through cross-validation, which involves partitioning the data set into multiple subsets and repeatedly training the model on some while validating it on the remaining ones (see more in Browne, 2000; Koul et al., 2018; Stone, 1974). For classification trees, a standard measure of predictive accuracy is the misclassification rate, which quantifies the proportion of incorrectly classified observations in the entire test set:

1 N test ∑ i = 1 N test I ( Y i ^ ≠ Y i ) ,

(3)

where Y ^ i is the final predicted class label for individual i and I ( Y ^ i = Y i ) is an indicator function that takes a value of 1 if the prediction is correct and 0 otherwise. N test is the total number of observations in the test set. Building increasingly complex trees is expected to effectively reduce misclassification in the training data. However, highly complex predictive models are also known to have large variance because they are sensitive to even small random fluctuations in the sample data. This is famously known as the bias-variance trade-off. Overly complex trees tend to capture not only the meaningful data patterns but also the random noises, leading to poor model performance when presented with new observations. Therefore, with a single tree, to reduce the risk of overfitting, it is important to assess a model’s predictive accuracy on data not used for growing the tree. With a single tree, it is a common strategy to first grow a large (i.e., bushy) tree and then prune off the branches, which helps address the issue of overfitting. A pruned tree is a subtree of the large bushy tree by collapsing its nodes (i.e., undo the splitting) from bottom up toward the root node. This is usually done through cost-complexity pruning, which aims to balance model fit against complexity by penalizing complex models, combined with cross-validation to select the optimal level of pruning.

Extending the discussion to regression trees follows a similar logic except that the outcome is continuous and the measures of goodness are naturally different. To grow a tree rather than using measures such as entropy or Gini index to evaluate the possible splits, regression trees instead minimize the variability of the responses in each node by computing the sum of squared errors (SSE):

SSE t = ∑ i ∈ ℜ t ( Y i − Y ¯ . t ) 2 ,

(4)

where Y ¯ . t = 1 n t ∑ i ∈ ℜ t Y i is the average of the observed outcome values within node t. SSE _t is thus the squared difference between the observed outcome and the predicted outcome, summed over all observations in each node t (n_t denotes the total number of observations falling into node t). The split that yields the minimum overall SSE is then selected. To evaluate the prediction accuracy with test data, the test SSE can be similarly computed:

SSE test = ∑ i N test ( Y i − Y ^ i ) 2 .

(5)

Ensemble methods: bagging and RFs

Individual decision trees usually suffer from instability, which means they are sensitive to even small changes or outliers in data, especially for the large complex models that are overfit to the training data. That is to say, partitioning based on a slightly different sample can yield an entirely different tree model and thus very different predicted results for the same observations. Therefore, a single tree is unstable with high variance. Such instability is generally undesirable for predictive tasks. To address this problem, ensemble methods have been developed, such as bagging (bootstrap aggregating; Breiman, 1996) and RFs (Breiman, 2001). ³ These approaches explicitly leverage the inherent instability of individual trees to construct a robust collective committee of models. By perturbing the training data, typically through bootstrap sampling, multiple predictive models are generated. These predictions are then aggregated, often through averaging or majority voting, to produce a more stable and accurate ensemble predictor. Ensemble methods are useful because they take the average across multiple trees, which can stabilize the predictions and reduce the impacts of random noises.

Let the training data be denoted as D train = { X i , Y i } of size n, i = 1, 2, . . ., n. The bagging procedure draws B bootstrap samples D * = { X i * b , Y i * b | b = 1 , 2 , … , B } from D_train, typically by repeated sampling with replacement. ⁴ Each bootstrap sample often contains n observations, thus preserving the structure of the original data. Suppose Y_i is a binary outcome (e.g., depression status). For each bootstrap sample D * b = { X i * b , Y i * b } , a classification tree T^*b is independently grown, thereby generating B classification trees in total. Unlike the single-tree scenario, each tree is usually grown without pruning. Theoretically, the bootstrapped sample is expected to reflect the same data-generating process as the original data. But the procedure of bootstrapping introduces random variations between each data set, thus leading to variations among the trees grown (Fig. 2). To predict the outcome for a new observation X_test, one needs to pass the observation through each classification tree T^*b, sending it down to a unique terminal node in each tree. Therefore, this new observation will obtain B predicted values regarding outcome Y. For instance, consider a hypothetical new observation with X_SleepQuality = 2, X_{PhysicalActivity} = 2, X_Stress = 3, and X_{SocialSupport} = 2. With the first bootstrap classification tree T^*1(Fig. 2a), the prediction is Y ^ * 1 = depressed ; the second bootstrap classification tree T^*2 (Fig. 2b) predicts Y ^ * 2 = not depressed , so on and so forth. After collecting the prediction results from all the B classification trees, one simply counts the number of “votes” each outcome category receives. The final predicted class for the test case X_test is determined by the majority vote across the B trees. For instance, if 100 bootstrapped samples and classification trees are grown and 73 of them predict Y ^ test = depressed , then the final prediction outcome from this bagging procedure is “depressed.”

OPEN IN VIEWER

Fig. 2.

A hypothetical conceptual example of bagging. (a) The first bootstrap classification tree T^*1; (b) The second bootstrap classification tree T^*2; (c) The third bootstrap classification tree T^*3.

Bagging with regression trees follows a similar logic to that of classification trees. Other than the difference of evaluation measures, a key distinction between bagging regression trees and classification trees lies in how predictions are synthesized. While classification trees rely on a majority vote to determine the final outcome, regression trees aggregate predictions by computing the average of the outcome values predicted by each of the B individual trees:

Y ^ bag ( X test ) = 1 B ∑ b = 1 B Y ^ * b ( X test ) .

(6)

A methodological extension of bagging is RFs (Breiman, 2001). In bagging, randomness is introduced through the process of bootstrap sampling from the original data. Ideally, the individual trees grown out of these bootstrapped samples should be independent. Thus, bagging can effectively reduce the variance by aggregating over independent predictions while lowering bias by growing large, bushy individual trees. In practice, however, these tree models can be correlated, undermining the variance-reduction effect of bagging. To further de-correlate the individual trees, RFs introduce another layer of randomness on top of bagging. This typically involves randomly selecting a subset of features at each split in the individual trees (Ho, 1998).

The algorithm of RFs closely resembles that of bagging. It begins by drawing B bootstrap samples from the original training set via repeated sampling. But in each bootstrap sample, the way to grow a classification tree or regression tree is slightly different than in bagging. More specifically, at each candidate split, rather than evaluating all p available covariates, the algorithm randomly selects a subset of size m (m < p; when m = p, RF is simply equivalent to bagging). Among these m covariates, the RF algorithm identifies the optimal split based on a node-impurity measure, such as entropy or Gini index for classification or SSE for regression task. The process of random feature selection and recursive partitioning continues until a predefined stopping criterion is reached. Repeating this tree-growing process for each of the B bootstrap samples yields an ensemble of B classification or regression trees, which collectively defines an RF. Predictions for new observations follow the same rules as in bagging: majority voting for classification and averaging predictions for regression. In contrast to a single regression tree that assigns identical predicted outcome to all observations in the same terminal node, the averaged predictions in bagging and RFs typically vary from one observation to another even when some trees can place them in the same leaf. This results in a more finely grained prediction surface.

Assessing the variable importance (VI) is an important and common consideration in implementing RFs. Beyond obtaining prediction results, researchers are often interested in knowing which variables contribute to the classification or which variables are most influential in explaining the response variable. Although tree-based algorithms are mainly used for prediction purpose and are data-driven and exploratory in nature, evaluating the VI can help with feature selection for subsequent analyses and can potentially contribute to future theory development (e.g., Brick et al., 2018). Nonetheless, the evaluation of VI is already challenging in a single-tree model given the nonlinearity and complex interactions; it becomes even less intuitive with ensembles of trees, particularly compared with parametric models that most researchers are more familiar with (e.g., linear regression). One widely used approach is briefly introduced here. For RFs, in each individual tree T^*b, there is an out-of-bag (OOB) sample, defined as the collection of observations that are not selected into the bth bootstrap sample: D OOB * b = { X * , Y * | X * , Y * ∉ D * b } . These OOB observations are dropped down the bth tree T^*b to obtain predicted values. The prediction error (e.g., misclassification rate for classification and SSE for regression) can thus be computed for this OOB sample, PE OOB * b . Next, for a given covariate X_j, its observed values are randomly permuted across all the observations in this OOB sample D OOB * b while the rest of the data remain unchanged. This updated OOB data set, in which X_j has been randomly permuted, is denoted as PD _j ^*b; it is then passed through the bth tree again to generate a new set of predictions. The new prediction error is computed as PE OOB j * b , which is expected to be larger than the original OOB error PE OOB * b without permutations. The rationale is if a variable is important for making correct predictions, randomly permutating its values in the sample will greatly reduce the prediction accuracy. A larger decrease in model accuracy thus indicates a more important predictor. Therefore, the importance measure of covariate X_j in the bth tree can be defined as

VI j * b = PE OOB j * b − PE OOB * b .

(7)

This process is repeated in each of the B trees for covariate X_j. Finally, aggregating these importance measures across the B trees yields a total importance measure for X_j:

V I j = 1 B ∑ b = 1 B V I j * b .

(8)

Although it seems computationally cumbersome, this method has been conventially integrated into the RF algorithm, and importance scores are often automatically computed for all variables in software programs. Figure 3 provides a hypothetical example of the permutation-based VI measures, computed as the increase in misclassification rate, with variables ordered from more important to less important.

OPEN IN VIEWER

Fig. 3.

Variable-importance plot from random forests for the hypothetical example predicting depression.

The RF algorithm involves three main tuning parameters (also referred to as “hyperparameters”) that researchers need to set values for growing a forest: the tree depth (i.e., the number of splits along the longest path from the root node to a terminal node), the number of randomly selected covariates at each potential split (m), and the number of bootstrap samples (or the number of trees grown), B. To ensure the effectiveness of RFs, the hyperparameters’ values have to be carefully specified. Some common suggestions recommend m = p for classification and m = p / 3 for regression (Hastie et al., 2009, p. 592), but the optimal values really depend on the specific context and data set. Hyperparameters can be empirically determined using cross-validation, although this process can become computationally expensive. A more detailed discussion on hyperparameter tuning is provided in the next section.

Practical issues in applying RFs

Although decision trees and RFs have gained an increasing popularity in psychological research, the effectiveness of RFs, which were originally developed outside psychological-research contexts, has not yet been thoroughly investigated. Important questions remain about its applicability and how it can be effectively used to address the practical challenges that are more specific in psychology studies. In applying RFs, psychologists need to navigate through a series of methodological considerations that can substantially influence the effectiveness of the RF models and the interpretation of results.

Goals of implementing RFs

To begin, the purpose for using RFs can vary widely (Probst et al., 2019; Shmueli, 2010; Yarkoni & Westfall, 2017). Many researchers employ RFs as a predictive tool aiming to enhance accuracy in predicting future outcomes. Some predominantly rely on RFs because of its strength and flexibility in handling high-dimensional data in which higher-order nonlinearity and complex interactions are anticipated. Alternatively, others may choose to include RFs as part of comparison between various predictive models, ranging from traditional parametric models (e.g., linear regression or logistic regression) to nonparametric models (e.g., bagging, boosting), to empirically determine which method offers the best predictive accuracy. A third application of RFs is in preliminary variable selection, in which they assess the relative importance of predictors in explaining the outcome variable rather than just serving as a predictive tool. It usually involves using RFs to identify important covariates that will be further examined in inferential statistical modeling (and sometimes predictive modeling as well). Regardless of whether to predict, to explain, to select features, or to achieve multiple goals, RFs’ performance and applicability need to be assessed accordingly depending on the specific research purposes (Shmueli, 2010; Yarkoni & Westfall, 2017).

Data dimensions

Given that psychological research often concerns a very specific target population and in-person data collection, the data conditions in psychological studies can be quite different from those in other fields. Some common challenges include small sample size, data dimensionality, and missing data. Psychologists often need to work with small data sets with limited number of observations, particularly when data collection is costly. With traditional inferential statistics, small sample sizes can be challenging because they can introduce larger bias and reduce estimation precision. Does it also affect the performance of RFs as a predictive model or as a variable-selection method? While RFs are believed to be applicable in “small n large p” conditions (Fife & D’Onofrio, 2023; Matsuki et al., 2016; Strobl et al., 2009), potential challenges may still exist, particularly in terms of model overfitting (e.g., Yarkoni & Westfall, 2017). Recall that RFs fundamentally rely on bootstrap sampling. Each individual tree is trained on a bootstrap sample, which, on average, is expected to contain approximately 63% of the original data (if sampling with replacement). For small data sets, this can substantially reduce the effective training-set size for individual trees. With a limited training set, individual trees can be overfitted, especially if the trees are allowed to grow deep and large. This overfitting can be further exacerbated if irrelevant predictors are included in the training process without prescreening. For example, consider a small study with 30 participants, 10 of whom are clinically depressed. Suppose by pure chance, all of the depressed participants use iPhones to fill out the survey and some of the nondepressed participants use Android phones. With such a limited sample size, the individual tree may try to capture this random data pattern (e.g., suggests using iPhone is an important predictor for depression), rendering the results less generalizable. Unfortunately, sample size is not typically carefully evaluated and justified when predictive models are employed (Dhiman et al., 2023), and there are no clear guidelines regarding what should be considered a good sample size or how to plan for a sufficient sample size with predictive modeling. Some research has shown that small sample size can lead to larger bias when certain a machine-learning algorithm is applied (e.g., support vector machine; Vabalas et al., 2019), but the methodological implications of small n on RFs in psychology studies remains unclear.

On the other hand, it is also not yet clear whether having a small number of predictors will affect the effectiveness of RFs. With high-dimensional data, the key methodological advantage of RFs over bagging is that they de-correlate the trees by randomly selecting a subset of predictors at each split. With the total number of predictors being limited, this advantage may diminish. The same predictors may frequently appear across multiple trees, making the ensemble less independent, potentially affecting the overall performance of this forest. This consideration can be illustrated through a simple example. To make it more intuitive, we ran a simple simulation as a proof as concept: From p predictors, we randomly drew two independent sets of m candidate variables (without replacement) and counted the number of shared predictors; this process was repeated 1,000 times. With p = 10, the average overlap across replications was about 0.38, 2.50, and 6.39 when m = 2, 5, and 8, respectively; that is, with m/p = 0.5, roughly half of the five candidate predictors are expected to be shared by two trees, and with m/p = 0.8, about 80% of them are expected to be identical. In contrast, when p = 50, the average overlap between two independent draws dropped to about 0.09, 0.49, and 1.29 for m = 2, 5, and 8, respectively. Thus, the potential overlap between each tree is substantially reduced. ⁵ This example shows that when the total number of available predictors p is small, the m/p ratio can grow big easily with even small changes in m, causing the trees to consider largely the same variables at each split and thus become more correlated, undermining the de-correlation property of RFs.

Imbalanced data

In psychology studies, imbalanced data sets are common, particularly when the outcome of interest is relatively rare or underrepresented in the population. For example, in mental-health research, the number of children with suicide attempts (e.g., see Harman et al., 2021) is substantially smaller than the number of healthy control subjects. With such class imbalance, the RF algorithm can be biased because it will prioritize the majority class, resulting in poor prediction of the minority class. This can be problematic because the minority class is more often the focal research interest (e.g., identifying at-risk children for potential suicide).

As a consequence, when the outcome is highly imbalanced, the overall misclassification rate, or accuracy, can be a very misleading indicator of model performance because it is dominated by the majority class. For example, if 95% of the training cases are individuals without suicide attempts, a model that predicts all cases as “no suicide attempts” can achieve a 95% accuracy (a low misclassification rate = 5%), but the model itself is completely useless for identifying individuals at risk (both sensitivity and precision are zero ⁶ ). In such scenarios, alternative metrics that assess the model performance in each specific outcome class are preferable to better align with the intended research purposes, such as sensitivity (i.e., recall), specificity, ⁷ precision, and F1 score. ⁸ In addition, visual inspections can be informative as well. The precision-recall curve, which is conceptually similar to the receiver operating characteristic curve but more appropriate with imbalanced data, can be used to assess the model performance across varying decision thresholds.

To address class imbalance, some resampling techniques can be applied (Chawla et al., 2002; Japkowicz, 2000; Ling & Li, 1998). One possible approach is to randomly oversample the minority class or randomly undersample the majority class (Japkowicz, 2000; Ling & Li, 1998). Random oversampling balances the data set by duplicating observations from the minority class until the resampled minority class consists of as many data points as the majority class; in contrast, random undersampling balances the data set by sampling only a smaller proportion from the majority class until their size matches the minority class. As a more sophisticated alternative, synthetic minority oversampling technique (SMOTE; Chawla et al., 2002) further perturbs the data by randomly generating synthetic minority class samples through interpolation between each minority unit and its nearest neighbors in the feature space, which can be identified using a distance metric, such as Euclidean distance, until the sample size of the minority class and majority class size are equalized.

Missing data

Missing data are also typical in psychological studies involving human subjects because of participant nonresponse, dropouts, or other practical issues. There are different strategies for handling missing data in RFs, each with its own advantages and disadvantages. Listwise deletion is the most widely used approach, which simply discards all observations with incomplete data before fitting RF models. However, listwise deletion can be a huge waste of available information and can be infeasible if the sample size is small to begin with.

Among other more modern missing-data-handling methods, by far the most popular approach in CART is surrogate splits (Hapfelmeier et al., 2014; Hothorn et al., 2006). With this approach, missing data are handled internally in each tree. Surrogates are defined locally at each split; when the best splitting variable is selected, other candidate splitting variables that best mimic the current optimal splitting result are ranked and labeled as the surrogates. When the primary splitting predictor is missing for a given observation, RF algorithm instead uses the best surrogate available to assign this observation further down through the tree. The methodological challenges of surrogate splits have been noted, however, given the computational burden, and thus, many other approaches have been proposed (Tang & Ishwaran, 2017). Instead of using surrogate splits, some software (e.g., the ranger R package) implements an algorithm in which missing values are initially ignored when evaluating potential splits. Once the best split is determined, observations with missing values on the splitting variable are temporarily assigned to both child nodes, and the assignment that optimizes the split results is chosen. This optimal assignment is then stored as the default assignment so that when new observations (i.e., future data points) have missing values for that variable come in, they are directed to the chosen child node during prediction. Alternatively, some other software (e.g., the partykit R package) offers an option to treat the missing value as a unique value/category in its own right. This approach is referred to as “missingness incorporated in attributes” (MIA; Twala et al., 2008), which treats missing status as informative and allows it to predict the outcome. All these approaches handle missingness internally without the need to remove or impute the missing data.

An alternative and common approach for handling missing data is to preimpute the missing data and fit the analytical RF model to the complete preimputed data. In some software implementation (e.g., the randomForest R package), the preimputation involves iterative steps that start with simple imputation in which the median or the mode of observed values is used to impute the missing variable. An RF is grown using this simple-imputed data, and a proximity matrix is calculated based on how frequently two cases share the same terminal node across the trees. The proximity matrix is then used to update the initially imputed values (e.g., the missing variable is imputed as a weighted mean of the observed values across all the other cases, weighted by the proximity measure). This process is iterated until convergence. Another available option for missing data before imputation is via the missForest algorithm (Stekhoven & Bühlmann, 2012), which treats missing data as a prediction problem. It grows an RF for each variable with missing values using all other variables as predictors; the trained forest is then used to predict the missing values. This process continues until convergence is achieved. A computationally faster version of missForest has also been proposed as mForest, which employs multivariate forests for imputation and thus reduces the number of RFs needed (Tang & Ishwaran, 2017).

Missing values can also be imputed adaptively during the analytic process (i.e., on the fly) rather than beforehand. For instance, the on-the-fly-imputation (Ishwaran et al., 2008; Tang & Ishwaran, 2017) algorithm uses only complete data to determine the best split at each step. Once a splitting variable is selected, for observations with missing data on this selected splitting variable, a random value is drawn from the nonmissing in-bag data to “impute” for this missing value. This observation is then assigned to a child node according to this temporarily imputed value. These temporary imputed values are discarded after this observation is passed down to a child node; thus, the missingness in this observation is preserved onward.

Although missing data have been studied extensively in conventional inferential-statistical modeling (e.g., Enders, 2022, 2025), it remains unclear which approach is best suitable for missing-data handling when implementing RFs in psychological studies. The optimal choice likely depends on the underlying missingness mechanism (missing completely at random, missing at random, missing not at random; Rubin, 1976), the missing data rate, and the specific structure of the data.

Software tools

Many common software, such as Python, R, SAS, MATLAB, and SPSS, can be used to implement the RF algorithms. Among them, R and Python are two of the most popular statistical tools and programming languages, both of which are open-source environments and offer various specialized packages for RFs.

In R, multiple stand-alone packages have been developed for different modeling engines. For instance, the classic RF algorithm (Breiman, 2001) that uses CART as the base learner is implemented in both the randomForest package (Liaw & Wiener, 2002) and the ranger package (Wright & Ziegler, 2017); ranger is a computationally faster option. Alternatively, conditional inference forests, which uses CTree as the base learner and select splits via statistical tests, are implemented in the partykit package and in its predecessor party package (Hothorn et al., 2006; Hothorn & Zeileis, 2015). On the other hand, there are also two wrapper packages, caret (Kuhn, 2008) and tidymodels (Kuhn & Wickham, 2020), that integrate a collection of R packages for machine-learning tasks. These wrapper packages in R do not have their own unique modeling engines for RFs; rather, they are designed to call up the external modeling engines in R ⁹ while facilitating a more unified workflow of model specification, data preprocessing, hyperparameter tuning, and model fitting. Parallelly, in Python, the Scikit- learn library implements the classic CART-based RF algorithm via the RandomForestClassifier and RandomForestRegressor functions (Pedregosa et al., 2011). For a summary of these popular packages, including their base learners, key function arguments, and default hyperparameter values, see Table 1.

OPEN IN VIEWER

Table 1.

Commonly Used Software Packages for Random Forests

Software package		R			Python
		randomForest	ranger	partykit	Scikit-learn
Key function name		randomForest ()	ranger()	cforest ()	Randomforestclassifier ()	RandomForestRegressor ()
Base learner algorithm		CART	CART	CTree	CART
Hyperparameter arguments and defaults	Number of trees	ntree = 500	num.trees = 500	ntree = 500	n_estimators = 100
	Number of variables randomly selected	mtry = sqrt(p); mtry = p/3	mtry = sqrt (p)	mtry = sqrt(p)	max_features = "sqrt"	max_features = n_features
	Sample with replacement	replace = TRUE	replace = TRUE	replace = FALSE
	Size of the bootstrap sample	sampsize = n	sample.fraction = 1	fraction = 0.632	max_samples = None
	Minimum size of terminal nodes	nodesize = 1; nodesize = 5	min.bucket = 1	minbucket = 7	min_samples_leaf = 1
	Maximum number of terminal nodes	maxnodes = NULL			max_leaf_nodes = None
	Minimal node size for split		min.node.size = 1; min.node.size = 5	minsplit = 20	min_samples_split = 2
	Maximal tree depth		max.depth = NULL	maxdepth = Inf	max_depth = None
Software package		R
		caret			tidymodels
Key function name		train(. . ., method = "rf", . . .)	train(. . ., method = "ranger", . . .)	train(. . ., method = "cforest", . . .)	rand_forest(. . .) %>%set_engine("randomForest")	rand_forest(. . .) %>% set_engine("ranger")	rand_forest(. . .) %>% set_engine("partykit")
Base learner algorithm		CART	CART	CTree	CART	CART	CTree
Hyperparameter arguments and defaults	Number of trees	ntree = 500	num.trees = 500	ntree = 500	trees = 500	trees = 500	trees = 500
	Number of variables randomly selected	By default, caret conducts a grid search with three levels for a better mtry value.			mtry = sqrt (p); mtry = p/3	mtry = sqrt (p)	mtry = sqrt (p)
	Sample with replacement	replace = TRUE	replace = TRUE	replace = TRUE	Can be set up in set_engine() accordingly. The argument names and default values are defined locally in each engine.
	Size of the bootstrap sample	sampsize = n	sample.fraction = 1	fraction = 0.632
	Minimum size of terminal nodes	nodesize = 1; nodesize = 5	min.bucket = 1	minbucket = 7
	Maximum number of terminal nodes	maxnodes = NULL
	Minimal node size for split		min.node.size = 1; min.node.size = 5	minsplit = 20	min_n = 10; min_n = 5	min_n = 10; min_n = 5	min_n = 20
	Maximal tree depth		max.depth = NULL	maxdepth = 0	Can be set up in set_engine() accordingly. The argument names and default values are defined locally within each engine.

Note: For cells with two default values, the one denoted with underline is the default for classification tasks, and the other is the default for regression tasks. CART = classification and regression tree; CTree = conditional inference tree.

Besides the differences outlined in Table 1 that are directly relevant to model building, these packages also vary in their approaches of handling data, which also affect the results. For example, the options available for dealing with imbalanced data can differ substantially from package to package. Typically, data processing and model fitting are separate steps performed sequentially. In the R environment, the two wrapper packages ( caret and tidymodels ) are designed to streamline the entire process; thereby, both provide convenient internal functions that can be used to oversample or undersample the training data to address the class imbalance (e.g., step_upsample , step_downsample , and step_smote in tidymodels; upSample and downSample in caret ). In contrast, other stand-alone R packages generally do not offer such a fully built-in functionality for handling imbalanced data explicitly before model fitting. Although some limited controls are possible in randomForest and ranger , such as stratified bootstrap sampling with equal sizes at the tree level (as an approximate of simple undersampling), these cannot sufficiently replace the resampling techniques described earlier. Therefore, researchers usually need to rely on external R packages to perform the data-resampling step. Examples of external packages that can be used for resampling include ROSE (Lunardon et al., 2014) and smotefamily (Siriseriwan, 2024). In Python, on the other hand, the imbalanced-learn library (Lemaître et al., 2017) provides useful functions for random oversampling ( RandomOverSampler ), random undersampling ( RandomUnderSampler ), and SMOTE ( SMOTE ) before the RF model is fit via Scikit-learn .

Packages also differ substantially in how they handle missing data. In terms of the stand-alone modeling engines, the randomForest R package cannot handle missing values during the analysis process, and therefore, users need to either remove missing data or impute missing data beforehand as a separate step before passing data to the model. It provides internal functions for this purpose, including the na.roughfix() function that does simple imputation for the outcome and predictors and the rfImpute() function that can iteratively impute the missing data in predictors based on proximity measures. But researchers can also employ any external software for missing-data imputation. For instance, a more advanced extension for RF-based missing-data imputation, the missForest algorithm, is implemented in the missForest R package (Stekhoven & Bühlmann, 2012). The ranger and partykit R packages, on the other hand, offer better internal functionalities for handling missing data on the fly without the need of preimputation or preprocessing. In ranger , this can be specified via the na.action = “na.learn” argument, which, as detailed in previous section, allows the RF algorithm to dynamically “learn” the optimal branch for observations to go down the tree when there are missing values in a splitting variable. In the partykit R package, surrogate splits and/or the MIA procedure can be enabled to handle missing data in predictor variables (by setting a positive integer for the maxsurrogate argument or setting MIA = TRUE , respectively). However, both ranger and partykit do not allow missing values in the outcome variable, so any observations with missing outcomes must be removed or imputed before model fitting.

When the wrapper R packages are used, researchers have three general options to handle missing data. They can choose to rely on the engine’s native functionalities as described above, use the wrapper’s preprocessing tools to process missing values directly in the pipeline, or impute the missing values using other preferred packages before model fitting. To begin with, in caret and tidymodels , missing data can be handled in the same way as in each specific modeling engine as long as the proper arguments aforementioned are passed to the corresponding modeling engine. ¹⁰ Alternatively, both caret and tidymodels provide their own built-in functionalities in preprocessing missing data (e.g., preProcess() in caret; step_impute_ and step_unknown() functions in tidymodels ).

Finally, in Python, starting from Version 1.40 of the Scikit-learn library, the RandomForestClassifier and RandomForestRegressor functions also support internal missing-data handling when the predictor variables contain missing data. The underlying procedure is similar to what is implemented in the ranger package with na.learn , in which missing data are handled on the fly. Likewise, missing data in outcome variables are not allowed.

Hyperparameter tuning

After the data are properly preprocessed and loaded into the software of choice, the next step in implementing the RF algorithm is to determine the values of the hyperparameters. The current development of software usually provides some default values of the hyperparameters. Table 1 provides a summary of the hyperparameter options and default values used in the stand-alone modeling engine packages ( randomForest , ranger , partykit ), the wrapper R packages ( caret and tidymodels ), and Python’s Scikit-learn library. If the users do not explicitly specify any hyperparameter values, it will automatically use the software default setting. However, the optimal hyperparameter values are inherently data dependent, and most likely, the default setup will not yield the best performance of RFs. A more robust approach is to tune the hyperparameters either by manual cross-validation or via the automated optimization procedures. But at the end of the day, the choice has to be carefully considered because the value of hyperparameters can largely influence the performance of RFs.

Regarding the number of randomly selected covariates for candidate splitting (m), although the default setting is a convenient option (typically m = p for classification and m = p / 3 for regression), it is by no means universally optimal. In certain cases, for example, if most of the variables in the data are relevant and important for predicting the outcome, setting a smaller m can benefit the model prediction such that not all trees are dominated by the strongest predictor; on the contrary, if there are many irrelevant variables in the data, using a larger m can improve the prediction results by ensuring that at least some important predictors are included in each tree (Probst et al., 2019).

Tree complexity is another important hyperparameter in decision trees and RFs, which is usually controlled by setting the terminal-node size. Setting a smaller terminal-node size produces a larger tree with more splits, whereas a larger terminal-node size will effectively limit the tree depth. In many software, the terminal-node size is set to 1 for classification and 5 for regression by default. In addition to terminal-node size, other hyperparameters can also be used to control the tree complexity. These typically include the smallest node size required for a possible further split ¹¹ (e.g., minsplit in partykit R package), the total number of terminal nodes in a tree (e.g., maxnodes in randomForest R package), or more straightforwardly, the maximum depth of a tree (e.g., maxdepth in partykit R package). All these hyperparameters can be specified in combination simultaneously for an RF model, and the tree will be stopped from further splitting if any of their conditions are met.

The number of trees is also a key factor to consider in growing an RF. Most software implementations default to 500 trees. But unlike other tuning parameters, increasing the number of trees does not generally lead to overfitting. In fact, research suggests that it is preferable to set it to a reasonably large value for the best predictive performance (Probst & Boulesteix, 2018) as long as the computational resource allows. However, once the prediction performance stabilizes, adding more trees will have diminishing returns in improving the predictive accuracy. ¹²

The performance of RF models is highly dependent on the specification of hyperparameters, and therefore, to achieve the optimal results, they must be carefully tuned to adapt to the data and the modeling context. For example, having a smaller data with smaller n may necessitate fitting fewer trees in a forest, enforcing a larger size of terminal nodes, and limiting the depth of each tree to avoid overfitting. On the other hand, having a smaller number of predictors may require the users to increase the number of randomly selected variables at each step when growing an RF. In the presence of missing data, it is also crucial to confirm whether the chosen software supports internal missing-data handling, and if so, researchers may need to specify additional tuning parameters as needed (e.g., the maximum number of surrogates). Finally, all tuning parameters can affect the computational time, which can be an important practical consideration for researchers. Employing a smaller m, limiting the tree depth, and growing fewer trees in an RF can reduce the computational burden.

With all that, are there any principled ways to find the optimal values of the hyperparameters? Several different strategies can be used for hyperparameter tuning (Owen, 2022). The most straightforward approach is manual tuning, in which researchers manually change the value of one or more tuning parameters and check which value results in the most accurate prediction, usually via cross-validation or OOB error. An extension of manual search is grid search, which automatically loops through all the possible hyperparameter-value combinations in a predefined search space set by the researcher (e.g., tree depths of 5, 10, 15 combined with 50, 100, 200 trees, yielding nine combinations to test). Essentially, grid search is an automated version of manual hyperparameter tuning by iterating over each combination using nested loops. Although both manual search and grid search are conceptually straightforward and easy to implement, they can be time-consuming and rely on researchers to provide reasonably good candidate values to do the comparison in the first place. As the number of hyperparameters and possible values increase, they can become computationally expensive very quickly. Extending the grid search, a more efficient option for automated hyperparameter tuning is random search, which randomly picks hyperparameter values from given probability distributions (e.g., a uniform distribution) rather than relying on any specific user-supplied values. At each iteration, a value for each tuned hyperparameter is randomly and independently picked. Random search thus requires less prior knowledge from the users about the hyperparameter values and is more computationally efficient because it does not test all possible combinations of the values. Random search is also reported to frequently outperform a basic grid search when dealing with many hyperparameters or wide ranges of values (Owen, 2022). However, researchers still have to manually define the total number of iterations (or the number of randomly sampled hyperparameter values) when implementing random search.

Manual search, grid search, and random search can all be considered as exhaustive search strategies. In addition to these exhaustive search approaches, another well-established and successful tunning strategy is sequential model-based optimization (SMBO; Jones et al., 1998), also referred to as Bayesian optimization (BO), which employs a more adaptive strategy through iterations such that the next iteration is informed by previous iterations. In its implementation, SMBO starts by drawing several random values from the hyperparameter space and evaluating the RF performance accordingly. A surrogate model (a probabilistic regression model, such as Gaussian process) is fit to these initial training results, roughly assessing how changes in hyperparameters affect the model’s prediction accuracy. It then proposes the next set of hyperparameter values within the predefined hyperparameter space, where the proposed values have the best expected prediction result under the current surrogate model. The proposed hyperparameter values are empirically evaluated for prediction accuracy on the training data, and this new training result is added to the previous training results to further update the surrogate model. This process repeats iteratively to find the optimal hyperparameter values (for a more detailed and intuitive illustration of SMBO, see Appendix A in the Supplemental Material available online). SMBO is more computationally efficient than the other approaches, making it particularly appealing when there are a large number of possible hyperparameter configurations. It does, however, require the researchers to have a better statistical knowledge of the process to properly implement the procedure.

For a summary of the hyperparameter-tuning options available across the various packages in R and Python, see Table 2. They differ in terms of the available tuning options and the flexibility of automated hyperparameter search. Of all the packages reviewed here, only tidymodels offers native SMBO/BO options. When using other packages, researchers have to employ external packages or libraries that specialize in these optimization techniques. For example, the randomForest , ranger , ¹³ party , and caret R packages can be used along with the mlrMBO package (Bischl et al., 2018) and mlr package (Bischl et al., 2016) to perform SMBO/BO. Likewise, the scikit-optimize ( skopt ; Head et al., 2021) and Optuna library (Akiba et al., 2019) in Python also offer automated hyperparameter optimization via SMBO that can apply to scikit-learn ’s RF models.

OPEN IN VIEWER

Table 2.

Hyperparameter Tuning for Random Forest in Commonly Used Software Packages

		Hyperparameter-tuning options
Software package		Manual search	Grid search	Random search	Internal SMBO/BO	External SMBO/BO
R	randomForest	Users can manually loop over different values of hyperparameters and compare the OOB error or other evaluation measures	The helper function tuneRF() can be used to find an optimal value for mtry.	x	x	tuneParams() in mlrMBO package; tuneRanger() in tuneRanger package (for ranger package only)
	ranger		x	x	x
	partykit		x	x	x
	caret		User can call up the grid-search procedure by specifying the tuneGrid or tuneLength argument in train() function; more detailed control can be specified in the trControl argument.	User can call up the random-search procedure by specifying search = “random” in the trControl argument.	x
	tidymodels		Users can set relevant hyperparameter arguments to =tune() in rand_forest(. . .) and conduct grid search using the tune_grid() function, with search grid defined by grid_regular().	Users can set relevant hyperparameter arguments to =tune() in rand_forest( . . . ) and conduct random search using the tune_grid()function, with search grid defined by grid_random().	Users can set relevant hyperparameter arguments to =tune() in rand_forest( . . . ) and conduct BO using the tune_bayes()function.
Python	Scikit-learn		GridSearchCV()	RandomizedSearchCV()	x	BayesSearchCV() from scikit-optimize library; study.optimize() from Optuna library.

Note: x = not internally supported; SMBO = sequential model-based optimization; BO = Bayesian optimization . OOB = out-of-bag.