Day 5: Tree-based models

0. Introduction

Microplastics (MP) can impact various soil functions, and understanding these effects is crucial. In this tutorial, we will analyze a dataset of 150 soil samples from Germany, each tested for the impact of microplastics on soil functions. The dataset contains several soil property features (e.g. bulk density, nutrient levels, etc.) and two outcome variables: ES_WSA and ES_Dec, which are effect sizes of microplastics on two soil functions (WSA indicates water-stable aggregates stability, and Dec decomposition rate). Our goal is to use tree-based machine learning models to explore how soil properties influence the direction and magnitude of microplastic effects on a soil function. We will focus on one response variable ES_WSA in the tutorial and leave the other ES_Dec for you to practice.

What we will cover

Decision Tree modeling: Building a regression tree to predict the microplastic effect size (ES_WSA) from soil properties, visualizing the tree, and interpreting the results.
Model evaluation: Using k-fold cross-validation to assess the decision tree’s performance and discussing the bias–variance trade-off (controlling overfitting via pruning).
Random Forest modeling: Building an ensemble of trees (random forest) for the same task, interpreting feature importance, and comparing its performance to the single decision tree.

Throughout the tutorial, we will explain key concepts and include commented R code for clarity. By the end, you should understand how to implement and interpret tree-based models for this soil microplastics dataset.

Structure of the `MP_Effects_150.csv` dataset

Import the dataset and run the str() command should show that we have 150 observations and 8 variables. The key columns include:

bulk_density: Bulk density of the soil (a measure of compaction, g/cm³).
Total_N: Total nitrogen content of the soil (ppm).
CN_ratio: Carbon-to-Nitrogen ratio of the soil.
Olsen_P: Olsen phosphorus (available P) level in soil
WPT: Water penetration time, evaluates soil water repellency.
LUI_15_20: Land Use Intensity index (from 2025 to 2020).
ES_WSA: Effect size of microplastics on the soil’s water-stable aggregates (WSA) – our target variable for this tutorial.
ES_Dec: Effect size of microplastics on decomposition – the other target.

1. Decision Tree Model

Building the Tree

We’ll use the rpart package to build a regression tree for ES_WSA using the soil properties as predictors. Let’s load the necessary library and train the model:

# Install rpart if not already installed
install.packages("rpart")    # run this line if rpart is not installed
install.packages("rpart.plot")  # for plotting the tree

library(rpart)
library(rpart.plot)

# Fit a regression tree model for ES_WSA using all soil property predictors
tree_model <- rpart(
  ES_WSA ~ bulk_density + Total_N + CN_ratio + Olsen_P + WPT + LUI_15_20,
  data = data,
  method = "anova"  # "anova" for regression trees (splits by minimizing squared error)
)

# Print a summary of the tree model
tree_model

## n= 150 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 150 0.118355800  0.0031723350  
##    2) LUI_15_20>=2.79 7 0.030130430 -0.0238406300 *
##    3) LUI_15_20< 2.79 143 0.082867420  0.0044946480  
##      6) WPT>=15.16667 117 0.048824900  0.0009078162  
##       12) bulk_density>=1.04 7 0.012321640 -0.0307762900 *
##       13) bulk_density< 1.04 110 0.029028890  0.0029240780  
##         26) Total_N< 20.975 102 0.022244350  0.0014496690  
##           52) Olsen_P< 35.365 95 0.017580610  0.0002286634 *
##           53) Olsen_P>=35.365 7 0.002599971  0.0180204600 *
##         27) Total_N>=20.975 8 0.003735677  0.0217227900 *
##      7) WPT< 15.16667 26 0.025763660  0.0206353900  
##       14) bulk_density>=0.84 19 0.015203410  0.0132458900 *
##       15) bulk_density< 0.84 7 0.006706717  0.0406926100 *

When you print the tree_model, by default it will show the splits the tree made. Each line of the output describes a node: which variable and cutoff it splits on, how many observations go left/right, and the predicted value in that node. It might look a bit cryptic textually, so it’s often easier to visualize the tree structure:

# Plot the decision tree
rpart.plot(tree_model)

This will draw the tree diagram for us. The root node (top) contains all data and shows the initial mean of ES_WSA. The tree then splits on the most informative variable. For our data, you might find that the first split is on LUI_15_20 (land use intensity). For example, the tree might discover a rule like:

If LUI_15_20 <= 2.0 (low-intensity land use sites), then follow the left branch.
If LUI_15_20 > 2.0 (higher-intensity use), follow the right branch.

This kind of split suggests that land use intensity has a significant effect on how microplastics impact WSA. Lower-intensity sites might respond differently (perhaps showing slight positive effects or less negative impact) whereas very high intensity might be associated with negative effects. In our fitted tree, indeed one of the terminal nodes corresponded to soils with very high LUI_15_20 and high nitrogen, which had a negative predicted effect size (indicating microplastics reduce WSA in highly managed soils). You can interpret each branch similarly. For instance, a branch might further split on bulk_density or Total_N, meaning those properties further distinguish the effect. Each leaf (terminal node) gives a predicted ES_WSA. Because this is a regression, that predicted value is simply the mean ES_WSA of training samples in that leaf.

Avoiding Overfitting: Pruning the Tree

One issue with decision trees is that they can grow very complex and overfit the training data (fitting noise rather than true signal). A fully grown tree might have many splits that only apply to a few data points, which doesn’t generalize well. To control this, we prune the tree by limiting how far it grows. In rpart, this is handled via the complexity parameter (cp). The cp essentially sets a threshold on the improvement needed for a split: a smaller cp allows the tree to keep splitting on smaller improvements (making a bigger tree), while a larger cp makes the tree stop earlier (simpler tree). Using a very low cp will lead to an overfit, overly complex tree

By default, rpart already uses a certain cp (often 0.01) and also performs 10-fold cross-validation internally to help pick an optimal tree size. We can examine the cross-validation results for different tree sizes using the printcp() function:

# Print the cross-validation results (complexity parameter table)
printcp(tree_model)

## 
## Regression tree:
## rpart(formula = ES_WSA ~ bulk_density + Total_N + CN_ratio + 
##     Olsen_P + WPT + LUI_15_20, data = data, method = "anova")
## 
## Variables actually used in tree construction:
## [1] bulk_density LUI_15_20    Olsen_P      Total_N      WPT         
## 
## Root node error: 0.11836/150 = 0.00078904
## 
## n= 150 
## 
##         CP nsplit rel error xerror    xstd
## 1 0.057609      0   1.00000 1.0085 0.30931
## 2 0.032559      3   0.82163 1.1326 0.32667
## 3 0.025760      4   0.78907 1.1346 0.32374
## 4 0.017437      5   0.76331 1.1272 0.32443
## 5 0.010000      6   0.74587 1.1326 0.32508

This will output a table with columns: CP, number of splits, relative error, xerror, and xstd. Here, xerror is the cross-validation error (as a fraction of the total sum of squares, since it’s relative to the root node error). The row with the smallest xerror indicates the optimal complexity. We can retrieve that and prune the tree:

# Find the optimal cp that minimizes cross-validated error
best_cp <- tree_model$cptable[which.min(tree_model$cptable[,"xerror"]), "CP"]

# Prune the tree using the optimal cp
pruned_tree <- prune(tree_model, cp = best_cp)

# Plot the pruned tree
rpart.plot(pruned_tree)

The pruned tree is simpler (some of the lower-importance splits removed) and likely generalizes better. In many cases, the pruning might remove splits that were based on only a few points (potentially noise).

Cross-Validation Performance: We can explicitly measure how well our tree is likely to perform on new data by using k-fold cross-validation. Let’s use 10-fold CV (k=10) to estimate the model error (e.g., using the caret package for convenience):

# Install caret if needed
install.packages("caret")

library(caret)

set.seed(123)  # for reproducibility

# Set up 10-fold cross-validation
train_ctrl <- trainControl(method = "cv", number = 10)

# Train a decision tree with cross-validation (caret will tune the cp parameter automatically)
tree_cv <- train(ES_WSA ~ bulk_density + Total_N + CN_ratio + Olsen_P + WPT + LUI_15_20,
                 data = data,
                 method = "rpart",
                 trControl = train_ctrl,
                 tuneLength = 10)
                 
# View cross-validation results
tree_cv$results

