1 Why a Part 2?

In Part 1 we saw the basic complementarity:

  • Machine learning (ML) is powerful for prediction and pattern discovery.
  • Bayesian statistics is powerful for uncertainty, generative thinking, and decision-making under incomplete information.

In Part 2 we go one level deeper. We ask:

When does a model become confidently wrong, and how can statistics and ML rescue each other?

Research themes include astrophysics, theoretical sciences, physics of complex systems, chemical and biological sciences, macromolecular sciences, condensed matter, and materials physics. Therefore, our toy examples below are chosen from domains that feel close to natural-science research: physical chemistry, geoscience, solar/astronomical time series, and biological morphology.

The aim is not to teach all of Bayesian statistics or all of ML. The aim is to give researchers a working instinct:

ML learns patterns; statistics asks whether those patterns are stable, interpretable, calibrated, and scientifically meaningful.

2 Packages and helper functions

The document is designed to run in RStudio and publish to RPubs. It uses light packages only. The Bayesian calculations are written mostly from scratch so that students can see what is happening.

pkgs <- c(
  "ggplot2", "MASS", "dplyr", "tidyr", "tibble",
  "rpart", "rpart.plot", "palmerpenguins", "knitr", "scales"
)

new_pkgs <- pkgs[!(pkgs %in% installed.packages()[, "Package"])]
if (length(new_pkgs) > 0) {
  install.packages(new_pkgs, repos = "https://cloud.r-project.org")
}

invisible(lapply(pkgs, library, character.only = TRUE))
# MASS also has a function named select(); use dplyr::select() explicitly below.

theme_set(theme_minimal(base_size = 12))
set.seed(20260708)

2.1 Small utility functions

rmse <- function(y, yhat) sqrt(mean((y - yhat)^2, na.rm = TRUE))
mae  <- function(y, yhat) mean(abs(y - yhat), na.rm = TRUE)

