Statistical Inference and Hypothesis Testing for Basic-Science Data Science

Part 1: From Probability Models to Evidence, Uncertainty, Simulation, and Real Data

1 Why this tutorial exists

Many students in physics, chemistry, astronomy, geoscience, biological sciences, and materials science already know probability, random variables, distributions, and perhaps stochastic processes. But statistical inference asks a different question:

Given finite, noisy, incomplete data, what can we say about the unknown scientific mechanism?

A physicist may write a model for a process. A statistician asks:

• What is the estimand?
• What is random?
• What is the sampling distribution of the estimator?
• How uncertain is the estimate?
• Is the observed pattern stronger than expected by chance?
• Does the model fit the data?
• Would the conclusion survive new data?

This Part 1 is not Bayesian yet. Bayesian inference will come later. Here we build the classical/frequentist statistical backbone: estimators, sampling distributions, confidence intervals, hypothesis tests, power, likelihood, regression inference, resampling, and model criticism.

2 Data science for basic science: Statistics and ML

In scientific data science, statistics and machine learning should not be rivals. They answer related but distinct questions.

stats_ml <- data.frame(
  Perspective = c("Statistical inference", "Machine learning", "Scientific data science"),
  Central_question = c(
    "What can we infer about a population/mechanism from finite noisy data?",
    "How accurately can we predict or classify new observations?",
    "How do we combine prediction, mechanism, uncertainty, and decision?"
  ),
  Typical_output = c(
    "Estimate, standard error, confidence interval, p-value, uncertainty statement",
    "Prediction, classifier, feature importance, validation error",
    "Model-based insight plus prediction plus uncertainty-aware scientific interpretation"
  ),
  Risk = c(
    "A beautiful test can be invalid if assumptions fail",
    "A high-accuracy model can exploit shortcuts or fail outside training distribution",
    "Scientific claims need both predictive validation and inferential discipline"
  )
)

kable(stats_ml, caption = "Statistics, ML, and scientific data science: complementary roles.")

Statistics, ML, and scientific data science: complementary roles.

Perspective	Central_question	Typical_output	Risk
Statistical inference	What can we infer about a population/mechanism from finite noisy data?	Estimate, standard error, confidence interval, p-value, uncertainty statement	A beautiful test can be invalid if assumptions fail
Machine learning	How accurately can we predict or classify new observations?	Prediction, classifier, feature importance, validation error	A high-accuracy model can exploit shortcuts or fail outside training distribution
Scientific data science	How do we combine prediction, mechanism, uncertainty, and decision?	Model-based insight plus prediction plus uncertainty-aware scientific interpretation	Scientific claims need both predictive validation and inferential discipline

Important philosophy. ML often optimizes predictive loss. Statistical inference asks whether the observed pattern is scientifically reliable, uncertain, and transportable. A basic-science researcher needs both.

3 Probability model, data, parameter, estimand

A statistical model is a family of probability laws:

\[ \mathcal{P}=\{P_\theta:\theta\in\Theta\}. \]

Data are usually written as

\[ X_1,\ldots,X_n \sim P_{\theta_0}, \]

where \(\theta_0\) is the unknown true parameter, if the model is correctly specified.

An estimand is the target quantity. Examples:

\[ \mu = E(X),\quad \sigma^2=\operatorname{Var}(X),\quad p=P(X>c),\quad \beta_1=\text{slope in a regression model}. \]

An estimator is a function of data:

\[ \widehat{\theta}=T(X_1,\ldots,X_n). \]

A central question is:

\[ \widehat{\theta}\quad\text{versus}\quad \theta_0. \]

We want small bias, small variance, and meaningful uncertainty quantification.

4 The inference triangle

The basic triangle is:

\[ \text{Data} \longrightarrow \text{Estimator/Test statistic} \longrightarrow \text{Sampling distribution}. \]

The sampling distribution is the distribution of the estimator under repeated sampling:

\[ T(X_1,\ldots,X_n). \]

Most confusion in statistics disappears when students understand this point:
a confidence interval is random before data are observed; the parameter is fixed in the frequentist framework.

5 Simulation 1: Law of Large Numbers and consistency

The sample mean

\[ \bar{X}_n=\frac{1}{n}\sum_{i=1}^n X_i \]

is a natural estimator of \(\mu=E(X)\). Under mild conditions,

\[ \bar{X}_n \to \mu \]

as \(n\to\infty\). This is the statistical meaning of “truth is approached as larger data are obtained” — but only when the model and sampling mechanism are appropriate.

set.seed(1)

mu_true <- 2
sigma_true <- 1.5
n_path <- 1000
x <- rnorm(n_path, mean = mu_true, sd = sigma_true)
running_mean <- cumsum(x) / seq_along(x)

df_lln <- data.frame(
  n = seq_along(x),
  running_mean = running_mean,
  truth = mu_true
)

ggplot(df_lln, aes(x = n, y = running_mean)) +
  geom_line(color = "#1f77b4", linewidth = 0.8) +
  geom_hline(aes(yintercept = truth), color = "red", linetype = "dashed", linewidth = 0.9) +
  labs(
    title = "Law of Large Numbers: sample mean approaches the true mean",
    subtitle = "One simulated experiment from N(2, 1.5^2)",
    x = "Sample size n",
    y = "Running sample mean"
  )

6 Simulation 2: Sampling distribution and Central Limit Theorem

The Central Limit Theorem says that for many distributions,

