# Create a dataframe
data <- data.frame(
Value = c(300000, 350000, 280000, 500000, 450000, 400000, 380000),
Style = c("Ranch", "Ranch", "Ranch", "Colonial", "Colonial", "CapeCod","Colonial"),
SqFt = c(1600, 1400, 1700, 1800, 2100, 1500, 1900),
Bedrooms = c(3, 2, 3, 4, 3, 2, 3))# Build the Tree
tree_model <- rpart(Value ~ Style + SqFt + Bedrooms, data = data, method = "anova",
control = rpart.control(cp = 0.02, minsplit = 2))# Plot the tree
rpart.plot(tree_model, type = 3, extra = 0, box.palette = "Blues",
under = TRUE, fallen.leaves = TRUE)| Argument | Meaning | Effect |
|---|---|---|
| tree_model | The fitted model created with rpart(). | Input tree to plot. |
| type = 3 | Controls layout of split and node labels. | Places split labels above branches and node labels inside boxes for readability. |
| extra = 0 | Controls what extra info is displayed in leaf nodes. | Shows only predicted value (for regression) or class (for classification). |
| box.palette = ‘Blues’ | Applies a color gradient from the RColorBrewer palette ‘Blues’. | Darker color = higher predicted value. |
| under = TRUE | Places split condition text under the node box. | Makes the diagram more compact. |
| fallen.leaves = TRUE | Arranges all terminal nodes at the bottom of the plot. | Makes the tree look like a traditional tree with leaves at the bottom. |
# prediction for a new house:
predict(tree_model, newdata = data.frame(Style = "CapeCod", SqFt = 1600, Bedrooms = 3))## 1
## 390000
##
## Regression tree:
## rpart(formula = Value ~ Style + SqFt + Bedrooms, data = data,
## method = "anova", control = rpart.control(cp = 0.02, minsplit = 2))
##
## Variables actually used in tree construction:
## [1] Bedrooms SqFt Style
##
## Root node error: 3.7e+10/7 = 5285714286
##
## n= 7
##
## CP nsplit rel error xerror xstd
## 1 0.695270 0 1.000000 1.36111 0.49042
## 2 0.164189 1 0.304730 1.04362 0.26396
## 3 0.064865 2 0.140541 0.93716 0.22770
## 4 0.020000 4 0.010811 1.25946 0.39740
The Iris dataset is a classic and widely used dataset in statistics and machine learning. It contains 150 observations of iris flowers, collected from three different species: Iris setosa, Iris versicolor, and Iris virginica — with 50 samples per species.
Each observation includes four numerical features describing the physical dimensions of the flowers:
Sepal length (cm)
Sepal width (cm)
Petal length (cm)
Petal width (cm)
The goal of analyses using this dataset is often to classify the species of iris based on these measurements or to explore patterns and relationships among the features. The dataset is well-balanced, easy to visualize, and frequently used for demonstrating data analysis, pattern recognition, and machine learning techniques such as classification and clustering.
# Build a decision tree to predict Species
tree_model <- rpart(Species ~ ., data = iris, method = "class")
# Visualize the tree
rpart.plot(tree_model)# Predict for a new observation
new_data <- data.frame(Sepal.Length=5.1, Sepal.Width=3.5, Petal.Length=1.4, Petal.Width=0.2)
predict(tree_model, new_data, type="class")## 1
## setosa
## Levels: setosa versicolor virginica
Recall: The Classification Error Rate measures the proportion of observations that are incorrectly classified in a node.
\[E=1-\max_k(p_k)\]
where \(p_k\) is the proportion of observations belonging to class \(k\) in the node.
A smaller classification error rate indicates a purer node. A value of zero means all observations belong to the same class.
# Calculate classification error rate at each node.
node_results <- data.frame(
Node = rownames(tree_model$frame),
Observations = tree_model$frame$n,
Misclassified = tree_model$frame$dev,
Classification_Error = tree_model$frame$dev / tree_model$frame$n
)
knitr::kable(node_results, digits = 4)| Node | Observations | Misclassified | Classification_Error |
|---|---|---|---|
| 1 | 150 | 100 | 0.6667 |
| 2 | 50 | 0 | 0.0000 |
| 3 | 100 | 50 | 0.5000 |
| 6 | 54 | 5 | 0.0926 |
| 7 | 46 | 1 | 0.0217 |
library(rpart)
library(rpart.plot)
data(iris)
# Classification tree using Gini index
tree_gini <- rpart(
Species ~ .,
data = iris,
method = "class",
parms = list(split = "gini"),
control = rpart.control(cp = 0.01)
)
# Classification tree using Cross-Entropy
tree_entropy <- rpart(
Species ~ .,
data = iris,
method = "class",
parms = list(split = "information"),
control = rpart.control(cp = 0.01)
)
# Plot the Gini tree
rpart.plot(tree_gini,
type = 2,
extra = 104,
fallen.leaves = TRUE,
main = "Classification Tree: Gini Index")# Plot the Cross-Entropy tree
rpart.plot(tree_entropy,
type = 2,
extra = 104,
fallen.leaves = TRUE,
main = "Classification Tree: Cross-Entropy")##
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
## [1] 0.6666667
## [1] 1.098612
Consider a simple split using Petal.Length.
# Now consider a simple split using Petal.Length.
# Divide observations into two groups
left_node <- subset(iris, Petal.Length < 2.5)
right_node <- subset(iris, Petal.Length >= 2.5)
# Examine species counts
table(left_node$Species)##
## setosa versicolor virginica
## 50 0 0
##
## setosa versicolor virginica
## 0 50 50
The first group contains all 50 Setosa flowers. The second group contains the remaining 100 flowers.
# Functions for Gini and cross-entropy
gini_index <- function(y) {
p <- prop.table(table(y))
1 - sum(p^2)
}
cross_entropy <- function(y) {
p <- prop.table(table(y))
-sum(p[p > 0] * log(p[p > 0]))
}
# Calculate impurity in each node
gini_index(left_node$Species)## [1] 0
## [1] 0.5
## [1] 0
## [1] 0.6931472
This split separates
Setosa perfectly, so the left node has zero impurity. The right node
still contains two species in equal proportions.
Why the tree considers this split useful?
# Weighted Gini after splitting
weighted_gini <-
(50/150) * gini_index(left_node$Species) +
(100/150) * gini_index(right_node$Species)
# Gini improvement
gini_improvement <-
gini_index(iris$Species) - weighted_gini
# Weighted cross-entropy after splitting
weighted_entropy <-
(50/150) * cross_entropy(left_node$Species) +
(100/150) * cross_entropy(right_node$Species)
# Information gain
information_gain <-
cross_entropy(iris$Species) - weighted_entropy
gini_improvement## [1] 0.3333333
## [1] 0.6365142
| Measure | Before Split | After Split | Reduction |
|---|---|---|---|
| Gini Index | 0.6667 | 0.3333 | 0.3333 |
| Cross-Entropy | 1.0986 | 0.4621 | 0.6365 |
A classification tree searches for splits that maximize the reduction in impurity. Both Gini and cross-entropy reward splits that produce purer child nodes, although they can sometimes prefer different splits.
We can use the same iris dataset and the same split (Petal.Length < 2.5) to calculate Residual Mean Deviance for classification trees.
For classification trees, the formula from our lecture is:
\[\text{Residual Mean Deviance} = -\frac{1}{n-|T|}\sum_{m=1}^{|T|}\sum_{k=1}^{K}n_{mk}\ln(p_{mk})\]
# Load iris data
data(iris)
# Use the same split
left_node <- subset(iris, Petal.Length < 2.5)
right_node <- subset(iris, Petal.Length >= 2.5)
# Function to calculate node deviance
node_deviance <- function(y) {
counts <- table(y)
p <- counts / sum(counts) # Only include nonzero probabilities
-sum(counts[p > 0] * log(p[p > 0]))
}
# Deviance before splitting
deviance_before <- node_deviance(iris$Species)
# Deviance after splitting
deviance_after <-
node_deviance(left_node$Species) +
node_deviance(right_node$Species)
# Number of observations
n <- nrow(iris)
# Residual Mean Deviance
rmd_before <- deviance_before / (n - 1)
rmd_after <- deviance_after / (n - 2)
# Create comparison table
results_rmd <- data.frame(
Model = c("Before Split", "After Split"),
Terminal_Nodes = c(1, 2),
Deviance = c(deviance_before, deviance_after),
Residual_Mean_Deviance = c(rmd_before, rmd_after)
)
# Display results
knitr::kable(
results_rmd,
digits = 4,
col.names = c("Model", "Terminal Nodes",
"Deviance", "Residual Mean Deviance")
)| Model | Terminal Nodes | Deviance | Residual Mean Deviance |
|---|---|---|---|
| Before Split | 1 | 164.7918 | 1.1060 |
| After Split | 2 | 69.3147 | 0.4683 |
Before splitting, all 150 observations are in one node, and the Residual Mean Deviance is approximately 1.1060.
After splitting at Petal.Length = 2.5, the observations
are divided into two groups. The Residual Mean Deviance decreases to
approximately 0.4683.
This indicates that the split produces more homogeneous groups and improves the fit according to this measure.
Connection to your previous example: Cross-entropy measures the impurity within a node, while Residual Mean Deviance combines the cross-entropy contributions across terminal nodes and adjusts by \(n-|T|\).
One important distinction: this is a manually calculated entropy-based residual mean deviance. The dev column in rpart classification-tree output uses a different definition (misclassification counts), so its values will not match these calculations.
Residual Mean Deviance can be used to compare different decision trees fitted to the same dataset.
The general idea is: - Smaller Residual Mean Deviance indicates a better fit. - Larger Residual Mean Deviance indicates a poorer fit.
For example, using our iris dataset:
| Tree | Terminal Nodes | Residual Mean Deviance |
|---|---|---|
| Tree 1 (No split) | 1 | 1.1060 |
| Tree 2 (One split) | 2 | 0.4683 |
Here, Tree 2 fits the training data better because it has a smaller Residual Mean Deviance.
One important distinction
Gini index and cross-entropy: Primarily used to decide where to split a classification tree.
Residual Mean Deviance: Used to assess and compare the fit of trees after they have been constructed.
However, a smaller training Residual Mean Deviance does not necessarily mean a better tree for prediction. A complex tree may fit the training data very well but overfit.
For choosing the best tree, we generally prefer cross-validation error or test-set performance, particularly when comparing trees with different numbers of terminal nodes.
Let’s compare two different splits of the iris dataset and see which one produces a smaller Residual Mean Deviance.
We’ll use:
Tree 2: Petal.Length < 2.5 (our previous
split)
Tree 3: Petal.Width < 1.5 (a different
split)
Both trees have two terminal nodes, so this is a straightforward comparison.
data(iris)
# Function to calculate node deviance
node_deviance <- function(y) {
counts <- table(y)
p <- counts / sum(counts)
-sum(counts[p > 0] * log(p[p > 0]))
}
# Function to calculate Residual Mean Deviance for a split
calculate_rmd <- function(data, variable, cutoff) {
left_node <- data[data[[variable]] < cutoff, ]
right_node <- data[data[[variable]] >= cutoff, ]
n <- nrow(data)
T <- 2 # Two terminal nodes
total_deviance <-
node_deviance(left_node$Species) +
node_deviance(right_node$Species)
total_deviance / (n - T)
}
# Tree 2: Split using Petal.Length
rmd_tree2 <- calculate_rmd(iris, "Petal.Length", 2.5)
# Tree 3: Split using Petal.Width
rmd_tree3 <- calculate_rmd(iris, "Petal.Width", 1.5)
# Compare the trees
comparison <- data.frame(
Tree = c("Tree 2", "Tree 3"),
Split = c("Petal.Length < 2.5",
"Petal.Width < 1.5"),
Terminal_Nodes = c(2, 2),
Residual_Mean_Deviance = c(rmd_tree2, rmd_tree3)
)
knitr::kable(comparison, digits = 4)| Tree | Split | Terminal_Nodes | Residual_Mean_Deviance |
|---|---|---|---|
| Tree 2 | Petal.Length < 2.5 | 2 | 0.4683 |
| Tree 3 | Petal.Width < 1.5 | 2 | 0.6614 |
Now, Examine the splits
It’s also helpful to see which species are in each terminal node.
## Species
## Node setosa versicolor virginica
## Left 50 0 0
## Right 0 50 50
## Species
## Node setosa versicolor virginica
## Left 50 35 1
## Right 0 15 49
The split with the smaller Residual Mean Deviance produces purer terminal nodes and fits these training observations better according to the entropy-based measure.
Because both trees have the same number of terminal nodes, the comparison is particularly easy to interpret.
Which split would you choose based on Residual Mean Deviance? Does the composition of the terminal nodes explain your choice?
We can use the ipred package for bagging classification
trees.
# Load required libraries
library(ipred)
library(ggplot2)
# Use the iris dataset
data(iris)
# Fit a bagging model
set.seed(123)
bag_model <- bagging(Species ~ ., data = iris, nbagg = 50)
# Predict and evaluate accuracy
pred_bag <- predict(bag_model, iris)
mean(pred_bag == iris$Species)## [1] 1
Use randomForest directly — it’s the most popular
ensemble method.
library(randomForest)
set.seed(123)
rf_model <- randomForest(Species ~ ., data = iris, ntree = 100, importance = TRUE)
# Model summary
print(rf_model)##
## Call:
## randomForest(formula = Species ~ ., data = iris, ntree = 100, importance = TRUE)
## Type of random forest: classification
## Number of trees: 100
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 4.67%
## Confusion matrix:
## setosa versicolor virginica class.error
## setosa 50 0 0 0.00
## versicolor 0 46 4 0.08
## virginica 0 3 47 0.06
# Fancy
importance_df <- data.frame(Variable = rownames(importance(rf_model)),
Importance = importance(rf_model)[,5])
# Plot variable importance
ggplot(importance_df, aes(x = reorder(Variable, Importance), y = Importance)) +
geom_bar(stat = "identity", fill = "skyblue") +
coord_flip() +
labs(title = "Variable Importance ",
x = "Variable", y = "MeanDecreaseGini")For classification problems (like Species in the iris dataset):
MeanDecreaseAccuracy – how much model accuracy drops when the variable’s values are randomly permuted.
MeanDecreaseGini – how much each variable reduces node impurity (Gini index) across the forest.
For regression problems:
%IncMSE – percent increase in mean squared error when the variable is permuted.
IncNodePurity – total decrease in residual sum of squares (RSS) from splits involving that variable.
We’ll use gbm (gradient boosting machine) for a simple
demonstration.
library(gbm)
set.seed(123)
# gbm requires numeric target, so we convert
iris_boost <- iris
iris_boost$Species <- as.numeric(iris$Species) - 1
boost_model <- gbm(Species ~ ., data = iris_boost,
distribution = "multinomial",
n.trees = 200,
interaction.depth = 3,
shrinkage = 0.1,
cv.folds = 5)
# Summary gives variable importance
summary(boost_model)# Plot training loss
plot(boost_model$train.error, type = "l",
main = "Boosting Training Error",
xlab = "Number of Trees", ylab = "Error")# =========================================
# Ensemble Methods Demo with Visualization
# =========================================
# Load required libraries
library(ipred)
library(randomForest)
library(gbm)
library(ggplot2)
library(dplyr)
# Prepare data
data(iris)
iris2 <- iris[, c("Petal.Length", "Petal.Width", "Species")]
# Create a grid for plotting decision boundaries
make_grid <- function(data, n = 200) {
x1 <- seq(min(data$Petal.Length), max(data$Petal.Length), length.out = n)
x2 <- seq(min(data$Petal.Width), max(data$Petal.Width), length.out = n)
grid <- expand.grid(Petal.Length = x1, Petal.Width = x2)
return(grid)
}
grid <- make_grid(iris2)
# ==============================
# 1️⃣ Bagging
# ==============================
set.seed(123)
bag_model <- bagging(Species ~ ., data = iris2, nbagg = 50)
grid$pred_bag <- predict(bag_model, grid)
# ==============================
# 2️⃣ Random Forest
# ==============================
set.seed(123)
rf_model <- randomForest(Species ~ ., data = iris2, ntree = 100, importance = TRUE)
grid$pred_rf <- predict(rf_model, grid)
# Variable importance plot
varImpPlot(rf_model, main = "Variable Importance - Random Forest")# ==============================
# 3️⃣ Boosting
# ==============================
iris_boost <- iris2
iris_boost$Species <- as.numeric(iris_boost$Species) - 1 # convert to numeric
boost_model <- gbm(Species ~ ., data = iris_boost,
distribution = "multinomial",
n.trees = 200,
interaction.depth = 3,
shrinkage = 0.1,
cv.folds = 5,
verbose = FALSE)
# Predict on grid (boosting)
boost_pred <- predict(boost_model, newdata = grid, n.trees = 200, type = "response")
grid$pred_boost <- factor(apply(boost_pred, 1, which.max) - 1,
labels = levels(iris$Species))
# ==============================
# ✨ Visualization
# ==============================
plot_decision <- function(grid, model_col, title) {
ggplot(grid, aes(Petal.Length, Petal.Width)) +
geom_tile(aes(fill = !!sym(model_col)), alpha = 0.3) +
geom_point(data = iris2, aes(color = Species), size = 2) +
scale_fill_manual(values = c("#B3E5FC", "#C8E6C9", "#FFCDD2")) +
scale_color_manual(values = c("#0288D1", "#388E3C", "#D32F2F")) +
labs(title = title, fill = "Predicted Class", color = "Actual Species") +
theme_minimal()
}
p1 <- plot_decision(grid, "pred_bag", "Bagging Decision Boundary")
p2 <- plot_decision(grid, "pred_rf", "Random Forest Decision Boundary")
p3 <- plot_decision(grid, "pred_boost", "Boosting Decision Boundary")
# Display plots
p1| Method | Key Idea | Visual Interpretation |
|---|---|---|
| Bagging | Averages predictions of many trees built on bootstrapped samples. | Boundaries smoother than a single tree but may still have some noise. |
| Random Forest | Adds random feature selection to decorrelate trees. | Boundaries even smoother and more generalized. |
| Boosting | Sequentially fits trees to correct previous errors. | Produces sharp, precise boundaries (can overfit with too many trees). |
Summary Points
Bagging reduces variance by averaging multiple trees.
Random Forest improves bagging by introducing randomness in feature selection.
Boosting builds trees sequentially, each correcting errors of the previous ones.