DEV Community: Isaac Land

Regression discontinuity designs break in six common ways

Isaac Land — Wed, 21 Oct 2026 11:33:56 +0000

Regression discontinuity (RD) is one of the most credible designs for causal inference from observational data. When a rule assigns treatment based on whether a score crosses a cutoff (a scholarship for test scores above 80, a programme for towns under a population threshold, a policy on one side of a border), people just above and just below the cutoff should be comparable. Any jump in outcomes at the cutoff can then be attributed to the treatment.

But a jump in a plot is not the same thing as a credible effect. This post goes through the mistakes that most often undermine RD analyses and the tools that help avoid them.

A quick refresher

The variable that determines treatment is the running variable (also called the forcing variable or score). RD comes in two forms:

Sharp RD: treatment switches on completely at the cutoff.
Fuzzy RD: crossing the cutoff only changes the probability of treatment. The effect is estimated as the jump in the outcome divided by the jump in treatment take-up, which is essentially an instrumental-variables estimate at the cutoff. It needs a non-zero first stage, an exclusion restriction, and monotonicity, and it identifies a local effect for compliers near the cutoff.

The key identifying assumption is not that the observed outcome is smooth. It is that the potential outcomes, meaning what would happen with and without treatment, change smoothly through the cutoff. Lee and Lemieux (2010) give an accessible treatment, and the design itself goes back to Thistlethwaite and Campbell (1960).

Mistake 1: fitting high-order global polynomials

A common older practice was to fit a 4th-, 5th-, or 6th-degree polynomial to the whole range of the running variable on each side and read off the gap at the cutoff. Gelman and Imbens (2019), in the Journal of Business & Economic Statistics, argue strongly against this, for three reasons:

Global polynomials put surprisingly large and erratic weight on observations far from the cutoff.
Estimates are sensitive to the polynomial order.
Confidence intervals tend to be too narrow, so real coverage falls short of the nominal level.

Better: use local linear (or local quadratic) regression within a bandwidth around the cutoff, with a triangular kernel that gives more weight to points near the threshold.

Mistake 2: choosing the bandwidth by eye (or by result)

How much data on either side of the cutoff should you use? A narrow window reduces bias but increases variance. A wide window does the opposite. Trying several bandwidths and reporting the one with the cleanest result is a garden of forking paths.

Better: use a data-driven, MSE-optimal bandwidth selector and report it up front. Then show sensitivity to, say, half and double that bandwidth as a robustness check, not as a menu to choose from.

Mistake 3: conventional confidence intervals with an optimal bandwidth

Here is a subtle one. MSE-optimal bandwidths are, by construction, wide enough that some bias remains. Conventional confidence intervals ignore that bias, so they under-cover. Calonico, Cattaneo and Titiunik (2014), in Econometrica, proposed robust bias-corrected confidence intervals that account for this. They are now standard, and implemented in the rdrobust package for R, Stata, and Python:

library(rdrobust)
out &lt;- rdrobust(y = outcome, x = score, c = 0)   # local linear, MSE-optimal bandwidth
summary(out)          # report the robust bias-corrected CI
rdplot(y = outcome, x = score, c = 0)            # binned scatter with local fits

Report the point estimate together with the robust confidence interval, not the conventional one.

Mistake 4: not checking for manipulation of the running variable

If people can precisely control their score (retaking a test until they pass, or misreporting income to stay under a threshold), those just above the cutoff are no longer comparable to those just below. A tell-tale sign is bunching: a spike in the density of the running variable on one side.

McCrary (2008) introduced the standard density test. A more recent local-polynomial version is available in the rddensity package. A histogram of the running variable with fine bins is a cheap first look that every RD paper should show.

Note the asymmetry. A clean density test does not prove the absence of sorting, but clear bunching is a serious warning.

Mistake 5: skipping covariate balance checks

