1 Purpose of this RPubs note

This note is designed as a foundation module before Bayesian inference. The aim is to motivate researchers in basic sciences—physics, astronomy, geoscience, biology, materials science, and environmental science—to see statistical inference and machine learning as complementary components of data science.

We assume basic probability, but we carefully explain the statistical notation, uncertainty language, stochastic-process viewpoint, and filtering idea. Bayesian inference is intentionally kept mostly outside this note; it can come as the next lecture. Here the focus is:

The real-data examples use small built-in R datasets so that students can run everything immediately.

2 Data science in basic science: the four-layer view

In many scientific problems, data science is not merely applying a software package. It has four conceptual layers.

Layer Scientific meaning Statistical/ML object
Data layer What was measured? Under what instrument, design, and noise? random variables, sampling scheme, missingness
Model layer What scientific structure is being represented? regression, classification, time series, stochastic process
Inference layer What can be learned, and with what uncertainty? estimators, confidence intervals, tests, residuals
Prediction/decision layer What will happen next, or what should be done? prediction error, cross-validation, risk, control/filtering

A scientific data analysis fails when these layers are confused. For example, a model may have high predictive accuracy but answer the wrong scientific question; or a beautiful mechanistic model may have invalid uncertainty because its assumptions are wrong.

3 Core notation for statistical inference

Let \[ D_n = \{(X_i,Y_i): i=1,\ldots,n\} \] be observed data. Here:

A statistical model usually says that the data distribution belongs to a family \[ P_\theta, \qquad \theta \in \Theta, \] where \(\theta\) is an unknown parameter or collection of parameters.

An estimand is the scientific target, for example:

An estimator is a data-dependent rule: \[ \widehat\theta = T(D_n). \]

A basic frequentist uncertainty statement studies the distribution of \(\widehat\theta\) under repeated sampling: \[ \widehat\theta - \theta. \]

Three important quantities are: \[ \operatorname{Bias}(\widehat\theta)=E(\widehat\theta)-\theta, \] \[ \operatorname{Var}(\widehat\theta)=E\{\widehat\theta-E(\widehat\theta)\}^2, \] \[ \operatorname{MSE}(\widehat\theta)=E(\widehat\theta-\theta)^2 =\operatorname{Bias}^2+\operatorname{Var}. \]

4 Statistics and ML: two philosophies inside data science

Statistics and machine learning overlap heavily, but their emphasis is different.

Question Statistical inference emphasis Machine learning emphasis
What is the target? estimand, parameter, scientific effect prediction function, classifier, score
What is uncertainty? sampling uncertainty, standard errors, confidence intervals generalization error, validation error, calibration
What is a good model? interpretable, identifiable, well-checked predictive, flexible, scalable
Main danger wrong assumptions, underfitting, misinterpreting p-values overfitting, shortcut learning, extrapolation failure
Scientific strength effect estimation and uncertainty pattern recognition and flexible prediction
Best use together estimand + uncertainty + prediction + validation ML learns structure; statistics audits uncertainty

A useful working principle is:

Statistics asks: what can we infer, with what uncertainty, under what assumptions? ML asks: what can we predict, and how well does it generalize? Data science needs both.

5 Simulation 1: inference, prediction, and uncertainty

We simulate a nonlinear scientific relationship. A linear model is interpretable but misspecified; a polynomial model is more flexible; a tree is a simple ML model. This illustrates that prediction and inference are related but not identical.

set.seed(101)
n <- 90
sim <- data.frame(x = runif(n, -3, 3))
sim$f_true <- with(sim, sin(1.4 * x) + 0.25 * x)
sim$y <- sim$f_true + rnorm(n, sd = 0.35)

train_id <- sample(seq_len(n), size = 65)
train <- sim[train_id, ]
test  <- sim[-train_id, ]

fit_lm   <- lm(y ~ x, data = train)
fit_poly <- lm(y ~ poly(x, 5, raw = TRUE), data = train)
fit_tree <- rpart::rpart(y ~ x, data = train, control = rpart::rpart.control(cp = 0.01))

model_score <- data.frame(
  model = c("Linear statistical model", "Polynomial regression", "ML regression tree"),
  test_RMSE = c(
    rmse(test$y, predict(fit_lm, newdata = test)),
    rmse(test$y, predict(fit_poly, newdata = test)),
    rmse(test$y, predict(fit_tree, newdata = test))
  ),
  test_MAE = c(
    mae(test$y, predict(fit_lm, newdata = test)),
    mae(test$y, predict(fit_poly, newdata = test)),
    mae(test$y, predict(fit_tree, newdata = test))
  )
)

