Learn to Build Random Forest Models in R: A Step-by-Step Tutorial

When data scientists encounter complex modeling challenges where the relationship between a set of predictor features and a response variable is highly non-linear and intricate, conventional statistical methods often prove insufficient. These demanding scenarios necessitate the deployment of advanced non-linear techniques capable of robustly capturing underlying data patterns and interactions.

A foundational technique in the realm of predictive modeling is the decision tree. While simple and intuitive, relying on a solitary decision tree model introduces a critical vulnerability: it is acutely sensitive to high variance. This inherent volatility implies that minimal perturbations in the training dataset can result in dramatically different tree structures and, consequently, highly unstable predictions.

To fully grasp this instability, one might consider training two separate decision trees on slightly varying subsets of the same data; the resulting models and their performance characteristics are likely to diverge substantially. This lack of predictive stability mandates a move toward more resilient methodologies designed to smooth out these fluctuations.

The standard solution is to leverage ensemble learning, utilizing the powerful structure of a random forest model. Random forests are engineered precisely to mitigate high variance by aggregating the results from numerous individual decision trees. This collective approach yields significantly more stable, accurate, and reliable predictions than any single tree could possibly achieve.

Understanding Random Forests: The Core Mechanism

A random forest is fundamentally an ensemble learning method that operates by systematically constructing a multitude of decorrelated decision trees during the training phase. For a classification task, the final prediction is determined by the mode (most frequent class) of the individual trees, while for regression tasks, the prediction is the mean of the individual tree outputs.

The remarkable predictive power of random forests stems from two distinct mechanisms designed specifically to ensure that the individual trees within the forest are statistically independent and decorrelated. This forced diversity prevents them from making the same errors:

1. Bagging (Bootstrap Aggregating): The algorithm initiates the process by utilizing bootstrapped samples. This means that each constituent tree is constructed using a random sample of the original dataset, drawn with replacement. This ensures that each tree sees a slightly different view of the training data.

2. Feature Randomness: Crucially, when determining the optimal split at each node of an individual decision tree, the algorithm considers only a random subset of predictor variables, rather than evaluating the full set of features. Typically, this subset size ($m$) is set to the square root of the total number of predictors ($p$), forcing diversification and preventing a few overwhelmingly strong predictors from dominating every tree structure.

By averaging the predictions derived from these numerous, yet independent, trees, the random forest significantly reduces the detrimental effects of overall variance. This strategic reduction in variance is what ultimately results in superior performance when compared to single decision trees or less sophisticated standard bagged models.

Step 1: Setting up the R Environment and Data Preparation

To commence the process of building our random forest model within the R programming language, the first essential step is ensuring that all required packages are correctly installed and loaded. For this practical demonstration, we will rely predominantly on the widely used and authoritative randomForest package.

The following command is used to load the necessary library into the current R session:

library(randomForest)

Our example will utilize the readily available R dataset, airquality, which contains 153 daily measurements of air quality in New York. Before any model fitting can commence, it is mandatory practice to inspect the dataset’s structure (using str()) and rigorously handle any data quality issues, most commonly involving missing values.

#view structure of airquality dataset
str(airquality)

'data.frame':	153 obs. of  6 variables:
 $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
 $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
 $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
 $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
 $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
 $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

#find number of rows with missing values
sum(!complete.cases(airquality))

[1] 42

The structural inspection confirms that 42 rows within the dataset contain missing values (represented as NAs). Since random forests, by default, cannot process missing data directly, we must perform imputation. A straightforward and commonly accepted approach involves replacing the missing entries in each column with the respective column’s median value, thereby preserving the distribution’s central tendency while ensuring the data is complete:

#replace NAs with column medians
for(i in 1:ncol(airquality)) {
  airquality[ , i][is.na(airquality[ , i])] <- median(airquality[ , i], na.rm=TRUE)
}

Related: How to Impute Missing Values in R

Step 2: Fitting the Initial Regression Model

With the data successfully prepared and cleaned, we are now ready to fit the initial random forest model using the powerful randomForest() function available in the randomForest package. Our objective is to predict the Ozone concentration level (our response variable) based on all other variables in the dataset, succinctly expressed using the formula notation: Ozone ~ ..

Prior to executing the fitting process, it is standard practice to set a random seed (e.g., set.seed(1)). This crucial step ensures that the results of the stochastic processes inherent in random forest construction—namely, the selection of bootstrap samples and features—are fully reproducible across different computational runs.

#make this example reproducible
set.seed(1)

#fit the random forest model
model <- randomForest(
  formula = Ozone ~ .,
  data = airquality
)

#display fitted model summary
model

Call:
 randomForest(formula = Ozone ~ ., data = airquality) 
               Type of random forest: regression
                     Number of trees: 500