If the design works, pre-treatment characteristics should not jump at the cutoff. Run the same RD estimator with baseline covariates (age, prior outcomes, demographics) as the outcome. A significant jump in a variable the treatment could not have affected suggests sorting or a second policy that shares the same cutoff.

Mistake 6: over-extrapolating the result

An RD estimate is a local effect at the cutoff. A scholarship's effect on students scoring exactly 80 says little about students scoring 50. Reporting it as "the effect of scholarships" stretches the design beyond what it identifies. Be explicit about the population the estimate applies to.

A cautionary case: the Huai River heating study

China historically provided free winter heating north of the Huai River and not south of it, which creates a geographic discontinuity. Chen, Ebenstein, Greenstone and Li (2013), in PNAS, used this border to estimate the effect of air pollution on life expectancy. The paper drew wide attention, and also methodological criticism, notably from Andrew Gelman, centring on the use of a global cubic polynomial in distance from the river and how sensitive the estimate was to that choice. A follow-up study, Ebenstein, Fan, Greenstone, He and Zhou (2017), also in PNAS, revisited the question with more data and different specifications.

Whatever your view of the substantive result, the episode is a useful teaching case. It shows why specification choices in RD (polynomial order, bandwidth, what counts as "near" the cutoff) need to be justified and stress-tested rather than picked once.

A practical RD checklist

Confirm the assignment rule exists independently of the outcome data. Don't go hunting for a kink first and call it a policy cutoff afterwards.
Plot the raw data with rdplot or a binned scatter.
Check the density of the running variable for bunching.
Check covariate balance at the cutoff.
Estimate with local linear regression, an MSE-optimal bandwidth, and robust bias-corrected CIs.
Show sensitivity to bandwidth and kernel; avoid high-order global polynomials.
For fuzzy RD, report the first stage.
Run placebo cutoffs where no effect should exist.
State clearly that the estimate is local.

I wrote a longer version of this discussion, with more on the Huai River debate, in this article on common regression discontinuity mistakes.

References

Calonico, S., Cattaneo, M. D., & Titiunik, R. (2014). Robust nonparametric confidence intervals for regression-discontinuity designs. Econometrica, 82(6), 2295–2326.
Chen, Y., Ebenstein, A., Greenstone, M., & Li, H. (2013). Evidence on the impact of sustained exposure to air pollution on life expectancy from China's Huai River policy. PNAS, 110(32), 12936–12941.
Ebenstein, A., Fan, M., Greenstone, M., He, G., & Zhou, M. (2017). New evidence on the impact of sustained exposure to air pollution on life expectancy from China's Huai River Policy. PNAS, 114(39), 10384–10389.
Gelman, A., & Imbens, G. (2019). Why high-order polynomials should not be used in regression discontinuity designs. Journal of Business & Economic Statistics, 37(3), 447–456.
Lee, D. S., & Lemieux, T. (2010). Regression discontinuity designs in economics. Journal of Economic Literature, 48(2), 281–355.
McCrary, J. (2008). Manipulation of the running variable in the regression discontinuity design: A density test. Journal of Econometrics, 142(2), 698–714.
Thistlethwaite, D. L., & Campbell, D. T. (1960). Regression-discontinuity analysis: An alternative to the ex post facto test. Journal of Educational Psychology, 51(6), 309–317.

EEG classifiers can score well by recognising the person

Isaac Land — Wed, 21 Oct 2026 02:00:22 +0000

Many EEG machine-learning papers report accuracies that look too good. One common reason has nothing to do with the brain signal you care about: the model has learned to recognise who a recording came from. If the same person's data sit in both the training and the test set, the test score answers an easier question than the one the paper claims to answer. This post walks through how that happens, what a recent small replication found, and how to set up the split so your number means what you think it means.

How windowing creates the problem

A typical EEG pipeline cuts each recording into many short windows (say, 4 seconds each), extracts features per window, and treats each window as a sample. A study with 36 participants can easily end up with a couple of thousand windows. That looks like plenty of data.