##             cp       RMSE   Rsquared        MAE     RMSESD RsquaredSD
## 1  0.000000000 0.03073906 0.06407568 0.02142924 0.01076780 0.08423074
## 2  0.006401039 0.03051363 0.06761290 0.02087703 0.01077592 0.10313306
## 3  0.012802077 0.03033920 0.06480009 0.02011338 0.01087363 0.11127202
## 4  0.019203116 0.02994271 0.07836617 0.01964024 0.01051672 0.13369534
## 5  0.025604155 0.02973896 0.08793589 0.01957754 0.01070569 0.15756587
## 6  0.032005193 0.02952987 0.09190881 0.01931088 0.01084119 0.15592387
## 7  0.038406232 0.02919186 0.09531429 0.01911680 0.01123274 0.15533176
## 8  0.044807271 0.02882921 0.09104206 0.01885248 0.01129399 0.15623109
## 9  0.051208309 0.02861034 0.10236141 0.01861955 0.01142156 0.16161117
## 10 0.057609348 0.02845764 0.11121938 0.01835801 0.01153557 0.16991056
##          MAESD
## 1  0.004834204
## 2  0.004975121
## 3  0.004969050
## 4  0.004605698
## 5  0.004582870
## 6  0.004919286
## 7  0.005196206
## 8  0.004863799
## 9  0.005240635
## 10 0.005421448

The output will show the model performance (e.g. RMSE(Root Mean Square Error)) for various values of cp tried. It will also indicate the best cp chosen. You should see that the cross-validated RMSE of the tree is relatively low (on the order of a few hundredths for ES_WSA). For example, in our case the 10-fold CV RMSE was around 0.03 (effect size units). This is comparable to the standard deviation of ES_WSA in the data, meaning the tree does have some predictive power but not very high (which is not surprising given the subtle effect sizes).

The takeaway is that cross-validation gives us an unbiased estimate of how the tree will do on unseen data. It helps ensure we aren’t too optimistic from just looking at training error. In this case, if we had grown a very complex tree, the cross-val error would have been high, signaling overfitting. By using the optimal cp (or pruning), we balance the bias–variance trade-off: not too simple (high bias) and not too complex (high variance). Cross-validation helps find the sweet spot that minimizes prediction error on unseen data

Interpretation: From the decision tree model, we learned which soil properties are most influential. In our example, Land Use Intensity (LUI_15_20) appeared at the top of the tree, suggesting it’s a key factor in how microplastics affect soil structure (WSA). Other properties like Total_N (soil nitrogen) and bulk_density also appeared in splits, indicating they contribute to differences in MP effects. For instance, one branch of the tree suggested that at sites with very high land use intensity and high nitrogen, microplastics had a negative effect on WSA (possibly due to already stressed soil conditions), whereas in low-intensity sites the effect was closer to neutral or slightly positive. Decision trees make it easy to explain such interactions as simple rules (e.g., “if land use is low and bulk density is medium, then slight positive effect on WSA”). However, a single tree may not capture all complex relationships, and small changes in data could alter the tree. To address that instability and potentially improve accuracy, we turn to an ensemble method: the random forest.

2. Random Forest Model (optional)

A random forest is essentially an ensemble of many decision trees. The idea is to build a large number of trees on random subsets of the data and predictors, and then aggregate their results. Each tree contributes a vote (for classification) or an estimate (for regression), and the forest’s prediction is the average of all those tree predictions. By averaging many trees, a random forest reduces variance compared to a single tree, making the predictions more stable and often more accurate. In fact, random forests are designed to combat overfitting by using the bagging (bootstrap aggregating) approach and random feature selection for splits – this de-correlates the trees so that their errors cancel out to some extent. The result is a robust model that usually has better generalization performance than a single decision tree, at the cost of interpretability (we no longer have one simple tree to visualize, but we can still interpret the model in other ways).

Unlike a single decision tree, a random forest cannot be neatly drawn as one tree – it’s the combination of hundreds of trees. Instead, we look at summary metrics, such as overall error rate and variable importance plots that show which features were most influential across all those trees.

Let’s build a random forest model for ES_WSA using the randomForest package:

# Install randomForest if not installed
install.packages("randomForest")

library(randomForest)

set.seed(123)

