15 Linear Models in rtemis

setup_GLM() fits a generalized linear model with stats::glm(): linear regression for a numeric outcome and logistic regression for a binary factor outcome. This chapter fits both, examines their coefficients, shows the effect of collinearity, and records their cross-validated performance as the reference for the more flexible models of later chapters.

15.1 Linear regression

A linear regression of disease progression on the ten baseline variables of the diabetes data, fitted to all cases:

db <- read_diabetes()
mod_lm <- train(
  db,
  hyperparameters = setup_GLM(),
  execution_config = setup_SerialExecution(seed = 2026L),
  question = diabetes_question,
  verbosity = 0L
)
coefs <- summary(mod_lm@model)[["coefficients"]]
knitr::kable(coefs, digits = 3)

	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	-357.427	67.058	-5.330	0.000
age	-0.036	0.217	-0.168	0.867
sex2	-22.860	5.836	-3.917	0.000
bmi	5.603	0.717	7.813	0.000
bp	1.117	0.225	4.958	0.000
s1	-1.090	0.573	-1.901	0.058
s2	0.746	0.531	1.406	0.160
s3	0.372	0.782	0.475	0.635
s4	6.534	5.959	1.097	0.273
s5	68.483	15.670	4.370	0.000
s6	0.280	0.273	1.025	0.306

Each coefficient is the expected difference in progression per unit difference in the feature, with the others fixed: for example, 5.60 units of progression per unit of BMI and 1.12 per unit of blood pressure. get_varimp() returns the coefficients as the model’s variable importance measure, and plot_varimp() displays them:

plot_varimp(mod_lm)

Raw coefficients are in the units of their features and are not comparable across features: the coefficient of s5, which ranges from 3.2581 to 6.107, is much larger than that of bp, which ranges from 62 to 133. Standardizing the numeric features within training (Section 13.2) expresses every coefficient per standard deviation:

mod_lm_std <- train(
  db,
  preprocessor_config = setup_SupervisedPreprocessor(scale = TRUE),
  hyperparameters = setup_GLM(),
  execution_config = setup_SerialExecution(seed = 2026L),
  verbosity = 0L
)
plot_varimp(mod_lm_std)

15.2 Collinearity

The six serum measurements are correlated, some strongly:

serum <- db[, .(s1, s2, s3, s4, s5, s6)]
draw_heatmap(cor(serum), zlim = c(-1, 1), show_values = TRUE)

s1 and s2 have a correlation of 0.90. Their coefficients in the standardized model have opposite signs. The coefficients of collinear features vary between training sets; their sum, which determines the prediction, varies less. Refitting the standardized model in each fold of a 10-fold cross-validation shows this:

mod_lm_cv <- train(
  db,
  preprocessor_config = setup_SupervisedPreprocessor(scale = TRUE),
  hyperparameters = setup_GLM(),
  outer_resampling_config = setup_KFold(n_resamples = 10L, seed = 2026L),
  execution_config = setup_SerialExecution(seed = 2026L),
  question = diabetes_question
)

<Resampled Regression Model>
GLM (GLM (generalized linear model))
▣ Preprocessed using centering, scaling.
⟳ Tested using 10 independent folds.

  <Resampled Regression Training Metrics>
  Showing mean (sd) across resamples.
   MAE: 43.239 (0.482)
   MSE: 2853.158 (57.301)
  RMSE: 53.413 (0.537)
    R²: 0.519 (0.009)

  <Resampled Regression Test Metrics>
  Showing mean (sd) across resamples.
   MAE: 44.131 (4.284)
   MSE: 2984.432 (539.155)
  RMSE: 54.432 (4.901)
    R²: 0.493 (0.089)

fold_coefs <- rbindlist(lapply(mod_lm_cv@models, function(m) {
  as.list(coef(m@model)[c("bmi", "s1", "s2")])
}))
fold_coefs[, s1_plus_s2 := s1 + s2]
fold_coefs[, lapply(.SD, sd)]

     bmi    s1    s2 s1_plus_s2
   <num> <num> <num>      <num>
1: 1.274 6.327 3.837      2.737

Across folds, the standard deviation of the s1 coefficient is 6.3 and that of s2 is 3.8, against 1.3 for bmi; the standard deviation of their sum is 2.7.

15.3 Logistic regression

For the heart failure data, setup_GLM() fits a logistic regression. The fitted glm object returns the linear predictor \(\eta = x^\top \hat\beta\) with type = "link" and the probability with type = "response":

hf <- read_heart_failure()
mod_logit <- train(
  hf,
  hyperparameters = setup_GLM(),
  execution_config = setup_SerialExecution(seed = 2026L),
  question = hf_question,
  verbosity = 0L
)
eta <- predict(mod_logit@model, type = "link")
prob <- predict(mod_logit@model, type = "response")
draw_scatter(
  eta, prob,
  group = hf[["death"]],
  xlab = "Linear predictor",
  ylab = "Predicted probability of death"
)

Every patient lies on the logistic curve of Figure 14.1. The model ranks patients by their linear predictor, and the coefficients determine how each feature moves a patient along the curve. predict() on the rtemis model returns the same probabilities:

all.equal(
  as.numeric(predict(mod_logit, features(hf))),
  as.numeric(prob)
)

[1] TRUE

15.4 Fixed transformations

Serum creatinine and creatinine phosphokinase are strongly right-skewed (Section 7.2.4). A logarithm is a fixed transformation (Section 12.1) and can be applied before resampling. The cross-validated performance with and without it:

hf_log <- copy(hf)
hf_log[, `:=`(
  serum_creatinine = log(serum_creatinine),
  creatinine_phosphokinase = log(creatinine_phosphokinase)
)]
cv <- setup_KFold(n_resamples = 10L, seed = 2026L)
compare_log <- rbindlist(lapply(list(original = hf, log = hf_log), function(dat) {
  mod <- train(
    dat,
    hyperparameters = setup_GLM(),
    outer_resampling_config = cv,
    execution_config = setup_SerialExecution(seed = 2026L),
    verbosity = 0L
  )
  mod@metrics_test@mean_metrics[c("auc", "brier_score")]
}), idcol = "features")
compare_log

   features    auc brier_score
     <char>  <num>       <num>
1: original 0.7085      0.1395
2:      log 0.7062      0.1382

The transformation changes the cross-validated AUC and Brier score by less than their variation across folds. A transformation that improves a model’s assumptions need not improve its predictions; the cross-validated estimate decides.

15.5 Reference performance

The linear models fitted with 10-fold cross-validation, using the resampling seed of every comparison in this book, provide the reference against which later models are compared:

mod_hf_cv <- train(
  hf,
  hyperparameters = setup_GLM(),
  outer_resampling_config = cv,
  execution_config = setup_SerialExecution(seed = 2026L),
  question = hf_question,
  verbosity = 0L
)
reference <- data.table(
  dataset = c("Heart failure", "Diabetes"),
  metric = c("AUC", "R-squared"),
  mean = c(
    mod_hf_cv@metrics_test@mean_metrics[["auc"]],
    mod_lm_cv@metrics_test@mean_metrics[["rsq"]]
  ),
  sd = c(
    mod_hf_cv@metrics_test@sd_metrics[["auc"]],
    mod_lm_cv@metrics_test@sd_metrics[["rsq"]]
  )
)
reference

         dataset    metric   mean      sd
          <char>    <char>  <num>   <num>
1: Heart failure       AUC 0.7085 0.11702
2:      Diabetes R-squared 0.4935 0.08864