auc_simple <- function(y01, score) {
  y01 <- as.integer(y01)
  if (length(unique(y01)) < 2) return(NA_real_)
  r <- rank(score, ties.method = "average")
  n1 <- sum(y01 == 1)
  n0 <- sum(y01 == 0)
  (sum(r[y01 == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}

brier <- function(y01, p) mean((as.integer(y01) - p)^2)

classification_table <- function(y01, p, cutoff = 0.5) {
  pred <- ifelse(p >= cutoff, 1, 0)
  y01 <- as.integer(y01)
  tibble(
    accuracy    = mean(pred == y01),
    sensitivity = ifelse(sum(y01 == 1) > 0, mean(pred[y01 == 1] == 1), NA_real_),
    specificity = ifelse(sum(y01 == 0) > 0, mean(pred[y01 == 0] == 0), NA_real_),
    auc         = auc_simple(y01, p),
    brier       = brier(y01, p)
  )
}

# Inverse-Gamma sampler with density proportional to x^(-shape-1) exp(-rate/x)
rinvgamma <- function(n, shape, rate) {
  1 / rgamma(n, shape = shape, rate = rate)
}

3 Core Bayesian engine: conjugate Bayesian linear regression

Many scientific models begin with a regression layer. We use a simple conjugate model:

\[ \mathbf y \mid \boldsymbol\beta,\sigma^2 \sim N(\mathbf X\boldsymbol\beta,\sigma^2 I), \]

\[ \boldsymbol\beta\mid \sigma^2 \sim N(\boldsymbol\beta_0,\sigma^2 V_0), \qquad \sigma^2\sim \mathrm{Inv\text{-}Gamma}(a_0,b_0). \]

This is not the most modern Bayesian model, but it is excellent pedagogically because the posterior can be computed analytically.

nig_posterior <- function(X, y, beta0 = NULL, V0 = NULL, a0 = 2, b0 = 1) {
  X <- as.matrix(X)
  y <- as.numeric(y)
  n <- nrow(X)
  p <- ncol(X)

  if (is.null(beta0)) beta0 <- rep(0, p)
  if (is.null(V0)) V0 <- diag(100, p)

  beta0 <- matrix(beta0, ncol = 1)
  V0_inv <- solve(V0)

  Vn <- solve(V0_inv + t(X) %*% X)
  beta_n <- Vn %*% (V0_inv %*% beta0 + t(X) %*% y)
  an <- a0 + n / 2

  quad0 <- as.numeric(t(beta0) %*% V0_inv %*% beta0)
  quadn <- as.numeric(t(beta_n) %*% solve(Vn) %*% beta_n)
  bn <- as.numeric(b0 + 0.5 * (sum(y^2) + quad0 - quadn))

  list(beta_n = as.vector(beta_n), Vn = Vn, an = an, bn = bn,
       colnames = colnames(X))
}

sample_lm_posterior <- function(post, ndraw = 4000) {
  p <- length(post$beta_n)
  sigma2 <- rinvgamma(ndraw, shape = post$an, rate = post$bn)
  Z <- matrix(rnorm(ndraw * p), nrow = ndraw, ncol = p) %*% chol(post$Vn)
  beta <- sweep(Z, 1, sqrt(sigma2), "*")
  beta <- sweep(beta, 2, post$beta_n, "+")
  colnames(beta) <- post$colnames
  list(beta = beta, sigma2 = sigma2)
}

predict_lm_posterior <- function(post, Xnew, ndraw = 4000, predictive = TRUE) {
  Xnew <- as.matrix(Xnew)
  draws <- sample_lm_posterior(post, ndraw = ndraw)
  mu <- Xnew %*% t(draws$beta)
  if (predictive) {
    eps <- matrix(rnorm(nrow(Xnew) * ndraw), nrow = nrow(Xnew), ncol = ndraw)
    ydraw <- mu + sweep(eps, 2, sqrt(draws$sigma2), "*")
  } else {
    ydraw <- mu
  }
  tibble(
    fit = rowMeans(ydraw),
    lo  = apply(ydraw, 1, quantile, 0.025),
    hi  = apply(ydraw, 1, quantile, 0.975)
  )
}

4 Core Bayesian engine: approximate Bayesian logistic regression

For binary scientific decisions, e.g. strong vs ordinary earthquake, species A vs not species A, stable vs unstable material, we often use logistic models. Exact Bayesian logistic regression usually needs MCMC. For teaching, we use a Laplace approximation around the posterior mode.

Prior:

\[ \boldsymbol\beta \sim N(0,s^2I). \]

Likelihood:

\[ Y_i\mid \boldsymbol\beta \sim \mathrm{Bernoulli}(p_i), \qquad \mathrm{logit}(p_i)=\mathbf x_i^\top\boldsymbol\beta. \]

log1pexp <- function(eta) {
  ifelse(eta > 0, eta + log1p(exp(-eta)), log1p(exp(eta)))
}

bayes_logistic_laplace <- function(X, y01, prior_sd = 2.5) {
  X <- as.matrix(X)
  y01 <- as.numeric(y01)
  p <- ncol(X)

  neg_log_post <- function(beta) {
    eta <- as.vector(X %*% beta)
    loglik <- sum(y01 * eta - log1pexp(eta))
    logprior <- sum(dnorm(beta, mean = 0, sd = prior_sd, log = TRUE))
    -(loglik + logprior)
  }

  opt <- optim(rep(0, p), neg_log_post, hessian = TRUE, method = "BFGS",
               control = list(maxit = 2000))

  H <- opt$hessian
  V <- tryCatch(solve(H + diag(1e-6, p)), error = function(e) MASS::ginv(H))

  list(mean = opt$par, cov = V, colnames = colnames(X), converged = opt$convergence == 0)
}

predict_bayes_logistic <- function(fit, Xnew, ndraw = 4000) {
  Xnew <- as.matrix(Xnew)
  beta_draws <- MASS::mvrnorm(ndraw, mu = fit$mean, Sigma = fit$cov)
  prob_draws <- plogis(Xnew %*% t(beta_draws))
  tibble(
    prob = rowMeans(prob_draws),
    lo = apply(prob_draws, 1, quantile, 0.025),
    hi = apply(prob_draws, 1, quantile, 0.975)
  )
}

5 Theory: why misspecification matters

A model is misspecified when the mathematical model does not match the scientific data-generating process. This is not a rare event; it is the normal state of applied science.

Common misspecifications:

Scientific setting Possible misspecification Why ML may fail Why classical statistics may fail What a joint approach does
Astronomy time series Random split ignores time dependence Overoptimistic validation Too-rigid AR model misses nonlinearity Chronological validation + probabilistic forecasting
Geoscience hazards Rare-event imbalance High accuracy but low sensitivity Wrong link or omitted spatial variables ML screening + Bayesian calibration
Physical chemistry Wrong scale/physics Tree extrapolates poorly Linear model on raw scale wrong Physics transform + uncertainty
Biology/ecology Small training sample Flexible classifier overfits Strong parametric assumptions wrong Regularized classifier + uncertainty flags

A useful slogan:

ML is excellent at finding signal. Statistics is essential for asking whether the signal is reliable, transportable, and scientifically interpretable.

6 Simulation 1: flexible ML interpolates, but extrapolation can fail

We simulate a nonlinear scientific response. We train in a limited experimental range and then test beyond it.

set.seed(1)

n <- 180
x <- runif(n, -2.5, 2.5)
y <- sin(1.5 * x) + 0.35 * x + rnorm(n, sd = 0.25)
train_sim <- tibble(x = x, y = y)

grid_sim <- tibble(x = seq(-5, 5, length.out = 400))

# A misspecified statistical model: linear regression
fit_lm_wrong <- lm(y ~ x, data = train_sim)

# A flexible ML model: regression tree
fit_tree <- rpart(y ~ x, data = train_sim, control = rpart.control(cp = 0.002, minsplit = 8))

# A Bayesian polynomial model: still imperfect, but with uncertainty
X_train_poly <- model.matrix(~ x + I(x^2) + I(x^3), data = train_sim)
X_grid_poly  <- model.matrix(~ x + I(x^2) + I(x^3), data = grid_sim)
post_poly <- nig_posterior(X_train_poly, train_sim$y, V0 = diag(20, ncol(X_train_poly)))
pred_poly <- predict_lm_posterior(post_poly, X_grid_poly, ndraw = 3000)

plot_sim <- grid_sim %>%
  mutate(
    true = sin(1.5 * x) + 0.35 * x,
    lm = predict(fit_lm_wrong, newdata = grid_sim),
    tree = predict(fit_tree, newdata = grid_sim),
    bayes_poly = pred_poly$fit,
    bayes_lo = pred_poly$lo,
    bayes_hi = pred_poly$hi,
    region = ifelse(x >= min(train_sim$x) & x <= max(train_sim$x), "inside training range", "extrapolation")
  )

ggplot() +
  geom_rect(data = tibble(xmin = min(train_sim$x), xmax = max(train_sim$x), ymin = -Inf, ymax = Inf),
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            fill = "grey90", alpha = 0.6) +
  geom_point(data = train_sim, aes(x, y), alpha = 0.45, size = 1.5) +
  geom_line(data = plot_sim, aes(x, true), linewidth = 1.1, color = "black") +
  geom_line(data = plot_sim, aes(x, lm), linewidth = 0.9, linetype = 2, color = "firebrick") +
  geom_line(data = plot_sim, aes(x, tree), linewidth = 0.9, color = "steelblue") +
  geom_ribbon(data = plot_sim, aes(x = x, ymin = bayes_lo, ymax = bayes_hi), alpha = 0.18) +
  geom_line(data = plot_sim, aes(x, bayes_poly), linewidth = 0.9, color = "darkgreen") +
  labs(
    title = "Simulation: interpolation success does not guarantee extrapolation success",
    subtitle = "Grey band = training range; black = true curve; blue = ML tree; green = Bayesian polynomial; red dashed = wrong linear model",
    x = "scientific input x", y = "response y"
  )

6.1 Lesson

  • The wrong linear statistical model fails because it is too rigid.
  • The ML tree adapts well inside the observed range, but extrapolates as a flat step function.
  • The Bayesian polynomial is not guaranteed true, but at least gives uncertainty that grows outside the central data region.

Scientific message: flexible ML needs careful extrapolation diagnostics; Bayesian modelling needs good scientific structure. Neither should be used blindly.

7 Simulation 2: confounding and distribution shift

Now we show a common data-science trap. A hidden group variable changes both the input distribution and the response baseline. If ignored, a model may learn the wrong relationship.

set.seed(2)

n_train <- 250
n_test <- 250

make_confounded <- function(n, group_prob) {
  g <- rbinom(n, 1, group_prob)
  x <- rnorm(n, mean = 1.2 * g, sd = 1)
  y <- 1 + 1.0 * x + 2.8 * g + rnorm(n, sd = 1)
  tibble(x = x, group = factor(g), y = y)
}

train_c <- make_confounded(n_train, group_prob = 0.25)
test_c  <- make_confounded(n_test,  group_prob = 0.75)

fit_ignore <- lm(y ~ x, data = train_c)
fit_oracle <- lm(y ~ x + group, data = train_c)
fit_tree_c <- rpart(y ~ x, data = train_c)

# Bayesian regression with group included and weakly informative prior
X_train_c <- model.matrix(~ x + group, data = train_c)
X_test_c  <- model.matrix(~ x + group, data = test_c)
post_c <- nig_posterior(X_train_c, train_c$y, V0 = diag(10, ncol(X_train_c)))
pred_c <- predict_lm_posterior(post_c, X_test_c, ndraw = 3000)

perf_c <- tibble(
  model = c("LM ignoring group", "LM with group", "ML tree ignoring group", "Bayesian with group"),
  RMSE = c(
    rmse(test_c$y, predict(fit_ignore, newdata = test_c)),
    rmse(test_c$y, predict(fit_oracle, newdata = test_c)),
    rmse(test_c$y, predict(fit_tree_c, newdata = test_c)),
    rmse(test_c$y, pred_c$fit)
  ),
  MAE = c(
    mae(test_c$y, predict(fit_ignore, newdata = test_c)),
    mae(test_c$y, predict(fit_oracle, newdata = test_c)),
    mae(test_c$y, predict(fit_tree_c, newdata = test_c)),
    mae(test_c$y, pred_c$fit)
  )
)

knitr::kable(perf_c, digits = 3, caption = "Distribution-shift test: group composition changes between training and test data.")
Distribution-shift test: group composition changes between training and test data.
model RMSE MAE
LM ignoring group 1.673 1.390
LM with group 0.958 0.781
ML tree ignoring group 1.798 1.497
Bayesian with group 0.960 0.784 
test_plot_c <- test_c %>%
  mutate(
    pred_ignore = predict(fit_ignore, newdata = test_c),
    pred_bayes = pred_c$fit,
    lo = pred_c$lo,
    hi = pred_c$hi
  )

ggplot(test_plot_c, aes(x, y, color = group)) +
  geom_point(alpha = 0.55) +
  geom_line(aes(y = pred_ignore), color = "firebrick", linewidth = 1, linetype = 2) +
  geom_line(aes(y = pred_bayes), color = "darkgreen", linewidth = 1) +
  labs(
    title = "Confounding simulation: ignoring group creates transport failure",
    subtitle = "Red dashed = model ignoring group; green = Bayesian regression with group adjustment",
    x = "x", y = "response"
  )

7.1 Lesson

The issue is not that ML is bad or statistics is bad. The issue is that the target was misspecified. If the model ignores the group variable, it learns a relationship that does not transport well when the group composition changes.

This is the same broad logic behind many failures in scientific ML: site effects in medical imaging, instrument effects in spectroscopy, survey selection effects in astronomy, and training-domain effects in materials prediction.

8 Real data 1: physical chemistry/materials — vapor pressure of mercury

The built-in R dataset pressure contains 19 observations on the relationship between temperature in degrees Celsius and vapor pressure of mercury in millimeters of mercury.

Scientific motivation: in physical chemistry, modelling on the wrong scale can create a bad model. Here the raw relation is highly nonlinear. A physics-aware transform, such as modelling log pressure against inverse absolute temperature, is closer to the Clausius–Clapeyron intuition.

data("pressure")
press <- as_tibble(pressure) %>%
  mutate(
    temp_K = temperature + 273.15,
    invK = 1 / temp_K,
    log_pressure = log(pressure)
  )

knitr::kable(head(press), digits = 4, caption = "Mercury vapor pressure data from R's built-in pressure dataset.")
Mercury vapor pressure data from R’s built-in pressure dataset.
temperature pressure temp_K invK log_pressure
0 0.0002 273.15 0.0037 -8.5172
20 0.0012 293.15 0.0034 -6.7254
40 0.0060 313.15 0.0032 -5.1160
60 0.0300 333.15 0.0030 -3.5066
80 0.0900 353.15 0.0028 -2.4079
100 0.2700 373.15 0.0027 -1.3093 
ggplot(press, aes(temperature, pressure)) +
  geom_point(size = 2.5) +
  geom_line() +
  scale_y_log10() +
  labs(title = "Mercury vapor pressure grows nonlinearly with temperature",
       subtitle = "A log scale reveals structure hidden on the raw scale.",
       x = "Temperature (degree C)", y = "Pressure (mm Hg, log scale)")

8.1 Train-low, predict-high: an extrapolation challenge

train_p <- press %>% filter(temperature <= 250)
test_p  <- press %>% filter(temperature > 250)

# Wrong statistical model: log pressure linear in Celsius temperature
fit_p_wrong <- lm(log_pressure ~ temperature, data = train_p)

# ML tree on the same training data
fit_p_tree <- rpart(log_pressure ~ temperature, data = train_p,
                    control = rpart.control(cp = 0.001, minsplit = 4))

# Physics-aware Bayesian regression: log pressure linear in inverse Kelvin
X_train_p <- model.matrix(~ invK, data = train_p)
X_test_p  <- model.matrix(~ invK, data = test_p)
post_p <- nig_posterior(X_train_p, train_p$log_pressure, V0 = diag(100, ncol(X_train_p)))
pred_p <- predict_lm_posterior(post_p, X_test_p, ndraw = 5000)

pred_table_p <- test_p %>%
  transmute(
    temperature,
    observed_log_pressure = log_pressure,
    wrong_linear = predict(fit_p_wrong, newdata = test_p),
    tree_ML = predict(fit_p_tree, newdata = test_p),
    bayes_physics = pred_p$fit,
    bayes_lo = pred_p$lo,
    bayes_hi = pred_p$hi
  )

knitr::kable(pred_table_p, digits = 3, caption = "High-temperature extrapolation on log-pressure scale.")
High-temperature extrapolation on log-pressure scale.
temperature observed_log_pressure wrong_linear tree_ML bayes_physics bayes_lo bayes_hi
260 4.564 6.088 3.454 -0.993 -8.295 6.388
280 5.056 7.104 3.454 -0.994 -8.436 6.472
300 5.509 8.119 3.454 -0.935 -8.325 6.353
320 5.930 9.135 3.454 -0.962 -8.400 6.519
340 6.324 10.151 3.454 -0.992 -8.546 6.196
360 6.692 11.166 3.454 -0.983 -8.120 6.505 
metrics_p <- tibble(
  model = c("Wrong linear in Celsius", "ML tree", "Bayesian physics transform"),
  RMSE_log_pressure = c(
    rmse(test_p$log_pressure, pred_table_p$wrong_linear),
    rmse(test_p$log_pressure, pred_table_p$tree_ML),
    rmse(test_p$log_pressure, pred_table_p$bayes_physics)
  )
)
knitr::kable(metrics_p, digits = 3, caption = "Extrapolation error on the held-out high-temperature region.")
Extrapolation error on the held-out high-temperature region.
model RMSE_log_pressure
Wrong linear in Celsius 3.116
ML tree 2.341
Bayesian physics transform 6.695 
grid_p <- tibble(
  temperature = seq(min(press$temperature), max(press$temperature), length.out = 200)
) %>%
  mutate(temp_K = temperature + 273.15,
         invK = 1 / temp_K)

X_grid_p <- model.matrix(~ invK, data = grid_p)
pred_grid_p <- predict_lm_posterior(post_p, X_grid_p, ndraw = 5000)

grid_p <- grid_p %>%
  mutate(
    wrong_linear = predict(fit_p_wrong, newdata = grid_p),
    tree_ML = predict(fit_p_tree, newdata = grid_p),
    bayes = pred_grid_p$fit,
    lo = pred_grid_p$lo,
    hi = pred_grid_p$hi
  )

ggplot() +
  geom_point(data = press, aes(temperature, log_pressure), size = 2.4) +
  geom_vline(xintercept = max(train_p$temperature), linetype = 2) +
  geom_line(data = grid_p, aes(temperature, wrong_linear), color = "firebrick", linetype = 2, linewidth = 0.9) +
  geom_line(data = grid_p, aes(temperature, tree_ML), color = "steelblue", linewidth = 0.9) +
  geom_ribbon(data = grid_p, aes(temperature, ymin = lo, ymax = hi), alpha = 0.18) +
  geom_line(data = grid_p, aes(temperature, bayes), color = "darkgreen", linewidth = 1) +
  labs(
    title = "Physics-aware Bayes vs ML tree under extrapolation",
    subtitle = "Vertical dashed line separates training and test region; green ribbon = Bayesian predictive interval",
    x = "Temperature (degree C)", y = "log pressure"
  )

8.2 Lesson

  • The tree can interpolate but cannot physically extrapolate.
  • The wrong statistical model can also fail if the scientific scale is wrong.
  • A simple Bayesian model with a meaningful physical transform gives prediction plus uncertainty.

9 Real data 2: geoscience — Fiji earthquakes

The built-in R dataset quakes gives locations of 1000 seismic events of magnitude greater than 4.0 near Fiji since 1964. Variables include latitude, longitude, depth, magnitude, and number of reporting stations.

Our toy question:

Can we classify relatively strong events, say magnitude at least 5.0, from location, depth, and station information?

This is a toy problem, not an operational earthquake warning model. The point is to illustrate rare-event classification and calibration.

data("quakes")
q <- as_tibble(quakes) %>%
  mutate(
    strong = ifelse(mag >= 5.0, 1, 0),
    strong_label = factor(ifelse(strong == 1, "strong", "ordinary")),
    depth_s = as.numeric(scale(depth)),
    stations_s = as.numeric(scale(stations)),
    lat_s = as.numeric(scale(lat)),
    long_s = as.numeric(scale(long))
  )

knitr::kable(q %>% count(strong_label), caption = "Class balance: strong vs ordinary seismic events.")
Class balance: strong vs ordinary seismic events.
strong_label n
ordinary 802
strong 198 
ggplot(q, aes(long, lat, color = mag, size = depth)) +
  geom_point(alpha = 0.7) +
  scale_color_viridis_c() +
  scale_size_continuous(range = c(1, 4)) +
  labs(title = "Fiji seismic events: location, magnitude, and depth",
       x = "Longitude", y = "Latitude", color = "Magnitude", size = "Depth")

9.1 ML tree vs Bayesian logistic vs hybrid calibration

We split into three parts:

  • training set: fit the ML tree;
  • calibration set: calibrate the ML probability using Bayesian logistic regression;
  • test set: evaluate honestly.
set.seed(3)
idx <- sample(seq_len(nrow(q)))
n_train <- floor(0.60 * nrow(q))
n_cal   <- floor(0.20 * nrow(q))

train_q <- q[idx[1:n_train], ]
cal_q   <- q[idx[(n_train + 1):(n_train + n_cal)], ]
test_q  <- q[idx[(n_train + n_cal + 1):nrow(q)], ]

# ML tree: flexible but not automatically calibrated
fit_q_tree <- rpart(strong_label ~ depth_s + stations_s + lat_s + long_s,
                    data = train_q, method = "class",
                    control = rpart.control(cp = 0.01, minsplit = 20))

prob_tree_cal <- predict(fit_q_tree, newdata = cal_q, type = "prob")[, "strong"]
prob_tree_test <- predict(fit_q_tree, newdata = test_q, type = "prob")[, "strong"]

# Bayesian logistic model using scientific covariates directly
X_train_q <- model.matrix(~ depth_s + stations_s + lat_s + long_s, data = train_q)
X_test_q  <- model.matrix(~ depth_s + stations_s + lat_s + long_s, data = test_q)
fit_q_bayes <- bayes_logistic_laplace(X_train_q, train_q$strong, prior_sd = 2.5)
pred_q_bayes <- predict_bayes_logistic(fit_q_bayes, X_test_q, ndraw = 4000)

# Hybrid: Bayesian calibration of ML probability using calibration set
# Add a small epsilon to avoid infinite logits.
eps <- 1e-4
cal_q <- cal_q %>% mutate(tree_logit = qlogis(pmin(pmax(prob_tree_cal, eps), 1 - eps)))
test_q <- test_q %>% mutate(tree_logit = qlogis(pmin(pmax(prob_tree_test, eps), 1 - eps)))

X_cal_hybrid  <- model.matrix(~ tree_logit + depth_s + stations_s, data = cal_q)
X_test_hybrid <- model.matrix(~ tree_logit + depth_s + stations_s, data = test_q)
fit_q_hybrid <- bayes_logistic_laplace(X_cal_hybrid, cal_q$strong, prior_sd = 2.5)
pred_q_hybrid <- predict_bayes_logistic(fit_q_hybrid, X_test_hybrid, ndraw = 4000)

metrics_q <- bind_rows(
  classification_table(test_q$strong, prob_tree_test) %>% mutate(model = "ML tree"),
  classification_table(test_q$strong, pred_q_bayes$prob) %>% mutate(model = "Bayesian logistic"),
  classification_table(test_q$strong, pred_q_hybrid$prob) %>% mutate(model = "Hybrid: ML + Bayesian calibration")
) %>% dplyr::select(model, dplyr::everything())

knitr::kable(metrics_q, digits = 3, caption = "Test-set classification metrics. Accuracy alone can be misleading in rare-event problems.")
Test-set classification metrics. Accuracy alone can be misleading in rare-event problems.
model accuracy sensitivity specificity auc brier
ML tree 0.950 0.788 0.982 0.917 0.045
Bayesian logistic 0.945 0.788 0.976 0.983 0.042
Hybrid: ML + Bayesian calibration 0.945 0.727 0.988 0.985 0.039 
plot_q <- test_q %>%
  mutate(
    p_tree = prob_tree_test,
    p_bayes = pred_q_bayes$prob,
    p_hybrid = pred_q_hybrid$prob,
    hybrid_lo = pred_q_hybrid$lo,
    hybrid_hi = pred_q_hybrid$hi
  )

p1 <- ggplot(plot_q, aes(p_tree, p_hybrid, color = factor(strong))) +
  geom_point(alpha = 0.7) +
  geom_abline(linetype = 2) +
  labs(title = "Hybrid calibration: ML probability adjusted by Bayesian layer",
       x = "ML tree probability", y = "Hybrid calibrated probability", color = "Strong event")

p2 <- ggplot(plot_q, aes(long, lat, color = p_hybrid, size = mag)) +
  geom_point(alpha = 0.75) +
  scale_color_viridis_c(labels = percent) +
  scale_size_continuous(range = c(1, 4)) +
  labs(title = "Hybrid predicted probability on seismic map",
       x = "Longitude", y = "Latitude", color = "P(strong)", size = "Magnitude")

p1

p2

9.2 Lesson

  • ML gives a useful nonlinear screening probability.
  • Bayesian logistic regression gives interpretable coefficients and posterior uncertainty.
  • The hybrid model treats ML output as a feature and then calibrates it statistically.

This is a practical workflow for scientific data science:

\[ \text{ML prediction} \quad + \quad \text{Bayesian calibration} \quad \Rightarrow \quad \text{usable uncertainty-aware risk score}. \]

10 Real data 3: astronomical/solar time series — sunspot numbers

The built-in R dataset sunspot.year contains yearly sunspot numbers from 1700 to 1988. Solar activity has cyclic structure, so this is a nice example where random train-test splits can be misleading.

data("sunspot.year")
sun <- tibble(
  year = as.numeric(time(sunspot.year)),
  sunspots = as.numeric(sunspot.year)
) %>%
  mutate(
    lag1 = lag(sunspots, 1),
    lag2 = lag(sunspots, 2),
    lag11 = lag(sunspots, 11)
  ) %>%
  drop_na()

ggplot(sun, aes(year, sunspots)) +
  geom_line(color = "darkorange", linewidth = 0.8) +
  labs(title = "Yearly sunspot numbers: temporal dependence is not optional",
       x = "Year", y = "Sunspot number")

10.1 The validation trap: random split vs chronological split

If we randomly split a time series, training data from the future may help predict the past. That is leakage.

set.seed(4)

# Chronological split
train_s_chrono <- sun %>% filter(year <= 1950)
test_s_chrono  <- sun %>% filter(year > 1950)

# Random split
id_s <- sample(seq_len(nrow(sun)), size = floor(0.7 * nrow(sun)))
train_s_random <- sun[id_s, ]
test_s_random  <- sun[-id_s, ]

fit_tree_chrono <- rpart(sunspots ~ lag1 + lag2 + lag11, data = train_s_chrono)
fit_tree_random <- rpart(sunspots ~ lag1 + lag2 + lag11, data = train_s_random)

rmse_chrono <- rmse(test_s_chrono$sunspots, predict(fit_tree_chrono, newdata = test_s_chrono))
rmse_random <- rmse(test_s_random$sunspots, predict(fit_tree_random, newdata = test_s_random))

knitr::kable(
  tibble(
    validation = c("Chronological split", "Random split"),
    RMSE = c(rmse_chrono, rmse_random)
  ),
  digits = 2,
  caption = "Random validation can be optimistic for dependent scientific time series."
)
Random validation can be optimistic for dependent scientific time series.
validation RMSE
Chronological split 35.42
Random split 20.44

10.2 Bayesian autoregression with predictive uncertainty

X_train_s <- model.matrix(~ lag1 + lag2 + lag11, data = train_s_chrono)
X_test_s  <- model.matrix(~ lag1 + lag2 + lag11, data = test_s_chrono)
post_s <- nig_posterior(X_train_s, train_s_chrono$sunspots,
                        V0 = diag(100, ncol(X_train_s)), a0 = 2, b0 = 50)
pred_s <- predict_lm_posterior(post_s, X_test_s, ndraw = 5000)

fit_lm_s <- lm(sunspots ~ lag1 + lag2 + lag11, data = train_s_chrono)

sun_pred <- test_s_chrono %>%
  mutate(
    lm_pred = predict(fit_lm_s, newdata = test_s_chrono),
    tree_pred = predict(fit_tree_chrono, newdata = test_s_chrono),
    bayes = pred_s$fit,
    lo = pred_s$lo,
    hi = pred_s$hi
  )

metrics_s <- tibble(
  model = c("Linear AR-style model", "ML tree", "Bayesian AR-style model"),
  RMSE = c(
    rmse(sun_pred$sunspots, sun_pred$lm_pred),
    rmse(sun_pred$sunspots, sun_pred$tree_pred),
    rmse(sun_pred$sunspots, sun_pred$bayes)
  )
)
knitr::kable(metrics_s, digits = 2, caption = "Chronological test performance for sunspot prediction.")
Chronological test performance for sunspot prediction.
model RMSE
Linear AR-style model 23.56
ML tree 35.42
Bayesian AR-style model 23.55 
ggplot(sun_pred, aes(year, sunspots)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), fill = "darkgreen", alpha = 0.18) +
  geom_line(color = "black", linewidth = 0.8) +
  geom_line(aes(y = tree_pred), color = "steelblue", linewidth = 0.8) +
  geom_line(aes(y = bayes), color = "darkgreen", linewidth = 0.8) +
  labs(title = "Sunspot prediction: ML point prediction and Bayesian uncertainty",
       subtitle = "Black = observed; blue = ML tree; green = Bayesian AR-style predictive mean and interval",
       x = "Year", y = "Sunspot number")

10.3 Lesson

  • Time series validation must respect time.
  • ML point predictions are useful, but uncertainty is essential for scientific forecasting.
  • A simple Bayesian AR-style model may be misspecified, but it gives predictive intervals and makes its uncertainty visible.

11 Real data 4: biological morphology — Palmer penguins

The palmerpenguins dataset contains size measurements for adult foraging penguins near Palmer Station, Antarctica. It has 344 rows and variables such as species, island, bill length, bill depth, flipper length, body mass, sex, and year.

We use it as a toy biological morphology problem:

Can we classify species from bill and flipper measurements? And can we identify uncertain cases?

data("penguins", package = "palmerpenguins")
pg <- penguins %>%
  drop_na(species, bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g) %>%
  mutate(species = droplevels(species))

knitr::kable(pg %>% count(species), caption = "Palmer penguins species counts.")
Palmer penguins species counts.
species n
Adelie 151
Chinstrap 68
Gentoo 123 
ggplot(pg, aes(bill_length_mm, bill_depth_mm, color = species)) +
  geom_point(alpha = 0.8, size = 2.3) +
  labs(title = "Morphological separation of penguin species",
       x = "Bill length (mm)", y = "Bill depth (mm)")

11.1 Small-data classification: ML tree vs regularized Bayesian Gaussian classifier

We deliberately use a small training set: 15 individuals per species. This mimics many natural-science projects where data collection is expensive.

log_dmvnorm <- function(X, mu, Sigma) {
  X <- as.matrix(X)
  p <- ncol(X)
  R <- chol(Sigma)
  z <- backsolve(R, t(X) - mu, transpose = TRUE)
  -0.5 * p * log(2 * pi) - sum(log(diag(R))) - 0.5 * colSums(z^2)
}

fit_eb_gaussian_classifier <- function(data, features, class_var, kappa = 5, ridge = 0.25) {
  X <- as.matrix(data[, features])
  y <- data[[class_var]]
  classes <- levels(y)
  global_mu <- colMeans(X)
  pooled_cov <- cov(X) + diag(ridge, ncol(X))

  params <- lapply(classes, function(cl) {
    Xc <- X[y == cl, , drop = FALSE]
    nc <- nrow(Xc)
    mu_raw <- colMeans(Xc)
    # Empirical-Bayes shrinkage of class mean toward global mean
    mu_shrunk <- (nc * mu_raw + kappa * global_mu) / (nc + kappa)
    Sigma <- if (nc > ncol(X) + 2) cov(Xc) else pooled_cov
    Sigma <- 0.7 * Sigma + 0.3 * pooled_cov + diag(ridge, ncol(X))
    list(class = cl, prior = nc / nrow(X), mu = mu_shrunk, Sigma = Sigma)
  })
  names(params) <- classes
  list(params = params, features = features, classes = classes)
}

predict_eb_gaussian <- function(fit, newdata) {
  X <- as.matrix(newdata[, fit$features])
  logprob <- sapply(fit$params, function(par) {
    log(par$prior) + log_dmvnorm(X, par$mu, par$Sigma)
  })
  maxlp <- apply(logprob, 1, max)
  prob <- exp(logprob - maxlp)
  prob <- prob / rowSums(prob)
  pred <- colnames(prob)[max.col(prob)]
  tibble(pred = factor(pred, levels = fit$classes), confidence = apply(prob, 1, max)) %>%
    bind_cols(as_tibble(prob, .name_repair = "minimal"))
}
set.seed(5)
features_pg <- c("bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g")

train_pg <- pg %>%
  group_by(species) %>%
  slice_sample(n = 15) %>%
  ungroup()

test_pg <- anti_join(pg, train_pg, by = colnames(pg))

fit_pg_tree <- rpart(species ~ bill_length_mm + bill_depth_mm + flipper_length_mm + body_mass_g,
                     data = train_pg, method = "class",
                     control = rpart.control(cp = 0.001, minsplit = 5))
pred_tree_pg <- predict(fit_pg_tree, newdata = test_pg, type = "class")
prob_tree_pg <- predict(fit_pg_tree, newdata = test_pg, type = "prob")
conf_tree_pg <- apply(prob_tree_pg, 1, max)

fit_pg_bayes <- fit_eb_gaussian_classifier(train_pg, features_pg, "species", kappa = 8, ridge = 0.5)
pred_bayes_pg <- predict_eb_gaussian(fit_pg_bayes, test_pg)

metrics_pg <- tibble(
  model = c("ML tree", "Regularized Bayesian Gaussian classifier"),
  accuracy = c(
    mean(pred_tree_pg == test_pg$species),
    mean(pred_bayes_pg$pred == test_pg$species)
  ),
  mean_confidence = c(mean(conf_tree_pg), mean(pred_bayes_pg$confidence)),
  uncertain_fraction_below_70pct = c(mean(conf_tree_pg < 0.70), mean(pred_bayes_pg$confidence < 0.70))
)
knitr::kable(metrics_pg, digits = 3, caption = "Small-data species classification: accuracy plus uncertainty/confidence.")
Small-data species classification: accuracy plus uncertainty/confidence.
model accuracy mean_confidence uncertain_fraction_below_70pct
ML tree 0.912 0.945 0.000
Regularized Bayesian Gaussian classifier 0.976 0.907 0.071 
plot_pg <- test_pg %>%
  mutate(
    pred_tree = pred_tree_pg,
    conf_tree = conf_tree_pg,
    pred_bayes = pred_bayes_pg$pred,
    conf_bayes = pred_bayes_pg$confidence,
    bayes_correct = pred_bayes == species
  )

# Show the most uncertain Bayesian cases
uncertain_pg <- plot_pg %>%
  arrange(conf_bayes) %>%
  dplyr::select(species, pred_bayes, conf_bayes, bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g) %>%
  head(10)

knitr::kable(uncertain_pg, digits = 2, caption = "Most uncertain Bayesian classifications: these are scientifically interesting cases to inspect.")
Most uncertain Bayesian classifications: these are scientifically interesting cases to inspect.
species pred_bayes conf_bayes bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
Adelie Chinstrap 0.46 44.1 18.0 210 4000
Chinstrap Adelie 0.47 45.2 17.8 198 3950
Gentoo Gentoo 0.49 44.5 14.3 216 4100
Chinstrap Adelie 0.50 42.4 17.3 181 3600
Chinstrap Adelie 0.54 42.5 16.7 187 3350
Chinstrap Adelie 0.55 46.0 18.9 195 4150
Chinstrap Chinstrap 0.55 45.6 19.4 194 3525
Chinstrap Chinstrap 0.58 45.7 17.0 195 3650
Adelie Adelie 0.58 45.8 18.9 197 4150
Chinstrap Adelie 0.59 42.5 17.3 187 3350 
ggplot(plot_pg, aes(bill_length_mm, bill_depth_mm)) +
  geom_point(aes(color = species, shape = bayes_correct, size = 1 - conf_bayes), alpha = 0.8) +
  scale_size_continuous(range = c(1.5, 5), name = "Bayesian uncertainty") +
  labs(title = "Bayesian classifier identifies uncertain biological specimens",
       subtitle = "Large points are less certain; shape indicates whether prediction was correct",
       x = "Bill length (mm)", y = "Bill depth (mm)", color = "True species")

11.2 Lesson

ML can give accurate classification. Bayesian/regularized statistical classifiers add a second layer: which cases are uncertain?

For research, uncertain cases are not merely errors. They may be:

  • boundary morphologies,
  • measurement anomalies,
  • biologically mixed signals,
  • or samples requiring expert review.

12 A counterexample where statistics alone fails

Let us revisit the pressure dataset. A linear statistical model on raw pressure vs temperature is interpretable but scientifically wrong for extrapolation.

fit_raw_lm <- lm(pressure ~ temperature, data = train_p)
raw_pred <- predict(fit_raw_lm, newdata = test_p)

knitr::kable(
  tibble(
    temperature = test_p$temperature,
    observed_pressure = test_p$pressure,
    raw_linear_prediction = raw_pred
  ),
  digits = 3,
  caption = "A simple statistical model can fail badly if the scale and physics are wrong."
)
A simple statistical model can fail badly if the scale and physics are wrong.
temperature observed_pressure raw_linear_prediction
260 96 32.791
280 157 36.131
300 247 39.470
320 376 42.810
340 558 46.149
360 806 49.489

This is a useful warning:

Statistics is not automatically good science. A clean formula with the wrong scientific scale can be worse than a flexible ML model.

13 A counterexample where ML alone fails

A regression tree cannot extrapolate beyond the observed response structure. In physical science, extrapolation matters: high temperature, rare extremes, new materials, higher energy regimes, unseen survey conditions.

knitr::kable(metrics_p, digits = 3,
             caption = "ML tree can be weaker in extrapolation when the scientific law has not been encoded.")
ML tree can be weaker in extrapolation when the scientific law has not been encoded.
model RMSE_log_pressure
Wrong linear in Celsius 3.116
ML tree 2.341
Bayesian physics transform 6.695

This is another useful warning:

ML is not automatically scientific. It may learn the training domain very well and still fail where the scientific question is most important.

14 How to combine statistics and ML wisely

A practical scientific workflow:

  1. Start with the scientific question. What is the estimand? Prediction, effect, mechanism, risk, or decision?
  2. Use ML for flexible pattern discovery. Trees, random forests, neural networks, kernels, embeddings.
  3. Use statistical modelling for uncertainty and structure. Priors, likelihoods, censoring, hierarchy, calibration, checking.
  4. Use domain knowledge. Transform variables using physics, geometry, conservation laws, symmetries, or pathway information.
  5. Validate honestly. Chronological split for time series, spatial split for maps, site split for biomedical imaging, out-of-distribution stress tests for physical extremes.
  6. Report uncertainty. Credible intervals, calibration curves, posterior predictive checks, sensitivity analysis.

15 Summary table: what each example taught

Part 2 take-home map.
Example Domain ML role Bayesian/statistical role Main warning
Nonlinear simulation Generic scientific response Flexible interpolation Predictive uncertainty Interpolation is not extrapolation
Confounding simulation Distribution shift / hidden structure Pattern learning but vulnerable to missing variables Adjustment and uncertainty under group shift Omitted structure breaks transportability
Mercury vapor pressure Physical chemistry / materials Nonlinear learner but poor extrapolation Physics-aware transform and interval prediction Wrong scale defeats clean modelling
Fiji earthquakes Geoscience hazard toy problem Hazard screening probability Calibration and interpretable uncertainty Accuracy alone hides rare-event failure
Sunspot time series Astronomy / solar time series Flexible lag-based prediction Forecast intervals and chronological validation Random split leaks future information
Palmer penguins Biology / ecological morphology Accurate small-data classification Uncertainty flags for ambiguous specimens Accuracy without uncertainty hides boundary cases

16 Final message for students

The future of data science in natural science is not “statistics versus ML”.

It is:

\[ \textbf{scientific structure} + \textbf{machine learning flexibility} + \textbf{Bayesian/statistical uncertainty} + \textbf{honest validation}. \]

That combination is what makes a model not only accurate, but scientifically credible.

17 References and data notes

  • Gelman, A. et al. (2013). Bayesian Data Analysis, 3rd ed. Chapman & Hall/CRC.
  • Jaynes, E. T. (2003). Probability Theory: The Logic of Science. Cambridge University Press.
  • Breiman, L. (2001). Statistical Modeling: The Two Cultures. Statistical Science, 16(3), 199–231.
  • Lazer, D. et al. (2014). The Parable of Google Flu: Traps in Big Data Analysis. Science, 343(6176), 1203–1205.
  • Carleo, G. et al. (2019). Machine learning and the physical sciences. Reviews of Modern Physics, 91, 045002.
  • R built-in pressure dataset: vapor pressure of mercury as a function of temperature.
  • R built-in quakes dataset: 1000 seismic events near Fiji with magnitude greater than 4.0.
  • R built-in sunspot.year dataset: yearly sunspot numbers from 1700 to 1988.
  • palmerpenguins package: penguin morphometric measurements near Palmer Station, Antarctica.