1 Why this note?

This short tutorial is designed for students who may know regression or machine learning, but may not yet see clearly what Bayesian inference adds.

The main message is not “Bayes versus ML”. The message is:

Machine learning is strong at prediction; Bayesian inference is strong at uncertainty, small-data reasoning, regularisation, and scientific interpretability. In natural-science applications, the most useful workflow often combines both.

We shall use small examples from natural-science style datasets already available in R:

  • Simulation: nonlinear signal under small sample size.
  • Environmental science: New York air-quality ozone data from datasets::airquality.
  • Botany / morphology: Fisher–Anderson iris data from datasets::iris.
  • Zoology / allometry: brain and body weights from MASS::Animals.
  • Biomedical survival: lung cancer survival data from survival::lung.

The examples are intentionally small, so that the statistical ideas are visible.

2 A minimal conceptual map

2.1 Frequentist / classical estimation

A classical estimator is usually a function of data:

\[ \widehat\theta = T(Y_1,\ldots,Y_n). \]

The unknown parameter \(\theta\) is treated as fixed, and uncertainty is described through repeated-sampling behaviour of \(\widehat\theta\).

2.2 Bayesian inference

Bayesian inference treats the unknown quantity as uncertain and updates belief using Bayes’ rule:

\[ \pi(\theta \mid y) \propto L(\theta;y)\,\pi(\theta), \]

where

  • \(L(\theta;y)\) is the likelihood,
  • \(\pi(\theta)\) is the prior,
  • \(\pi(\theta\mid y)\) is the posterior.

A Bayesian answer is often a distribution, not only one number. From the posterior we can compute:

  • posterior mean or median,
  • credible interval,
  • posterior probability of a scientifically meaningful event, e.g. \(P(\beta>0\mid y)\),
  • predictive distribution for a future observation.

2.3 AI/ML perspective

Many ML methods optimise predictive performance:

\[ \widehat f = \arg\min_f \sum_{i=1}^n \ell(y_i,f(x_i)). \]

This is extremely powerful when the main aim is prediction. But in science, we often also need:

  • uncertainty intervals,
  • small-data inference,
  • mechanistic interpretation,
  • robustness to misspecification,
  • careful extrapolation.

That is where Bayesian thinking becomes complementary.

3 Utility functions used throughout

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

train_test_split <- function(n, train_prop = 0.7) {
  sample(seq_len(n), size = floor(train_prop * n))
}

# Conjugate Bayesian Gaussian linear regression:
# y | beta, sigma2 ~ N(X beta, sigma2 I)
# beta | sigma2 ~ N(beta0, sigma2 V0)
# sigma2 ~ Inv-Gamma(a0, b0)
# Here Inv-Gamma is parameterised through 1/sigma2 ~ Gamma(shape=a, rate=b).
bayes_lm_nig <- 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)
  V0_inv <- solve(V0)
  Vn <- solve(V0_inv + crossprod(X))
  betan <- Vn %*% (V0_inv %*% beta0 + crossprod(X, y))
  an <- a0 + n / 2
  bn <- as.numeric(b0 + 0.5 * (crossprod(y) + t(beta0) %*% V0_inv %*% beta0 - t(betan) %*% solve(Vn) %*% betan))
  list(beta_mean = as.numeric(betan), Vn = Vn, an = an, bn = bn, X = X, y = y)
}

draw_bayes_lm <- function(fit, Xnew, S = 4000, include_noise = TRUE) {
  Xnew <- as.matrix(Xnew)
  p <- length(fit$beta_mean)
  sigma2 <- 1 / rgamma(S, shape = fit$an, rate = fit$bn)
  beta_draws <- matrix(NA_real_, nrow = S, ncol = p)
  for (s in seq_len(S)) {
    beta_draws[s, ] <- MASS::mvrnorm(1, mu = fit$beta_mean, Sigma = sigma2[s] * fit$Vn)
  }
  mu_draws <- beta_draws %*% t(Xnew)
  if (include_noise) {
    y_draws <- mu_draws + matrix(rnorm(length(mu_draws), sd = sqrt(rep(sigma2, ncol(Xnew)))),
                                 nrow = S, ncol = nrow(Xnew))
    return(y_draws)
  }
  mu_draws
}

