Model Selection on EEG Signals: Least Squares, AIC/BIC and Rejection ABC by Hand

An MSc statistics assignment handed me 201 samples of simulated EEG data: four input channels x1–x4, one output channel y, sampled every 2 ms over 0.4 seconds. Five candidate polynomial regression architectures were given. Pick the right one, justify the choice, validate it, then estimate the two most influential parameters with Approximate Bayesian Computation.

The interesting constraint was that everything had to be implemented from the mathematics — no lm(), no AIC(), no abc package. Normal equations, log-likelihood, information criteria and rejection sampling all written out in base R matrix algebra. Source and rendered notebook: github.com/rustamniraula90/eeg-signal-analysis.

That constraint turned out to matter. Twice, writing a formula out by hand made visible something a library call would have handled silently — and in both cases the pipeline produced numbers that looked fine but weren’t. Those are the parts of this worth reading.

---

The data, before any model

Before fitting anything, four questions about the data are worth answering, because each one rules something in or out later.

Are the inputs independent? They are emphatically not:

      x1    x2    x3    x4     y
x1 1.000 0.259 0.584 0.805 0.857
x2 0.259 1.000 0.523 0.364 0.194
x3 0.584 0.523 1.000 0.708 0.526
x4 0.805 0.364 0.708 1.000 0.775
y  0.857 0.194 0.526 0.775 1.000

x1 and x4 correlate at 0.805 with each other, and 0.857 / 0.775 with the output. x2 is nearly useless linearly (0.194 with y). This predicts something concrete: any model containing both x1 and x4 terms will have a nearly singular XᵀX, and their individual coefficients will be unstable even if the model’s predictions are fine. Collinearity attacks interpretation before it attacks fit.

Is anything non-Gaussian? The per-channel skewness says yes, and only in one place. The inputs are all near-symmetric (|skew| ≤ 0.32), but the output has skewness −1.79 — a long left tail, driven by a minimum of −9.60 against a median of 0.41. The inputs are drawn symmetric; the output isn’t. Something in the generating process is asymmetric, and a symmetric polynomial in the inputs can only reproduce that through odd powers.

What shape are the relationships? Loess curves through y against each input come back visibly bent, not straight. That’s the whole justification for a polynomial basis rather than plain multiple linear regression.

How much data is there? n = 201, and the candidate models want 3–6 parameters. That ratio is comfortable, but it’s also why BIC (whose penalty is k·log n = 5.30k here) will bite noticeably harder than AIC (2k).

Least squares written out is a conditioning problem

The estimator is one line of linear algebra:

θ̂ = (Xᵀ X)⁻¹ Xᵀ y

least_squares <- function(Xmat, yvec) {
  XtX <- t(Xmat) %*% Xmat
  theta_hat <- solve(XtX) %*% t(Xmat) %*% yvec
  return(theta_hat)
}

That’s the textbook form, and calling it is where the polynomial basis starts to hurt. The design matrix for Model 1 holds a column of x4 (range ±8.7) next to a column of x1⁴ (range up to 5,800). Squaring the matrix to form XᵀX squares that disparity too. The condition number of XᵀX:

Model	κ(XᵀX)
Model 1 (with x1², x1³, x1⁴, x2⁴)	2.8 × 10⁶
Model 2 (with x4, x1³, x3⁴)	7.4 × 10⁴

Double precision carries about 16 significant digits, so at κ ≈ 10⁶ you’re spending six of them on the inversion. Model 1 survives it. A model with x1⁵ or x1⁶ in the same basis would not, and solve() would return numbers rather than an error.

What I’d change. solve(t(X) %*% X) is the formula, not the algorithm. qr.solve(X, y) reaches the same estimate through a QR decomposition without ever forming XᵀX, so it operates at κ(X) rather than κ(X)² — roughly half the digits lost. Writing the normal equations out by hand is the right call for an assignment about understanding least squares; it’s the wrong call for one about computing it. That the two differ, and why, was the first thing this exercise taught me.

RSS ranks fit; it does not rank models

With five architectures fitted on the full dataset:

Model	Terms	k	RSS
Model 1	x4, x1², x1³, x2⁴, x1⁴	6	35.397
Model 2	x4, x1³, x3⁴	4	2.140
Model 3	x3³, x3⁴	3	463.312
Model 4	x2, x1³, x3⁴	4	20.259
Model 5	x4, x1², x1³, x3⁴	5	2.136

The first thing to notice looks like a bug: Model 1 has the most parameters and nearly the worst fit. Six parameters, RSS 35.4 — sixteen times worse than Model 2’s four parameters.