\[ \sqrt{n}\frac{\bar{X}_n-\mu}{\sigma}\Rightarrow N(0,1). \]

This is why normal approximations, standard errors, z-tests and t-tests appear everywhere.

set.seed(2)

B <- 5000
n <- 40

# Use an exponential distribution, which is skewed, to show CLT robustness
mu_exp <- 1
sd_exp <- 1

z <- replicate(B, {
  sample_x <- rexp(n, rate = 1)
  sqrt(n) * (mean(sample_x) - mu_exp) / sd_exp
})

df_clt <- data.frame(z = z)

ggplot(df_clt, aes(x = z)) +
  geom_histogram(aes(y = after_stat(density)), bins = 45, fill = "#9ecae1", color = "white") +
  stat_function(fun = dnorm, color = "red", linewidth = 1.1) +
  labs(
    title = "Central Limit Theorem from skewed data",
    subtitle = "Even exponential observations give an approximately normal standardized mean when n is moderate",
    x = expression(sqrt(n) * (bar(X) - mu) / sigma),
    y = "Density"
  )

7 Estimator quality: bias, variance, MSE

For an estimator \(\widehat{\theta}\),

\[ \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}. \]

set.seed(3)

B <- 4000
theta <- 5
n_values <- c(5, 10, 20, 50, 100, 300)

est_summary <- do.call(rbind, lapply(n_values, function(n) {
  est <- replicate(B, mean(rnorm(n, mean = theta, sd = 2)))
  data.frame(
    n = n,
    mean_estimate = mean(est),
    bias = mean(est) - theta,
    variance = var(est),
    mse = mean((est - theta)^2)
  )
}))

kable(round_df(est_summary, 5),
      caption = "Simulation: sample mean becomes more stable as sample size increases.")

Simulation: sample mean becomes more stable as sample size increases.

n	mean_estimate	bias	variance	mse
5	4.97400	-0.02600	0.83580	0.83626
10	5.01550	0.01550	0.40773	0.40787
20	4.99823	-0.00177	0.19958	0.19953
50	4.99641	-0.00359	0.08337	0.08337
100	5.00395	0.00395	0.04070	0.04070
300	5.00295	0.00295	0.01356	0.01357

ggplot(est_summary, aes(x = n, y = variance)) +
  geom_point(size = 3, color = "#2ca25f") +
  geom_line(color = "#2ca25f", linewidth = 0.9) +
  scale_x_continuous(breaks = n_values) +
  labs(
    title = "Estimator variance decreases as sample size increases",
    subtitle = "For the sample mean of Normal data, variance is sigma^2/n",
    x = "Sample size",
    y = "Empirical variance of estimator"
  )

8 Confidence intervals

A \(95\%\) confidence interval procedure means:

\[ P_{\theta_0}\{\theta_0\in C(X_1,\ldots,X_n)\}\approx 0.95. \]

Interpretation:

If the same procedure is repeated over many independent datasets, around 95% of the constructed intervals should cover the true value.

It does not mean that after seeing the data, the fixed parameter has 95% probability of being in this particular interval. That Bayesian interpretation comes later.

set.seed(4)

B <- 3000
mu <- 10
sigma <- 3
n_grid <- c(5, 10, 20, 50, 100)

ci_summary <- do.call(rbind, lapply(n_grid, function(n) {
  out <- replicate(B, {
    x <- rnorm(n, mu, sigma)
    se <- sd(x) / sqrt(n)
    ci <- mean(x) + qt(c(0.025, 0.975), df = n - 1) * se
    c(estimate = mean(x), lower = ci[1], upper = ci[2], width = diff(ci), cover = ci[1] <= mu && mu <= ci[2])
  })
  out <- t(out)
  data.frame(
    n = n,
    mean_estimate = mean(out[, "estimate"]),
    mean_width = mean(out[, "width"]),
    empirical_coverage = mean(out[, "cover"])
  )
}))

kable(round_df(ci_summary, 4),
      caption = "Repeated-sampling verification of 95% confidence intervals.")

Repeated-sampling verification of 95% confidence intervals.

n	mean_estimate	mean_width	empirical_coverage
5	9.9797	7.0202	0.9563
10	10.0093	4.1532	0.9557
20	10.0001	2.7714	0.9543
50	9.9954	1.6978	0.9497
100	9.9967	1.1862	0.9460

ggplot(ci_summary, aes(x = n, y = mean_width)) +
  geom_point(size = 3, color = "#d95f0e") +
  geom_line(color = "#d95f0e", linewidth = 0.9) +
  scale_x_continuous(breaks = n_grid) +
  labs(
    title = "Confidence intervals become narrower as data increase",
    subtitle = "The target coverage remains approximately 95%, but precision improves",
    x = "Sample size",
    y = "Average width of 95% CI"
  )

9 Hypothesis testing: the logic

A hypothesis test begins with two statements:

\[ H_0:\theta=\theta_0 \quad\text{versus}\quad H_1:\theta\neq \theta_0. \]

A test statistic \(T(X)\) is chosen. A p-value is

\[ p=P_{H_0}\{\text{test statistic at least as extreme as observed}\}. \]

Important:

• A small p-value is evidence against \(H_0\), relative to the model.
• A p-value is not the probability that \(H_0\) is true.
• A large p-value is not proof of no effect.
• With very large \(n\), tiny scientifically irrelevant effects can become statistically significant.
• With small \(n\), scientifically important effects can be missed.