summarise_draws <- function(draw_mat) {
  data.frame(
    mean = colMeans(draw_mat),
    lwr = apply(draw_mat, 2, quantile, 0.025),
    upr = apply(draw_mat, 2, quantile, 0.975)
  )
}

coverage95 <- function(y, lwr, upr) mean(y >= lwr & y <= upr, na.rm = TRUE)

4 Simulation: when ML predicts better, but Bayes tells uncertainty

We simulate a small nonlinear relationship:

\[ y = \sin(2x)+0.5x+\epsilon. \]

A simple Bayesian linear regression is deliberately misspecified here because the true pattern is nonlinear. A decision tree can adapt to nonlinearity, but does not automatically give a calibrated posterior uncertainty interval. Then we combine them by using the tree prediction as an input to a Bayesian calibration model.

set.seed(1)
n <- 80
sim <- data.frame(
  x = runif(n, -2.5, 2.5)
)
sim$truth <- sin(2 * sim$x) + 0.5 * sim$x
sim$y <- sim$truth + rnorm(n, sd = 0.45)

idx <- train_test_split(nrow(sim), 0.7)
train <- sim[idx, ]
test  <- sim[-idx, ]

p <- ggplot(sim, aes(x, y)) +
  geom_point(alpha = 0.75) +
  stat_function(fun = function(z) sin(2*z) + 0.5*z, linewidth = 1) +
  labs(title = "Simulated nonlinear natural-science style signal",
       subtitle = "Black curve = true mean function",
       x = "x", y = "response") +
  theme_minimal()
p

4.1 Model A: Bayesian linear regression

Xtr <- model.matrix(~ x, data = train)
Xte <- model.matrix(~ x, data = test)
fit_bayes_lin <- bayes_lm_nig(Xtr, train$y)
draw_lin <- draw_bayes_lm(fit_bayes_lin, Xte, S = 3000, include_noise = TRUE)
pred_lin <- summarise_draws(draw_lin)

rmse_lin <- rmse(test$y, pred_lin$mean)
cov_lin <- coverage95(test$y, pred_lin$lwr, pred_lin$upr)

c(RMSE = rmse_lin, Coverage95 = cov_lin)
##       RMSE Coverage95 
##  0.9092169  0.9583333

4.2 Model B: ML decision tree

tree_sim <- rpart::rpart(y ~ x, data = train, method = "anova",
                         control = rpart.control(cp = 0.005, minsplit = 8))
pred_tree <- predict(tree_sim, newdata = test)
rmse_tree <- rmse(test$y, pred_tree)
rmse_tree
## [1] 0.4877032

4.3 Model C: ML + Bayesian calibration

Here the ML tree discovers nonlinear structure. We then fit a Bayesian calibration layer:

\[ y_i = \alpha + \gamma\,\widehat f_{ML}(x_i)+e_i. \]

This is a simple version of Bayesian post-processing of an ML model.

train$tree_hat <- predict(tree_sim, newdata = train)
test$tree_hat  <- predict(tree_sim, newdata = test)

Xtr_h <- model.matrix(~ tree_hat, data = train)
Xte_h <- model.matrix(~ tree_hat, data = test)
fit_hybrid <- bayes_lm_nig(Xtr_h, train$y)
draw_hybrid <- draw_bayes_lm(fit_hybrid, Xte_h, S = 3000, include_noise = TRUE)
pred_hybrid <- summarise_draws(draw_hybrid)

rmse_hybrid <- rmse(test$y, pred_hybrid$mean)
cov_hybrid <- coverage95(test$y, pred_hybrid$lwr, pred_hybrid$upr)

