Classification: From Logistic Regression to Neural Networks

In classification, every training observation has an established category. Binary has two possible classes; multiclass has three or more mutually exclusive classes; multilabel allows more than one true class per row. “Not diagnosed” is not a valid negative label merely because a code is absent.

1. Binary classification

Our synthetic diabetes target is binary. Each method takes the same predictor table, but learns a different rule.

No model family always wins. Model quality depends on labels, missingness, sample size, deployment shift, class prevalence, and clinical consequences, not only algorithm sophistication.

Example: load the data and count the two classes of the binary target. Every method in the table has to learn from these same counts.

library(tidymodels)
ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE) |>
  mutate(diabetes = factor(diabetes, levels = c("Yes", "No")))
count(ehr, diabetes)

# A tibble: 2 x 2
  diabetes     n
  <fct>    <int>
1 Yes         54
2 No         546

round(mean(ehr$diabetes == "Yes"), 3)

[1] 0.09

There are 54 “Yes” and 546 “No” rows, so 9% of the rows are positive. We list “Yes” first so that tidymodels treats it as the event of interest. That 9% is the rate any useful model has to beat.

2. A fair, runnable comparison on fake EHR data

Instead of fitting algorithms to different splits and declaring a winner, compare logistic regression, a decision tree, a random forest, and kNN with the same training folds. The included script uses training-fold preprocessing and saves a comparison table.

source("scripts/05_supervised_comparison.R")

Look in outputs/supervised_classification_cv.csv. No p-values appear there; compare ROC AUC and PR AUC (and the std_err column) and remember that only about 54 of the 600 synthetic records are positive, so the scores are noisy. The code does not inspect the locked test set to choose an algorithm.

Example: the same idea on a small scale. Two methods, logistic regression and a decision tree, are scored on the same five folds, with the same preprocessing learned inside each fold.

rec <- recipe(diabetes ~ age + bmi + sbp + pulse, data = ehr) |>
  step_impute_median(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())
set.seed(500)
folds <- vfold_cv(ehr, v = 5, strata = diabetes)
cv_auc <- function(spec) {
  workflow() |> add_recipe(rec) |> add_model(spec) |>
    fit_resamples(folds, metrics = metric_set(roc_auc)) |>
    collect_metrics()
}
bind_rows(logistic = cv_auc(logistic_reg()),
          tree = cv_auc(decision_tree(mode = "classification")),
          .id = "model") |>
  select(model, .metric, mean, std_err, n)

# A tibble: 2 x 5
  model    .metric  mean std_err     n
  <chr>    <chr>   <dbl>   <dbl> <int>
1 logistic roc_auc 0.828  0.0289     5
2 tree     roc_auc 0.701  0.0204     5

Logistic regression has a mean ROC AUC of about 0.83 and the tree about 0.70, each over 5 folds. The std_err column (0.03 and 0.02) shows how much the score moves from fold to fold, so only a gap clearly larger than that is worth reading as a difference.

3. What is a decision boundary?

A decision boundary separates regions assigned to categories. Logistic regression creates a straight line in a simple two-feature plane (on log odds); a decision tree creates axis-aligned steps; an RBF SVM or neural network can produce curved boundaries. More curved does not mean more correct. Increasing complexity can fit noise.

Example: a logistic regression on age and bmi gives a straight boundary. Suppose we flag a patient when the predicted probability of “Yes” reaches 0.2. The boundary is the set of (age, BMI) pairs where the log odds equal qlogis(0.2).

b <- coef(glm(diabetes == "Yes" ~ age + bmi, data = ehr, family = binomial))
round(b, 3)

(Intercept)         age         bmi 
     -9.675       0.068       0.103 

ages <- c(30, 50, 70)
bmi_at_boundary <- (qlogis(0.2) - b["(Intercept)"] - b["age"] * ages) / b["bmi"]
round(setNames(bmi_at_boundary, paste("age", ages)), 1)

age 30 age 50 age 70 
  60.5   47.4   34.3 

At age 30 a BMI above about 60.5 is needed to reach the flag, at age 50 about 47.4, and at age 70 about 34.3. The boundary slopes down as age rises, and because the model is linear it is one straight line in the age-BMI plane.

4. How does naive Bayes work?

Bayes’ rule relates a prior probability, feature evidence, and a posterior probability. The naive assumption is that predictors are conditionally independent given the class. BMI and SBP are often related, so that assumption can be violated. Even so, naive Bayes can be a useful benchmark. In R it is available in the optional e1071 package (e1071::naiveBayes()), which this book does not demonstrate; never treat its scores as automatically calibrated disease probabilities.

Example: the naive step is multiplication. Suppose 10% of patients have the condition, and a high BMI is seen in 30% of patients with it and 10% of those without, while a high SBP is seen in 40% and 20%. Treating the two as independent given the class, the evidence for each class is the prior times the two rates.