knitr::kable(model_score, digits = 3, caption = "Prediction performance on a held-out test set.")
Prediction performance on a held-out test set.
model test_RMSE test_MAE
Linear statistical model 0.736 0.653
Polynomial regression 0.346 0.261
ML regression tree 0.368 0.280 
grid <- data.frame(x = seq(min(sim$x), max(sim$x), length.out = 300))
grid$f_true <- with(grid, sin(1.4 * x) + 0.25 * x)

lm_ci <- as.data.frame(predict(fit_lm, newdata = grid, interval = "confidence"))
poly_pred <- predict(fit_poly, newdata = grid)
tree_pred <- predict(fit_tree, newdata = grid)

grid$lm_fit <- lm_ci$fit
grid$lm_lwr <- lm_ci$lwr
grid$lm_upr <- lm_ci$upr
grid$poly <- poly_pred
grid$tree <- tree_pred

ggplot2::ggplot(sim, ggplot2::aes(x, y)) +
  ggplot2::geom_point(alpha = 0.65) +
  ggplot2::geom_line(data = grid, ggplot2::aes(x, f_true), linewidth = 1.1, linetype = 2) +
  ggplot2::geom_ribbon(data = grid, ggplot2::aes(x = x, ymin = lm_lwr, ymax = lm_upr), alpha = 0.15, inherit.aes = FALSE) +
  ggplot2::geom_line(data = grid, ggplot2::aes(x, lm_fit), linewidth = 1) +
  ggplot2::geom_line(data = grid, ggplot2::aes(x, poly), linewidth = 1, linetype = 3) +
  ggplot2::geom_line(data = grid, ggplot2::aes(x, tree), linewidth = 1, linetype = 4) +
  ggplot2::labs(
    title = "Statistical model, flexible regression, and ML tree on a nonlinear scientific relationship",
    subtitle = "Dashed curve = true mean function; ribbon = linear-model confidence band",
    y = "response", x = "scientific covariate"
  )

Interpretation. The linear model gives a clean slope and confidence band, but it is misspecified. The tree can adapt locally but gives no natural uncertainty interval. The polynomial model is flexible but can become unstable if the degree is too high. The lesson is not that one family wins always: the analyst must match the model to the scientific purpose.

6 ML + inference: split conformal uncertainty for an ML predictor

Many ML methods produce point predictions. Researchers, however, often need uncertainty. A simple distribution-free idea is split conformal prediction. We train on one part of data, calibrate prediction residuals on a second part, and build an interval around future predictions.

set.seed(202)
idx <- sample(seq_len(nrow(sim)))
tr  <- sim[idx[1:45], ]
cal <- sim[idx[46:65], ]
te  <- sim[idx[66:nrow(sim)], ]

ml_tree <- rpart::rpart(y ~ x, data = tr, control = rpart::rpart.control(cp = 0.01))
cal_resid <- abs(cal$y - predict(ml_tree, newdata = cal))
alpha <- 0.10
q_level <- min(1, ceiling((length(cal_resid) + 1) * (1 - alpha)) / length(cal_resid))
qhat <- as.numeric(quantile(cal_resid, probs = q_level, type = 1))

te$pred <- predict(ml_tree, newdata = te)
te$lwr <- te$pred - qhat
te$upr <- te$pred + qhat
coverage <- mean(te$y >= te$lwr & te$y <= te$upr)

knitr::kable(
  data.frame(n_train = nrow(tr), n_calibration = nrow(cal), n_test = nrow(te), qhat = qhat, empirical_coverage = coverage),
  digits = 3,
  caption = "Split conformal interval around a tree predictor."
)
Split conformal interval around a tree predictor.
n_train n_calibration n_test qhat empirical_coverage
45 20 25 0.806 0.88 
ggplot2::ggplot(te, ggplot2::aes(x, y)) +
  ggplot2::geom_point() +
  ggplot2::geom_errorbar(ggplot2::aes(ymin = lwr, ymax = upr), width = 0.03, alpha = 0.55) +
  ggplot2::geom_point(ggplot2::aes(y = pred), shape = 4, size = 3) +
  ggplot2::labs(
    title = "ML prediction plus statistical calibration",
    subtitle = "Crosses are ML predictions; vertical bars are conformal predictive intervals",
    x = "x", y = "observed/predicted response"
  )