data.frame(
  Model = c("Bayesian linear", "ML tree", "ML tree + Bayesian calibration"),
  RMSE = c(rmse_lin, rmse_tree, rmse_hybrid),
  Coverage95 = c(cov_lin, NA, cov_hybrid)
)
##                            Model      RMSE Coverage95
## 1                Bayesian linear 0.9092169  0.9583333
## 2                        ML tree 0.4877032         NA
## 3 ML tree + Bayesian calibration 0.4859872  0.8750000
plot_df <- test %>%
  mutate(
    bayes_mean = pred_lin$mean,
    bayes_lwr = pred_lin$lwr,
    bayes_upr = pred_lin$upr,
    tree_pred = pred_tree,
    hybrid_mean = pred_hybrid$mean,
    hybrid_lwr = pred_hybrid$lwr,
    hybrid_upr = pred_hybrid$upr
  ) %>%
  arrange(x)

ggplot(plot_df, aes(x, y)) +
  geom_point(size = 2) +
  geom_line(aes(y = bayes_mean), linewidth = 1, linetype = 2) +
  geom_ribbon(aes(ymin = bayes_lwr, ymax = bayes_upr), alpha = 0.15) +
  geom_line(aes(y = hybrid_mean), linewidth = 1) +
  geom_ribbon(aes(ymin = hybrid_lwr, ymax = hybrid_upr), alpha = 0.15) +
  labs(title = "Bayesian linear vs ML+Bayesian calibration",
       subtitle = "Dashed = Bayesian linear; solid = Bayesian calibration of tree predictions",
       x = "x", y = "y") +
  theme_minimal()

Lesson: ML captures flexible pattern; Bayes reports uncertainty and calibrates predictions. Together, we can often obtain both adaptability and uncertainty.

5 Real data 1: Environmental science using airquality

The airquality dataset contains daily air-quality measurements in New York, May–September 1973. We model ozone using temperature, wind, and solar radiation.

data("airquality")
air <- airquality %>%
  dplyr::select(Ozone, Solar.R, Wind, Temp, Month) %>%
  tidyr::drop_na()

summary(air)
##      Ozone          Solar.R           Wind            Temp      
##  Min.   :  1.0   Min.   :  7.0   Min.   : 2.30   Min.   :57.00  
##  1st Qu.: 18.0   1st Qu.:113.5   1st Qu.: 7.40   1st Qu.:71.00  
##  Median : 31.0   Median :207.0   Median : 9.70   Median :79.00  
##  Mean   : 42.1   Mean   :184.8   Mean   : 9.94   Mean   :77.79  
##  3rd Qu.: 62.0   3rd Qu.:255.5   3rd Qu.:11.50   3rd Qu.:84.50  
##  Max.   :168.0   Max.   :334.0   Max.   :20.70   Max.   :97.00  
##      Month      
##  Min.   :5.000  
##  1st Qu.:6.000  
##  Median :7.000  
##  Mean   :7.216  
##  3rd Qu.:9.000  
##  Max.   :9.000
ggplot(air, aes(Temp, Ozone, colour = Wind)) +
  geom_point(size = 2) +
  scale_colour_viridis_c() +
  labs(title = "Environmental data: ozone, temperature and wind",
       x = "Temperature (F)", y = "Ozone (ppb)") +
  theme_minimal()

5.1 Bayesian regression for ozone

set.seed(2)
idx <- train_test_split(nrow(air), 0.72)
train <- air[idx, ]
test  <- air[-idx, ]

# Standardise predictors using training data only.
vars <- c("Solar.R", "Wind", "Temp")
mu <- sapply(train[vars], mean)
sdval <- sapply(train[vars], sd)
scale_df <- function(dat) {
  out <- dat
  out[vars] <- sweep(sweep(out[vars], 2, mu, "-"), 2, sdval, "/")
  out
}
train_s <- scale_df(train)
test_s <- scale_df(test)

Xtr <- model.matrix(~ Solar.R + Wind + Temp, data = train_s)
Xte <- model.matrix(~ Solar.R + Wind + Temp, data = test_s)
fit_air_bayes <- bayes_lm_nig(Xtr, train_s$Ozone)
draw_air <- draw_bayes_lm(fit_air_bayes, Xte, S = 4000, include_noise = TRUE)
pred_air_bayes <- summarise_draws(draw_air)