It isn’t a bug, and understanding why sharpened my intuition more than the rest of the fitting did. “More parameters never increase RSS” is true within a nested family — where the bigger model can always reproduce the smaller one by zeroing a coefficient. These five models are not nested. Model 1 contains x2⁴ and x1⁴ but has no x3⁴ term, and Model 2 has one. There is no setting of Model 1’s six coefficients that reproduces Model 2. The count of parameters is irrelevant; what matters is whether the span of the columns contains the truth.

And x3⁴ is clearly part of the truth. Every model containing it lands at RSS ≈ 2; every model without it lands between 20 and 463. Model 3, which has x3³ and x3⁴ but nothing else, is the worst of the five — so x3⁴ is necessary but nowhere near sufficient. The signal needs a term from x1 as well.

The comparison that actually decides the assignment is Model 2 against Model 5, because those two are nested: Model 5 is Model 2 plus an x1² term. It fits better, as it must — but by 0.004 in RSS, a 0.2% improvement. And its estimated coefficient on x1² is 0.00033, which is indistinguishable from zero.

Log-likelihood, and the σ² that got assumed

Turning RSS into a likelihood requires an assumption that RSS itself doesn’t make: that the errors are i.i.d. Gaussian. Under that assumption,

ln L = −(n/2)·ln(2π) − (n/2)·ln σ̂²  −  RSS / (2σ̂²)        where  σ̂² = RSS / (n − 1)

compute_loglik <- function(rss, n) {
  sigma2_hat <- rss / (n - 1)
  loglik <- -(n/2)*log(2*pi) - (n/2)*log(sigma2_hat) - (1/(2*sigma2_hat))*rss
  return(list(loglik = loglik, sigma2 = sigma2_hat))
}

Two things about that n - 1 are worth being explicit about, since it’s the specification I was given rather than the one I’d choose.

The unbiased residual variance for a regression divides by n - k, not n - 1, because fitting k coefficients already consumed k degrees of freedom. With n = 201 and k ≤ 6 the difference is under 3%, and — critically — it’s a monotone function of RSS applied identically to every model, so the AIC/BIC ranking is untouched. It’s the right shortcut for a comparison and the wrong one for reporting σ² as a quantity. Later in the pipeline, that distinction stops being cosmetic.

The Gaussian assumption is testable, and this is where the Q-Q plots earn their place in the pipeline rather than being a box to tick. For Model 2 the residuals sit on the reference line across the full range: skewness −0.04, excess kurtosis 0.15, Shapiro-Wilk W = 0.9975, p = 0.99. The likelihood is licensed. For Models 1, 3 and 4 the same plots bend away at both tails — their residuals still contain structure the model failed to absorb, which means their reported log-likelihoods are computed under an assumption their own residuals contradict.

That reframed the diagnostic for me. A Q-Q plot isn’t a check performed after selection. It’s a check on whether the quantity you’re selecting with is meaningful, and it should be read before the AIC table, not after.

AIC and BIC: the same answer, and a useful disagreement in margin

AIC = 2k − 2·ln L          BIC = k·ln(n) − 2·ln L

Model	k	RSS	ln L	AIC	BIC
Model 2	4	2.140	171.32	−334.65	−321.44
Model 5	5	2.136	171.52	−333.05	−316.53
Model 4	4	20.259	−54.59	117.18	130.39
Model 1	6	35.397	−110.67	233.34	253.16
Model 3	3	463.312	−369.14	744.27	754.18

Both criteria pick Model 2, which is the answer I want — agreement between two criteria with different penalty philosophies is much stronger evidence than either alone.

The Model 2 / Model 5 pair is where the arithmetic becomes legible. Model 5’s extra x1² term buys +0.20 of log-likelihood. AIC charges 1 unit of log-likelihood per parameter, so the trade is 0.20 earned against 1.00 charged, and Model 5 loses by 1.60. BIC charges log(201)/2 = 2.65 per parameter, so it loses by 4.11 — the same verdict, 2.6× more emphatic.

This is the parsimony penalty doing exactly the job it exists for. On RSS alone Model 5 wins every time, and would keep winning as I bolted on more powers of x1, each one soaking up a little noise. The criteria convert “does this term buy more fit than it costs in flexibility?” into a number, and the answer for x1² is no. Its coefficient of 0.00033 was already whispering that; AIC and BIC just make it a decision rather than a judgement call.

So Model 2 is:

y = 0.1436·x4  +  0.01004·x1³  −  0.001913·x3⁴  +  0.4831

Reading the coefficients requires care. x3⁴’s coefficient is 75× smaller than x4’s, which says nothing about importance — the regressors live on wildly different scales. Multiply each coefficient by the standard deviation of its own column to get comparable contributions:

Term	|θ|	sd(column)	sd(θ · column)
x4	0.1436	3.15	0.452
x1³	0.01004	144.4	1.450
x3⁴	0.001913	232.9	0.446

x1³ is the dominant driver by a factor of three, despite having the second-smallest coefficient. It also explains the output’s −1.79 skewness: a cubic term is the odd-powered, asymmetry-producing part of the model, and −0.0019·x3⁴ pulls hard in one direction only. The structure the exploratory analysis flagged is visible in the fitted parameters. Hold on to this table — it comes back to bite in the ABC section.

Validation: 95% intervals that covered 80%

Splitting 70/30 (140 train, 61 test), re-estimating on the training slice, and predicting the held-out points with 95% intervals:

pred_var <- sapply(seq_len(nrow(X_test)), function(i) {
  xi <- X_test[i, , drop = FALSE]
  sigma2_test * (1 + xi %*% XtX_inv_train %*% t(xi))
})

The (1 + xᵢ(XᵀX)⁻¹xᵢᵀ) factor is right, and worth pausing on: the 1 is the irreducible noise on a new observation and the quadratic form is the uncertainty in the fitted surface itself, which grows as you move away from the centre of the training data. Including both makes this a prediction interval rather than a confidence interval on the mean — the correct choice when the question is “where will the next observation fall.”

The training coefficients came back at (0.1404, 0.01017, −0.00188, 0.48595) against full-data (0.1436, 0.01004, −0.00191, 0.48307) — stable to three significant figures on 70% of the data, which is the real evidence that Model 2 isn’t overfitting. Validation RSS was 0.650 on 61 points, versus 2.140 on 201; per-point, 0.0107 against 0.0107. Out-of-sample error is identical to in-sample error.

And then the coverage came back at 80.3%.

That number is wrong in an interesting way. Under-coverage normally means overfitting — but the per-point RSS just ruled overfitting out. Something else was shrinking the intervals, and the culprit was the variance estimate:

sigma2_test <- rss_test / (nrow(X_train) - 1)   # 0.650 / 139

That divides the test RSS by the training sample count. It’s a units mismatch: 61 points’ worth of squared error spread across 140 degrees of freedom, giving σ² = 0.00468 where the model’s actual noise level is 0.0107 — a factor of 2.4 too small, so every interval is 1.56× too narrow. Rerunning with the denominators aligned:

σ̂² definition	σ̂²	Mean half-width	Coverage
rss_test / (n_train − 1)	0.00468	0.136	80.3%
rss_test / (n_test − k)	0.01141	0.213	93.4%
rss_train / (n_train − k)	0.01133	0.212	93.4%

93.4% against a nominal 95%, on 61 points — that’s four misses where you’d expect three, comfortably within binomial noise. The model was calibrated the whole time; the reported interval wasn’t.

This is the failure mode I’d never have found with predict(fit, interval = "prediction"), because that function’s whole job is knowing which n and which k belong in the denominator. Two symbols that both spell “sample size” got crossed, no error was raised, and the pipeline produced a plausible-looking plot of intervals that were half again too tight. The check that caught it was empirical coverage — counting how many held-out points actually landed inside — and it’s now the first thing I compute after any interval, precisely because it tests the arithmetic rather than trusting it.

Rejection ABC, and a tolerance that could never accept

The last phase estimates a posterior over two parameters without writing down a likelihood. Rejection ABC is almost embarrassingly simple:

1. Draw candidate parameters from a prior.
2. Simulate a dataset from them.
3. Keep the draw if the simulated data is within ε of the observed data.
4. The kept draws approximate the posterior.

Step 3 is the entire method, and step 3 is where I got it wrong.

The obvious tolerance is a fraction of the reference fit — 5% of Model 2’s RSS, ε = 0.107. Run 20,000 simulations against that and zero are accepted. Not “few.” Zero.

Writing the simulator by hand is what showed me why, because the noise term is right there in the line:

y_sim <- X_full %*% theta_candidate + rnorm(n, mean = 0, sd = sqrt(sigma2_ref))

Every simulated dataset gets a fresh draw of noise at all 201 points. So even with candidate parameters exactly equal to the truth, the distance Σ(y_obs − y_sim)² is a sum of squared differences between two independent noise realisations. Its expectation is roughly n·σ² = 201 × 0.0107 = 2.15:

Naive tolerance (5% of RSS):                  0.107
Expected distance from noise injection alone: 2.150   (20.1× larger)
Smallest distance seen across 20,000 draws:   3.231

