Non Linear Regression
Candace Grant
2025-11-10
Simulate data
library(mlbench)
library(caret)
set.seed(200)
trainingData <- mlbench.friedman1(200, sd = 1)
trainingData$x <- data.frame(trainingData$x)
featurePlot(trainingData$x, trainingData$y)
testData <- mlbench.friedman1(5000, sd = 1)
testData$x <- data.frame(testData$x)Tuning other models on the data
knnModel <- train(x = trainingData$x,
y = trainingData$y,
method = "knn",
preProc = c("center", "scale"),
tuneLength = 10)
knnModel## k-Nearest Neighbors
##
## 200 samples
## 10 predictor
##
## Pre-processing: centered (10), scaled (10)
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 200, 200, 200, 200, 200, 200, ...
## Resampling results across tuning parameters:
##
## k RMSE Rsquared MAE
## 5 3.466085 0.5121775 2.816838
## 7 3.349428 0.5452823 2.727410
## 9 3.264276 0.5785990 2.660026
## 11 3.214216 0.6024244 2.603767
## 13 3.196510 0.6176570 2.591935
## 15 3.184173 0.6305506 2.577482
## 17 3.183130 0.6425367 2.567787
## 19 3.198752 0.6483184 2.592683
## 21 3.188993 0.6611428 2.588787
## 23 3.200458 0.6638353 2.604529
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was k = 17.knnPred <- predict(knnModel, newdata = testData$x)
postResample(pred = knnPred, obs = testData$y)## RMSE Rsquared MAE
## 3.2040595 0.6819919 2.5683461Train Models
library(caret)
library(earth)
library(tidyverse)
# Set up cross-validation
ctrl <- trainControl(method = "cv", number = 10)
set.seed(200)
# Train models
models <- list(
Linear = train(x = trainingData$x, y = trainingData$y, method = "lm", trControl = ctrl),
KNN = knnModel,
MARS = train(x = trainingData$x, y = trainingData$y, method = "earth",
tuneGrid = expand.grid(degree = 1:2, nprune = seq(5, 15, by = 2)),
trControl = ctrl),
RF = train(x = trainingData$x, y = trainingData$y, method = "rf",
tuneLength = 3, trControl = ctrl),
SVM = train(x = trainingData$x, y = trainingData$y, method = "svmRadial",
preProc = c("center", "scale"), tuneLength = 5, trControl = ctrl)
)Question 7.2a: Which Model Performs Best?
# Test set performance
test_results <- sapply(models, function(m) {
pred <- predict(m, testData$x)
c(RMSE = RMSE(pred, testData$y),
Rsquared = cor(pred, testData$y)^2)
}) %>% t() %>% as.data.frame() %>%
rownames_to_column("Model") %>%
arrange(RMSE)
#TEST SET PERFORMANCE (Best to Worst)
print(test_results, row.names = FALSE)## Model RMSE Rsquared
## MARS 1.158995 0.9460418
## SVM 2.059224 0.8287291
## RF 2.409470 0.7875506
## Linear 2.697068 0.7084666
## KNN 3.204059 0.6819919# Visualization
test_results %>%
ggplot(aes(x = reorder(Model, -RMSE), y = RMSE, fill = Model)) +
geom_col(show.legend = FALSE) +
geom_text(aes(label = round(RMSE, 2)), vjust = -0.5) +
labs(title = "Model Performance Comparison",
subtitle = "Lower RMSE = Better Performance",
x = "Model", y = "Test RMSE") +
theme_minimal()
Question 7.2a Solution
cat("Answer: The", test_results$Model[1], "model achieves the best performance with RMSE =",
round(test_results$RMSE[1], 3), "and R² =", round(test_results$Rsquared[1], 3), "\n")## Answer: The MARS model achieves the best performance with RMSE = 1.159 and R² = 0.946Question 7.2b: Does MARS Select the Informative Predictors?
# Variable importance
mars_imp <- varImp(models$MARS)$importance %>%
rownames_to_column("Variable") %>%
mutate(
Type = ifelse(Variable %in% paste0("X", 1:5), "Informative (X1-X5)", "Noise (X6-X10)"),
Selected = Overall > 0
)
# Summary
informative <- mars_imp %>% filter(grepl("X[1-5]", Variable))
noise <- mars_imp %>% filter(grepl("X[6-9]|X10", Variable))
#MARS FEATURE SELECTION RESULTS
cat("Informative predictors selected:", sum(informative$Selected), "out of 5\n")## Informative predictors selected: 4 out of 5cat("Noise predictors excluded:", sum(!noise$Selected), "out of 5\n")## Noise predictors excluded: 0 out of 5cat("Selection accuracy:",
round(100 * (sum(informative$Selected) + sum(!noise$Selected)) / 10, 1), "%\n\n")## Selection accuracy: 40 %#Selected variables:
print(mars_imp %>% filter(Selected) %>% select(Variable, Overall, Type))## Variable Overall Type
## 1 X1 100.00000 Informative (X1-X5)
## 2 X4 75.23592 Informative (X1-X5)
## 3 X2 48.72974 Informative (X1-X5)
## 4 X5 15.51884 Informative (X1-X5)# Visualization
ggplot(mars_imp, aes(x = factor(Variable, levels = paste0("X", 1:10)),
y = Overall, fill = Type)) +
geom_col() +
scale_fill_manual(values = c("Informative (X1-X5)" = "cornflowerblue",
"Noise (X6-X10)" = "coral")) +
labs(title = "MARS Variable Importance: Feature Selection",
subtitle = "True function: y = 10·sin(π·X1·X2) + 20·(X3-0.5)² + 10·X4 + 5·X5 + noise",
x = "Predictor", y = "Importance", fill = NULL) +
theme_minimal() +
theme(legend.position = "top")
Question 7.2b Solution
cat("Answer: MARS correctly identified", sum(informative$Selected),
"out of 5 informative predictors and excluded", sum(!noise$Selected),
"out of 5 noise variables (",
round(100 * (sum(informative$Selected) + sum(!noise$Selected)) / 10, 1),
"% accuracy).\n", sep = "")## Answer: MARS correctly identified4out of 5 informative predictors and excluded0out of 5 noise variables (40% accuracy).Question 7.5 Data Preparation
library(caret)
library(tidyverse)
library(mlbench)
library(earth)
library(kernlab)
library(randomForest)
# Generate data (like you did for 7.2)
set.seed(200)
trainingData <- mlbench.friedman1(200, sd = 1)
trainingData$x <- data.frame(trainingData$x)
testData <- mlbench.friedman1(5000, sd = 1)
testData$x <- data.frame(testData$x)
# Setup (same as 6.3)
ctrl <- trainControl(method = "cv", number = 10)
pp <- c("center", "scale") # No missing data, so no medianImpute neededQuestion 7.5 Train Linear and Non Linear Models
cat("Training models...\n")## Training models...set.seed(200)
# Train all models
models <- list(
Linear = train(x = trainingData$x, y = trainingData$y,
method = "lm", trControl = ctrl),
GLMNET = train(x = trainingData$x, y = trainingData$y,
method = "glmnet", trControl = ctrl, preProc = pp),
MARS = train(x = trainingData$x, y = trainingData$y,
method = "earth", trControl = ctrl,
tuneGrid = expand.grid(degree = 1:2, nprune = seq(5, 15, by = 2))),
SVM = train(x = trainingData$x, y = trainingData$y,
method = "svmRadial", trControl = ctrl,
preProc = pp, tuneLength = 5),
RF = train(x = trainingData$x, y = trainingData$y,
method = "rf", trControl = ctrl,
tuneLength = 3, importance = TRUE),
KNN = train(x = trainingData$x, y = trainingData$y,
method = "knn", trControl = ctrl,
preProc = pp, tuneLength = 10)
)
cat(" All models trained\n")## All models trainedQuestion 7.5a Which Nonlinear Model Performs Best?
# Test set performance
test_results <- sapply(models, function(m) {
pred <- predict(m, testData$x)
c(RMSE = RMSE(pred, testData$y),
Rsquared = cor(pred, testData$y)^2,
MAE = MAE(pred, testData$y))
}) %>% t() %>% as.data.frame() %>%
rownames_to_column("Model") %>%
mutate(Type = ifelse(Model %in% c("Linear", "GLMNET"), "Linear", "Nonlinear")) %>%
arrange(RMSE)
#TEST SET PERFORMANCE (Best to Worst)
print(test_results, row.names = FALSE)## Model RMSE Rsquared MAE Type
## MARS 1.158995 0.9460418 0.925023 Nonlinear
## SVM 2.066014 0.8277434 1.561776 Nonlinear
## RF 2.377852 0.7941973 1.877786 Nonlinear
## GLMNET 2.673271 0.7120164 2.048565 Linear
## Linear 2.697068 0.7084666 2.060054 Linear
## KNN 3.122264 0.6690472 2.496365 Nonlinear# Visualization
test_results %>%
ggplot(aes(x = reorder(Model, RMSE), y = RMSE, fill = Type)) +
geom_col() +
coord_flip() +
geom_text(aes(label = round(RMSE, 2)), hjust = -0.2) +
scale_fill_manual(values = c("Linear" = "steelblue", "Nonlinear" = "coral")) +
labs(title = "Model Performance Comparison - Exercise 7.5",
subtitle = "Friedman1 Dataset: Lower RMSE = Better Performance",
x = "Model", y = "Test RMSE", fill = "Model Type") +
theme_minimal() +
theme(legend.position = "top")
# Identify best models
best_overall <- test_results$Model[1]
best_nonlinear <- test_results %>% filter(Type == "Nonlinear") %>% slice(1)
best_linear <- test_results %>% filter(Type == "Linear") %>% slice(1)Question 7.5b Variable Importance Comparison
# Get best nonlinear and best linear models
best_nl_model <- models[[best_nonlinear$Model]]
best_lin_model <- models[[best_linear$Model]]
# Variable importance for best nonlinear
imp_nonlinear <- varImp(best_nl_model, scale = TRUE)$importance %>%
rownames_to_column("Variable") %>%
arrange(desc(Overall)) %>%
slice(1:10) %>%
mutate(
Type = ifelse(Variable %in% paste0("X", 1:5),
"Informative (X1-X5)",
"Noise (X6-X10)"),
Selected = Overall > 0,
Rank = row_number()
)
# Variable importance for linear
imp_linear <- varImp(best_lin_model, scale = TRUE)$importance %>%
rownames_to_column("Variable") %>%
arrange(desc(Overall)) %>%
slice(1:10) %>%
mutate(
Type = ifelse(Variable %in% paste0("X", 1:5),
"Informative (X1-X5)",
"Noise (X6-X10)"),
Rank = row_number()
)
# Count informative vs noise
informative_nl <- imp_nonlinear %>% filter(grepl("X[1-5]", Variable))
noise_nl <- imp_nonlinear %>% filter(grepl("X[6-9]|X10", Variable))
cat("Top 10 Predictors in", best_nonlinear$Model, "(NONLINEAR):\n")## Top 10 Predictors in MARS (NONLINEAR):print(imp_nonlinear %>% select(Rank, Variable, Overall, Type), row.names = FALSE)## Rank Variable Overall Type
## 1 X1 100.00000 Informative (X1-X5)
## 2 X4 75.23592 Informative (X1-X5)
## 3 X2 48.72974 Informative (X1-X5)
## 4 X5 15.51884 Informative (X1-X5)
## 5 X3 0.00000 Informative (X1-X5)cat("\nInformative predictors (X1-X5):", nrow(informative_nl), "out of 10\n")##
## Informative predictors (X1-X5): 5 out of 10cat("Noise predictors (X6-X10):", nrow(noise_nl), "out of 10\n")## Noise predictors (X6-X10): 0 out of 10cat("Selection accuracy:",
round(100 * (nrow(informative_nl) + (5 - nrow(noise_nl))) / 10, 1), "%\n\n")## Selection accuracy: 100 %# Comparison with linear model
cat("Top 10 Predictors in", best_linear$Model, "(LINEAR):\n")## Top 10 Predictors in GLMNET (LINEAR):print(imp_linear %>% select(Rank, Variable, Overall, Type), row.names = FALSE)## Rank Variable Overall Type
## 1 X4 1.000000e+02 Informative (X1-X5)
## 2 X1 9.016243e+01 Informative (X1-X5)
## 3 X2 7.766616e+01 Informative (X1-X5)
## 4 X5 5.901595e+01 Informative (X1-X5)
## 5 X6 1.538641e+01 Noise (X6-X10)
## 6 X7 4.960110e+00 Noise (X6-X10)
## 7 X9 1.658735e+00 Noise (X6-X10)
## 8 X10 1.068399e-01 Noise (X6-X10)
## 9 X3 3.070241e-03 Informative (X1-X5)
## 10 X8 0.000000e+00 Noise (X6-X10)informative_lin <- imp_linear %>% filter(grepl("X[1-5]", Variable))
noise_lin <- imp_linear %>% filter(grepl("X[6-9]|X10", Variable))
cat("\nInformative predictors:", nrow(informative_lin), "out of 10\n")##
## Informative predictors: 6 out of 10cat("Noise predictors:", nrow(noise_lin), "out of 10\n\n")## Noise predictors: 5 out of 10# Compare overlap
overlap <- intersect(imp_nonlinear$Variable, imp_linear$Variable)
unique_nl <- setdiff(imp_nonlinear$Variable, imp_linear$Variable)
unique_lin <- setdiff(imp_linear$Variable, imp_nonlinear$Variable)
#COMPARISON: Nonlinear vs Linear Top 10
cat("Overlap:", length(overlap), "out of 10\n")## Overlap: 5 out of 10if (length(overlap) > 0) {
cat(" ", paste(overlap, collapse = ", "), "\n")
}## X1, X4, X2, X5, X3cat("\nUnique to", best_nonlinear$Model, ":", length(unique_nl), "\n")##
## Unique to MARS : 0if (length(unique_nl) > 0) {
cat(" ", paste(unique_nl, collapse = ", "), "\n")
}
cat("\nUnique to", best_linear$Model, ":", length(unique_lin), "\n")##
## Unique to GLMNET : 5if (length(unique_lin) > 0) {
cat(" ", paste(unique_lin, collapse = ", "), "\n")
}## X6, X7, X9, X10, X8# Visualization
ggplot(imp_nonlinear, aes(x = reorder(Variable, Overall),
y = Overall, fill = Type)) +
geom_col() +
coord_flip() +
scale_fill_manual(values = c("Informative (X1-X5)" = "darkgreen",
"Noise (X6-X10)" = "coral")) +
labs(title = paste("Variable Importance:", best_nonlinear$Model),
subtitle = paste("Informative:", nrow(informative_nl), "| Noise:", nrow(noise_nl)),
x = "Predictor", y = "Importance", fill = NULL) +
theme_minimal() +
theme(legend.position = "top")
Question 7.5c Unique Predictors
if (length(unique_nl) == 0) {
cat("No unique predictors - both models identified the same top 10.\n")
} else {
cat("Exploring", length(unique_nl), "predictor(s) unique to", best_nonlinear$Model, ":\n\n")
for (var in unique_nl) {
var_type <- imp_nonlinear$Type[imp_nonlinear$Variable == var]
var_rank <- imp_nonlinear$Rank[imp_nonlinear$Variable == var]
cat("---", var, "(", var_type, ") - Rank:", var_rank, "---\n")
# Create comparison plots
plot_data <- data.frame(
x = trainingData$x[[var]],
y = trainingData$y
)
# Linear fit
p1 <- ggplot(plot_data, aes(x = x, y = y)) +
geom_point(alpha = 0.5, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "red") +
labs(title = paste(var, "- Linear Fit"),
subtitle = var_type,
x = var, y = "Response (y)") +
theme_minimal()
# Nonlinear fit (LOESS)
p2 <- ggplot(plot_data, aes(x = x, y = y)) +
geom_point(alpha = 0.5, color = "coral") +
geom_smooth(method = "loess", se = TRUE, color = "darkgreen") +
labs(title = paste(var, "- Nonlinear Fit (LOESS)"),
subtitle = var_type,
x = var, y = "Response (y)") +
theme_minimal()
print(gridExtra::grid.arrange(p1, p2, ncol = 2))
}
}## No unique predictors - both models identified the same top 10.Interpretation of Results
The absence of unique predictors in MARS’s top 10 list is a positive finding that demonstrates superior feature selection. MARS achieved perfect accuracy (100%) by identifying all five truly informative predictors (X1-X5) while excluding all noise variables (X6-X10). In contrast, the linear model’s top 10 included five noise variables alongside the five informative ones. The 5-out-of-10 overlap represents agreement on the core informative variables, with the difference being that MARS appropriately discarded irrelevant predictors that the linear model retained.
This result exemplifies MARS’s strength: automatic feature selection through adaptive basis functions that focus modeling effort on genuine signal. The 56.6% RMSE improvement (1.159 vs 2.670) confirms that MARS not only selected the correct predictors but also effectively captured their complex nonlinear relationships from the true Friedman function. The finding that MARS identified exactly five important variables—matching the known structure of the data-generating process—validates its capability to distinguish signal from noise in high-dimensional settings.
Summary
# Calculate improvement here
if (!exists("improvement")) {
improvement <- ((best_linear$RMSE - best_nonlinear$RMSE) / best_linear$RMSE) * 100
}
cat("(a) Best nonlinear model:", best_nonlinear$Model, "\n")## (a) Best nonlinear model: MARScat(" Test RMSE =", round(best_nonlinear$RMSE, 4), "\n")## Test RMSE = 1.159cat(" Test R² =", round(best_nonlinear$Rsquared, 4), "\n")## Test R² = 0.946cat(" Improvement over linear:", round(improvement, 1), "%\n\n")## Improvement over linear: 56.6 %cat("(b) Variable dominance: INFORMATIVE (X1-X5)\n")## (b) Variable dominance: INFORMATIVE (X1-X5)cat(" Informative in top 10:", nrow(informative_nl), "out of 5 true informatives\n")## Informative in top 10: 5 out of 5 true informativescat(" Noise in top 10:", nrow(noise_nl), "out of 5 noise variables\n")## Noise in top 10: 0 out of 5 noise variablescat(" Selection accuracy:",
round(100 * (nrow(informative_nl) + (5 - nrow(noise_nl))) / 10, 1), "%\n")## Selection accuracy: 100 %cat(" Overlap with linear:", length(overlap), "out of 10\n\n")## Overlap with linear: 5 out of 10cat("(c) Unique predictors:", length(unique_nl), "\n")## (c) Unique predictors: 0if (length(unique_nl) > 0) {
cat(" List:", paste(unique_nl, collapse = ", "), "\n")
unique_informative <- sum(grepl("X[1-5]", unique_nl))
cat(" Informative:", unique_informative, "| Noise:", length(unique_nl) - unique_informative, "\n")
}