Supervised Regression: From Lines to Ensembles

This chapter opens Part III. The Part II capstone (chapter 18) ran one supervised workflow on a Yes/No outcome; Part III keeps that same workflow and adds more methods, starting with numeric outcomes. Regression is supervised learning because training rows contain the answer. Regression is not limited to fitting straight lines; tree ensembles, nearest neighbors, support vector methods, and neural networks can also estimate numerical outcomes.

1. The common workflow

For any regression approach: define the numeric outcome → choose inputs available at prediction time → split by patient/time → preprocess within training folds → compare to a training-mean or other baseline → evaluate with out-of-fold MAE/RMSE (and often R²) → select → assess once on locked test.

RMSE: square root of mean squared errors, sensitive to large mistakes. MAE: mean absolute error, easier to interpret in outcome units. R²: fraction of outcome variation explained relative to a constant prediction, which can be negative out of sample. None of them demonstrates causation. The script scripts/05_supervised_comparison.R reports this residual-based R² as rsq_trad; the metric named rsq is a different, squared-correlation measure that cannot be negative.

Example: five invented patients with observed SBP, a training-mean baseline that predicts 130 for everyone, and a model’s predictions. We compute each metric by hand.

obs   <- c(120, 135, 150, 128, 142)
base  <- rep(130, 5)                         # training-mean baseline
model <- c(124, 133, 146, 131, 138)
metrics <- function(pred) c(
  RMSE = sqrt(mean((obs - pred)^2)),
  MAE  = mean(abs(obs - pred)),
  R2   = 1 - sum((obs - pred)^2) / sum((obs - mean(obs))^2))
round(rbind(baseline = metrics(base), model = metrics(model)), 2)

          RMSE MAE    R2
baseline 11.60 9.8 -0.23
model     3.49 3.4  0.89

The model is off by 3.4 mmHg on average (MAE) against 9.8 for the baseline, with RMSE 3.49 against 11.60 and R-squared 0.89 against -0.23. The baseline’s R-squared is negative because its constant, 130, is not the mean of these five values (135): a prediction that is worse than predicting the true mean scores below 0.

2. Linear regression

library(tidymodels)
ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE)
numeric_rows <- ehr |>
  filter(!is.na(sbp)) # outcome must be known for supervised learning
set.seed(500)
split <- initial_split(numeric_rows, prop = .8)
tr <- training(split)

linear_model <- lm(sbp ~ age + bmi, data = tr, na.action = na.omit)
summary(linear_model)

Call:
lm(formula = sbp ~ age + bmi, data = tr, na.action = na.omit)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9798  -7.1824  -0.0525   7.1510  10.5110 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 75.56386    1.83527   41.17   <2e-16 ***
age          0.49739    0.01795   27.70   <2e-16 ***
bmi          0.70406    0.04856   14.50   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.06 on 430 degrees of freedom
  (31 observations deleted due to missingness)
Multiple R-squared:  0.6813,    Adjusted R-squared:  0.6798 
F-statistic: 459.6 on 2 and 430 DF,  p-value: < 2.2e-16

Both age and bmi are positive: each extra year of age goes with about 0.50 mmHg higher SBP and each extra BMI point with about 0.70 mmHg, holding the other fixed. The multiple R-squared of 0.68 says the two variables explain about 68% of the variation in SBP in these training rows, and the residual standard error of 7.06 means predictions are typically off by about 7 mmHg. The note “31 observations deleted due to missingness” is the missing bmi rows that lm() dropped.

This demonstration uses observed training rows for the lm fit. For general evaluation, do not globally delete missing values or impute before splitting; use a recipe within CV. A coefficient is a conditional association, not a causal effect.

3. Regularized regression: ridge, LASSO, elastic net

Ordinary regression can become unstable with many correlated inputs. Regularization penalizes large coefficients.

• Ridge (L2) shrinks many coefficients, usually not exactly to zero.
• LASSO (L1) can set coefficients exactly to zero, acting as feature selection.
• Elastic net combines L1 and L2.
• Tune penalty strength and mixture inside training CV; standardize numeric inputs inside training folds.

Optional example (requires glmnet; not part of the required course installation):

# install.packages("glmnet")
library(glmnet)
# x must be a numeric predictor matrix built using training-only preprocessing.
# y must be the training numeric outcome.
# cv.glmnet(x, y, family = "gaussian", alpha = 1) # LASSO

Example: ridge regression needs no package. For standardized inputs the ridge coefficients have a closed form, solve(t(X) %*% X + lambda * diag(p), t(X) %*% y), and we can watch them shrink as the penalty lambda grows.

cc <- na.omit(tr[, c("sbp", "age", "bmi", "pulse")])
X <- scale(as.matrix(cc[, c("age", "bmi", "pulse")]))
y <- cc$sbp - mean(cc$sbp)
ridge <- function(lambda) drop(solve(t(X) %*% X + lambda * diag(3), t(X) %*% y))
round(sapply(c(none = 0, small = 50, large = 500, huge = 5000), ridge), 2)

      none small large huge