No. of variables tried at each split: 1

          Mean of squared residuals: 327.0914
                    % Var explained: 61

#find number of trees that produce lowest test MSE
which.min(model$mse)

[1] 82

#find RMSE of best model
sqrt(model$mse[which.min(model$mse)]) 

[1] 17.64392

By default, the function constructs 500 individual trees. The summary output provides essential diagnostic metrics, including the final out-of-bag (OOB) estimate of the Mean of squared residuals (327.0914) and the percentage of variance in the Ozone level explained by the model (61%). Further analysis of the model object reveals the specific number of trees required to achieve the minimum test error. In this case, the lowest test Mean Squared Error (MSE) was achieved with only 82 trees. This optimal configuration corresponds to a Root Mean Squared Error (RMSE) of 17.64392, which quantifies the average magnitude of the error in our Ozone predictions.

Step 3: Interpreting Model Diagnostics and Variable Importance

Visualizing the model’s performance trajectory is crucial for gaining insight into its stability and determining the appropriate size of the ensemble. We can plot the out-of-bag test MSE against the cumulative number of trees used to confirm the point at which the error stabilizes, indicating that adding more trees offers diminishing returns.

#plot the test MSE by number of trees
plot(model)

Next, interpreting which predictor variables contribute most significantly to the model’s overall predictive accuracy is vital for scientific interpretation and drawing actionable conclusions. The varImpPlot() function is used to generate a graphic display of variable importance, typically based on metrics like the mean decrease in accuracy or, for regression problems, the average increase in node purity.

#produce variable importance plot
varImpPlot(model) 

The resulting plot’s x-axis often represents the average increase in node purity, measured by the collective improvement in the sum of squared errors achieved when splitting on a specific variable across all trees. A close inspection of this visualization clearly identifies Wind as the most critical predictor for estimating Ozone levels, followed closely by Temp. This objective assessment of feature contribution is indispensable for understanding the underlying environmental factors influencing the response variable.

Step 4: Optimizing Model Parameters (Tuning)

Although the initial model employed the default settings (500 trees and a starting value of $m=1$ variable tried at each split), optimal predictive performance frequently requires fine-tuning the key model hyperparameters. The most critical hyperparameter in random forest construction is mtry, which dictates the number of variables randomly sampled and considered at each split point.

We can systematically search for the optimal mtry value using the tuneRF() function. This function employs an iterative process, testing different numbers of predictors until the out-of-bag estimated error rate demonstrates no further significant improvement, effectively automating the tuning process.

The primary parameters necessary for effective use of tuneRF() are defined as follows:

• ntreeTry: Defines the number of trees constructed at each step during the search process.

• mtryStart: Specifies the initial number of predictor variables to be evaluated when the tuning begins.

• stepFactor: Determines the factor by which the mtry value is multiplied in subsequent iterative steps.

• improve: Sets the minimum required percentage improvement in the OOB error necessary to justify continuing the tuning search to the next step.

model_tuned <- tuneRF(
               x=airquality[,-1], #define predictor variables (all columns except Ozone)
               y=airquality$Ozone, #define response variable
               ntreeTry=500,
               mtryStart=4, 
               stepFactor=1.5,
               improve=0.01,
               trace=FALSE #don't show real-time progress
               )

The execution of this function generates a diagnostic plot that maps the number of predictors considered (mtry) against the corresponding Out-of-Bag (OOB) estimated error:

The resulting graph clearly illustrates that the lowest OOB error is attained when the algorithm utilizes 2 randomly chosen predictors at each split. Importantly, this optimized value closely aligns with the default calculated parameter for this specific dataset, confirming the fundamental robustness of the initial model structure, even before explicit tuning.

Step 5: Applying the Final Model for Prediction

Once the random forest model has been thoroughly fitted, diagnosed, and optimized (if necessary), the final step involves confidently utilizing it to predict the response variable for new, unseen observations. This crucial application validates the model’s practical utility in a real-world context.

We first define a new data frame that represents a hypothetical day characterized by specific atmospheric conditions that the model has not yet encountered:

#define new observation
new <- data.frame(Solar.R=150, Wind=8, Temp=70, Month=5, Day=5)

#use fitted bagged model to predict Ozone value of new observation
predict(model, newdata=new)

27.19442

Based on the input values provided for solar radiation, wind speed, temperature, month, and day, the final fitted random forest model predicts that the Ozone value for this specific environmental scenario will be approximately 27.19442. This result clearly demonstrates how a properly constructed and tuned random forest serves as an accurate and reliable predictive tool for complex regression tasks in data science.

The complete, executable R code used throughout this comprehensive example can be accessed and reviewed here.