db <- read_diabetes()
hf <- read_heart_failure()
cv <- setup_KFold(n_resamples = 10L, seed = 2026L)23 Additive Models and Basis Expansions in rtemis
rtemis fits generalized additive models with setup_GAM(), using mgcv::gam(); multivariate adaptive regression splines with setup_MARS(), using the earth package; the highly adaptive lasso with setup_HAL(), using hal9001; and isotonic regression with setup_Isotonic(). This chapter fits them to the diabetes and heart failure data and compares their cross-validated performance with that of the linear models of Section 15.5.
23.1 Generalized additive models
setup_GAM(k = 5L) fits a penalized spline with a basis of dimension 5 for each numeric feature and a linear term for each factor:
mod_gam <- train(
db,
hyperparameters = setup_GAM(k = 5L),
execution_config = setup_SerialExecution(seed = 2026L),
question = diabetes_question,
verbosity = 0L
)
summary(mod_gam@model)
Family: gaussian
Link function: identity
Formula:
progression ~ s(age, k = 5) + s(bmi, k = 5) + s(bp, k = 5) +
s(s1, k = 5) + s(s2, k = 5) + s(s3, k = 5) + s(s4, k = 5) +
s(s5, k = 5) + s(s6, k = 5) + sex
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 163.48 3.68 44.43 < 2e-16 ***
sex2 -24.23 5.75 -4.21 3.1e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(age) 2.68 3.21 2.38 0.064 .
s(bmi) 1.76 2.21 29.22 < 2e-16 ***
s(bp) 1.00 1.00 26.54 8.8e-07 ***
s(s1) 1.00 1.00 0.86 0.354
s(s2) 1.00 1.00 0.57 0.449
s(s3) 1.00 1.00 0.27 0.601
s(s4) 1.00 1.00 1.03 0.311
s(s5) 2.39 2.83 12.35 2.8e-05 ***
s(s6) 2.10 2.61 2.67 0.050 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.533 Deviance explained = 54.9%
GCV = 2880.4 Scale est. = 2776.5 n = 442The edf column gives the effective degrees of freedom of each smooth. A value of 1 means that the penalty reduced the smooth to a straight line: in this model, the smooths of bp, s1, s2, s3, s4 are linear, and those of age, bmi, s5, s6 are curved. The partial effect of each feature, its contribution \(f_j(x_j)\) to the prediction, is returned by predict(type = "terms") on the mgcv model:
terms_gam <- predict(mod_gam@model, type = "terms")
partial <- function(v) {
o <- order(db[[v]])
list(x = db[[v]][o], y = terms_gam[o, paste0("s(", v, ")")])
}
draw_panels(
lapply(c("age", "bmi", "s5", "s6"), function(v) {
pe <- partial(v)
draw_line(
pe[["x"]], list(`Partial effect` = pe[["y"]]),
points = FALSE, xlab = v, ylab = "Contribution to progression"
)
}),
ncol = 2L,
height = 640
)Each curve shows how the predicted progression changes with one feature, holding the others fixed, on the scale of the outcome. Cross-validation assesses the model:
mod_gam_cv <- train(
db,
hyperparameters = setup_GAM(k = 5L),
outer_resampling_config = cv,
execution_config = setup_SerialExecution(seed = 2026L),
question = diabetes_question,
verbosity = 0L
)23.2 Multivariate adaptive regression splines
setup_MARS() fits an additive MARS model by default (degree = 1L), and allows interactions between pairs of hinge functions with degree = 2L:
mod_mars <- train(
db,
hyperparameters = setup_MARS(),
execution_config = setup_SerialExecution(seed = 2026L),
verbosity = 0L
)
summary(mod_mars@model)Call: earth(x=matrix[442,11], y=c(151,75,141,20...), pmethod="backward",
trace=0, degree=1, nfold=0, ncross=1, stratify=TRUE, thresh=0.001,
minspan=0, endspan=0, newvar.penalty=0, fast.k=20, fast.beta=1)
coefficients
(Intercept) 64.248
sex_2 -24.666
h(38-age) 1.927
h(age-38) 0.662
h(bmi-22.2) 5.925
h(99-bp) -0.753
h(bp-99) 1.550
h(249-s1) 1.211
h(181.2-s2) -0.998
h(s5-4.0254) 77.420
h(s6-94) 1.365
Selected 11 of 20 terms, and 8 of 11 predictors
Termination condition: Reached nk 23
Importance: bmi, s5, bp, s1, s2, sex_2, s6, age, sex_1-unused, s3-unused, ...
Number of terms at each degree of interaction: 1 10 (additive model)
GCV 2951 RSS 1183375 GRSq 0.5046 RSq 0.5485Each term h(x - t) or h(t - x) is a hinge function with knot t. The summary reports how many of the terms generated by the forward pass were kept after pruning, and the generalized cross-validation estimate of \(R^2\) (GRSq). get_varimp() returns earth’s importance measures:
plot_varimp(mod_mars)mod_mars_cv <- train(
db,
hyperparameters = setup_MARS(),
outer_resampling_config = cv,
execution_config = setup_SerialExecution(seed = 2026L),
question = diabetes_question,
verbosity = 0L
)
mod_mars2_cv <- train(
db,
hyperparameters = setup_MARS(degree = 2L),
outer_resampling_config = cv,
execution_config = setup_SerialExecution(seed = 2026L),
question = diabetes_question,
verbosity = 0L
)23.3 The highly adaptive lasso
setup_HAL(max_degree = 1L) restricts the basis to functions of single features, an additive model; smoothness_orders = 1L, the default, uses hinge functions, and 0L would use indicator functions. The penalty is selected by internal cross-validation, whose fold assignment is fixed by seed:
mod_hal_cv <- train(
db,
hyperparameters = setup_HAL(max_degree = 1L, seed = 2026L),
outer_resampling_config = cv,
execution_config = setup_SerialExecution(seed = 2026L),
question = diabetes_question,
verbosity = 0L
)HAL is the slowest of the methods in this chapter: its basis contains one function per observed value of each feature, and its fit is a lasso over all of them, cross-validated within each training fold.
23.4 Comparison
mod_glm_cv <- train(
db,
hyperparameters = setup_GLM(),
outer_resampling_config = cv,
execution_config = setup_SerialExecution(seed = 2026L),
verbosity = 0L
)
compare_db <- rbindlist(lapply(
list(
GLM = mod_glm_cv,
GAM = mod_gam_cv,
`MARS, degree 1` = mod_mars_cv,
`MARS, degree 2` = mod_mars2_cv,
`HAL, degree 1` = mod_hal_cv
),
function(m) {
data.table(
rsq = m@metrics_test@mean_metrics[["rsq"]],
rsq_sd = m@metrics_test@sd_metrics[["rsq"]],
rmse = m@metrics_test@mean_metrics[["rmse"]]
)
}
), idcol = "model")
compare_db model rsq rsq_sd rmse
<char> <num> <num> <num>
1: GLM 0.4935 0.08864 54.43
2: GAM 0.5112 0.07796 53.51
3: MARS, degree 1 0.4721 0.08536 55.58
4: MARS, degree 2 0.4654 0.07196 56.01
5: HAL, degree 1 0.4660 0.08097 55.97The GAM has the highest mean test \(R^2\), 0.511 against 0.493 for the linear model, a difference much smaller than the standard deviation across folds. Since all models were assessed on the same folds, their per-fold results can be compared directly:
Fold_1 Fold_2 Fold_3 Fold_4 Fold_5 Fold_6 Fold_7 Fold_8
0.047113 -0.003861 -0.018685 0.032624 0.041075 0.011495 -0.001762 0.008602
Fold_9 Fold_10
0.042019 0.018976 The GAM’s \(R^2\) is higher than the linear model’s in 7 of the 10 folds. The standard deviation of the paired differences, 0.022, is smaller than that of either model’s fold results, 0.078 for the GAM, because both models face the same test cases in each fold. The differences are still not independent across folds, for the reasons of Section 10.9. The interaction terms of the degree-2 MARS model do not improve on the additive model, consistent with Section 19.3.
23.5 Classification
MARS fits binary outcomes with a logistic link. For the heart failure data:
mod_mars_hf <- train(
hf,
hyperparameters = setup_MARS(),
outer_resampling_config = cv,
execution_config = setup_SerialExecution(seed = 2026L),
question = hf_question
)<Resampled Classification Model>
MARS (Multivariate Adaptive Regression Splines (earth))
⟳ Tested using 10 independent folds.
<Resampled Classification Training Metrics>
Aggregate Confusion Matrix across resamples.
Predicted
Reference Yes No
Yes 146 340
No 68 2056
Showing mean (sd) across resamples.
Sensitivity: 0.300 (0.071)
Specificity: 0.968 (0.009)
Balanced Accuracy: 0.634 (0.035)
PPV: 0.682 (0.073)
NPV: 0.858 (0.012)
F1: 0.414 (0.075)
Accuracy: 0.844 (0.014)
AUC: 0.807 (0.021)
Brier Score: 0.115 (0.007)
<Resampled Classification Test Metrics>
Aggregate Confusion Matrix across resamples.
Predicted
Reference Yes No
Yes 11 43
No 11 225
Showing mean (sd) across resamples.
Sensitivity: 0.200 (0.132)
Specificity: 0.953 (0.047)
Balanced Accuracy: 0.577 (0.070)
PPV: NA (NA)
NPV: 0.840 (0.024)
F1: NA (NA)
Accuracy: 0.813 (0.045)
AUC: 0.714 (0.146)
Brier Score: 0.144 (0.032)
The mean test AUC, 0.714, is similar to that of the logistic regression with the same folds.
23.6 Isotonic regression
setup_Isotonic() fits a non-decreasing function of a single feature. Its use for recalibrating probabilities is shown in Chapter 43. As a regression of progression on s5 alone:
mod_iso <- train(
db[, .(s5, progression)],
hyperparameters = setup_Isotonic(),
execution_config = setup_SerialExecution(seed = 2026L),
verbosity = 0L
)
s5_grid <- data.table(s5 = sort(unique(db[["s5"]])))
draw_fits(
db[["s5"]], db[["progression"]],
x_fit = s5_grid[["s5"]],
fits = list(Isotonic = as.numeric(predict(mod_iso, s5_grid))),
xlab = "s5", ylab = "Disease progression",
point_alpha = 0.4
)