age   9.55  8.45  4.16 0.69
bmi   4.88  4.26  1.97 0.31
pulse 0.07  0.18  0.33 0.09

With lambda = 0 the coefficients are the ordinary least-squares ones (age 9.55, bmi 4.88 mmHg per standard deviation). As the penalty grows, age and bmi shrink toward 0 (0.69 and 0.31 at huge) without reaching exactly 0, which is the ridge behavior. pulse has almost no effect to begin with (0.07).

4. Regression trees and random forests

A regression tree splits the feature space into regions and predicts a number within each region. A random forest averages predictions from many trees.

library(rpart)
tree <- rpart(
  sbp ~ age + bmi + pulse,
  data = tr, method = "anova",
  control = rpart.control(maxdepth = 4, minsplit = 15)
)
head(predict(tree, newdata = tr))

       1        2        3        4        5        6 
140.9200 116.7903 110.7879 105.2308 133.4062 133.4043 

Each number is the tree’s predicted SBP for one of the first six training patients, such as 140.92 for patient 1. Rows that land in the same leaf get the same prediction (patients 5 and 6 are close but not identical here because they fall in different leaves). Only 10 distinct values exist in all, because this depth-limited tree has 10 leaves.

This is a fit demonstration, not a valid test score. Trees can learn interactions and nonlinearities but overfit when allowed to grow too deeply. A random forest provides smoother ensemble predictions, often sacrificing simplicity.

5. k-nearest-neighbor regression

For a new row, choose the k training rows closest in feature space and average their numerical outcomes. Small k can be noisy; large k oversmooths. Distances depend on the units, so learn scaling parameters from training folds.

Use the included comparison script to fit linear regression, tree, random forest, and kNN under matched resampling:

source("scripts/05_supervised_comparison.R")

The script writes training-fold CV comparisons into outputs/ and does not use final test scores for choosing the winner. It is the practical workbook for Chapters 19–21: its Part B compares regression models (this chapter), and its Part A compares binary classifiers (chapters 20 and 21).

Example: kNN regression by hand for one patient. We scale age and BMI, find the 5 training rows closest to patient 1 of the training table, and average their SBP.

knn_data <- na.omit(tr[, c("sbp", "age", "bmi")])
Z <- scale(knn_data[, c("age", "bmi")])
d <- sqrt((Z[, 1] - Z[1, 1])^2 + (Z[, 2] - Z[1, 2])^2)    # distance to row 1
nearest <- order(d)[2:6]                                    # skip row 1 itself
c(actual_sbp = knn_data$sbp[1], knn_prediction = mean(knn_data$sbp[nearest]))

    actual_sbp knn_prediction 
           151            143 

The five nearest training patients to patient 1 have an average SBP of 143, against an actual value of 151. kNN is a local average, so it lands near the SBP of similar patients but not exactly on this one.

6. Generalized linear models (GLMs)

Linear regression assumes an unbounded continuous response with a specific error structure. Other outcomes call for other models:

A tiny numeric example (unrelated to patient care):

set.seed(50)
visits <- data.frame(staff = rep(1:4, each = 8),
                     week = rep(1:8, 4))
visits$count <- rpois(nrow(visits), lambda = 2 + visits$staff)
count_model <- glm(count ~ staff, family = poisson(), data = visits)
coef(count_model)

(Intercept)       staff 
  0.7570478   0.2320399 

Poisson regression models the conditional mean count; it assumes a particular mean–variance relationship. Check dispersion, exposure time, and whether an offset is needed.

The coefficients are on the log scale. The slope 0.23 means each extra staff member multiplies the expected count by exp(0.23), about 1.26. The data were simulated with a mean count of 2 + staff, which goes from 3 at staff level 1 to 6 at level 4, about a 26% rise per level on average, and the fitted slope agrees.

7. Flexible smooth and robust regression

• GAM (generalized additive model): sums smooth functions of predictors, capturing bends while retaining interpretable individual effects. R: mgcv::gam() (often included with R).
• Quantile regression: predicts a conditional quantile, such as the 90th percentile, rather than the mean. R: quantreg.
• Robust regression: reduces sensitivity to extreme observations; example MASS::rlm().
• Mixed-effects/hierarchical models: account for patients nested within sites or repeated visits; examples lme4 and nlme. Still design the test split to match the prediction target.
• Support-vector regression (SVR): margin-based regression, sometimes with nonlinear kernels. R: kernlab or e1071.
• Gradient boosting: iteratively adds weak predictive trees; examples xgboost, lightgbm, catboost with different ecosystem support and installation requirements.
• Neural-network regression: learns flexible nonlinear functions, but needs sufficient data and careful regularization; R: nnet, torch, or keras.