rmse_air_bayes <- rmse(test_s$Ozone, pred_air_bayes$mean)
cov_air_bayes <- coverage95(test_s$Ozone, pred_air_bayes$lwr, pred_air_bayes$upr)

# Posterior coefficient summaries
coef_draws <- draw_bayes_lm(fit_air_bayes, Xnew = diag(ncol(Xtr)), S = 1, include_noise = FALSE)
# Better: draw beta directly for coefficient summaries
S <- 5000
sigma2 <- 1 / rgamma(S, fit_air_bayes$an, rate = fit_air_bayes$bn)
beta_draws <- matrix(NA_real_, S, ncol(Xtr))
for (s in seq_len(S)) {
  beta_draws[s, ] <- MASS::mvrnorm(1, fit_air_bayes$beta_mean, sigma2[s] * fit_air_bayes$Vn)
}
coef_tab <- data.frame(
  Term = colnames(Xtr),
  Median = apply(beta_draws, 2, median),
  L95 = apply(beta_draws, 2, quantile, 0.025),
  U95 = apply(beta_draws, 2, quantile, 0.975),
  Prob_Positive = colMeans(beta_draws > 0)
)
coef_tab
##          Term    Median         L95       U95 Prob_Positive
## 1 (Intercept) 39.864623  36.2417994 43.547716        1.0000
## 2     Solar.R  3.535845  -0.3419495  7.211755        0.9622
## 3        Wind -9.672489 -14.1211400 -5.263183        0.0000
## 4        Temp 17.718192  13.0640154 22.225853        1.0000

5.2 ML tree for ozone

tree_air <- rpart::rpart(Ozone ~ Solar.R + Wind + Temp, data = train_s,
                         method = "anova",
                         control = rpart.control(cp = 0.01, minsplit = 10))
pred_air_tree <- predict(tree_air, newdata = test_s)
rmse_air_tree <- rmse(test_s$Ozone, pred_air_tree)

5.3 Hybrid: Bayesian calibration of the ML tree

train_s$tree_hat <- predict(tree_air, newdata = train_s)
test_s$tree_hat <- predict(tree_air, newdata = test_s)

Xtr_h <- model.matrix(~ tree_hat, data = train_s)
Xte_h <- model.matrix(~ tree_hat, data = test_s)
fit_air_hybrid <- bayes_lm_nig(Xtr_h, train_s$Ozone)
draw_air_hybrid <- draw_bayes_lm(fit_air_hybrid, Xte_h, S = 4000, include_noise = TRUE)
pred_air_hybrid <- summarise_draws(draw_air_hybrid)

rmse_air_hybrid <- rmse(test_s$Ozone, pred_air_hybrid$mean)
cov_air_hybrid <- coverage95(test_s$Ozone, pred_air_hybrid$lwr, pred_air_hybrid$upr)

score_air <- data.frame(
  Model = c("Bayesian linear", "ML tree", "ML tree + Bayesian calibration"),
  RMSE = c(rmse_air_bayes, rmse_air_tree, rmse_air_hybrid),
  Coverage95 = c(cov_air_bayes, NA, cov_air_hybrid)
)
score_air
##                            Model     RMSE Coverage95
## 1                Bayesian linear 28.71827     0.8125
## 2                        ML tree 31.44423         NA
## 3 ML tree + Bayesian calibration 31.45280     0.6250
pred_df <- test_s %>%
  mutate(
    Bayesian = pred_air_bayes$mean,
    Hybrid = pred_air_hybrid$mean,
    Hybrid_L = pred_air_hybrid$lwr,
    Hybrid_U = pred_air_hybrid$upr
  ) %>%
  arrange(Ozone)