This is a good example of hand-in-hand data science: ML learns a flexible predictor, while statistical calibration turns the prediction into an uncertainty-aware statement.

7 Counterexample: when ML fails by extrapolation, and statistics helps only if the model is scientifically meaningful

We simulate a process that grows approximately exponentially. We train only on a limited experimental range and test outside that range. A tree cannot extrapolate beyond the training range. A physically motivated log-linear model can extrapolate better, but only because the transformation matches the scientific mechanism.

set.seed(303)
tr_ex <- data.frame(x = runif(80, 0, 3))
tr_ex$truth <- exp(0.55 * tr_ex$x)
tr_ex$y <- tr_ex$truth + rnorm(nrow(tr_ex), sd = 0.20)
tr_ex$y <- pmax(tr_ex$y, 0.05)

te_ex <- data.frame(x = seq(3.05, 5, length.out = 80))
te_ex$truth <- exp(0.55 * te_ex$x)
te_ex$y <- te_ex$truth + rnorm(nrow(te_ex), sd = 0.20)
te_ex$y <- pmax(te_ex$y, 0.05)

fit_naive_lm <- lm(y ~ x, data = tr_ex)
fit_log_lm   <- lm(log(y) ~ x, data = tr_ex)
fit_ex_tree  <- rpart::rpart(y ~ x, data = tr_ex)

te_ex$pred_linear <- predict(fit_naive_lm, newdata = te_ex)
te_ex$pred_loglinear <- exp(predict(fit_log_lm, newdata = te_ex))
te_ex$pred_tree <- predict(fit_ex_tree, newdata = te_ex)

ex_score <- data.frame(
  model = c("Naive linear", "Physics-inspired log-linear", "ML tree"),
  extrapolation_RMSE = c(
    rmse(te_ex$y, te_ex$pred_linear),
    rmse(te_ex$y, te_ex$pred_loglinear),
    rmse(te_ex$y, te_ex$pred_tree)
  )
)
knitr::kable(ex_score, digits = 3, caption = "Extrapolation error outside the training range.")
Extrapolation error outside the training range.
model extrapolation_RMSE
Naive linear 4.204
Physics-inspired log-linear 0.392
ML tree 6.036 
grid_ex <- rbind(
  dplyr::mutate(tr_ex, region = "training"),
  dplyr::mutate(te_ex[, c("x", "truth", "y")], region = "extrapolation")
)

ggplot2::ggplot() +
  ggplot2::geom_point(data = grid_ex, ggplot2::aes(x, y, shape = region), alpha = 0.55) +
  ggplot2::geom_line(data = te_ex, ggplot2::aes(x, truth), linewidth = 1.1, linetype = 2) +
  ggplot2::geom_line(data = te_ex, ggplot2::aes(x, pred_linear), linewidth = 1) +
  ggplot2::geom_line(data = te_ex, ggplot2::aes(x, pred_loglinear), linewidth = 1, linetype = 3) +
  ggplot2::geom_line(data = te_ex, ggplot2::aes(x, pred_tree), linewidth = 1, linetype = 4) +
  ggplot2::labs(
    title = "Extrapolation is a scientific assumption, not just an ML problem",
    subtitle = "The tree is flexible inside the observed range but cannot infer the exponential law outside it",
    x = "experimental condition", y = "response"
  )

Message. ML often excels at interpolation. Extrapolation requires structure. Statistical modelling can help when the structure is scientifically justified. But a wrong parametric structure can also fail. This is why scientific modelling is not reducible to software choice.

8 Real data 1: Old Faithful geyser as a geophysical process

The faithful dataset contains waiting times between eruptions and eruption durations for the Old Faithful geyser. It has 272 observations on two variables. The scientific question is: can eruption duration predict the waiting time to the next eruption?

data(faithful)
faith <- faithful

set.seed(404)
train_id <- sample(seq_len(nrow(faith)), size = round(0.75 * nrow(faith)))
faith_train <- faith[train_id, ]
faith_test  <- faith[-train_id, ]

faith_lm <- lm(waiting ~ eruptions, data = faith_train)
faith_lo <- loess(waiting ~ eruptions, data = faith_train, span = 0.75)
faith_tree <- rpart::rpart(waiting ~ eruptions, data = faith_train)