Example: robust regression with MASS::rlm() (MASS is part of R’s standard install list for this course). Ten typed points follow a line, and one data-entry error is added.

xx <- 1:10
yy <- c(2.1, 4.0, 6.2, 7.9, 10.1, 12.0, 13.8, 16.1, 18.0, 20.2)
yy_bad <- replace(yy, 10, 60)                     # a typing error at the last point
round(rbind(lm_clean = coef(lm(yy ~ xx)),
            lm_bad   = coef(lm(yy_bad ~ xx)),
            rlm_bad  = coef(MASS::rlm(yy_bad ~ xx))), 2)

         (Intercept)   xx
lm_clean        0.04 2.00
lm_bad         -7.92 4.17
rlm_bad         0.04 2.00

The straight-line fit on the clean data has slope 2.00. One typing error at the last point drags the least-squares slope to 4.17 and the intercept to -7.92, while rlm() stays at 2.00 and 0.04, the clean-data fit. Robust fitting is useful when a few extreme values should not steer the whole line.

8. When would each be sensible?

If a simple linear model meets the study goal, keep it. Add trees when meaningful nonlinear interactions seem plausible; consider forests or boosting if repeated CV shows material improvement; consider a GAM when smooth nonlinear effects are needed; select count models for event counts, and survival models for censored times rather than forcing them into an ordinary regression.

Example: a small cross-validated comparison of the four models in the script, with the same recipe and folds as scripts/05_supervised_comparison.R. We keep the result in reg_cv and use it again below.

set.seed(500)
folds <- vfold_cv(tr, v = 5)
rec <- recipe(sbp ~ age + bmi + pulse + family_history + sex + clinic_site, data = tr) |>
  step_impute_median(all_numeric_predictors()) |>
  step_impute_mode(all_nominal_predictors()) |>
  step_dummy(all_nominal_predictors()) |>
  step_zv(all_predictors()) |>
  step_normalize(all_numeric_predictors())
specs <- list(
  linear = linear_reg() |> set_engine("lm"),
  tree = decision_tree(tree_depth = 4, min_n = 10) |> set_engine("rpart") |> set_mode("regression"),
  forest = rand_forest(trees = 150, mtry = 3, min_n = 5) |> set_engine("ranger") |> set_mode("regression"),
  knn = nearest_neighbor(neighbors = 15, weight_func = "rectangular") |>
    set_engine("kknn") |> set_mode("regression"))
reg_cv <- purrr::imap_dfr(specs, function(spec, name) {
  set.seed(500)
  fit <- fit_resamples(workflow() |> add_recipe(rec) |> add_model(spec), folds,
                       metrics = metric_set(rmse, mae))
  collect_metrics(fit) |> mutate(model = name)
})
reg_cv |> select(model, .metric, mean, std_err) |> arrange(.metric, mean)

# A tibble: 8 x 4
  model  .metric  mean std_err
  <chr>  <chr>   <dbl>   <dbl>
1 linear mae      6.45  0.112 
2 tree   mae      7.19  0.289 
3 knn    mae      7.29  0.105 
4 forest mae      7.41  0.147 
5 linear rmse     7.26  0.0636
6 forest rmse     8.45  0.142 
7 tree   rmse     8.48  0.260 
8 knn    rmse     8.75  0.107 

Linear regression has the smallest error on both metrics (RMSE 7.26, MAE 6.45), ahead of the forest (8.45), tree (8.48) and kNN (8.75) on RMSE. The standard errors of the RMSE range from about 0.06 to 0.26, so linear’s lead of about 1.2 RMSE points over the next model is several times larger than any of them. With a gap that large, the simple model is the one to keep.

Clinical caution: Choosing a flexible model does not repair a missing or poorly defined target. Variable importance is not causal evidence.

reg_cv |> filter(.metric == "rmse") |> select(model, mean, std_err) |> arrange(mean)

# A tibble: 4 x 3
  model   mean std_err
  <chr>  <dbl>   <dbl>
1 linear  7.26  0.0636
2 forest  8.45  0.142 
3 tree    8.48  0.260 
4 knn     8.75  0.107 

c(sbp_missing = sum(is.na(ehr$sbp)), bmi_missing = sum(is.na(ehr$bmi)),
  bmi_missing_sbp_known = sum(is.na(ehr$bmi) & !is.na(ehr$sbp)))

          sbp_missing           bmi_missing bmi_missing_sbp_known 
                   20                    35                    34 

raw <- knn_data[, c("age", "bmi")]
nn <- function(m) order(sqrt((m[, 1] - m[1, 1])^2 + (m[, 2] - m[1, 2])^2))[2]
c(years = nn(as.matrix(raw)),
  months = nn(cbind(age = raw$age * 12, bmi = raw$bmi)),
  scaled = nn(scale(as.matrix(raw))))

 years months scaled 
   235    126    246