10 Simulation 3: p-values under the null and alternative

Under a valid continuous null hypothesis, p-values are approximately uniform on \([0,1]\). Under an alternative, p-values concentrate near zero.

set.seed(5)

B <- 5000
n <- 30

p_null <- replicate(B, {
  x <- rnorm(n, mean = 0, sd = 1)
  t.test(x, mu = 0)$p.value
})

p_alt <- replicate(B, {
  x <- rnorm(n, mean = 0.45, sd = 1)
  t.test(x, mu = 0)$p.value
})

df_p <- rbind(
  data.frame(p_value = p_null, setting = "Null true"),
  data.frame(p_value = p_alt, setting = "Alternative true")
)

ggplot(df_p, aes(x = p_value, fill = setting)) +
  geom_histogram(bins = 40, color = "white", alpha = 0.75, position = "identity") +
  facet_wrap(~ setting, ncol = 1) +
  scale_fill_manual(values = c("Null true" = "#a6bddb", "Alternative true" = "#fb9a99")) +
  labs(
    title = "What p-values look like under null and alternative",
    x = "p-value",
    y = "Count"
  ) +
  theme(legend.position = "none")

11 Type I error, Type II error, and power

A test has:

\[ \alpha=P(\text{reject }H_0\mid H_0\text{ true}), \]

\[ \beta=P(\text{fail to reject }H_0\mid H_1\text{ true}), \]

\[ \text{Power}=1-\beta. \]

set.seed(6)

B <- 2000
n <- 30
delta_grid <- seq(0, 1.2, by = 0.1)

power_df <- do.call(rbind, lapply(delta_grid, function(delta) {
  pvals <- replicate(B, {
    x <- rnorm(n, mean = delta, sd = 1)
    t.test(x, mu = 0)$p.value
  })
  data.frame(delta = delta, power = mean(pvals < 0.05))
}))

ggplot(power_df, aes(x = delta, y = power)) +
  geom_line(color = "#756bb1", linewidth = 1) +
  geom_point(color = "#756bb1", size = 2.5) +
  geom_hline(yintercept = 0.05, linetype = "dashed", color = "gray40") +
  labs(
    title = "Power curve for a one-sample t-test",
    subtitle = "Larger true effects are easier to detect",
    x = "True standardized effect size",
    y = "Power at alpha = 0.05"
  )

12 Likelihood and maximum likelihood

For independent observations with density \(f(x;\theta)\),

\[ L(\theta)=\prod_{i=1}^n f(X_i;\theta). \]

The maximum likelihood estimator is

\[ \widehat{\theta}_{MLE}=\arg\max_{\theta} L(\theta). \]

Likelihood is not a probability distribution over \(\theta\). It is a function of \(\theta\) after data are observed.

set.seed(7)

x <- rnorm(25, mean = 3, sd = 1.2)
sigma_known <- 1.2
theta_grid <- seq(1, 5, length.out = 400)

loglik <- sapply(theta_grid, function(theta) {
  sum(dnorm(x, mean = theta, sd = sigma_known, log = TRUE))
})

df_lik <- data.frame(theta = theta_grid, loglik = loglik)
theta_hat <- theta_grid[which.max(loglik)]

ggplot(df_lik, aes(x = theta, y = loglik)) +
  geom_line(color = "#08519c", linewidth = 1) +
  geom_vline(xintercept = mean(x), color = "red", linetype = "dashed") +
  labs(
    title = "Likelihood as a function of the unknown mean",
    subtitle = paste("MLE =", round(theta_hat, 3), "; sample mean =", round(mean(x), 3)),
    x = expression(theta),
    y = "Log-likelihood"
  )

13 Real-data example 1: Mercury vapor pressure and scientific transformation

The pressure dataset contains 19 observations on temperature in degrees Celsius and vapor pressure of mercury in millimeters of mercury. The relationship is strongly nonlinear on the raw scale.

Scientific lesson: a transformation can convert a nonlinear physical relation into a nearly linear inference problem.

data(pressure)

kable(head(pressure, 8), caption = "First rows of the mercury vapor-pressure dataset.")

First rows of the mercury vapor-pressure dataset.

temperature	pressure
0	0.0002
20	0.0012
40	0.0060
60	0.0300
80	0.0900
100	0.2700
120	0.7500
140	1.8500

ggplot(pressure, aes(x = temperature, y = pressure)) +
  geom_point(size = 3, color = "#d94801") +
  geom_smooth(method = "lm", se = TRUE, color = "gray30") +
  labs(
    title = "Raw scale: vapor pressure rises nonlinearly with temperature",
    x = "Temperature (degree Celsius)",
    y = "Pressure (mm Hg)"
  )

pressure$log_pressure <- log(pressure$pressure)
fit_pressure <- lm(log_pressure ~ temperature, data = pressure)

coef_table <- summary(fit_pressure)$coefficients
kable(round_df(as.data.frame(coef_table), 5),
      caption = "Linear model on log pressure: coefficient estimates and hypothesis tests.")

Linear model on log pressure: coefficient estimates and hypothesis tests.

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	-6.06814	0.48383	-12.54188	0
temperature	0.03979	0.00230	17.32976	0

new_temp <- data.frame(temperature = seq(min(pressure$temperature), max(pressure$temperature), length.out = 200))
pred_conf <- predict(fit_pressure, newdata = new_temp, interval = "confidence")
pred_pred <- predict(fit_pressure, newdata = new_temp, interval = "prediction")
plot_df <- cbind(new_temp, as.data.frame(pred_conf), pred_lwr = pred_pred[, "lwr"], pred_upr = pred_pred[, "upr"])