prior   <- c(Yes = 0.10, No = 0.90)
lik_bmi <- c(Yes = 0.30, No = 0.10)
lik_sbp <- c(Yes = 0.40, No = 0.20)
evidence <- prior * lik_bmi * lik_sbp
round(evidence / sum(evidence), 3)

Yes  No 
0.4 0.6 

Dividing by the total turns the two scores into posterior probabilities: 0.4 for Yes and 0.6 for No. The 10% prior rose to 40% after seeing both findings. If BMI and SBP move together, counting them as two separate pieces of evidence overstates the case.

5. Multiclass classification

Suppose a fictional clinic records referral categories A, B, or C, with exactly one true category per example. Methods include:

• Multinomial logistic regression (one coefficient set for each category except a reference category).
• One-vs-rest / one-vs-one strategies with binary classifiers.
• Multiclass trees, random forests, gradient boosting, and SVM.
• LDA/QDA, kNN, naive Bayes, and neural-network softmax.

The built-in iris dataset (flower measurements, not healthcare) provides a safe, established three-class illustration. Run:

source("scripts/07_multiclass_demo.R")

It performs fold-wise assessment for LDA, QDA, kNN, tree, and multinomial regression. A high iris score is not an EHR performance claim.

Example: a multinomial logistic regression on iris gives one probability per class for each flower, and the predicted class is the largest one.

m <- nnet::multinom(Species ~ ., data = iris, trace = FALSE)
round(head(predict(m, type = "probs"), 3), 3)

  setosa versicolor virginica
1      1          0         0
2      1          0         0
3      1          0         0

table(predicted = predict(m), truth = iris$Species)

            truth
predicted    setosa versicolor virginica
  setosa         50          0         0
  versicolor      0         49         1
  virginica       0          1        49

The first three flowers get probability 1 for setosa, and each row of probabilities sums to 1 because the classes are mutually exclusive. The table has 148 of 150 flowers on the diagonal; one versicolor and one virginica are swapped. This is accuracy on the training rows, so it is optimistic.

6. Ordinal and multilabel outcomes

Ordinal categories have an order (e.g. level 1 < 2 < 3). Ordinary multiclass models ignore this order, while cumulative link / proportional odds models and ordinal-specific algorithms exploit it under assumptions. R packages include ordinal and MASS::polr().

Multilabel means one row can simultaneously have “A” and “C”. That requires one binary indicator per label, classifier chains, multivariate approaches, or specialized neural outputs. Multiclass softmax, which forces one of A/B/C, is wrong for multilabel.

In real EHR data, diagnoses can co-occur, may be recorded at multiple encounters, and can be missing for administrative reasons. Define the observation and label process before modeling.

Example: an ordinal outcome and a multilabel outcome look different in the data. First, systolic blood pressure cut into three ordered levels, modeled with MASS::polr().

ord <- ehr |>
  filter(!is.na(sbp)) |>
  mutate(sbp_level = cut(sbp, c(0, 120, 140, Inf), right = FALSE,
                         labels = c("normal", "raised", "high"),
                         ordered_result = TRUE))
pol <- MASS::polr(sbp_level ~ age, data = ord)
round(coef(pol), 3)

  age 
0.103 

round(predict(pol, data.frame(age = c(30, 70)), type = "probs"), 2)

  normal raised high
1   0.82   0.18 0.00
2   0.07   0.74 0.19

The age coefficient is positive (0.103), so an older age moves the probability toward the higher levels: a 30-year-old is most likely “normal” (0.82) and a 70-year-old most likely “raised” (0.74). One slope serves both cut points, which is the proportional-odds assumption. A multilabel outcome instead needs one yes/no column per label, and a row can have several 1s.

label_table <- tibble(row = 1:3, A = c(1, 0, 1), B = c(0, 1, 1), C = c(1, 0, 0))
label_table

# A tibble: 3 x 4
    row     A     B     C
  <int> <dbl> <dbl> <dbl>
1     1     1     0     1
2     2     0     1     0
3     3     1     1     0

rowSums(label_table[, c("A", "B", "C")])

[1] 2 1 2

Rows 1 and 3 have two true labels each, so a one-of-three softmax output could never be fully right for them.

7. Imbalance and cost-sensitive decisions

When positives are rare, accuracy can be inflated by predicting the negative class for everyone. Use:

• A prevalence / majority baseline and class-support counts.
• PR AUC, sensitivity, specificity, precision, calibration and decision-relevant costs.
• Class weighting or resampling inside training folds only if supported by the method and justified.
• A threshold chosen from out-of-fold training predictions (not in-sample fits) based on a stated objective.
• An independent test with confirmed truth labels.

SMOTE creates synthetic feature-space points. It must never be applied before the split or to a held-out test set, and it is not a replacement for verifying outcomes.

Example: a model that predicts “No” for every patient is the majority baseline. Score it on the 600 rows.

all_no <- tibble(truth = ehr$diabetes,
                 pred  = factor(rep("No", nrow(ehr)), levels = c("Yes", "No")))
c(accuracy = accuracy(all_no, truth, pred)$.estimate,
  sensitivity = sens(all_no, truth, pred)$.estimate)

   accuracy sensitivity 
       0.91        0.00 