ggplot(pred_df, aes(Ozone, Hybrid)) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = Hybrid_L, ymax = Hybrid_U), alpha = 0.25) +
  geom_abline(slope = 1, intercept = 0, linetype = 2) +
  labs(title = "Ozone prediction: ML tree with Bayesian uncertainty calibration",
       x = "Observed ozone", y = "Predicted ozone") +
  theme_minimal()

Scientific interpretation: a tree may capture nonlinear environmental effects, while Bayesian calibration gives a predictive interval. For policy or monitoring, an interval is often more useful than a single point forecast.

6 Real data 2: Botany and morphology using iris

The iris data contain measurements for 150 flowers from three species. To make the problem nontrivial, we classify only versicolor vs virginica.

data("iris")
iris2 <- iris %>%
  filter(Species != "setosa") %>%
  droplevels()

summary(iris2)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.900   Min.   :2.000   Min.   :3.000   Min.   :1.000  
##  1st Qu.:5.800   1st Qu.:2.700   1st Qu.:4.375   1st Qu.:1.300  
##  Median :6.300   Median :2.900   Median :4.900   Median :1.600  
##  Mean   :6.262   Mean   :2.872   Mean   :4.906   Mean   :1.676  
##  3rd Qu.:6.700   3rd Qu.:3.025   3rd Qu.:5.525   3rd Qu.:2.000  
##  Max.   :7.900   Max.   :3.800   Max.   :6.900   Max.   :2.500  
##        Species  
##  versicolor:50  
##  virginica :50  
##                 
##                 
##                 
##

6.1 ML classification tree

set.seed(3)
idx <- train_test_split(nrow(iris2), 0.7)
train <- iris2[idx, ]
test <- iris2[-idx, ]

tree_iris <- rpart::rpart(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
                          data = train,
                          method = "class",
                          control = rpart.control(cp = 0.01, minsplit = 8))
prob_tree <- predict(tree_iris, newdata = test, type = "prob")
pred_tree <- colnames(prob_tree)[max.col(prob_tree)]
acc_tree <- mean(pred_tree == test$Species)
acc_tree
## [1] 0.9333333

6.2 Bayesian Gaussian naive Bayes from scratch

This is a small Bayesian classifier. For each class and each feature, we use a Normal–Inverse-Gamma model, yielding a Student-\(t\) posterior predictive density. Class probabilities are combined using Bayes’ rule.

fit_one_feature <- function(x, m0 = 0, k0 = 0.01, a0 = 2, b0 = 1) {
  x <- as.numeric(x)
  n <- length(x)
  xbar <- mean(x)
  ss <- sum((x - xbar)^2)
  kn <- k0 + n
  mn <- (k0 * m0 + n * xbar) / kn
  an <- a0 + n / 2
  bn <- b0 + 0.5 * ss + (k0 * n * (xbar - m0)^2) / (2 * kn)
  list(mn = mn, kn = kn, an = an, bn = bn)
}

log_pred_t <- function(xnew, fit) {
  df <- 2 * fit$an
  scale <- sqrt(fit$bn * (fit$kn + 1) / (fit$an * fit$kn))
  dt((xnew - fit$mn) / scale, df = df, log = TRUE) - log(scale)
}

bayes_nb_fit <- function(train, yvar, xvars) {
  y <- train[[yvar]]
  classes <- levels(y)
  alpha <- rep(1, length(classes))
  names(alpha) <- classes
  class_counts <- table(y)
  class_prior <- (as.numeric(class_counts[classes]) + alpha) /
    (sum(class_counts) + sum(alpha))
  names(class_prior) <- classes
  fits <- list()
  for (cl in classes) {
    fits[[cl]] <- lapply(xvars, function(v) fit_one_feature(train[train[[yvar]] == cl, v]))
    names(fits[[cl]]) <- xvars
  }
  list(classes = classes, prior = class_prior, fits = fits, xvars = xvars)
}