The trap is the next step. If you shuffle all the windows and do a random 80/20 split, windows from the same person land on both sides. Consecutive windows from one recording are highly similar, and so are different recordings from the same head. A flexible classifier can then learn person-specific boundaries ("for this participant, the alpha power shifts like this between rest and task") that do not carry over to anyone new.

The randomisation unit decides what you are estimating:

Window-level split: how well the model classifies new windows from people it has probably already seen.
Subject-level split: how well the model transfers to people it has never seen.

Neither is wrong in itself. The problem is reporting the first and claiming the second.

What a minimal replication showed

In October 2026, Renish Joshi published a deliberately simple replication in Towards Data Science. He used the public EEGMAT dataset on PhysioNet (EEG during mental arithmetic): 36 healthy students, each with a resting recording and a serial-subtraction recording, sampled at 500 Hz with 10/20 electrode placement (Zyma et al., 2019).

The pipeline was kept fixed:

19 scalp channels, 1–40 Hz band-pass filter, non-overlapping 4-second windows (2,134 windows in total).
Welch power spectra, relative power in theta, alpha, beta and low-gamma bands, giving 76 features per window.
Three classifiers: class-weighted logistic regression, an RBF-kernel SVM, and a 300-tree class-weighted random forest. Scaling sat inside a scikit-learn pipeline, so it was fit on training data only.

Only the evaluation changed. The window-level protocol used 10 stratified random 80/20 splits. The subject-level protocol used 6-fold GroupKFold with participant ID as the group, so no test participant contributed a single training window.

Balanced accuracy dropped for every model:

Classifier	Random windows	Held-out subjects
Logistic regression	0.770	0.648
RBF SVM	0.864	0.602
Random forest	0.664	0.532

ROC-AUC moved the same way. Because a linear model, a kernel method and a tree ensemble all showed the pattern, it is hard to blame one model family.

Do the features know who you are?

A score drop alone doesn't explain the mechanism, so Joshi ran a second test: use the same 76 spectral features to predict which of the 36 participants a window came from. Chance is about 2.8%. The identity classifier reached 95.6%. On a class-balanced subset (15 rest and 15 task windows per person, 1,080 windows), it still reached 94.0%, with 99.5% top-3 accuracy.

He is careful about what this shows. It does not prove the task classifier used identity as its decision rule. It shows the features carry a lot of person-specific information, which is enough to make overlapping splits optimistic. The window-versus-subject gap also persisted on the balanced subset, so class imbalance alone doesn't explain it.

The averages also hide a wide spread. Per-participant balanced accuracy for the subject-level logistic regression ranged from 0.289 to 0.956, and 5 of 36 participants scored below chance.

This isn't a one-off

The replication follows Brookshire et al. (2024) in Frontiers in Neuroscience, who reran two published convolutional networks. For Alzheimer's disease versus subjective cognitive decline, test accuracy was 99.8% with segment-level splits and 53.0% with subject-level splits, where the confidence interval included chance. For seizure detection on the Siena Scalp EEG database, accuracy fell from 79.1% to 65.1%. In a non-exhaustive review of 63 deep-learning EEG papers on disease identification, they found only 17 (27.0%) that clearly avoided this leakage.

Setting up the split correctly

The rule is simple to state: match the split to the generalisation claim. If you claim the model works on new people, no data from a test person may reach training, and that includes anything learned from data along the way.

In scikit-learn, that means grouped cross-validation, with every learned step inside the pipeline:

from sklearn.model_selection import GroupKFold, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression

# X: (n_windows, n_features), y: labels, groups: participant ID per window
pipe = make_pipeline(
    StandardScaler(),
    LogisticRegression(class_weight="balanced", max_iter=1000),
)

cv = GroupKFold(n_splits=6)
scores = cross_val_score(pipe, X, y, groups=groups, cv=cv,
                         scoring="balanced_accuracy")
print(scores.mean(), scores.std())

A few things to check beyond the outer split:

Tuning: if you tune hyperparameters, do it with an inner grouped split inside the training participants only (nested CV), not on the outer test folds.
Feature selection, PCA, ICA, normalisation: anything fit on data belongs inside the pipeline so it only sees training participants.
Report per-subject results: show the spread across held-out people, not just the mean.
Run the identity check: if your features can name the participant, window-level results are likely to be inflated.

Window-level evaluation is still fine when the question really is "new sessions from people already enrolled", for example a personalised model that only ever serves the users it was trained on. Just say so.

Limits worth keeping in mind

Joshi's study is one 36-person dataset with hand-crafted spectral features. The two protocols also differ in resampling structure (10 random splits versus 6 grouped folds), so the gap is descriptive, not a formal causal estimate of leakage. Where the person-specific signal comes from (physiology, electrode placement, session effects or residual artefacts) isn't resolved. And as of mid-October 2026, the article page did not link a public code repository, so the numbers can't yet be rerun independently. Brookshire et al. did publish code in their supplementary material, and EEGMAT is freely downloadable.

I first wrote up this case in Chinese, with the figures from Brookshire et al., on Scholay's research-methods section. The same logic applies to any data with repeated samples per person: patients, users, sensors, or schools.

References

Brookshire, G., Kasper, J., Blauch, N. M., Wu, Y. C., Glatt, R., Merrill, D. A., Gerrol, S., Yoder, K. J., Quirk, C., & Lucero, C. (2024). Data leakage in deep learning studies of translational EEG. Frontiers in Neuroscience, 18, 1373515. https://doi.org/10.3389/fnins.2024.1373515
Joshi, R. (2026, October 17). When EEG models see the same people in training and testing. Towards Data Science. https://towardsdatascience.com/when-eeg-models-see-the-same-people-in-training-and-testing/
Zyma, I., Seleznov, I., Popov, A., Chernykh, M., & Shpenkov, O. (2018). EEG during mental arithmetic tasks (version 1.0.0). PhysioNet. https://doi.org/10.13026/C2JQ1P
Zyma, I., Tukaev, S., Seleznov, I., Kiyono, K., Popov, A., Chernykh, M., & Shpenkov, O. (2019). Electroencephalograms during mental arithmetic task performance. Data, 4(1), 14. https://doi.org/10.3390/data4010014
scikit-learn developers. Cross-validation: evaluating estimator performance (section on cross-validation iterators for grouped data). https://scikit-learn.org/stable/modules/cross_validation.html

Cross-validation fails when data leaks between folds

Isaac Land — Tue, 20 Oct 2026 12:44:31 +0000

A model that scores well on the data it was trained on has shown that it can memorise. Cross-validation is the standard way to estimate how it will do on data it has never seen. It only works if one rule holds, though: no information from the validation data may reach the model-building process by any route. Break that rule, which is easy to do without noticing, and your estimate turns optimistic. This post covers k-fold and leave-one-out cross-validation and the leakage patterns that most often break them.

Why not just score on the training set?

The scikit-learn user guide gives the extreme case: a model that simply repeats the labels of samples it has already seen gets a perfect training score and is useless on new data. That is overfitting in its purest form.

The obvious fix is a held-out test set. With small datasets, though, holding data out shrinks the training set further, and the score depends heavily on which rows happened to land in the test split. Cross-validation reuses the data more efficiently and gives a more stable estimate.

The idea is old. Stone (1974), in the Journal of the Royal Statistical Society, Series B, used cross-validation for both choosing a prediction method and assessing it. The Elements of Statistical Learning (Hastie, Tibshirani & Friedman, chapter 7) keeps those two jobs separate:

Model selection: using CV to pick a model or tune hyperparameters.
Model assessment: using CV to report the final model's error on new data.

Mixing the two in a single CV loop is one of the leakage sources covered below.

k-fold and leave-one-out

k-fold CV randomly splits the data into k roughly equal folds. In each round, one fold is held out for validation and the model is trained on the other k − 1. Average the k validation scores and you have the CV estimate. Every observation is used for validation exactly once.