The tolerance was set 20× below the floor of the distance distribution. No parameter values could have cleared it, so the acceptance rate was zero for reasons that had nothing to do with parameters at all.

The fix is to stop guessing an absolute distance and let the simulations calibrate it — run all 20,000, then set ε to the 1st percentile of the observed distances:

acceptance_quantile <- 0.01
tolerance_epsilon <- as.numeric(quantile(distances, probs = acceptance_quantile))

ε lands at 4.239 and 200 draws are accepted, by construction. This is standard practice for rejection ABC and it’s better than a guess for a structural reason: the scale of the distance metric depends on n, on σ², and on how the simulator injects noise, none of which you know before running it. A quantile is scale-free. It also makes the honest trade explicit — ε and the acceptance rate are the same knob, and you’re choosing between posterior sharpness and effective sample size, not discovering a tolerance that’s objectively correct.

The generalisable lesson: when a rejection sampler accepts nothing, check the floor of your distance metric before touching the priors. I lost time assuming my prior windows were badly placed. They were fine. The threshold was below the minimum achievable value, and one line of arithmetic — n × σ² — would have said so immediately.

What the posterior says, and what selecting by magnitude cost me

The assignment asks for the two “most influential” parameters, operationalised as the two largest by absolute value. For Model 2 that’s bias (0.4831) and x4 (0.1436), with x1³ and x3⁴ held fixed. Priors were uniform, centred on the least-squares estimates, ±50%.

The 200 accepted draws:

Parameter	LS estimate	Posterior mean	Prior sd	Posterior sd	Shrinkage
bias	0.48307	0.48291	0.1401	0.0238	5.9×
x4	0.14358	0.14328	0.0417	0.0081	5.1×

Posterior means recover the least-squares estimates to four significant figures. That agreement is the real result: a likelihood-free method that only ever simulates forward, never inverting anything, lands on the same answer as the closed-form estimator. When an analytic solution exists, ABC should reproduce it — and here it does, which is what makes it trustworthy on problems where no closed form is available.

Shrinkage of ~5-6× from prior to posterior is the data doing work; the joint posterior is a compact blob rather than the ridge that collinearity would produce, and the posterior correlation between the two parameters is −0.08, effectively independent. That last part makes sense: x4 is mean-centred (mean 0.015), so shifting the intercept barely trades off against its slope.

But the parameter selection is where I’d push back on the specification. Ranking by |θ| is ranking by coefficient magnitude, and coefficient magnitude is an artifact of the regressor’s units. From the contribution table earlier:

Term	|θ| rank	Contribution rank
bias	1	— (constant, no variance)
x4	2	2 (sd 0.452)
x1³	3	1 (sd 1.450)
x3⁴	4	3 (sd 0.446)

The |θ| criterion selected the intercept — which contributes no variance at all, being constant — and skipped x1³, the term that drives three times more variation in y than anything else in the model. The ABC posterior is over the two parameters that happen to be attached to the smallest-scale columns.

The estimates are still correct; the question they answer is just less interesting than the one I’d have chosen. Ranking by |θⱼ| · sd(xⱼ) would have put x1³ and x4 under the microscope — the two terms actually shaping the signal — and the resulting joint posterior would likely have shown the collinearity structure the exploratory correlation matrix predicted, since x1 and x4 correlate at 0.805. That’s the follow-up run worth doing, and it costs nothing but changing one order() call.

What’s worth keeping from this

• RSS ranks models only within a nested family. Across non-nested architectures, six parameters can fit sixteen times worse than four; what matters is whether the column span contains the truth, and here that meant containing x3⁴.
• AIC and BIC agreeing is worth more than either alone, and the margin between them is informative — BIC’s log(n)/2 = 2.65 per parameter versus AIC’s 1.0 turned the same verdict on x1² into a 2.6× stronger one.
• Q-Q plots belong before the AIC table, not after. They test whether the Gaussian likelihood you’re ranking with is licensed at all; three of five models here failed that test while still reporting a log-likelihood.
• Coefficient magnitude is not importance. Multiply by the regressor’s standard deviation before ranking anything, or the intercept wins a contest it isn’t in.
• Empirical coverage is the check that catches interval bugs. Nominal 95% covering 80.3% traced back to one crossed denominator — test RSS over training count — that no error message would ever have flagged.
• When a rejection sampler accepts nothing, check the floor of your distance metric first. Re-simulating noise puts an n·σ² lower bound on the distance, and no tolerance below it can accept regardless of the parameters.
• Writing least squares as solve(t(X) %*% X) teaches you the estimator and hides the conditioning; QR gets the same answer at half the digit loss. The formula and the algorithm are different things, and for polynomial bases that difference is measured in lost significant figures.