faith_score <- data.frame(
  model = c("Linear inference model", "LOESS smoother", "ML tree"),
  test_RMSE = c(
    rmse(faith_test$waiting, predict(faith_lm, newdata = faith_test)),
    rmse(faith_test$waiting, predict(faith_lo, newdata = faith_test)),
    rmse(faith_test$waiting, predict(faith_tree, newdata = faith_test))
  )
)
knitr::kable(faith_score, digits = 3, caption = "Prediction of waiting time for Old Faithful eruptions.")
Prediction of waiting time for Old Faithful eruptions.
model test_RMSE
Linear inference model 5.872
LOESS smoother 5.447
ML tree 5.542 
faith_grid <- data.frame(eruptions = seq(min(faith$eruptions), max(faith$eruptions), length.out = 250))
faith_ci <- as.data.frame(predict(faith_lm, newdata = faith_grid, interval = "confidence"))
faith_grid$lm_fit <- faith_ci$fit
faith_grid$lm_lwr <- faith_ci$lwr
faith_grid$lm_upr <- faith_ci$upr
faith_grid$loess <- predict(faith_lo, newdata = faith_grid)
faith_grid$tree <- predict(faith_tree, newdata = faith_grid)

ggplot2::ggplot(faith, ggplot2::aes(eruptions, waiting)) +
  ggplot2::geom_point(alpha = 0.55) +
  ggplot2::geom_ribbon(data = faith_grid, ggplot2::aes(x = eruptions, ymin = lm_lwr, ymax = lm_upr), inherit.aes = FALSE, alpha = 0.15) +
  ggplot2::geom_line(data = faith_grid, ggplot2::aes(eruptions, lm_fit), linewidth = 1) +
  ggplot2::geom_line(data = faith_grid, ggplot2::aes(eruptions, loess), linewidth = 1, linetype = 3) +
  ggplot2::geom_line(data = faith_grid, ggplot2::aes(eruptions, tree), linewidth = 1, linetype = 4) +
  ggplot2::labs(
    title = "Old Faithful geyser: statistical inference and ML-style flexible prediction",
    subtitle = "Longer eruptions are associated with longer waiting times, but the relationship is not purely linear",
    x = "eruption duration (minutes)", y = "waiting time to next eruption (minutes)"
  )

knitr::kable(
  data.frame(term = names(coef(faith_lm)), estimate = coef(faith_lm), confint(faith_lm)),
  digits = 3,
  caption = "Linear model estimates with 95% confidence intervals."
)
Linear model estimates with 95% confidence intervals.
term estimate X2.5.. X97.5..
(Intercept) (Intercept) 33.247 30.565 35.930
eruptions eruptions 10.794 10.071 11.518

Interpretation. The linear model gives an interpretable effect: a longer eruption is associated with a longer next waiting time. The smoother and tree capture nonlinearity. The statistical model gives uncertainty around the scientific effect; ML-style methods help discover flexible structure.

9 Stochastic processes: a scientific language for time, noise, and dynamics

A stochastic process is a family of random variables indexed by time or space: \[ \{X_t: t\in T\}. \] Examples:

9.1 Three simulated processes

set.seed(505)
Tn <- 250
process_df <- data.frame(t = 1:Tn)
process_df$random_walk <- cumsum(rnorm(Tn, sd = 1))
process_df$ar1 <- as.numeric(arima.sim(list(ar = 0.85), n = Tn, sd = 1))
process_df$poisson_counts <- rpois(Tn, lambda = 4 + 2 * sin(2 * pi * (1:Tn) / 50))

long_proc <- tidyr::pivot_longer(process_df, cols = -t, names_to = "process", values_to = "value")

ggplot2::ggplot(long_proc, ggplot2::aes(t, value)) +
  ggplot2::geom_line() +
  ggplot2::facet_wrap(~process, scales = "free_y", ncol = 1) +
  ggplot2::labs(
    title = "Three basic stochastic processes",
    subtitle = "Random walk, AR(1), and time-varying Poisson counts",
    x = "time", y = "process value"
  )

Scientific translation. A time series is not merely a list of numbers. It may represent a dynamic system with memory, forcing, feedback, shocks, and measurement noise.

10 Filtering: inference when the true state is hidden

Many scientific systems have an unobserved state \(X_t\) and noisy observations \(Y_t\): \[ X_t = aX_{t-1}+\eta_t, \qquad \eta_t\sim N(0,q), \] \[ Y_t = X_t+\epsilon_t, \qquad \epsilon_t\sim N(0,r). \]