bayes_nb_predict <- function(object, newdata) {
  logp <- matrix(NA_real_, nrow = nrow(newdata), ncol = length(object$classes))
  colnames(logp) <- object$classes
  for (cl in object$classes) {
    lp <- rep(log(object$prior[cl]), nrow(newdata))
    for (v in object$xvars) {
      lp <- lp + log_pred_t(newdata[[v]], object$fits[[cl]][[v]])
    }
    logp[, cl] <- lp
  }
  # stable softmax
  logp_centered <- logp - apply(logp, 1, max)
  prob <- exp(logp_centered) / rowSums(exp(logp_centered))
  as.data.frame(prob)
}

xvars <- c("Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width")
bnb <- bayes_nb_fit(train, "Species", xvars)
prob_bayes <- bayes_nb_predict(bnb, test)
pred_bayes <- colnames(prob_bayes)[max.col(prob_bayes)]
acc_bayes <- mean(pred_bayes == test$Species)

# Hybrid ensemble: average probabilities from ML tree and Bayesian classifier
common_cols <- intersect(colnames(prob_tree), colnames(prob_bayes))
prob_hybrid <- 0.5 * prob_tree[, common_cols] + 0.5 * as.matrix(prob_bayes[, common_cols])
pred_hybrid <- common_cols[max.col(prob_hybrid)]
acc_hybrid <- mean(pred_hybrid == test$Species)

data.frame(
  Model = c("ML decision tree", "Bayesian Gaussian classifier", "Average ML + Bayesian probabilities"),
  Accuracy = c(acc_tree, acc_bayes, acc_hybrid),
  Mean_Max_Probability = c(mean(apply(prob_tree, 1, max)),
                           mean(apply(prob_bayes, 1, max)),
                           mean(apply(prob_hybrid, 1, max)))
)
##                                 Model  Accuracy Mean_Max_Probability
## 1                    ML decision tree 0.9333333            0.9916667
## 2        Bayesian Gaussian classifier 0.9000000            0.9166456
## 3 Average ML + Bayesian probabilities 0.9333333            0.9526973
iris_pred <- test %>%
  mutate(
    Bayes_MaxProb = apply(prob_bayes, 1, max),
    Bayes_Pred = pred_bayes,
    Correct = Bayes_Pred == Species
  )

ggplot(iris_pred, aes(Petal.Length, Petal.Width, colour = Bayes_MaxProb, shape = Correct)) +
  geom_point(size = 3) +
  scale_colour_viridis_c(limits = c(0.5, 1)) +
  labs(title = "Bayesian classifier uncertainty on iris species",
       subtitle = "Lower maximum posterior class probability indicates uncertainty",
       x = "Petal length", y = "Petal width", colour = "Max posterior\nclass probability") +
  theme_minimal()

Lesson: ML gives a decision boundary; Bayesian classification gives probability and uncertainty. The hybrid probability average can be a simple way to stabilise predictions.

7 Real data 3: Zoology and allometry using MASS::Animals

Allometry asks how biological quantities scale with body size. Here we study the relationship between body weight and brain weight in animals. With only a small number of species, a Bayesian regression is valuable because uncertainty in the allometric exponent matters.

data("Animals", package = "MASS")
animals <- MASS::Animals %>%
  tibble::rownames_to_column("Animal") %>%
  mutate(log_body = log(body), log_brain = log(brain))

head(animals)
##            Animal     body brain   log_body log_brain
## 1 Mountain beaver     1.35   8.1 0.30010459  2.091864
## 2             Cow   465.00 423.0 6.14203741  6.047372
## 3       Grey wolf    36.33 119.5 3.59264385  4.783316
## 4            Goat    27.66 115.0 3.31998733  4.744932
## 5      Guinea pig     1.04   5.5 0.03922071  1.704748
## 6     Dipliodocus 11700.00  50.0 9.36734412  3.912023
ggplot(animals, aes(log_body, log_brain, label = Animal)) +
  geom_point(size = 2) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Zoological allometry: brain weight vs body weight",
       x = "log body weight", y = "log brain weight") +
  theme_minimal()

X <- model.matrix(~ log_body, data = animals)
fit_animals <- bayes_lm_nig(X, animals$log_brain, a0 = 2, b0 = 1)