Leave-one-out CV (LOOCV) is the special case k = n: each round holds out a single observation.

Choosing k is a bias–variance trade-off. ESL notes that LOOCV is approximately unbiased for expected prediction error, but its n training sets are almost identical, so the estimate can have high variance. It can also be expensive. Five- or ten-fold CV is the common default and usually a sensible compromise.

In scikit-learn:

from sklearn.model_selection import cross_val_score, KFold
from sklearn.linear_model import LogisticRegression

cv = KFold(n_splits=5, shuffle=True, random_state=0)
scores = cross_val_score(LogisticRegression(max_iter=1000), X, y, cv=cv)
print(scores.mean(), scores.std())

Leakage pattern 1: preprocessing before splitting

This is the classic mistake. You standardise features, impute missing values, or run feature selection on the whole dataset, and only then run cross-validation. Every validation fold has now influenced the means, the imputation values, or the selected features.

Feature selection is the worst offender. ESL section 7.10.2, "The Wrong and Right Way to Do Cross-validation", describes screening thousands of predictors for those most correlated with the label and then cross-validating a classifier on the survivors. Even when labels are pure noise, this procedure can report error rates far below chance. Ambroise and McLachlan (2002) documented the same selection bias in gene-expression studies.

The fix is to put every data-dependent step inside the CV loop. In scikit-learn, that is what Pipeline is for:

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.feature_selection import SelectKBest, f_classif

pipe = make_pipeline(
    StandardScaler(),
    SelectKBest(f_classif, k=20),
    LogisticRegression(max_iter=1000),
)
scores = cross_val_score(pipe, X, y, cv=cv)  # scaler and selector refit per fold

Leakage pattern 2: tuning and reporting on the same folds

If you try 50 hyperparameter settings and report the best CV score, that score is biased upwards: you picked the winner partly because it got lucky on those folds. Varma and Simon (2006) showed how large this bias can be in small-sample settings.

The remedy is nested cross-validation. An inner loop tunes hyperparameters and an outer loop estimates the performance of the whole tuning procedure:

from sklearn.model_selection import GridSearchCV

inner = KFold(5, shuffle=True, random_state=1)
outer = KFold(5, shuffle=True, random_state=2)
search = GridSearchCV(pipe, {"logisticregression__C": [0.01, 0.1, 1, 10]}, cv=inner)
nested_scores = cross_val_score(search, X, y, cv=outer)

Leakage pattern 3: grouped data split row by row

When several rows come from the same unit (multiple scans per patient, many EEG windows per subject, several photos of the same object), a random row-level split puts the same unit in both training and validation. The model can learn to recognise the individual rather than the pattern you care about. The score then reflects identity matching, not generalisation to new people.

Split by group instead:

from sklearn.model_selection import GroupKFold

scores = cross_val_score(pipe, X, y, groups=subject_ids, cv=GroupKFold(n_splits=5))

A useful sanity check: if the group-wise score falls sharply compared with the row-wise score, your earlier number was mostly measuring leakage.

Leakage pattern 4: time

For time series, a random split lets the model train on the future and predict the past. Features built with rolling windows or "future-aware" aggregates make this worse. Use forward-chaining splits such as TimeSeriesSplit, which only ever validate on data that comes after the training window, and add a gap between them if features look back over several periods.

Leakage pattern 5: features that encode the target

Kaufman, Rosset, Perlich and Stitelman (2012) describe leakage as including information that would not be legitimately available at prediction time. Typical examples are an ID field that correlates with outcome because of how records were created, a "date of treatment" column when predicting diagnosis, or a field that is filled in only after the event you are predicting. No splitting scheme protects you from this. You have to ask, for each feature, "would I actually have this value at the moment of prediction?"

It is not just a Kaggle problem