Here:

The Kalman filter recursively updates the estimate of \(X_t\) after each new observation. It has two steps:

  1. Prediction: use the state equation to predict the next state.
  2. Update: use the observation to correct the prediction.
kalman_local_level <- function(y, q, r, m0 = y[1], C0 = 10 * r) {
  n <- length(y)
  m <- C <- a <- Rv <- K <- rep(NA_real_, n)
  m_prev <- m0
  C_prev <- C0
  for (t in seq_len(n)) {
    a[t] <- m_prev
    Rv[t] <- C_prev + q
    K[t] <- Rv[t] / (Rv[t] + r)
    m[t] <- a[t] + K[t] * (y[t] - a[t])
    C[t] <- (1 - K[t]) * Rv[t]
    m_prev <- m[t]
    C_prev <- C[t]
  }
  data.frame(time = seq_along(y), obs = y, pred = a, filtered = m,
             se = sqrt(C), lower = m - 1.96 * sqrt(C), upper = m + 1.96 * sqrt(C), K = K)
}

10.1 Simulation: hidden state plus noisy observations

set.seed(606)
Tn <- 120
x_true <- numeric(Tn)
x_true[1] <- 0
q_true <- 0.08
r_true <- 0.90
for (t in 2:Tn) x_true[t] <- x_true[t - 1] + rnorm(1, sd = sqrt(q_true))
y_obs <- x_true + rnorm(Tn, sd = sqrt(r_true))

kf_sim <- kalman_local_level(y_obs, q = q_true, r = r_true)
kf_sim$truth <- x_true

ggplot2::ggplot(kf_sim, ggplot2::aes(time)) +
  ggplot2::geom_point(ggplot2::aes(y = obs), alpha = 0.35) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = lower, ymax = upper), alpha = 0.15) +
  ggplot2::geom_line(ggplot2::aes(y = filtered), linewidth = 1) +
  ggplot2::geom_line(ggplot2::aes(y = truth), linewidth = 1, linetype = 2) +
  ggplot2::labs(
    title = "Kalman filtering: extracting a hidden state from noisy observations",
    subtitle = "Dashed line = true hidden state; solid line = filtered estimate; ribbon = uncertainty band",
    y = "state/observation", x = "time"
  )

11 Real data 2: Nile flow as hydrology and state-space inference

The Nile time series contains annual flow of the river Nile at Aswan from 1871 to 1970, length 100, with an apparent change point near 1898. We use it to show time-series inference, bad validation practice, and filtering.

data(Nile)
nile_df <- data.frame(year = as.integer(time(Nile)), flow = as.numeric(Nile))

ggplot2::ggplot(nile_df, ggplot2::aes(year, flow)) +
  ggplot2::geom_line() +
  ggplot2::geom_vline(xintercept = 1898, linetype = 2) +
  ggplot2::labs(
    title = "River Nile annual flow",
    subtitle = "A classical hydrological time series with apparent level shift near 1898",
    x = "year", y = expression(paste("flow (", 10^8, " m"^3, ")"))
  )

11.1 Simple change-point inference by comparing periods

nile_df$period <- ifelse(nile_df$year <= 1898, "1871--1898", "1899--1970")
period_summary <- nile_df |>
  dplyr::group_by(period) |>
  dplyr::summarise(n = dplyr::n(), mean_flow = mean(flow), sd_flow = sd(flow), .groups = "drop")
knitr::kable(period_summary, digits = 2, caption = "Nile flow before and after the apparent change point.")
Nile flow before and after the apparent change point.
period n mean_flow sd_flow
1871–1898 28 1097.75 135.00
1899–1970 72 849.97 124.78 
ttest_nile <- t.test(flow ~ period, data = nile_df)
knitr::kable(
  data.frame(
    estimate_difference = diff(ttest_nile$estimate),
    conf_low = ttest_nile$conf.int[1],
    conf_high = ttest_nile$conf.int[2],
    p_value = ttest_nile$p.value
  ),
  digits = 4,
  caption = "Two-sample comparison of mean flow across periods."
)
Two-sample comparison of mean flow across periods.
estimate_difference conf_low conf_high p_value
mean in group 1899–1970 -247.7778 188.5048 307.0508 0

11.2 Validation counterexample: random split is misleading for time series

nile_lag <- nile_df |>
  dplyr::mutate(lag1 = dplyr::lag(flow)) |>
  stats::na.omit()