S <- 6000
sigma2 <- 1 / rgamma(S, fit_animals$an, rate = fit_animals$bn)
beta_draws <- matrix(NA_real_, S, ncol(X))
for (s in seq_len(S)) {
  beta_draws[s, ] <- MASS::mvrnorm(1, fit_animals$beta_mean, sigma2[s] * fit_animals$Vn)
}
colnames(beta_draws) <- colnames(X)

allometry_tab <- data.frame(
  Parameter = colnames(beta_draws),
  Median = apply(beta_draws, 2, median),
  L95 = apply(beta_draws, 2, quantile, 0.025),
  U95 = apply(beta_draws, 2, quantile, 0.975),
  Prob_Positive = colMeans(beta_draws > 0)
)
allometry_tab
##               Parameter   Median       L95       U95 Prob_Positive
## (Intercept) (Intercept) 2.557552 1.7844757 3.3322046             1
## log_body       log_body 0.497074 0.3482656 0.6420763             1
beta_df <- data.frame(slope = beta_draws[, "log_body"])
ggplot(beta_df, aes(slope)) +
  geom_density(fill = "steelblue", alpha = 0.25) +
  geom_vline(xintercept = median(beta_df$slope), linetype = 2) +
  labs(title = "Posterior uncertainty in the allometric exponent",
       x = "Slope of log(brain) on log(body)", y = "Posterior density") +
  theme_minimal()

Lesson: when the scientific target is an interpretable exponent, a Bayesian posterior is more informative than only a fitted line.

8 Real data 4: Biomedical survival using survival::lung

Survival data are different from ordinary regression data because some patients are censored. The event time is not always observed exactly. We use a simple exponential survival model with censoring to demonstrate Bayesian small-data inference.

For an exponential event time with rate \(\lambda\), censored survival likelihood contributes

\[ L(\lambda) \propto \lambda^d \exp(-\lambda T), \]

where \(d\) is the number of observed events and \(T\) is total observed follow-up time. With prior

\[ \lambda \sim \mathrm{Gamma}(a_0,b_0), \]

the posterior is

\[ \lambda\mid data \sim \mathrm{Gamma}(a_0+d, b_0+T). \]

data("lung", package = "survival")
lung2 <- survival::lung %>%
  mutate(
    event = ifelse(status == 2, 1, 0),
    sex_label = ifelse(sex == 1, "Male", "Female")
  ) %>%
  filter(!is.na(time), !is.na(event), !is.na(sex_label))

table(lung2$sex_label, lung2$event)
##         
##            0   1
##   Female  37  53
##   Male    26 112
km <- survival::survfit(survival::Surv(time, event) ~ sex_label, data = lung2)
plot(km, col = c("firebrick", "navy"), lwd = 2,
     xlab = "Days", ylab = "Survival probability",
     main = "Kaplan--Meier curves by sex")
legend("topright", legend = names(km$strata), col = c("firebrick", "navy"), lwd = 2, bty = "n")

bayes_exp_rate <- function(time, event, a0 = 0.01, b0 = 0.01, S = 10000) {
  d <- sum(event == 1)
  Ttot <- sum(time)
  shape <- a0 + d
  rate <- b0 + Ttot
  lambda <- rgamma(S, shape = shape, rate = rate)
  med_surv <- log(2) / lambda
  data.frame(lambda = lambda, median_survival = med_surv)
}

post_male <- bayes_exp_rate(lung2$time[lung2$sex_label == "Male"],
                            lung2$event[lung2$sex_label == "Male"])
post_female <- bayes_exp_rate(lung2$time[lung2$sex_label == "Female"],
                              lung2$event[lung2$sex_label == "Female"])

post_surv <- bind_rows(
  post_male %>% mutate(Group = "Male"),
  post_female %>% mutate(Group = "Female")
)