# Fit a random forest model to predict ES_WSA
rf_model <- randomForest(
  ES_WSA ~ bulk_density + Total_N + CN_ratio + Olsen_P + WPT + LUI_15_20,
  data = data,
  ntree = 500,      # number of trees in the forest
  mtry = 3,         # number of variables sampled at each split (default is sqrt of total predictors)
  importance = TRUE # whether to calculate variable importance
)

# Print the random forest model summary
print(rf_model)

## 
## Call:
##  randomForest(formula = ES_WSA ~ bulk_density + Total_N + CN_ratio +      Olsen_P + WPT + LUI_15_20, data = data, ntree = 500, mtry = 3,      importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##           Mean of squared residuals: 0.000770949
##                     % Var explained: 2.29

When you print the rf_model, you’ll see output like:

Type of random forest: regression.
Number of trees grown (500).
Number of variables tried at each split (mtry).
Out-of-bag (OOB) error estimate: an MSE value.
Possibly an pseudo-R² (called % Var explained).

The OOB error is the error measured on samples not used in each tree’s bootstrap training set (on average, each tree is built on ~2/3 of the data, leaving 1/3 out-of-bag). This OOB error serves a similar role to a cross-validation error – it is an unbiased estimate of model performance without needing a separate test set. For our model, the OOB estimate might show an MSE around, say, 0.0007 (in effect size units²), which corresponds to an RMSE of ~0.026. That is a bit better than the single decision tree’s ~0.03 RMSE, indicating the random forest is indeed capturing more signal. The output also gives a “% Variance Explained” – don’t be alarmed if this is low (the effect size variance is very small to begin with), but any improvement is useful.

Feature Importance in Random Forest

One of the advantages of random forests is that we can examine variable importance. Because the model is a combination of many trees, we cannot just read one set of splits to see importance. Instead, randomForest computes metrics like %IncMSE: the increase in prediction error (MSE) if a given feature’s values are randomly permuted. A larger increase means the model relied more on that feature for accurate predictions.

Let’s extract and plot the importance:

# Get variable importance scores
importance(rf_model)

##                %IncMSE IncNodePurity
## bulk_density 10.704457   0.019829126
## Total_N       9.941886   0.021538126
## CN_ratio      3.217187   0.009469874
## Olsen_P      -1.134471   0.008599999
## WPT           7.965759   0.010922097
## LUI_15_20     6.486016   0.035641988

# Plot variable importance
varImpPlot(rf_model)

The importance scores will show up in a table (and the plot will visualize them).

(This is an illustrative example). The %IncMSE column is the most relevant for us – it tells how much the mean squared error increases when that variable is left out. A higher value means more importance. According to the random forest on our data, Bulk density tends to be the top contributor (it has the highest %IncMSE, meaning if we scramble soil bulk density data, the model error shoots up the most). The next most important is often Total_N (soil nitrogen), followed by WPT (soil water penetration time). This not align with what we saw in the single tree, because single tree model has very high variance for capturing the signals. The other variables (LUI_15_20, C:N ratio, Olsen P) show lower importance scores, suggesting that they don’t greatly influence the MP effect on WSA compared to Bulk density or Total_N.

The varImpPlot will display a bar chart of %IncMSE for each variable, making it easy to see the ranking.

3. Conclusion and Next steps

In this tutorial, we demonstrated how to use tree-based models to analyze the effects of soil properties on microplastic impact (effect size) in soils. The decision tree provided a straightforward model that we could visualize and interpret as rules. We took care to prune the tree and use cross-validation to avoid overfitting, illustrating the importance of model complexity control (a very low complexity parameter would overfit the data, so we chose an optimal value via cross-validation).

We then used a random forest, which by aggregating many trees achieved a more robust model (lower error) and gave us a ranked importance of predictors. The ensemble model is more resistant to overfitting due to its averaging of many de-correlated trees. We discussed how the bias-variance trade-off is handled: a single tree can have low bias but high variance (unless pruned), whereas a random forest reduces variance without increasing bias too much by averaging many low-bias tree

Next steps

Now that you have seen the workflow for ES_WSA, you can repeat a similar analysis for microplastic effects on soil decomposition rate ES_Dec, and try to answer the following questions:

How good is the model in predicting microplastic effects by using soil characteristics?
Which soil characteristic(s) is more important in predicting microplastic effects on soil decomposition rate?
What is cross-validation and Why do we use it?

Day 5: Tree-based models