set.seed(707)
rand_train <- sample(seq_len(nrow(nile_lag)), size = round(0.7 * nrow(nile_lag)))
train_random <- nile_lag[rand_train, ]
test_random  <- nile_lag[-rand_train, ]

time_cut <- sort(nile_lag$year)[round(0.7 * nrow(nile_lag))]
train_time <- nile_lag[nile_lag$year <= time_cut, ]
test_time  <- nile_lag[nile_lag$year > time_cut, ]

fit_random <- rpart::rpart(flow ~ year + lag1, data = train_random)
fit_time   <- rpart::rpart(flow ~ year + lag1, data = train_time)

validation_table <- data.frame(
  validation_design = c("Random split", "Chronological split"),
  test_RMSE = c(
    rmse(test_random$flow, predict(fit_random, newdata = test_random)),
    rmse(test_time$flow, predict(fit_time, newdata = test_time))
  ),
  test_MAE = c(
    mae(test_random$flow, predict(fit_random, newdata = test_random)),
    mae(test_time$flow, predict(fit_time, newdata = test_time))
  )
)
knitr::kable(validation_table, digits = 2, caption = "Random validation can be optimistic for time-dependent scientific data.")
Random validation can be optimistic for time-dependent scientific data.
validation_design test_RMSE test_MAE
Random split 132.22 98.38
Chronological split 133.44 103.52 
kf_nile <- kalman_local_level(as.numeric(Nile), q = 1500, r = 12000)
kf_nile$year <- as.integer(time(Nile))

ggplot2::ggplot(kf_nile, ggplot2::aes(year)) +
  ggplot2::geom_point(ggplot2::aes(y = obs), alpha = 0.45) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = lower, ymax = upper), alpha = 0.15) +
  ggplot2::geom_line(ggplot2::aes(y = filtered), linewidth = 1) +
  ggplot2::labs(
    title = "Local-level Kalman filter for Nile flow",
    subtitle = "A toy state-space smoother: hidden hydrological level plus noisy annual observations",
    x = "year", y = expression(paste("flow (", 10^8, " m"^3, ")"))
  )

Interpretation. Time-series data need time-aware validation. Filtering gives a dynamic inference tool: the scientific state is updated sequentially as new measurements arrive.

12 Real data 3: plant physiology with the CO2 dataset

The CO2 dataset concerns carbon dioxide uptake in grass plants under different concentrations, treatments, and plant types. It is useful for explaining regression, interaction, nonlinear dose-response, and ML prediction.

data(CO2)
co2_df <- as.data.frame(CO2)
co2_df$log_conc <- log(co2_df$conc)

co2_lm <- lm(uptake ~ log_conc * Treatment + Type, data = co2_df)

coef_table <- data.frame(
  term = names(coef(co2_lm)),
  estimate = coef(co2_lm),
  confint(co2_lm)
)
knitr::kable(coef_table, digits = 3, caption = "Regression effects for CO2 uptake.")
Regression effects for CO2 uptake.
term estimate X2.5.. X97.5..
(Intercept) (Intercept) -19.163 -30.929 -7.396
log_conc log_conc 9.646 7.649 11.644
Treatmentchilled Treatmentchilled 6.670 -9.902 23.243
TypeMississippi TypeMississippi -12.660 -14.780 -10.539
log_conc:Treatmentchilled log_conc:Treatmentchilled -2.325 -5.150 0.499 
ggplot2::ggplot(co2_df, ggplot2::aes(conc, uptake, color = Treatment, shape = Type)) +
  ggplot2::geom_point(size = 2) +
  ggplot2::scale_x_log10() +
  ggplot2::geom_smooth(method = "lm", formula = y ~ log(x), se = TRUE) +
  ggplot2::labs(
    title = "Plant physiology: CO2 uptake as a function of concentration",
    subtitle = "Statistical regression estimates treatment and type effects; ML can capture nonlinear patterns",
    x = "ambient CO2 concentration (log scale)", y = "CO2 uptake"
  )

set.seed(808)
tr_id <- sample(seq_len(nrow(co2_df)), size = round(0.75 * nrow(co2_df)))
co2_train <- co2_df[tr_id, ]
co2_test  <- co2_df[-tr_id, ]

co2_lm2 <- lm(uptake ~ log_conc * Treatment + Type, data = co2_train)
co2_tree <- rpart::rpart(uptake ~ conc + Treatment + Type, data = co2_train)