ggplot(pressure, aes(x = temperature, y = log_pressure)) +
  geom_point(size = 3, color = "#d94801") +
  geom_ribbon(data = plot_df, aes(x = temperature, ymin = pred_lwr, ymax = pred_upr), alpha = 0.12, fill = "#6baed6", inherit.aes = FALSE) +
  geom_ribbon(data = plot_df, aes(x = temperature, ymin = lwr, ymax = upr), alpha = 0.25, fill = "#3182bd", inherit.aes = FALSE) +
  geom_line(data = plot_df, aes(x = temperature, y = fit), color = "#08519c", linewidth = 1, inherit.aes = FALSE) +
  labs(
    title = "Log-transform: physical relation becomes much more linear",
    subtitle = "Dark band: confidence interval for mean response; light band: prediction interval",
    x = "Temperature (degree Celsius)",
    y = "log pressure"
  )

Inference statement. The slope test asks whether increasing temperature is associated with higher log vapor pressure:

\[ H_0:\beta_1=0,\quad H_1:\beta_1\neq 0. \]

The result is not merely a fitted curve. It gives an estimate, uncertainty, and evidence against no temperature effect.

14 Real-data example 2: Old Faithful geyser and regression inference

The faithful dataset contains waiting time between eruptions and eruption duration for Old Faithful geyser. This is a physical process with a clear pattern and uncertainty.

data(faithful)

fit_faithful <- lm(waiting ~ eruptions, data = faithful)
faithful_coef <- summary(fit_faithful)$coefficients

kable(round_df(as.data.frame(faithful_coef), 5),
      caption = "Regression inference: waiting time versus eruption duration.")

Regression inference: waiting time versus eruption duration.

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	33.47440	1.15487	28.98534	0
eruptions	10.72964	0.31475	34.08904	0

new_erupt <- data.frame(eruptions = seq(min(faithful$eruptions), max(faithful$eruptions), length.out = 200))
pred_faith <- predict(fit_faithful, newdata = new_erupt, interval = "prediction")
pred_faith_df <- cbind(new_erupt, as.data.frame(pred_faith))

ggplot(faithful, aes(x = eruptions, y = waiting)) +
  geom_point(alpha = 0.65, color = "#2b8cbe") +
  geom_ribbon(data = pred_faith_df, aes(x = eruptions, ymin = lwr, ymax = upr), fill = "#a6bddb", alpha = 0.25, inherit.aes = FALSE) +
  geom_line(data = pred_faith_df, aes(x = eruptions, y = fit), color = "#045a8d", linewidth = 1, inherit.aes = FALSE) +
  labs(
    title = "Old Faithful: regression as scientific prediction with uncertainty",
    subtitle = "Longer eruptions tend to be followed by longer waiting times",
    x = "Eruption duration (minutes)",
    y = "Waiting time to next eruption (minutes)"
  )

15 Hypothesis test versus scientific effect size

A p-value can be small because the effect is large, or because the sample size is large. Always ask:

\[ \text{Is the effect scientifically meaningful?} \]

r_faith <- cor(faithful$eruptions, faithful$waiting)
r2_faith <- summary(fit_faithful)$r.squared
effect_table <- data.frame(
  Quantity = c("Correlation", "R-squared", "Slope estimate", "Slope p-value"),
  Value = c(r_faith, r2_faith, coef(fit_faithful)[2], summary(fit_faithful)$coefficients[2,4])
)
kable(round_df(effect_table, 5),
      caption = "Effect size and statistical evidence should be read together.")

Effect size and statistical evidence should be read together.

Quantity	Value
Correlation	0.90081
R-squared	0.81146
Slope estimate	10.72964
Slope p-value	0.00000

16 Real-data example 3: Fiji earthquakes and geophysical inference

The quakes dataset contains 1000 seismic events near Fiji with latitude, longitude, depth, magnitude and number of reporting stations.

Scientific question examples:

• Are deeper events systematically different in magnitude?
• Is reported magnitude associated with number of stations?
• Can a simple model detect risk-related patterns?

data(quakes)

quakes_summary <- data.frame(
  variable = names(quakes),
  mean = sapply(quakes, mean),
  sd = sapply(quakes, sd),
  min = sapply(quakes, min),
  max = sapply(quakes, max)
)

kable(round_df(quakes_summary, 3), caption = "Summary of Fiji earthquake variables.")

Summary of Fiji earthquake variables.

	variable	mean	sd	min	max
lat	lat	-20.643	5.029	-38.59	-10.72
long	long	179.462	6.069	165.67	188.13
depth	depth	311.371	215.535	40.00	680.00
mag	mag	4.620	0.403	4.00	6.40
stations	stations	33.418	21.900	10.00	132.00

ggplot(quakes, aes(x = long, y = lat, color = depth, size = mag)) +
  geom_point(alpha = 0.65) +
  scale_color_viridis_c() +
  labs(
    title = "Fiji seismic events: location, depth, and magnitude",
    x = "Longitude",
    y = "Latitude",
    color = "Depth (km)",
    size = "Magnitude"
  )

fit_quake <- lm(mag ~ depth + stations, data = quakes)
kable(round_df(as.data.frame(summary(fit_quake)$coefficients), 5),
      caption = "Linear inference for earthquake magnitude.")