Accuracy is 0.91, which sounds good, but sensitivity is 0: this model finds none of the 54 positives. That is why the list above asks for the baseline, class counts, and metrics other than accuracy.

8. Probabilities and calibration

A classifier’s ranking score is not necessarily a probability. SVM margin, class-weighted or oversampled model, and ensembles may require calibration via Platt scaling, isotonic regression, or another properly validated method. Calibration itself is a learned transformation and belongs in training/validation, not fitted on the locked test.

Example: a calibration check groups the predicted probabilities into bins and compares each bin’s average prediction with the share of positives actually seen there.

ehr_bmi <- filter(ehr, !is.na(bmi))
p <- predict(glm(diabetes == "Yes" ~ age + bmi, data = ehr_bmi, family = binomial),
             type = "response")
tibble(p = p, y = ehr_bmi$diabetes == "Yes") |>
  mutate(bin = cut(p, c(0, 0.05, 0.1, 0.2, 1))) |>
  group_by(bin) |>
  summarise(n = n(), mean_pred = mean(p), actual = mean(y))

# A tibble: 4 x 4
  bin            n mean_pred actual
  <fct>      <int>     <dbl>  <dbl>
1 (0,0.05]     306    0.0186 0.0163
2 (0.05,0.1]    97    0.0722 0.0825
3 (0.1,0.2]     84    0.139  0.155 
4 (0.2,1]       78    0.316  0.295 

In each bin the average prediction is close to the actual rate (for example 0.316 predicted against 0.295 observed in the top bin), so this logistic model is reasonably calibrated on these rows. These are the same rows the model was fitted on, so a real check needs validation data.

9. Which method would you start with?

For modest tabular EHR-like data, start with majority baseline → logistic regression → a tree/forest. Add kNN, penalized models, boosting, SVM, or neural networks only when there is a scientific reason and enough credible validation data. Complexity is a cost as well as a capability.

Example: start with the baseline, then logistic regression, on a split with the test rows left alone until the end.

set.seed(500)
sp <- initial_split(ehr, prop = 0.75, strata = diabetes)
train <- training(sp)
test <- testing(sp)
fit_lr <- workflow() |> add_recipe(rec) |> add_model(logistic_reg()) |> fit(train)
res <- predict(fit_lr, test) |> bind_cols(test["diabetes"])
c(baseline_accuracy = mean(test$diabetes == "No"),
  logistic_accuracy = accuracy(res, diabetes, .pred_class)$.estimate,
  logistic_sensitivity = sens(res, diabetes, .pred_class)$.estimate)

   baseline_accuracy    logistic_accuracy logistic_sensitivity 
          0.92000000           0.92000000           0.08333333 

The baseline and the logistic model have the same accuracy (0.92), so accuracy does not show any gain. The logistic sensitivity of 0.083 shows it finds one in twelve positives, a small gain that the accuracy number hid. Only if this gain matters would a forest or boosting be worth the added cost.

fit1 <- glm(diabetes == "Yes" ~ age + bmi, data = ehr, family = binomial)
prob <- head(predict(fit1, type = "response"), 4)
round(prob, 3)

    1     2     3     4 
0.021 0.004 0.055 0.009 

ifelse(prob >= 0.2, "Yes", "No")

   1    2    3    4 
"No" "No" "No" "No" 

toy <- tibble(
  truth = factor(c("Yes","Yes","No","No","No","Yes"), levels = c("Yes","No")),
  pred  = factor(c("Yes","No","No","No","Yes","No"), levels = c("Yes","No")))
c(precision = precision(toy, truth, pred)$.estimate,
  sensitivity = sens(toy, truth, pred)$.estimate)

  precision sensitivity 
  0.5000000   0.3333333 

boosted <- bind_rows(ehr_bmi, filter(ehr_bmi, diabetes == "Yes") |> slice(rep(row_number(), 4)))
fit_os <- glm(diabetes == "Yes" ~ age + bmi, data = boosted, family = binomial)
c(true_rate = mean(ehr_bmi$diabetes == "Yes"),
  mean_prediction = mean(predict(fit_os, ehr_bmi, type = "response")))

      true_rate mean_prediction 
     0.08672566      0.25402787 

set.seed(500)
iris_folds <- vfold_cv(iris, v = 5, strata = Species)
fold_acc <- purrr::map_dfr(iris_folds$splits, function(s) {
  tr <- analysis(s); va <- assessment(s)
  tibble(LDA = mean(predict(MASS::lda(Species ~ ., tr), va)$class == va$Species),
         QDA = mean(predict(MASS::qda(Species ~ ., tr), va)$class == va$Species))
})
fold_acc

# A tibble: 5 x 2
    LDA   QDA
  <dbl> <dbl>
1 1     1    
2 0.967 0.967
3 0.967 0.967
4 0.967 0.967
5 1     1    

colMeans(fold_acc)

 LDA  QDA 
0.98 0.98