Statistics for Neuroscience
A working tour of the statistics you will actually use: describing a distribution, putting an honest error bar on it, choosing a test that respects your data, and using resampling when the textbook formulas run out.
Distributions
Here are four samples of “firing rate” measurements. For each, would you summarize it with the mean or the median — and which choice would mislead a reader?
From left to right these are Gaussian, right-skewed, bimodal, and long-tailed. The mean and median agree only for the first. On the skewed sample the mean is dragged toward the long tail; on the bimodal sample the mean lands in the empty valley between two clusters, describing no real cell. The fourth looks like a narrow spike with empty space to the right — its long tail is invisible on a linear axis. The next figure uses a log scale to reveal what is hiding there, and introduces the closely related log-normal distribution.
The two summaries are defined as
\[\bar{x} = \frac{1}{n}\sum_{i=1}^{n} x_i \qquad\qquad \text{median} = \text{middle value of the sorted } x_i.\]
The mean uses every value (so outliers move it); the median only cares about rank (so outliers do not).
When you need a log axis
The two distributions on the left each look like a clump near zero with “nothing” to the right. Is there really nothing out there — and what happens if we plot the same data on a log axis?
Top — a log-normal distribution is strongly right-skewed on a linear axis (left) but becomes a symmetric bell once the x-axis is logarithmic (right). That is exactly what “log-normal” means: the logarithm* of the values is normally distributed (firing rates, synaptic weights, and many biological quantities are log-normal). Bottom — a long-tailed distribution looks like a single spike on a linear axis, but the log axis reveals a tail spanning more than two orders of magnitude. Whenever values are positive and span orders of magnitude, plot them on a log axis — and label it clearly.*
The Sampling Distribution
Almost everything that follows — error bars, confidence intervals, p-values, the bootstrap — is one idea wearing different hats. This section is that idea.
You record 20 neurons and compute their mean firing rate. Tomorrow you run the identical experiment on 20 different neurons from the same population. Write down your prediction: will you get the same number? If not, how far off?
Left is the truth you never observe: every neuron in the population. Right is what you actually get — 20 draws from it, and one sample mean that lands near the true mean but not on it. That gap is not a mistake. It is sampling variability, and quantifying it is what statistics is for.
Repeating the experiment
Suppose you could run the experiment hundreds of times. What would the collection of sample means look like?
Population SD (σ): 6.00 Hz
Observed SD of the 5000 sample means: 1.36 Hz
Formula prediction σ/√n = 6.0/√20 = 1.34 HzThis histogram of sample means is the sampling distribution. Three things to notice: it is centered on the true mean, it is far narrower than the raw data, and its standard deviation matches \(\sigma/\sqrt{n}\) almost exactly. That is what the SEM is — not a property of your neurons, but the width of this distribution: how much your estimate of the mean would bounce around if you repeated the experiment.
Why \(\sqrt{n}\)
More neurons per experiment ⇒ a tighter sampling distribution ⇒ a more precise estimate. But the width falls as \(\sqrt{n}\), not \(n\): quadrupling your sample size only halves your error bar. This is the arithmetic behind every “just add more cells” conversation — going from 20 to 80 neurons buys you a factor of two, and no more.
The Central Limit Theorem
Everything above used a Gaussian population. Real data usually is not Gaussian. Does any of this survive?
The population (left) is emphatically not normal. Yet the distribution of its sample means (right) is very nearly a bell curve. This is the Central Limit Theorem: averages tend toward normality even when the underlying data do not, and faster as \(n\) grows. It is why t-tests are more forgiving than their assumptions suggest — and why that forgiveness runs out when \(n\) is small or the skew is severe, which is exactly when we reach for nonparametric tests later.
The catch: you only ever run the experiment once. You never see the sampling distribution. So how can you know its width?
Two answers, and the rest of this page is built on them:
- Assume a form and use a formula. If the sampling distribution is normal, its width is \(s/\sqrt{n}\) — the SEM — estimated from your single sample. Confidence intervals and t-tests follow directly.
- Simulate it from your own data. Treat your sample as a stand-in for the population and resample it thousands of times. That is the bootstrap, and it works for statistics that have no formula at all.
Everything below is one of these two moves.
Error Bars and Effect Size
Error bars: SD vs. SEM vs. CI
These three bars answer three different questions and are constantly confused. With the sampling distribution in hand, the difference is easy to state:
- Standard Deviation (SD) — the spread of your data. A property of the neurons. Does not shrink as you collect more. \[s = \sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2}\]
- Standard Error of the Mean (SEM) — the width of the sampling distribution you just saw. A property of your estimate, not your neurons. Shrinks as \(\sqrt{n}\). \[\text{SEM} = \frac{s}{\sqrt{n}}\]
- 95% Confidence Interval (CI) — built by walking out from your estimate far enough to cover 95% of that sampling distribution. \[\bar{x} \pm t_{0.975,\,n-1}\cdot \text{SEM}\]
Same 30 neurons, three error bars. Which one looks the most impressive — and is that a good reason to plot it?
SEM is the smallest bar — but only because it is the SD divided by \(\sqrt{30} \approx 5.5\), not because the data are tight. Nothing about the neurons changed between the three bars. Use SD to describe your sample, CI to express uncertainty about the mean, and always say in the caption which one you plotted.
What a confidence interval actually means
A statement you will hear constantly: “the 95% CI is [18.2, 22.4], so there’s a 95% chance the true mean lies in that range.” Is that correct? Commit to yes or no before scrolling.
Studies run: 100
CIs containing the true mean: 97
CIs that miss it entirely (red): 3Each horizontal line is one lab running the same experiment and reporting its 95% CI. The vertical line is the truth. Close to 95 of the 100 intervals catch it — the exact count wobbles from run to run, because coverage is itself a random quantity — and the red ones miss completely. Those labs have no way of knowing they are the red ones.
So the statement above is wrong, subtly. The 95% describes the procedure, not any particular interval: intervals built this way capture the truth 95% of the time. Your one interval either contains the true mean or it does not — the probability is 0 or 1, you just don’t know which. The honest reading of a CI is “values in this range are compatible with my data,” and its width is the useful part: a wide CI is a confession that the experiment could not pin the answer down.
Effect size
A small p-value says an effect is detectable, not large. The standardized effect size, Cohen’s d, reports magnitude independent of sample size:
\[d = \frac{\bar{x}_1 - \bar{x}_2}{s_\text{pooled}}, \qquad s_\text{pooled} = \sqrt{\frac{(n_1-1)s_1^2 + (n_2-1)s_2^2}{n_1+n_2-2}}.\]
Conventionally \(d \approx 0.2\) is small, \(0.5\) medium, \(0.8\) large.
Both panels show the same effect size. Only \(n\) differs. Why does the p-value collapse?
Cohen’s d is identical by construction in both panels — the groups overlap exactly as much, and a randomly chosen neuron from B beats one from A just as often. Only \(n\) differs. Yet the p-value falls from “not significant” to vanishingly small. “Significant” answers can we detect it?; effect size answers does it matter?* Report both, always — a p-value alone tells you as much about your sample size as about your biology.*
p-values
Every test in the next section ends in a p-value. Before running any of them, it is worth being precise about what that number is.
What a p-value is
You compare two groups and get p = 0.03. In one sentence, what does that number mean? Write it down — most published descriptions of it are wrong.
The definition has a specific and narrow shape:
A p-value is the probability of observing data at least as extreme as yours, if the null hypothesis were true.
It is a statement about data given a hypothesis, never about the hypothesis given the data. We can build one directly: simulate a world where the null is true (both groups drawn from the same population), collect the differences that arise from noise alone, and see where the real result falls.
Observed difference: 2.29 Hz
Fraction of null world at least
this extreme (red tails) = p = 0.1246
t-test p for comparison: 0.1290The grey histogram is a world with no real effect — every difference in it is pure noise. The red tails are the outcomes at least as extreme as yours. The p-value is simply their share of the total. Note what this does not tell you: nothing here is a probability that the null is true, and nothing is a measure of how big the effect is.
Three things a p-value is not
| Claim | Why it’s wrong |
|---|---|
| “p = 0.03, so there’s a 3% chance the null is true.” | Backwards. p is P(data | null), not P(null | data). Getting the second requires a prior. |
| “p = 0.03, so there’s a 97% chance it replicates.” | Unrelated. Replication probability depends on true effect size and power, and is usually far lower than people assume. |
| “p = 0.001 means a bigger effect than p = 0.04.” | No. p mixes effect size with sample size. A huge n makes a trivial effect tiny in p — see the effect-size figure above. |
The most useful fact about p-values
No real effect: 5.2% of experiments give p < 0.05
Real effect present: 56.1% of experiments give p < 0.05 (this is POWER)When there is no effect, p-values are uniformly distributed — every value from 0 to 1 is equally likely, so p < 0.05 happens exactly 5% of the time. That flat histogram is the entire reason α = 0.05 means what it means, and the reason 100 tests on noise yield about 5 “discoveries.” When a real effect exists, the distribution piles up near zero, and the fraction below 0.05 is the test’s power.
Choosing a Test
Three questions decide your test: Is the data roughly normal? Are the measurements paired? Are you running many tests at once?
Start with the distribution: is it normal?
A parametric test (t-test, Pearson correlation) assumes the data come from a known distribution family — almost always the normal, defined by two parameters, mean \(\mu\) and SD \(\sigma\). If that assumption is wrong, the test’s p-values are wrong.
Before running a t-test: how would you check whether these data are normal? What tools do you have?
Shapiro–Wilk: W = 0.897, p = 0.0001
D'Agostino–Pearson: K2 = 11.85, p = 0.0027
Kolmogorov–Smirnov: D = 0.163, p = 0.0730
Anderson–Darling: A2 = 2.03 (5% crit = 0.74)The histogram and QQ plot are your first line of defense — the QQ points bend away from the line whenever data are skewed. Formal tests put a number on it:
| Test | Best at | Caution |
|---|---|---|
| Shapiro–Wilk | General-purpose, strong for small n | The usual default |
| D’Agostino–Pearson | Detecting skew & kurtosis | Needs n ≳ 20 |
| Kolmogorov–Smirnov | Comparing to any reference distribution | Weak unless parameters are known; low power |
| Anderson–Darling | Sensitive to the tails | Compare statistic to critical value, not a p |
A caveat: with large n these tests flag trivial, harmless deviations; with small n they miss real ones. Use them to support what the QQ plot already shows, not to replace your eyes.
Paired vs. unpaired
Data are paired when the same unit (neuron, animal, subject) is measured twice. Pairing removes between-unit variability and dramatically increases power.
The two clouds overlap almost completely. Is there no effect here?
Every neuron went up — the change is small but perfectly consistent. The big between-neuron baseline spread is what makes the unpaired view look uninformative. The paired test removes that spread and finds the effect easily. If the design is within-subject, the test must be paired.
Parametric vs. nonparametric
When normality fails (small n, skew, outliers), use a rank-based test: Mann–Whitney U (unpaired), Wilcoxon signed-rank (paired), Spearman (correlation).
The two tests below give different p-values on the same data. Which do you trust?
With skewed data and n < ~15, the nonparametric p-value is the trustworthy one. Rank tests cost a little power when data really are normal, but protect you when they are not — a good default for small neuroscience samples.
Comparing more than two groups
Experiments rarely stop at two conditions — three drug doses, four genotypes, five cortical layers.
With three groups, why not just run three t-tests (A vs B, A vs C, B vs C) and report whichever is significant?
Three groups drawn from the SAME population — no real differences:
'Any of 3 pairwise t-tests significant': 12.9%
One-way ANOVA significant: 5.2%Three pairwise tests give you three chances to be wrong, so the false-positive rate climbs to roughly 1 in 8. The ANOVA asks one question — “is there any difference among these groups?” — and holds the error rate at 5%.
What the F statistic is
ANOVA compares two estimates of variance:
\[F = \frac{\text{variance between group means}}{\text{variance within groups}}\]
If the groups are truly identical, both numerator and denominator estimate the same noise and \(F \approx 1\). A real difference inflates the numerator. Crucially, the group means alone cannot tell you — the same means can be convincing or meaningless depending on the scatter around them.
Both panels have the same three group means (18, 22, 26 Hz). On the left the within-group scatter is small, so the separation between means is large relative to the noise — F is big and the result is clear. On the right the identical means are buried in scatter, F drops toward 1, and the same biological difference is undetectable. F is a signal-to-noise ratio, which is why reporting group means without their spread tells a reader nothing.
After a significant ANOVA
A significant ANOVA says something differs — not which groups differ. That question needs a follow-up:
| Situation | Use |
|---|---|
| All pairwise comparisons after ANOVA | Tukey’s HSD (corrects for all pairs) |
| Comparing every group to one control | Dunnett’s test (fewer comparisons, more power) |
| A few specific comparisons planned in advance | Planned contrasts, corrected as in the next section |
| Non-normal data / small n | Kruskal–Wallis (the rank-based analog), then Dunn’s test |
| Repeated measures on the same animals | Repeated-measures ANOVA or a mixed-effects model |
One-way ANOVA: F(2, 39) = 21.17, p = 5.9e-07
Kruskal–Wallis: H = 22.27, p = 1.5e-05
Post-hoc pairwise (uncorrected -> Bonferroni x3):
Low vs Mid p = 0.0020 -> 0.0059
Low vs High p = 0.0000 -> 0.0000
Mid vs High p = 0.0067 -> 0.0200The ANOVA is significant, and the post-hoc comparisons localize it. Note the Bonferroni column: multiplying each p by the number of comparisons is the simplest honest correction, and it is exactly the idea the next section generalizes.
Multiple comparisons
Every test at \(\alpha = 0.05\) has a 5% false-positive rate. Run \(m\) of them and false positives accumulate.
- Bonferroni controls the chance of any false positive (family-wise error). Test each at \[\alpha_\text{corrected} = \frac{\alpha}{m}.\]
- Benjamini–Hochberg (FDR) controls the expected proportion of false positives among your hits. Sort the p-values \(p_{(1)} \le \dots \le p_{(m)}\) and find the largest \(k\) with \[p_{(k)} \le \frac{k}{m}\,\alpha;\] reject all hypotheses up to that \(k\). Less conservative, more power.
These 100 tests were run on pure noise — there is no real effect anywhere. How many came out “significant,” and what should we do about it?
Uncorrected (p<0.05): 9 'significant' — all false positives
Bonferroni: 0
Benjamini–Hochberg: 0About five tests cross the uncorrected line by chance alone — exactly the 5% you signed up for. Bonferroni and FDR both clean this up. Use Bonferroni when any single false positive is costly; use FDR when you can tolerate a known fraction and want more power (typical for screens across many ROIs, channels, or cells).
Correlation and Regression
So far every question compared groups. Just as often you want to relate two continuous measurements — firing rate against contrast, spine density against age, performance against dose.
What a correlation actually looks like
Below are four relationships. Write down your estimate of \(r\) for each before scrolling. Most people badly overestimate the weak ones.
Calibrate your eye on these. An \(r\) of 0.3 — which is common in published neuroscience and often described as “correlated” — is a cloud with barely visible structure. Even \(r = 0.6\) leaves enormous scatter. The reason is \(r^2\): at \(r = 0.3\) the predictor accounts for 9% of the variance, leaving 91% unexplained.
Fitting a line
Regression asks a directional question that correlation does not: given \(x\), what is my best guess for \(y\)? It answers with a line, \(\hat{y} = \beta_0 + \beta_1 x\), chosen to minimize the squared vertical distances (residuals) from the points to the line.
Left: the red segments are the residuals the fit minimizes. Two numbers describe it, and they answer different questions — the slope is biological (“each 1% of contrast adds ~0.18 Hz”, and it carries units), while \(r\) is a unitless measure of scatter about the line. You can have a steep slope with weak \(r\), or a shallow slope with near-perfect \(r\); never use one as a stand-in for the other.
Right: always plot residuals against fitted values. It should look like a structureless band, as it does here. Any curve means you fitted a line to something that bends; any funnel shape means the variance grows with the mean, and the p-value on that slope is not trustworthy.
Regression to the mean
You measure every neuron on two days with no treatment whatsoever. You take the 10% most active neurons from day 1 and re-examine them on day 2. What happens to their average — and why?
Population mean (all cells, both days): 20.3 Hz
Selected top 10% — day 1 mean: 31.3 Hz
Same cells — day 2 mean: 25.5 Hz
Apparent 'decline': -5.8 Hz
Paired t-test on that decline: p = 1.1e-09The selected cells dropped by several Hz, and a paired t-test calls it highly significant — yet nothing was done to them and their underlying rates never changed. Cells land in the top 10% on day 1 partly because they are genuinely active and partly because that day’s noise happened to favor them. The noise does not repeat, so on day 2 they fall back. This is regression to the mean, and it is pure arithmetic, not biology.
It is the reason any design that selects on an extreme and then re-measures is booby-trapped: responders to a treatment, the most impaired animals, the strongest-tuned cells, the worst-performing sessions. All of them will “improve” on their own. The defenses are a control group selected the same way (both regress; only the difference between them is real), or selecting on one dataset and measuring on independent data.
Resampling
When a quantity has no clean formula for its error bar, or the data break parametric assumptions, let the data generate the answer by resampling.
Bootstrap
The measurement. You record one visual-cortex neuron while showing drifting gratings at eight directions, 12 trials each, and count spikes on every trial. Averaging across trials gives the direction-tuning curve below. What is this neuron telling us?
Preferred direction: 90°, R_pref = 10.8 Hz
Null direction: 270°, R_null = 5.6 Hz
Direction-selectivity index (DSI) = 0.32This tuning curve is the raw measurement: each point is the mean firing rate at one direction over 12 trials, with SEM error bars. The neuron clearly fires more to its preferred direction than to the opposite (“null”) direction, and we compress that selectivity into a single number — the direction-selectivity index \(\text{DSI} = (R_\text{pref}-R_\text{null})/(R_\text{pref}+R_\text{null})\). But each point rests on only 12 noisy trials, and the DSI is a ratio with no simple formula for its error bar.
The DSI came out positive. With only 12 trials, is it reliably greater than zero, or could noise alone produce it?
Resample the 12 trials with replacement, recompute the DSI each time, and repeat thousands of times. The fraction of resampled DSIs above zero is your confidence that the true value is positive.
Observed DSI: 0.317
P(DSI > 0): 100.0%
95% CI: [0.216, 0.405]Each resample is a “parallel universe” version of your 12-trial experiment, drawn from your own data. Collecting thousands of resampled DSIs builds the sampling distribution; the fraction above zero estimates your confidence that the true DSI is positive, and the 2.5th–97.5th percentiles give a 95% CI. If that CI comfortably excludes zero, the selectivity is reliable; if it straddles zero, 12 trials simply were not enough. The bootstrap delivers this for any metric, with no formula required — its great appeal for messy neuroscience quantities.
Jackknife
The measurement. Across 20 recording sessions you measure how much a drug changes each session’s mean firing rate. Every session contributes one number, and your headline result is the average change across sessions.
This is the raw measurement: one drug-induced change per session, with the reported population mean drawn through them. Most sessions cluster near +6 Hz, but one (red) sits far above the rest. Is the headline mean a genuine population effect, or is a single session dragging it upward?
How much does any single session move the mean? The jackknife recomputes the mean while leaving out one session at a time.
The jackknife is older and simpler than the bootstrap — still useful as a fast influence check, and common in legacy literature.
Watch the right panel: most leave-one-out means cluster tightly, but dropping the single high point (the red, circled observation) visibly pulls the mean down. That lone outlier is driving the estimate — exactly the kind of influence the jackknife is built to expose.
Permutation (Monte Carlo) test
The measurement. In a behavioral task you record response latency on 16 control trials and 16 stimulation trials. Latencies are skewed and the samples are small, so you don’t trust a t-test. Here are the raw data.
The raw measurement: 16 latencies per condition, with the group means marked. Stim trials look slower on average (the observed mean difference is annotated on the next figure), but the data are skewed and n is small. Is that gap larger than we’d expect if the condition labels meant nothing?
Is the observed difference more than chance — without assuming any distribution?
The test statistic is the difference between the two group means — not a comparison of individual samples. The logic: pool all 32 latencies, randomly reassign them into a fake “control” and “stim” group of the original sizes, and recompute the difference of means. Repeat thousands of times to build the null distribution of “no difference,” then see how extreme the real difference is against it.
The animation above used 60 shuffles to show the mechanism. The figure below runs 10,000 of them for a stable p-value.
The shuffled null says “if the condition labels were meaningless, here’s how big a difference between the group means we’d see by chance.” Our observed difference of means sits far in the tail, so it is unlikely under chance. The second panel shows the estimate stabilizing — a few thousand permutations is plenty.
A caution: resampling does not fix pseudoreplication
You recorded 30 cells from each of 3 mice — 90 cells. Is your sample size 90?
Resampling assumes independent observations, but 90 cells from 3 mice are nested, not independent — bootstrapping the cells gives a falsely tight CI. For questions at the animal level, your effective n is 3. Analyze at the level of the independent unit (animal means) or use a hierarchical / mixed-effects model that respects the nesting. More cells per mouse does not buy more independent samples.