Linear inference for earthquake magnitude.

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	4.20322	0.01522	276.10291	0
depth	-0.00032	0.00003	-10.70456	0
stations	0.01543	0.00029	53.13506	0

quakes$depth_group <- ifelse(quakes$depth > median(quakes$depth), "Deep", "Shallow")
quake_test <- t.test(mag ~ depth_group, data = quakes)

quake_group_table <- aggregate(mag ~ depth_group, data = quakes, function(v) c(mean = mean(v), sd = sd(v), n = length(v)))
quake_group_table <- data.frame(
  depth_group = quake_group_table$depth_group,
  mean = quake_group_table$mag[, "mean"],
  sd = quake_group_table$mag[, "sd"],
  n = quake_group_table$mag[, "n"]
)

kable(round_df(quake_group_table, 4), caption = "Magnitude by depth group.")

Magnitude by depth group.

depth_group	mean	sd	n
Deep	4.5232	0.3748	500
Shallow	4.7176	0.4067	500

quake_test

## 
##  Welch Two Sample t-test
## 
## data:  mag by depth_group
## t = -7.8602, df = 991.42, p-value = 9.963e-15
## alternative hypothesis: true difference in means between group Deep and group Shallow is not equal to 0
## 95 percent confidence interval:
##  -0.2429337 -0.1458663
## sample estimates:
##    mean in group Deep mean in group Shallow 
##                4.5232                4.7176

ggplot(quakes, aes(x = depth_group, y = mag, fill = depth_group)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.15, alpha = 0.25) +
  scale_fill_manual(values = c("Deep" = "#fb6a4a", "Shallow" = "#6baed6")) +
  labs(
    title = "Hypothesis test as a scientific comparison",
    subtitle = "Magnitude distribution for shallow versus deep events",
    x = "Depth group",
    y = "Magnitude"
  ) +
  theme(legend.position = "none")

Warning. Earthquakes are spatially dependent. A simple t-test is only an introductory demonstration; serious geophysical inference needs spatial, temporal, and physical modelling.

17 Real-data example 4: Nile flow and changepoint thinking

The Nile dataset records annual Nile flow at Aswan from 1871 to 1970 and has an apparent changepoint near 1898.

A basic inference question:

\[ H_0:\mu_{\text{before}}=\mu_{\text{after}} \]

versus

\[ H_1:\mu_{\text{before}}\neq\mu_{\text{after}}. \]

data(Nile)

nile_df <- data.frame(
  year = as.numeric(time(Nile)),
  flow = as.numeric(Nile),
  period = ifelse(as.numeric(time(Nile)) <= 1898, "1871-1898", "1899-1970")
)

nile_table <- aggregate(flow ~ period, data = nile_df, function(v) c(mean = mean(v), sd = sd(v), n = length(v)))
nile_table <- data.frame(
  period = nile_table$period,
  mean = nile_table$flow[, "mean"],
  sd = nile_table$flow[, "sd"],
  n = nile_table$flow[, "n"]
)

kable(round_df(nile_table, 3), caption = "Nile flow before and after 1898.")

Nile flow before and after 1898.

period	mean	sd	n
1871-1898	1097.750	134.996	28
1899-1970	849.972	124.776	72

ggplot(nile_df, aes(x = year, y = flow)) +
  geom_line(color = "#225ea8", linewidth = 0.8) +
  geom_point(color = "#225ea8", size = 1.8) +
  geom_vline(xintercept = 1898, color = "red", linetype = "dashed") +
  labs(
    title = "Nile annual flow at Aswan",
    subtitle = "A simple changepoint-style scientific question",
    x = "Year",
    y = expression(paste("Flow (", 10^8, " ", m^3, ")"))
  )

t.test(flow ~ period, data = nile_df)

## 
##  Welch Two Sample t-test
## 
## data:  flow by period
## t = 8.4145, df = 45.991, p-value = 7.308e-11
## alternative hypothesis: true difference in means between group 1871-1898 and group 1899-1970 is not equal to 0
## 95 percent confidence interval:
##  188.5048 307.0508
## sample estimates:
## mean in group 1871-1898 mean in group 1899-1970 
##               1097.7500                849.9722

Interpretation. This is an introductory test for a mean shift. It is not a complete time-series analysis. Time dependence and physical hydrology should be considered for serious modelling.

18 Real-data example 5: Solar activity and dependence

The yearly sunspot series is a classic solar-activity time series. It is not independent. It has cyclic structure.

A naive i.i.d. test can miss the scientific structure. Statistical inference must respect dependence.

data(sunspot.year)

sun_df <- data.frame(
  year = as.numeric(time(sunspot.year)),
  sunspots = as.numeric(sunspot.year)
)

ggplot(sun_df, aes(x = year, y = sunspots)) +
  geom_line(color = "#e6550d", linewidth = 0.8) +
  labs(
    title = "Yearly sunspot numbers",
    subtitle = "A scientific time series: cyclic, dependent, and non-i.i.d.",
    x = "Year",
    y = "Sunspot number"
  )

acf_vals <- acf(sun_df$sunspots, plot = FALSE, lag.max = 40)
acf_df <- data.frame(lag = as.numeric(acf_vals$lag), acf = as.numeric(acf_vals$acf))