surv_summary <- post_surv %>%
  group_by(Group) %>%
  summarise(
    median_rate = median(lambda),
    L95_rate = quantile(lambda, 0.025),
    U95_rate = quantile(lambda, 0.975),
    median_survival_days = median(median_survival),
    L95_median_survival = quantile(median_survival, 0.025),
    U95_median_survival = quantile(median_survival, 0.975),
    .groups = "drop"
  )
surv_summary
## # A tibble: 2 × 7
##   Group  median_rate L95_rate U95_rate median_survival_days L95_median_survival
##   <chr>        <dbl>    <dbl>    <dbl>                <dbl>               <dbl>
## 1 Female     0.00173  0.00129  0.00223                 401.                311.
## 2 Male       0.00285  0.00236  0.00341                 243.                203.
## # ℹ 1 more variable: U95_median_survival <dbl>
ggplot(post_surv, aes(median_survival, fill = Group)) +
  geom_density(alpha = 0.35) +
  labs(title = "Bayesian posterior for median survival under exponential model",
       x = "Median survival in days", y = "Posterior density") +
  theme_minimal()

# A direct Bayesian probability question:
# What is the posterior probability that median survival is larger for females than males?
prob_female_better <- mean(post_female$median_survival > post_male$median_survival)
prob_female_better
## [1] 0.9988

Lesson: Bayesian survival inference lets us answer direct probability questions such as:
“What is the posterior probability that one group has larger median survival than another?”

9 What did we learn?

Question ML/AI answer Bayesian answer Combined answer
Is prediction possible? Often yes, using flexible learners. Yes, especially with probabilistic models. Use ML predictors inside Bayesian calibration.
How uncertain is the answer? Not always automatic. Natural through posterior/predictive intervals. Bayesian post-processing of ML outputs.
What happens in small data? Risk of overfitting. Priors regularise and stabilise. ML features + Bayesian shrinkage.
Is the model interpretable? Depends on the learner. Often clearer through parameters. ML for pattern, Bayes for explanation.
Natural-science decision? Point prediction. Uncertainty-aware decision. Predictive + calibrated + interpretable.

10 A small concluding message

Bayesian inference and AI/ML should not be presented as enemies. In natural-science data analysis:

  • ML can discover complex predictive structure.
  • Bayesian inference can quantify uncertainty, encode scientific prior information, and handle small-data regimes.
  • Hybrid workflows can take an ML output and place a Bayesian uncertainty layer on top.

For scientific students, the key habit is:

Do not ask only: what is the prediction? Also ask: how uncertain is it, why should I believe it, and what scientific decision depends on it?

11 Reproducibility and package citations

citation()
citation("MASS")
citation("survival")
citation("rpart")
citation("ggplot2")

Useful references:

  1. Gelman, A. et al. (2013). Bayesian Data Analysis, 3rd ed. Chapman and Hall/CRC.
  2. Berger, J. O. (1985). Statistical Decision Theory and Bayesian Analysis, 2nd ed. Springer.
  3. Breiman, L. (2001). Random Forests. Machine Learning, 45, 5–32.
  4. Cox, D. R. (1972). Regression Models and Life-Tables. JRSS Series B, 34(2), 187–220.
  5. Anderson, E. (1935). The Irises of the Gaspe Peninsula. Bulletin of the American Iris Society, 59, 2–5.
  6. Fisher, R. A. (1936). The use of multiple measurements in taxonomic problems. Annals of Eugenics, 7(2), 179–188.
  7. Chambers, J. M., Cleveland, W. S., Kleiner, B., and Tukey, P. A. (1983). Graphical Methods for Data Analysis. Wadsworth & Brooks/Cole.
  8. Venables, W. N. and Ripley, B. D. (2002). Modern Applied Statistics with S, 4th ed. Springer.
  9. Therneau, T. M. (survival package documentation and vignettes).

12 How to publish this on RPubs

  1. Save this file as Bayes_ML_Natural_Science_RPubs.Rmd.
  2. Open it in RStudio.
  3. Click Knit to produce HTML.
  4. Click Publish and choose RPubs.
  5. Add a short title such as: Bayesian Inference Meets AI/ML in Natural-Science Data.