co2_scores <- data.frame(
  model = c("Interpretable statistical regression", "ML tree"),
  test_RMSE = c(
    rmse(co2_test$uptake, predict(co2_lm2, newdata = co2_test)),
    rmse(co2_test$uptake, predict(co2_tree, newdata = co2_test))
  )
)
knitr::kable(co2_scores, digits = 3, caption = "Predicting plant CO2 uptake.")
Predicting plant CO2 uptake.
model test_RMSE
Interpretable statistical regression 3.696
ML tree 6.027

Interpretation. Regression gives interpretable treatment effects and confidence intervals. The tree can capture simple nonlinear thresholds and interactions. A biological scientist often needs both: explanation and prediction.

13 Real data 4: topographic surface of Maunga Whau volcano

The volcano dataset gives topographic information for Maunga Whau on a 10m by 10m grid. It is a toy example for spatial data, surface modelling, residual maps, and geoscience visualization.

data(volcano)
volcano_df <- data.frame(
  row = rep(seq_len(nrow(volcano)), times = ncol(volcano)),
  col = rep(seq_len(ncol(volcano)), each = nrow(volcano)),
  elevation = as.vector(volcano)
)

ggplot2::ggplot(volcano_df, ggplot2::aes(col, row, fill = elevation)) +
  ggplot2::geom_raster() +
  ggplot2::coord_equal() +
  ggplot2::labs(
    title = "Maunga Whau volcano topography",
    subtitle = "A spatial scientific surface: every pixel has location and elevation",
    x = "south-north grid", y = "east-west grid", fill = "elevation"
  )

set.seed(909)
vid <- sample(seq_len(nrow(volcano_df)), size = round(0.7 * nrow(volcano_df)))
vol_train <- volcano_df[vid, ]
vol_test  <- volcano_df[-vid, ]

vol_lm <- lm(elevation ~ row + col + I(row^2) + I(col^2) + row:col, data = vol_train)
vol_tree <- rpart::rpart(elevation ~ row + col, data = vol_train, control = rpart::rpart.control(cp = 0.001))

vol_scores <- data.frame(
  model = c("Quadratic statistical surface", "ML regression tree"),
  test_RMSE = c(
    rmse(vol_test$elevation, predict(vol_lm, newdata = vol_test)),
    rmse(vol_test$elevation, predict(vol_tree, newdata = vol_test))
  )
)
knitr::kable(vol_scores, digits = 3, caption = "Predicting elevation from grid location.")
Predicting elevation from grid location.
model test_RMSE
Quadratic statistical surface 13.675
ML regression tree 6.666 
volcano_df$pred_lm <- predict(vol_lm, newdata = volcano_df)
volcano_df$pred_tree <- predict(vol_tree, newdata = volcano_df)
volcano_df$resid_lm <- volcano_df$elevation - volcano_df$pred_lm
volcano_df$resid_tree <- volcano_df$elevation - volcano_df$pred_tree

resid_long <- tidyr::pivot_longer(
  volcano_df,
  cols = c(resid_lm, resid_tree),
  names_to = "model",
  values_to = "residual"
)

ggplot2::ggplot(resid_long, ggplot2::aes(col, row, fill = residual)) +
  ggplot2::geom_raster() +
  ggplot2::coord_equal() +
  ggplot2::facet_wrap(~model) +
  ggplot2::labs(
    title = "Residual maps: what the models fail to explain",
    subtitle = "Residual structure is scientific information, not merely error",
    x = "south-north grid", y = "east-west grid", fill = "residual"
  )

Interpretation. Spatial residual maps are a very useful diagnostic. If residuals show structure, the model is missing spatial geometry. ML may reduce prediction error, but a statistical residual map explains where the model is scientifically incomplete.

14 Real data 5: petroleum rock samples and small-data inference

The rock dataset contains measurements on 48 petroleum rock samples. Variables include pore area, perimeter, shape, and permeability. This is a good example of small-data scientific inference: ML may be unstable, and uncertainty matters.

data(rock)
rock_df <- rock |>
  dplyr::mutate(
    log_perm = log(perm),
    log_area = log(area),
    log_peri = log(peri)
  )

rock_lm <- lm(log_perm ~ log_area + log_peri + shape, data = rock_df)
rock_tree <- rpart::rpart(log_perm ~ log_area + log_peri + shape, data = rock_df,
                          control = rpart::rpart.control(minsplit = 10, cp = 0.01))