ggplot(acf_df, aes(x = lag, y = acf)) +
  geom_col(fill = "#fb6a4a") +
  geom_hline(yintercept = 0, color = "gray30") +
  labs(
    title = "Autocorrelation of yearly sunspot numbers",
    subtitle = "Dependence violates the simplest independent-sampling assumptions",
    x = "Lag",
    y = "Autocorrelation"
  )

19 Lesson: inference depends on the data-generating structure

For i.i.d. data:

\[ \operatorname{SE}(\bar{X})\approx s/\sqrt{n}. \]

For dependent data, this can be too optimistic. In time series, nearby observations carry partially repeated information.

set.seed(8)

x <- sun_df$sunspots
n <- length(x)

iid_se <- sd(x) / sqrt(n)

# Simple block bootstrap for the mean
B <- 1000
block_len <- 11
starts <- 1:(n - block_len + 1)
block_means <- replicate(B, {
  idx <- integer(0)
  while (length(idx) < n) {
    s <- sample(starts, 1)
    idx <- c(idx, s:(s + block_len - 1))
  }
  idx <- idx[1:n]
  mean(x[idx])
})

block_se <- sd(block_means)

se_table <- data.frame(
  Method = c("Naive i.i.d. SE", "Block bootstrap SE"),
  Standard_error = c(iid_se, block_se)
)

kable(round_df(se_table, 4),
      caption = "Time dependence can change uncertainty quantification.")

Time dependence can change uncertainty quantification.

Method	Standard_error
Naive i.i.d. SE	2.3220
Block bootstrap SE	3.5967

20 Real-data example 6: Galaxy velocities and multimodality

The MASS::galaxies dataset contains velocities of 82 galaxies from the Corona Borealis region. Multimodality in such surveys is relevant to voids and superclusters.

A simple mean and standard deviation may hide multiple structures.

data(galaxies, package = "MASS")
gal_df <- data.frame(velocity = as.numeric(galaxies) / 1000)

ggplot(gal_df, aes(x = velocity)) +
  geom_histogram(aes(y = after_stat(density)), bins = 18, fill = "#9ecae1", color = "white") +
  geom_density(color = "#08519c", linewidth = 1) +
  labs(
    title = "Galaxy velocities: a distribution may contain structure",
    subtitle = "A single normal distribution may not capture multimodality",
    x = "Velocity (1000 km/sec)",
    y = "Density"
  )

mu_gal <- mean(gal_df$velocity)
sd_gal <- sd(gal_df$velocity)

ggplot(gal_df, aes(x = velocity)) +
  geom_histogram(aes(y = after_stat(density)), bins = 18, fill = "#bdbdbd", color = "white", alpha = 0.7) +
  geom_density(color = "#08519c", linewidth = 1) +
  stat_function(fun = dnorm, args = list(mean = mu_gal, sd = sd_gal), color = "red", linewidth = 1, linetype = "dashed") +
  labs(
    title = "Model criticism: normal curve versus empirical density",
    subtitle = "Inference requires checking whether the model family is scientifically plausible",
    x = "Velocity (1000 km/sec)",
    y = "Density"
  )

21 A bridge to ML: clustering suggests structure, inference asks uncertainty

A simple k-means clustering can reveal groups, but by itself it does not give uncertainty about the number of groups or the physical mechanism.

set.seed(9)

k_grid <- 1:6
wss <- sapply(k_grid, function(k) {
  kmeans(gal_df$velocity, centers = k, nstart = 25)$tot.withinss
})

wss_df <- data.frame(k = k_grid, withinss = wss)

ggplot(wss_df, aes(x = k, y = withinss)) +
  geom_point(size = 3, color = "#756bb1") +
  geom_line(color = "#756bb1", linewidth = 1) +
  labs(
    title = "ML-style exploratory clustering",
    subtitle = "K-means can suggest structure; inference must still ask uncertainty and physical interpretation",
    x = "Number of clusters",
    y = "Total within-cluster sum of squares"
  )

22 Simulation 4: larger data reveal truth — if the model is right

Let the true model be:

\[ Y_i=\beta_0+\beta_1 x_i+\epsilon_i,\quad \epsilon_i\sim N(0,\sigma^2). \]

As \(n\) increases, the regression slope estimate concentrates near the truth.

set.seed(10)

B <- 1000
beta0 <- 1
beta1 <- 2
sigma <- 1
n_grid <- c(20, 50, 100, 300, 1000)

reg_cons <- do.call(rbind, lapply(n_grid, function(n) {
  slopes <- replicate(B, {
    x <- runif(n, -1, 1)
    y <- beta0 + beta1 * x + rnorm(n, 0, sigma)
    coef(lm(y ~ x))[2]
  })
  data.frame(
    n = n,
    mean_slope = mean(slopes),
    sd_slope = sd(slopes),
    q025 = quantile(slopes, 0.025),
    q975 = quantile(slopes, 0.975)
  )
}))

kable(round_df(reg_cons, 4),
      caption = "Regression slope estimate concentrates near the truth as n increases.")

Regression slope estimate concentrates near the truth as n increases.

	n	mean_slope	sd_slope	q025	q975
2.5%	20	2.0024	0.3994	1.2418	2.7896
2.5%1	50	1.9947	0.2497	1.5299	2.5042
2.5%2	100	2.0056	0.1780	1.6632	2.3624
2.5%3	300	2.0047	0.0999	1.8121	2.1930
2.5%4	1000	2.0028	0.0529	1.8994	2.1030