Kapoor and Narayanan (2023), in Patterns, surveyed published research that used machine learning across many scientific fields. They found leakage to be a widespread cause of overstated results, and proposed "model info sheets" so authors could document how they had ruled it out. If you are using ML in a research setting, it is worth reading their taxonomy before submitting anything.

A short checklist

Every step that learns from data (scaling, imputation, encoding, feature selection, resampling such as SMOTE) sits inside the CV loop.
Hyperparameter tuning happens in an inner loop, and the reported score comes from an outer loop or a truly untouched test set.
If rows share a subject, site, or source, split by group.
If time matters, validate only on the future.
Each feature is available at prediction time.
The test set is used once, at the end.

For a broader look at how validation choices affect what a study can claim, I keep notes on research-method topics like this one on Scholay.

References

Ambroise, C., & McLachlan, G. J. (2002). Selection bias in gene extraction on the basis of microarray gene-expression data. PNAS, 99(10), 6562–6566.
Hastie, T., Tibshirani, R., & Friedman, J. (2009). The Elements of Statistical Learning (2nd ed.), chapter 7. Springer.
Kapoor, S., & Narayanan, A. (2023). Leakage and the reproducibility crisis in machine-learning-based science. Patterns, 4(9).
Kaufman, S., Rosset, S., Perlich, C., & Stitelman, O. (2012). Leakage in data mining: Formulation, detection, and avoidance. ACM Transactions on Knowledge Discovery from Data, 6(4).
Stone, M. (1974). Cross-validatory choice and assessment of statistical predictions. Journal of the Royal Statistical Society: Series B, 36(2), 111–147.
Varma, S., & Simon, R. (2006). Bias in error estimation when using cross-validation for model selection. BMC Bioinformatics, 7, 91.

Nested data needs a multilevel model

Isaac Land — Tue, 20 Oct 2026 12:43:16 +0000

When your data come in groups (students inside classrooms, patients inside hospitals, repeated measurements inside the same person), ordinary regression quietly assumes something false: that every row is independent. This post walks through why that matters, how a multilevel (mixed-effects) model fixes it, and how to fit one in R with lme4.

The problem: rows that know each other

Imagine 3,000 students spread across 100 classrooms. Students in the same classroom share a teacher, a timetable, and often a neighbourhood. Their test scores will be more alike than scores drawn at random from the whole population.

Ordinary least squares (OLS) does not know that. It treats 3,000 rows as 3,000 independent pieces of information. In reality, part of the information is duplicated within each classroom, so OLS standard errors come out too small, confidence intervals too narrow, and p-values too optimistic. Classroom-level predictors such as "teacher has a master's degree" are hit hardest, because there are really only 100 independent units for that variable, not 3,000.

Measuring the similarity: the intraclass correlation

The standard way to quantify "how alike are members of the same group" is the intraclass correlation (ICC). Start with an empty (intercept-only) multilevel model:

y_ij = γ00 + u_j + e_ij

u_j is classroom j's deviation from the grand mean, with variance τ²
e_ij is student i's deviation from their classroom mean, with variance σ²

Then

ICC = τ² / (τ² + σ²)

It is the share of total variance that lies between groups, and it also equals the expected correlation between two observations from the same group.

From ICC to the design effect

Survey statisticians have long used the design effect to translate the ICC into lost precision. Following Kish (1965), for groups of average size m:

deff ≈ 1 + (m − 1) × ICC

A purely hypothetical illustration: 100 classrooms of 30 students and an ICC of 0.2 give deff = 1 + 29 × 0.2 = 6.8. The effective sample size is about 3000 / 6.8 ≈ 441, and naive standard errors would be roughly 1/√6.8 ≈ 0.38 of what they should be. Even a modest ICC of 0.05 gives a design effect of 2.45 with classes of 30. Small ICCs matter when groups are large.

Three ways to handle groups

Gelman and Hill (2007) frame the choice as three options:

Complete pooling. Ignore groups and fit one regression. Simple, but it erases between-group variation (which is often the very thing you care about) and gives overconfident standard errors.
No pooling. Fit each group separately, or add a dummy variable per group. Small groups get noisy estimates, and group dummies are perfectly collinear with any group-level predictor, so you cannot estimate the effect of, say, school type at the same time.
Partial pooling (multilevel model). Each group gets its own parameters, but those parameters are drawn from a common distribution. Groups with little data are "shrunk" towards the overall mean; groups with lots of data mostly keep their own estimates.

Partial pooling is the core idea. It is a data-driven compromise between "all groups are the same" and "all groups are unrelated".

Random intercepts and random slopes

A random-intercept model lets each group have its own baseline:

y_ij = γ00 + γ10 x_ij + u0j + e_ij

A random-slope model also lets the effect of a predictor vary by group:

y_ij = γ00 + (γ10 + u1j) x_ij + u0j + e_ij

For example, the link between hours of homework and test scores might be steeper in some classrooms than in others. The variance of u1j tells you how much that slope varies, and the covariance between u0j and u1j tells you whether high-baseline groups also have steeper slopes.

Fitting it in R with lme4

The lme4 package (Bates, Mächler, Bolker & Walker, 2015) is the workhorse for this in R. It ships with the sleepstudy dataset: reaction times for 18 participants measured over 10 days of sleep restriction, so measurements are nested in people.

library(lme4)

# Complete pooling: one line for everyone
m_ols &lt;- lm(Reaction ~ Days, data = sleepstudy)

# Random intercept per participant
m_ri &lt;- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)

# Random intercept and random slope for Days
m_rs &lt;- lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)

summary(m_rs)

Things worth comparing yourself:

The standard error of Days in m_ols versus m_rs. The OLS version treats 180 rows as independent.
The Subject variance components in summary(m_ri). Plug them into the ICC formula above.
coef(m_rs)$Subject against per-subject lm() fits. You will see shrinkage: participants with extreme individual slopes are pulled towards the average.

lmer() uses REML (restricted maximum likelihood) by default, which gives less biased variance estimates than plain ML. If you compare models that differ in their fixed effects with a likelihood-ratio test, refit with REML = FALSE first. anova() in lme4 does this automatically for merMod objects.

Common mistakes

Too few groups. Variance components are estimated from the number of groups, not the number of rows. Maas and Hox (2005) found that with very few groups the standard errors of the group-level variances were noticeably biased. Five schools with 1,000 students each is a weak basis for a random effect.

Ignoring the design and just clustering standard errors. Cluster-robust standard errors can fix inference for fixed effects, which is a reasonable choice in some settings. But they do not give you estimates of between-group variation, group-specific predictions, or random slopes. Pick the tool that answers your question.

Adding random slopes for everything. Models with many random slopes and correlations often fail to converge or produce singular fits (isSingular() returns TRUE). Start from a theoretically motivated structure and simplify when the data cannot support it.

Forgetting centring. With a random slope on x, the intercept variance depends on where x = 0. Centring predictors, either at the grand mean or within each group, changes what the intercept means. Within-group centring also separates within-group and between-group effects, which are often different things entirely.

Treating the ICC as a test of whether you "need" a multilevel model. A small ICC can still produce a large design effect when groups are big, and the grouping is part of the design whatever the ICC turns out to be.

When to reach for a multilevel model

Use one when:

observations are grouped, and the grouping comes from how the data were collected (sampling, randomisation, or repeated measurement);
you care about variation between groups, or about predictors measured at the group level;
some groups are small and you want stable estimates for them.

If you want a longer, education-focused walk-through of the same ideas, including a worked example with classroom data, the original article on nested data and multilevel models goes into more detail.

References

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48.
Gelman, A., & Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.
Kish, L. (1965). Survey Sampling. Wiley.
Maas, C. J. M., & Hox, J. J. (2005). Sufficient sample sizes for multilevel modeling. Methodology, 1(3), 86–92.
Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical Linear Models: Applications and Data Analysis Methods (2nd ed.). Sage.