knitr::kable(
  data.frame(term = names(coef(rock_lm)), estimate = coef(rock_lm), confint(rock_lm)),
  digits = 3,
  caption = "Small-data petrophysics: regression effects for log permeability."
)
Small-data petrophysics: regression effects for log permeability.
term estimate X2.5.. X97.5..
(Intercept) (Intercept) 4.848 0.125 9.572
log_area log_area 3.297 2.250 4.344
log_peri log_peri -3.747 -4.634 -2.860
shape shape 0.907 -2.836 4.650 
ggplot2::ggplot(rock_df, ggplot2::aes(log_area, log_perm, size = shape)) +
  ggplot2::geom_point(alpha = 0.65) +
  ggplot2::geom_smooth(method = "lm", se = TRUE) +
  ggplot2::labs(
    title = "Petroleum rock samples: pore area and permeability",
    subtitle = "Small data: uncertainty and scientific interpretability matter strongly",
    x = "log pore area", y = "log permeability", size = "shape"
  )

14.1 Leave-one-out comparison: small data and model stability

loocv_rmse_lm <- function(data) {
  pred <- numeric(nrow(data))
  for (i in seq_len(nrow(data))) {
    fit <- lm(log_perm ~ log_area + log_peri + shape, data = data[-i, ])
    pred[i] <- predict(fit, newdata = data[i, ])
  }
  rmse(data$log_perm, pred)
}

loocv_rmse_tree <- function(data) {
  pred <- numeric(nrow(data))
  for (i in seq_len(nrow(data))) {
    fit <- rpart::rpart(log_perm ~ log_area + log_peri + shape, data = data[-i, ],
                        control = rpart::rpart.control(minsplit = 10, cp = 0.01))
    pred[i] <- predict(fit, newdata = data[i, ])
  }
  rmse(data$log_perm, pred)
}

rock_cv <- data.frame(
  model = c("Linear small-data inference model", "ML tree"),
  LOOCV_RMSE = c(loocv_rmse_lm(rock_df), loocv_rmse_tree(rock_df))
)
knitr::kable(rock_cv, digits = 3, caption = "Leave-one-out cross-validation in a small scientific dataset.")
Leave-one-out cross-validation in a small scientific dataset.
model LOOCV_RMSE
Linear small-data inference model 0.948
ML tree 1.091

Interpretation. In small scientific datasets, a transparent statistical model may be preferable to a highly flexible learner. ML is not wrong, but flexibility must be paid for by data volume and validation strength.

15 Where statistics fails, where ML fails, and how they cooperate

Situation Statistical danger ML danger Good data-science response
Nonlinear relation underfitting by linear model overfitting if too flexible compare interpretable and flexible models
Time series false independent-sample uncertainty random-split leakage chronological validation, stochastic-process models
Hidden state ignoring measurement noise treating observations as truth state-space modelling and filtering
Small sample unstable p-values if assumptions fail unstable flexible learners regularization, resampling, transparent modelling
Extrapolation wrong parametric law no extrapolation mechanism scientific structure plus validation
Spatial surface residual spatial dependence ignored good prediction but little explanation residual maps, spatial models, uncertainty

16 A practical workflow for basic-science researchers

  1. Define the scientific question. Prediction, effect estimation, mechanism, monitoring, control, or decision?
  2. Define the observational unit. Patient, star, specimen, time point, pixel, sensor, or molecule?
  3. Write the estimand. Mean, slope, classification probability, hidden state, change point, hazard, or risk?
  4. Plot first. Scatterplots, time plots, residual maps, density plots.
  5. Fit a transparent baseline. Linear/logistic model, ANOVA, AR model, simple smoother.
  6. Fit a flexible ML model. Tree, random forest, boosting, kernel, neural net, etc.
  7. Validate honestly. Use test sets, cross-validation, chronological split, spatial blocking, or leave-one-out.
  8. Quantify uncertainty. Confidence intervals, prediction intervals, bootstrap, conformal intervals, filters.
  9. Diagnose misspecification. Residuals, calibration, subgroup errors, out-of-distribution checks.
  10. Report scientific meaning, not only accuracy. What was learned about the phenomenon?

17 Closing message

Statistical inference and ML are not enemies. They are two languages for scientific data science.

For basic sciences, the most powerful workflow is not statistics versus ML, but:

\[ \text{scientific structure} + \text{statistical inference} + \text{machine learning} + \text{uncertainty quantification}. \]

18 References and further reading