ggplot(reg_cons, aes(x = n, y = mean_slope)) +
  geom_ribbon(aes(ymin = q025, ymax = q975), fill = "#c7e9c0", alpha = 0.6) +
  geom_line(color = "#238b45", linewidth = 1) +
  geom_point(color = "#238b45", size = 3) +
  geom_hline(yintercept = beta1, color = "red", linetype = "dashed") +
  scale_x_log10(breaks = n_grid) +
  labs(
    title = "As data increase, a correct estimator finds the truth",
    subtitle = "Simulation distribution of regression slope estimates",
    x = "Sample size (log scale)",
    y = "Estimated slope"
  )

23 Simulation 5: larger data do not save a wrong model

Now suppose the true relation is nonlinear:

\[ Y_i = 1 + 2X_i - 3X_i^2 + \epsilon_i. \]

If we fit a linear model, larger \(n\) makes us confidently estimate the wrong linear projection.

set.seed(11)

n <- 600
x <- runif(n, -2, 2)
y <- 1 + 2 * x - 3 * x^2 + rnorm(n, 0, 1)

df_mis <- data.frame(x = x, y = y)
fit_wrong <- lm(y ~ x, data = df_mis)
fit_right <- lm(y ~ x + I(x^2), data = df_mis)

newx <- data.frame(x = seq(-2, 2, length.out = 300))
newx$wrong <- predict(fit_wrong, newdata = newx)
newx$right <- predict(fit_right, newdata = newx)

ggplot(df_mis, aes(x = x, y = y)) +
  geom_point(alpha = 0.35, color = "gray40") +
  geom_line(data = newx, aes(y = wrong), color = "red", linewidth = 1) +
  geom_line(data = newx, aes(y = right), color = "#08519c", linewidth = 1) +
  labs(
    title = "A wrong model becomes confidently wrong with larger data",
    subtitle = "Red: linear fit; Blue: quadratic fit",
    x = "x",
    y = "y"
  )

mis_table <- rbind(
  Linear_model = c(AIC = AIC(fit_wrong), R2 = summary(fit_wrong)$r.squared),
  Quadratic_model = c(AIC = AIC(fit_right), R2 = summary(fit_right)$r.squared)
)
kable(round_df(as.data.frame(mis_table), 4),
      caption = "Model comparison: misspecification can dominate sample size.")

Model comparison: misspecification can dominate sample size.

	AIC	R2
Linear_model	3224.261	0.2771
Quadratic_model	1721.034	0.9412

24 Hypothesis testing under model misspecification

The test for the slope in the wrong linear model may be highly significant, but the scientific conclusion “linear law” is wrong.

kable(round_df(as.data.frame(summary(fit_wrong)$coefficients), 5),
      caption = "The wrong model can give precise-looking coefficients.")

The wrong model can give precise-looking coefficients.

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	-2.70512	0.14467	-18.69851	0
x	1.96006	0.12947	15.13925	0

Lesson:

Statistical significance is conditional on the model. If the model is wrong, the test can answer the wrong question very precisely.

25 Bootstrap: uncertainty without closed-form formulas

Bootstrap approximates the sampling distribution by resampling the observed data.

For an estimator \(T\),

\[ \widehat{T}^{*b}=T(X_1^{*b},\ldots,X_n^{*b}). \]

set.seed(12)

B <- 1500
n <- nrow(faithful)
boot_slopes <- replicate(B, {
  idx <- sample(seq_len(n), size = n, replace = TRUE)
  coef(lm(waiting ~ eruptions, data = faithful[idx, ]))[2]
})

boot_ci <- quantile(boot_slopes, c(0.025, 0.975))
lm_ci <- confint(fit_faithful)[2, ]

boot_table <- data.frame(
  Method = c("Bootstrap percentile CI", "Classical lm CI"),
  lower = c(boot_ci[1], lm_ci[1]),
  upper = c(boot_ci[2], lm_ci[2])
)

kable(round_df(boot_table, 4), caption = "Bootstrap versus classical confidence interval for the geyser slope.")

Bootstrap versus classical confidence interval for the geyser slope.

	Method	lower	upper
2.5%	Bootstrap percentile CI	10.1648	11.3383
2.5 %	Classical lm CI	10.1100	11.3493

ggplot(data.frame(slope = boot_slopes), aes(x = slope)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40, fill = "#c7e9c0", color = "white") +
  geom_vline(xintercept = coef(fit_faithful)[2], color = "red", linewidth = 1) +
  labs(
    title = "Bootstrap sampling distribution of a regression slope",
    subtitle = "Resampling-based uncertainty for the Old Faithful regression",
    x = "Bootstrap slope estimate",
    y = "Density"
  )

26 Permutation test: evidence without normality assumptions

A permutation test asks: if group labels were irrelevant, how extreme would the observed difference be?

We demonstrate with the Nile before/after comparison. This simple permutation ignores time ordering, so it is mainly pedagogical.

set.seed(13)

obs_diff <- with(nile_df, mean(flow[period == "1871-1898"]) - mean(flow[period == "1899-1970"]))
B <- 3000

perm_diff <- replicate(B, {
  perm_period <- sample(nile_df$period)
  mean(nile_df$flow[perm_period == "1871-1898"]) - mean(nile_df$flow[perm_period == "1899-1970"])
})

perm_p <- mean(abs(perm_diff) >= abs(obs_diff))

data.frame(
  observed_difference = obs_diff,
  permutation_p_value = perm_p
) |> round_df(4) |> kable(caption = "Permutation test for the Nile before/after contrast.")

Permutation test for the Nile before/after contrast.

observed_difference	permutation_p_value
247.7778	0

ggplot(data.frame(diff = perm_diff), aes(x = diff)) +
  geom_histogram(bins = 45, fill = "#9ecae1", color = "white") +
  geom_vline(xintercept = obs_diff, color = "red", linewidth = 1) +
  geom_vline(xintercept = -obs_diff, color = "red", linewidth = 1, linetype = "dashed") +
  labs(
    title = "Permutation distribution under 'no group effect'",
    subtitle = "Pedagogical: time dependence should be respected in serious hydrology",
    x = "Difference in group means under random label permutations",
    y = "Count"
  )

27 Multiple testing: why scientific screening needs correction

High-dimensional science often tests many features. If all null hypotheses are true and we test \(m=1000\) features at level 0.05, we expect about 50 false positives.

set.seed(14)

m <- 1000
p_null <- runif(m)
num_sig <- sum(p_null < 0.05)

mt_table <- data.frame(
  number_of_tests = m,
  nominal_alpha = 0.05,
  expected_false_positives = m * 0.05,
  observed_false_positives_in_one_simulation = num_sig
)

kable(mt_table, caption = "Multiple testing under all null hypotheses.")

Multiple testing under all null hypotheses.

number_of_tests	nominal_alpha	expected_false_positives	observed_false_positives_in_one_simulation
1000	0.05	50	54

p_adj <- p.adjust(p_null, method = "BH")
mt_table2 <- data.frame(
  Rule = c("Raw p < 0.05", "BH-adjusted p < 0.05"),
  Discoveries = c(sum(p_null < 0.05), sum(p_adj < 0.05))
)

kable(mt_table2, caption = "False discoveries are controlled by multiplicity-aware procedures.")

False discoveries are controlled by multiplicity-aware procedures.

Rule	Discoveries
Raw p < 0.05	54
BH-adjusted p < 0.05	0

df_mt <- data.frame(p_value = p_null)

ggplot(df_mt, aes(x = p_value)) +
  geom_histogram(bins = 40, fill = "#bdbdbd", color = "white") +
  geom_vline(xintercept = 0.05, color = "red", linetype = "dashed", linewidth = 1) +
  labs(
    title = "Under the global null, some small p-values occur by chance",
    subtitle = "This is why high-dimensional science needs multiple-testing correction",
    x = "p-value",
    y = "Count"
  )

28 Summary of core ideas

summary_core <- data.frame(
  Concept = c(
    "Estimator",
    "Sampling distribution",
    "Standard error",
    "Confidence interval",
    "p-value",
    "Power",
    "Likelihood",
    "Bootstrap",
    "Permutation test",
    "Model checking"
  ),
  Meaning = c(
    "A statistic used to estimate an unknown scientific quantity",
    "Distribution of the estimator under repeated sampling",
    "Scale of estimator uncertainty",
    "Procedure with long-run coverage",
    "Extremeness of observed statistic under the null model",
    "Chance of detecting a real effect",
    "Data-supported score over parameter values",
    "Resampling approximation to sampling uncertainty",
    "Randomization-based evidence against exchangeability/no group effect",
    "Assess whether the assumed model is scientifically adequate"
  ),
  Warning = c(
    "Can be biased or unstable",
    "Depends on model and sampling design",
    "Can be wrong under dependence/misspecification",
    "Not posterior probability",
    "Not probability that the null is true",
    "Small samples can miss important effects",
    "High likelihood does not prove model truth",
    "Bootstrap also assumes data are representative",
    "Exchangeability assumption matters",
    "Significance without checking can mislead"
  )
)

kable(summary_core, caption = "General statistical inference: the essential vocabulary.")

General statistical inference: the essential vocabulary.

Concept	Meaning	Warning
Estimator	A statistic used to estimate an unknown scientific quantity	Can be biased or unstable
Sampling distribution	Distribution of the estimator under repeated sampling	Depends on model and sampling design
Standard error	Scale of estimator uncertainty	Can be wrong under dependence/misspecification
Confidence interval	Procedure with long-run coverage	Not posterior probability
p-value	Extremeness of observed statistic under the null model	Not probability that the null is true
Power	Chance of detecting a real effect	Small samples can miss important effects
Likelihood	Data-supported score over parameter values	High likelihood does not prove model truth
Bootstrap	Resampling approximation to sampling uncertainty	Bootstrap also assumes data are representative
Permutation test	Randomization-based evidence against exchangeability/no group effect	Exchangeability assumption matters
Model checking	Assess whether the assumed model is scientifically adequate	Significance without checking can mislead

29 Final message for basic-science researchers

Statistical inference is the bridge from probability to scientific evidence. Probability tells us how randomness behaves under a model. Statistical inference asks how finite data constrain unknown scientific quantities.

For researchers, the most important lessons are:

1. Define the estimand before computing.
2. Always separate estimation, uncertainty, and testing.
3. Larger data improve inference when the model and sampling design are appropriate.
4. Larger data can make wrong models confidently wrong.
5. A p-value is not a scientific conclusion.
6. ML and statistics are complementary: ML helps predict; statistics helps quantify evidence, uncertainty, and reliability.
7. Model criticism is not optional in science.

Bayesian inference will continue from here by treating unknown parameters themselves probabilistically through prior, likelihood, posterior, and posterior predictive distributions.