Bundle review table

Eutrophication is a switch? A test of the posted phosphorus-threshold claim

README.md

# Eutrophication is a switch? A test of the posted phosphorus-threshold claim

**Hypothesis:** "Eutrophication Is a Switch" — Original Theories, Volume XVI #10
**Posted:** <https://incoprea.com/research-library/original-theories-vol16-20261030/>
**Test:** `tests/05-lake-switch`, run 30 October 2026
**Data:** EPA National Lakes Assessment 2007 and 2017 (1,146 and 993 target lakes)
**Status:** working result, not peer reviewed, AI-assisted

## The result

A sharp change point is statistically demonstrable on the full 2017 sample — a
continuous two-segment fit beats a straight line by ΔBIC 63.6 and a quadratic by
ΔBIC 14.5, with a break at 145 µg/L — but it is carried by the top 5% of the
phosphorus distribution: removing the 50 lakes above 448 µg/L removes the break
by the pre-specified BIC rule (ΔBIC −0.79 against a line, 0.18 against a
quadratic), and where it appears the break is a ceiling (slope 1.03 below, 0.11
above), not an accelerating switch. The proposal's threshold is therefore not
established as a general property of the gradient; the two-state signature is
weak (residual valley 0.616 against a null 0.117) while the persistence
signature is strong (2007→2017 chlorophyll coefficient 0.502 across 216 lakes).

The proposal's directional picture — a lake that is clear below the threshold
and then crosses a boundary into algal dominance above it — is the opposite of
what the full sample shows. Chlorophyll rises in proportion to phosphorus
across most of the gradient (slope 1.03 on log–log axes) and then stops
responding (slope 0.11). That is a saturation curve, not an accelerating switch,
and the flat upper segment is a property of a few dozen hypereutrophic lakes.

## What was tested

The claim has three testable implications, and all were fixed in
`hypothesis.md` before any model was fitted:

- **P1 — a sharp break.** log10 chlorophyll-a against log10 total phosphorus
  should be better described by a continuous two-segment change point (M1) than
  by a straight line (M0) or a quadratic (M2).
- **P2 — two states.** Residuals from a smooth fit should be bimodal, not
  merely skewed.
- **P3 — path dependence.** Among lakes sampled in both 2007 and 2017, 2017
  chlorophyll should depend on 2007 chlorophyll after conditioning on 2017
  phosphorus.

The pre-specified support rule for P1 was: ΔBIC(M1 vs M0) > 6 **and** bootstrap
likelihood-ratio p < 0.05 **and** ΔBIC(M1 vs M2) > 0. All three were required;
none was added after the fact.

## Data and sample

Two EPA National Lakes Assessment surveys, one index visit per lake
(`VISIT_NO == 1`), target-sampled lakes only (2007 `LAKE_SAMP ==
"Target_Sampled"`; 2017 `SITETYPE == "PROB"`). Total phosphorus and chlorophyll-a
are in µg/L and both are logged. Full provenance and row counts are in
`SOURCE.md`; the working CSVs are byte-identical members of the official EPA
archives.

| | 2007 | 2017 |
|---|---|---|
| Target lakes with TP > 0 and CHLA > 0 | 1,146 | 993 |
| TP range (µg/L) | 1–4,679 | 1.36–11,238.8 |
| CHLA range (µg/L) | 0.067–936 | 0.1–3,299.18 |
| TP values flagged below reporting limit | 116 | 25 |
| Lakes with maximum depth | 1,145 | 988 |

Below-limit values are kept at their reported numeric value in the primary
analysis, because both surveys report a number rather than only a "less-than"
token and dropping them would preferentially delete the low-phosphorus end. A
sensitivity analysis drops them instead (S4 below). The 2017 sample is the
primary sample; 2007 is a replication in a different survey year, though not an
independent one — the 2017 revisit panel re-sampled 2007 sites, and 216 lakes
match one-to-one in the path-dependence sample. The two surveys are not pooled
in the decisive test.

## Result 1 — a tail-driven ceiling

All models are ordinary least squares on log10 scales, n = 993.

| Model | Parameters | BIC | RSS |
|---|---|---|---|
| M0 — linear | 2 | −1609.05 | 193.72 |
| M1 — continuous two-segment, break at τ | 4 | −1672.62 | 179.20 |
| M2 — quadratic | 3 | −1658.15 | 183.10 |

ΔBIC(M1 vs M0) = 63.56 and ΔBIC(M1 vs M2) = 14.47, both favouring the break.
The bootstrap likelihood-ratio statistic is 77.36 with p = 0.0005 (2,000
parametric draws; 0.0005 is the smallest value 2,000 draws can resolve). The
break is at **τ = 144.8 µg/L**, with a 95% percentile bootstrap interval of
**69.8–274.2 µg/L** from 10,000 paired resamples; 10.1% of resamples put τ on a
grid edge, so the location is bounded rather than pinned down.

By the pre-specified rule P1 is supported on the full sample: the relation is
not a single gradual curve. The support is for the existence of a break, not
for the direction the proposal predicts — and, as the next subsection shows, not
for the break as a general property of the gradient.

The shape is the important part. Below the break the slope is **1.026**; above
it the slope is **0.105**. The slope change is **−0.921** (95% CI −1.674 to
−0.561). A lake at 50 µg/L has predicted chlorophyll of about 14 µg/L; at
145 µg/L it is about 41 µg/L; by 1,000 µg/L it is only about 50 µg/L. The
break marks where chlorophyll stops tracking phosphorus, not where it starts
to accelerate.

The 2007 survey gives the same qualitative answer: n = 1,146, ΔBIC(M1 vs M0) =
89.10, ΔBIC(M1 vs M2) = 53.39, LRT 103.19, p = 0.0005, τ = 237.9 µg/L (95% CI
95.9–255.0), slopes 0.826 → −0.392 (change −1.218). Two qualifications: 30.5%
of the 2007 bootstrap resamples put τ at a grid edge, because the 2007
break sits at the top of the search window, so only the existence of a break —
not its location — replicates; and the 2007 upper-segment slope is actually
negative, i.e. chlorophyll is flat or declining at the highest phosphorus
concentrations.

Spearman correlations on log–log axes are 0.749 in 2017 and 0.762 in 2007;
M0 explains R² = 0.512 of the 2017 variance with a slope of 0.808. The break is
a departure from a single power law on the full sample, but the departure is
confined to the hypereutrophic tail.

### The break disappears when the top 5% of TP is removed

The 95th percentile of 2017 TP is 447.9 µg/L. Removing the 50 lakes above it
leaves 943. On that sample ΔBIC(M1 vs M0) = −0.79 and ΔBIC(M1 vs M2) = 0.18: the
break no longer earns its two extra parameters. The likelihood-ratio statistic
is still nominally significant (p = 0.001), which is why the two criteria can
disagree — the break improves the fit slightly but not enough to pay the BIC
penalty. The same trim on 2007 leaves 1,088 lakes and a marginal break (ΔBIC
6.45 and 13.37, τ = 101.4 µg/L, slopes 0.869 → 0.316); the 2007 break survives
the BIC bar, but it is a slope reduction, not a ceiling. On the 2017 deep
stratum, trimming removes the break (n = 613, ΔBIC 3.58 and 6.26). These trims
were added after the full-sample fit and are labelled post-hoc; the
pre-specified full-sample result is unchanged, but its interpretation is not:
the break is a property of the hypereutrophic tail, not of the main gradient.

## Result 2 — conditional bimodality, but a small second mode

Residuals from M2 have a Gaussian-KDE valley depth of **0.616** on the
`dip_like` statistic. Against 2,000 Gaussian residual resamples the null median
is 0.0 and the 95th percentile is 0.117; the raw p is 0.001 and the
Benjamini–Hochberg adjusted p is 0.002 across the two bimodality tests. That
passes the pre-specified rule (adjusted p < 0.05 and depth ≥ 0.10).

The second mode is small. Twenty-seven lakes (2.7%) have M2 residuals below
−1 log10 unit, and the unconditional log-chlorophyll distribution is not
bimodal: its valley depth is 0.026 with p = 0.106. So the evidence is for a
long left tail of unusually clear lakes conditional on phosphorus, not for two
comparable algal states. Two post-hoc checks, not part of the pre-specified
decision, probe this: residuals from a 20-bin flexible smoother still give a
valley depth of 0.373 against a null 95th percentile of 0.091 (p = 0.0055), and
dropping the two most extreme M2 residuals reduces the depth from 0.616 to
0.322. The statistical valley survives, but at roughly half its original depth,
and its substantive interpretation as a distinct clear state is weak.

`dip_like` is a descriptive Gaussian-KDE valley-depth statistic, **not
Hartigan's dip test**, and should be read as such.

## Result 3 — path dependence is present, but it is not a test of reversal

The surveys do resample the same lakes: matching on the NHD COMID with one site
per survey per waterbody gives **216 matched lakes** (median relative area
difference 10⁻⁶, confirming the same waterbody). This is above the 50-lake
threshold, so the clause is testable rather than marked UNTESTED.

The pre-specified model is

  log10 CHLA₂₀₁₇ = b₀ + b₁ log10 TP₂₀₁₇ + b₂ log10 CHLA₂₀₀₇ + ε.

The estimate is **b₂ = 0.502** (se 0.054; t = 9.23; 95% paired-bootstrap CI
0.358–0.640; one-sided log10 p = −16.87; 10,000-label permutation p = 0.0001).
The model R² is 0.592, b₁ = 0.332, the partial correlation of 2007 and 2017
chlorophyll given 2017 phosphorus is 0.535, and the simple Spearman correlation
between the two chlorophyll measurements is 0.745. Adding 2007 phosphorus
leaves the result essentially unchanged (b₂ = 0.539, b₃ = −0.122).

This is persistence conditional on current phosphorus. It is consistent with
the proposal's "not reversed" clause, but it is equally consistent with fixed
lake-level differences — colour, morphometry, nitrogen, grazing — that make a
lake consistently clearer or greener than its phosphorus alone would predict.
No lake in the data is observed returning to its earlier phosphorus
concentration, so the strong form of the clause is not directly tested.

## Depth and alternative definitions

The depth split was fixed at 3 m before fitting. The proposal's internal-loading
mechanism implies the break should matter more in shallow, polymictic lakes; the
data do not follow that pattern.

| Survey | Stratum | n | ΔBIC(M1 vs M0) | ΔBIC(M1 vs M2) | τ (µg/L) | P1 verdict |
|---|---|---|---|---|---|---|
| 2017 | shallow (< 3 m) | 342 | 34.99 | −3.00 | 176 | not supported |
| 2017 | deep (≥ 3 m) | 646 | 34.13 | 14.65 | 71.6 | supported |
| 2007 | shallow (< 3 m) | 288 | 10.92 | −0.44 | 246 | not supported |
| 2007 | deep (≥ 3 m) | 857 | 100.21 | 77.49 | 109 | supported |

In shallow lakes the smooth quadratic is at least as good as the break (it wins
by 3.00 BIC in 2017 and by 0.44 in 2007), so the pre-specified rule does not
support a threshold there. On the full samples the break is a deep-lake feature
in both surveys, but it is itself tail-dependent: after trimming the top 5% of
TP within the 2017 deep stratum, ΔBIC(M1 vs M0) falls to 3.58 (n = 613), below
the pre-specified bar.

The verdict does not move under S1–S6, but it does move under the post-hoc tail
trims S7–S9:

| Variant | n | ΔBIC(M1 vs M0) | ΔBIC(M1 vs M2) | τ (µg/L) | Same P1 verdict? |
|---|---|---|---|---|---|
| S1 — 2007 replication | 1,146 | 89.10 | 53.39 | 237.9 | yes |
| S2 — pooled 2007 + 2017 | 2,139 | 139.59 | 87.78 | 200.3 | yes |
| S3 — 2017 including hand-picked sites | 1,097 | 68.49 | 18.45 | 146.1 | yes |
| S4 — drop below-reporting-limit values | 968 | 63.62 | 8.06 | 143.7 | yes |
| S5 — τ window central 60% | 993 | 58.15 | 9.05 | 75.3 | yes |
| S6 — 2007 excluding reference lakes | 1,023 | 81.70 | 52.35 | 241.4 | yes |
| S7 — 2017, trim top 5% TP (post-hoc) | 943 | −0.79 | 0.18 | 76.1 | **no** |
| S8 — 2007, trim top 5% TP (post-hoc) | 1,088 | 6.45 | 13.37 | 101.4 | yes, marginal |
| S9 — 2017 deep, trim top 5% TP (post-hoc) | 613 | 3.58 | 6.26 | 77.8 | **no** |

S2 double-counts the 216 matched lakes and is reported only as a stress test.
S5 shows that narrowing the search window pulls τ down to 75 µg/L without
changing the verdict, which is a reminder that the break location is
model-dependent even where its existence is not. S7–S9 were added after the
full-sample fit, when the upper segment was seen to be dominated by
hypereutrophic lakes; each removes the lakes above its own 95th percentile of
TP (50 lakes above 447.9 µg/L in S7, 58 above 449.25 µg/L in S8, 33 above
242.2 µg/L in the 2017 deep stratum in S9). They are labelled post-hoc because
they were not in `hypothesis.md`, and they are reported because they change the
verdict.

## What this test cannot see

A single visit per lake measures a snapshot, not a trajectory. P1 and P2 test a
cross-sectional signature that a threshold mechanism would produce, but the same
signature — a kink and a long left tail — can arise without bistability: from a
smooth but concave relation, from a change in the limiting nutrient along the
gradient, from light limitation at high phosphorus, or from the way lakes are
distributed across the landscape. The fitted break is continuous by
construction; a genuinely discontinuous jump would be smoothed into a steep
segment, so this design cannot separate a stability boundary from a steep but
continuous response.

P3 is the closest available approach to the irreversibility clause, and it is
weaker than the clause. Two snapshots ten years apart can show that history
matters after conditioning on current phosphorus; they cannot show a lake
crossing a boundary and returning, and they cannot separate hysteresis from
fixed lake-level typology. The 216 matched lakes are the 2017 revisit panel, a
selected subset of the 2007 sample, not a fresh random draw from it.

Other limits, specific to this design:

- **Unweighted fits.** The NLA is a probability sample, but the models are
  unweighted OLS. The estimand is the relation among sampled lakes, not a
  population total; design weights would change the target, not just the
  precision.
- **Chlorophyll is a proxy for algal dominance.** Non-algal turbidity and
  light limitation can both flatten the chlorophyll–phosphorus relation, and
  this test does not use turbidity or light measurements, so it cannot separate
  them.
- **Index-season snapshot.** Within-lake seasonal variation, which is large for
  chlorophyll, is invisible.
- **Revisit-panel composition.** The 2017 sample mixes new probability sites with
  a revisit panel drawn from the 2007 sample. Unweighted fits give the revisit
  lakes their sample share rather than their population weight, so the 2017
  cross-section is not a clean probability sample of the national population.
- **Survey population.** Very small and non-probability water bodies are outside
  the target population, and the 2007 and 2017 laboratory programmes are not
  identical.
- **No manipulation.** Every statement here is an association; the data cannot
  identify the mechanism the proposal invokes, and nothing in this file should
  be read as a causal claim or as management advice.

## What would have changed the answer

- A positive slope change would have matched the proposal's directional
  prediction. The observed change is negative and its 95% CI excludes zero
  (−1.674 to −0.561), so the data are not merely underpowered on this point.
- If M2 had beaten M1 by more than 0 BIC in 2017, P1 would have failed by the
  pre-specified rule. It did not (M1 wins by 14.47); in shallow lakes it does.
- If the matched-lake count had been below 50, P3 would have been marked
  UNTESTED. At 216 it is testable, but the estimate is persistence, not an
  observed reversal.
- Removing the top 5% of TP changes the P1 verdict from supported to not
  supported under the BIC rule (ΔBIC −0.79 and 0.18). That is the most
  consequential result in this file: the threshold is a property of the
  hypereutrophic tail, not of the main gradient.
- If the bimodality rule had been applied to raw log chlorophyll rather than
  M2 residuals, P2 would have failed (p = 0.106). The pre-specified residual
  version passes, but the post-hoc binned-smoother check gives a valley depth
  of 0.373 rather than 0.616, so the two-state reading is the least secure of
  the three results.

## Every number quoted, with its `results.json` key

| Where | Numbers as quoted | `results.json` keys, in the same order |
|---|---|---|
| Primary sample | 993; 1.36–11,238.8; 0.1–3,299.18; 25; 988 | `primary_2017.n`; `samples.ranges_2017.tp_ugl`; `samples.ranges_2017.chla_ugl`; `samples.n_2017_below_rl_tp`; `samples.n_2017_depth_available` |
| 2007 replication sample | 1,146; 1–4,679; 0.067–936; 116; 1,145 | `replication_2007.n`; `samples.ranges_2007.tp_ugl`; `samples.ranges_2007.chla_ugl`; `samples.n_2007_below_rl_tp`; `samples.n_2007_depth_available` |
| Primary association | 0.749; 0.762; 0.512; 0.808 | `primary_2017.spearman_loglog`; `replication_2007.spearman_loglog`; `primary_2017.M0.r2`; `primary_2017.M0.slope` |
| Primary BIC | −1609.05; −1672.62; −1658.15; 193.72; 179.20; 183.10; 63.56; 14.47 | `primary_2017.M0.bic`; `primary_2017.M1.bic`; `primary_2017.M2.bic`; `primary_2017.M0.rss`; `primary_2017.M1.rss`; `primary_2017.M2.rss`; `primary_2017.delta_bic_M1_M0`; `primary_2017.delta_bic_M1_M2` |
| Primary LRT | 77.36; 0.0005 | `primary_2017.lrt.stat`; `primary_2017.lrt.p_boot` |
| Primary break | 144.8; 69.8–274.2; 10.1% | `primary_2017.M1.tau_ugl`; `primary_2017.M1.tau_ci95_ugl`; `primary_2017.M1.frac_tau_at_grid_edge` |
| Primary slopes | 1.026; 0.105; −0.921; −1.674 to −0.561 | `primary_2017.M1.slope_left`; `primary_2017.M1.slope_right`; `primary_2017.M1.slope_change`; `primary_2017.M1.slope_change_ci95`; `primary_2017.M1.predicted_chla_ugl.tp_50` = 13.7; `.tp_145` = 40.8; `.tp_1000` = 50.0 |
| 2007 replication | 89.10; 53.39; 103.19; 0.0005; 237.9; 95.9–255.0; 30.5%; 0.826; −0.392; −1.218 | `replication_2007.delta_bic_M1_M0`; `replication_2007.delta_bic_M1_M2`; `replication_2007.lrt.stat`; `replication_2007.lrt.p_boot`; `replication_2007.M1.tau_ugl`; `replication_2007.M1.tau_ci95_ugl`; `replication_2007.M1.frac_tau_at_grid_edge`; `replication_2007.M1.slope_left`; `replication_2007.M1.slope_right`; `replication_2007.M1.slope_change` |
| Bimodality (pre-specified) | 0.616; 0.117; 0.001; 0.002; 0.026; 0.106 | `bimodality.residual.observed_depth`; `bimodality.residual.null_q95`; `bimodality.residual.p`; `bimodality.residual.p_bh`; `bimodality.raw_chla.observed_depth`; `bimodality.raw_chla.p` |
| Bimodality (post-hoc) | 0.373; 0.091; 0.0055; 0.322; 27; 2.7% | `bimodality.posthoc.binned_smoother_20.bins` = 20; `bimodality.posthoc.binned_smoother_20.observed_depth`; `bimodality.posthoc.binned_smoother_20.null_q95`; `bimodality.posthoc.binned_smoother_20.p`; `bimodality.posthoc.m2_depth_drop_two_most_negative`; `bimodality.posthoc.n_residual_below_minus1`; `bimodality.posthoc.share_residual_below_minus1` |
| Path dependence | 216; 10⁻⁶; 0.502; 0.054; 9.23; 0.358–0.640; −16.87; 0.0001; 0.592; 0.332; 0.535; 0.745; 0.539; −0.122 | `path_dependence.n_matched`; `path_dependence.area_median_rel_diff`; `path_dependence.model.b2`; `path_dependence.model.se_b2`; `path_dependence.model.t_b2`; `path_dependence.model.b2_ci95`; `path_dependence.model.p_one_sided_log10`; `path_dependence.permutation.p_two_sided`; `path_dependence.model.r2`; `path_dependence.model.b1`; `path_dependence.partial_corr_chla_given_tp17`; `path_dependence.spearman_chla_07_17`; `path_dependence.extended_model.b2`; `path_dependence.extended_model.b3` |
| Depth table | 342; 34.99; −3.00; 176; 646; 34.13; 14.65; 71.6; 288; 10.92; −0.44; 246; 857; 100.21; 77.49; 109 | `depth.2017.shallow.n`, `.delta_bic_M1_M0`, `.delta_bic_M1_M2`, `.M1.tau_ugl`; `depth.2017.deep.*`; `depth.2007.shallow.*`; `depth.2007.deep.*` (same four fields in each block) |
| Tail trims (post-hoc) | 50; 448; 943; −0.79; 0.18; 76.1; 58; 449; 1,088; 6.45; 13.37; 101.4; 33; 242; 613; 3.58; 6.26; 77.8 | `sensitivity.S7_2017_trim_top5pct_TP.n_removed`; `.trim_threshold_ugl`; `.n`; `.delta_bic_M1_M0`; `.delta_bic_M1_M2`; `.M1.tau_ugl`; `sensitivity.S8_2007_trim_top5pct_TP.n_removed`; `.trim_threshold_ugl`; `.n`; `.delta_bic_M1_M0`; `.delta_bic_M1_M2`; `.M1.tau_ugl`; `sensitivity.S9_2017_deep_trim_top5pct_TP.n_removed`; `.trim_threshold_ugl`; `.n`; `.delta_bic_M1_M0`; `.delta_bic_M1_M2`; `.M1.tau_ugl` |
| Pre-specified thresholds | 6; 0.05; 0; 0.05; 0.10; 0.025; 50; 3 | `settings.criteria.delta_bic_m1_m0_min`; `settings.criteria.lrt_alpha`; `settings.criteria.delta_bic_m1_m2_min`; `settings.criteria.p2_alpha`; `settings.criteria.dip_depth_min`; `settings.criteria.p3_alpha`; `settings.criteria.p3_min_matched_lakes`; `settings.depth_cut_m` |
| Sensitivity table | 1,146; 2,139; 1,097; 968; 993; 1,023; 89.10; 139.59; 68.49; 63.62; 58.15; 81.70; 53.39; 87.78; 18.45; 8.06; 9.05; 52.35; 237.9; 200.3; 146.1; 143.7; 75.3; 241.4 | `sensitivity.S1_2007_replication.*`; `S2_pooled`; `S3_2017_including_HAND`; `S4_drop_below_rl`; `S5_tau_window_60pct`; `S6_2007_excluding_reference` — fields `.n`, `.delta_bic_M1_M0`, `.delta_bic_M1_M2`, `.M1.tau_ugl` in each block |

The overall coding of the verdict is in `results.json` under `verdict`
(`P1_supported`, `P1_robust_to_tail_trim`, `P1_direction`,
`claim_direction_met`, `P2_supported`, `P3_supported`, `overall`).

## Reproduction

```
cd tests/05-lake-switch
export MPLCONFIGDIR=[local provenance path]
[local provenance path] run.py
```

`run.py` reads only `data/raw/`, writes `data/lake_table.csv`, `results.json`,
`figures/fig1_change_point.png` and `figures/fig2_bootstrap_and_null.png`, and
takes about 2–4 minutes. Every random draw uses seed **20261030**. The two
figures show the 2017 scatter with the M0 and M1 fits (Figure 1) and the
10,000-draw bootstrap distribution of τ beside the 2,000-draw null distribution
of the bimodality statistic (Figure 2).

## Disclosure

Designed and executed with AI assistance (DeepSeek via its coding agent), not
peer reviewed. The pre-specified criteria are in `hypothesis.md`; the analysis
code is `run.py`; every quoted number is in `results.json`; source URLs and
checksums are in `SOURCE.md`.

hypothesis.md

# Eutrophication is a switch? — pre-specified test

**Proposal tested (verbatim, Incoprea Research, Original Theories Volume XVI #10,
posted 30 October 2026):**

> "A lake switches to algal dominance above a phosphorus threshold rather than
> degrading gradually, and the switch is not reversed by returning to the same
> concentration."

**Mechanism (verbatim):**

> "Nutrient loading crosses a stability boundary between clear and turbid states;
> internal loading from sediment maintains the turbid state after external
> loading falls."

**Prediction (verbatim):**

> "Recovery should require lower phosphorus than the level that caused the
> switch."

**Refute criterion (verbatim):**

> "Reversible, gradual response would remove the bistability."

Posted at: <https://incoprea.com/research-library/original-theories-vol16-20261030/>

This file was written before any model was fitted. It fixes the samples, the
models, the decisive comparison, the support/refute rules and the sensitivity
re-runs. Results are in `README.md` and `results.json`; nothing below was
changed after the first fit.

---

## 1. What the claim implies, operationally

The claim contains two separable assertions:

- **P1 — a switch.** Along a phosphorus gradient, chlorophyll-a should change
  slope sharply at some threshold, not smoothly.
- **P2 — two states.** At a given phosphorus concentration, lakes should cluster
  into a clear and a turbid state, so the scatter around a smooth fit should be
  bimodal rather than unimodal.
- **P3 — hysteresis.** A lake's current chlorophyll should carry information
  about its past chlorophyll beyond what its current phosphorus implies:
  lakes that were turbid at a given phosphorus should remain turbid when
  phosphorus returns to that same value.

P1 and P2 are cross-sectional signatures. P3 needs repeated measurement of the
same lakes.

## 2. Data inspected before this file was written

Only structure, never a model fit. Specifically: column names, units and flag
semantics in the five metadata files; the number of rows, sites and repeat
visits per survey; counts of missing and below-reporting-limit values; the
observed ranges of total phosphorus, chlorophyll-a, area and maximum depth;
the number of 2007 sites that also appear in the 2017 survey, matched on the
National Hydrography Dataset COMID; and the availability of a 2017 maximum
depth field. No correlation, change-point fit, ΔBIC or bimodality statistic had
been computed when the criteria in §5–§8 were fixed.

## 3. Data sources and lake-level table

Two EPA National Lakes Assessment probability surveys, treated as two
cross-sections of the same national population:

- **NLA 2007** — `NLA2007_WaterQuality_20091123.csv` (one row per site-visit;
  `PTL` total phosphorus µg/L, `CHLA` chlorophyll-a µg/L) joined to
  `NLA2007_SampledLakeInformation_20091113.csv` on `SITE_ID` and `VISIT_NO`
  (`DEPTHMAX` m, `AREA_HA`, `WSA_ECO9`, `COM_ID`). Both files are byte-identical
  members of the EPA "NLA 2007 All Data" archive (sha256 in `SOURCE.md`).
- **NLA 2017** — `nla2017_all_sites-visits_data_for_population_estimates_220705_0.csv`
  (`PTL_RESULT` µg/L, `CHLA_RESULT` µg/L, `SITETYPE`, `WGT_TP`), joined to
  `nla_2017_site_information-data.csv` (`COMID`, `AREA_HA`, `US_L3NAME`,
  `AG_ECO9`). Maximum depth is not in those two files; it is taken from
  `nla_2017_profile-data.csv` inside `NLA2017_Report_Data_Files.zip` as the
  maximum `INDEX_SITE_DEPTH` for the visit. All three are byte-identical
  members of the EPA 2017 report-data archive.

`data/lake_table.csv` is one row per lake per survey: survey year, site id,
COMID, state, ecoregion, area, maximum depth, TP, CHLA, the raw quality flags,
and a below-reporting-limit indicator.

**Inclusion rules (fixed in advance).**

- One row per lake per survey: `VISIT_NO == 1`. Repeat visits are excluded so
  a lake cannot enter twice within a survey.
- 2007: `LAKE_SAMP == "Target_Sampled"`. 2017: `SITETYPE == "PROB"`
  (probability sample). Hand-picked 2017 sites (`HAND`) are excluded from the
  primary sample and used only in sensitivity S3.
- Both TP and CHLA present and strictly positive (a log is undefined at zero).

**Below-reporting-limit rule (fixed in advance).** Both surveys report a
*numeric* value alongside the below-limit flag (`PTL_FLAG` / `PTL_NARS_FLAG` /
`CHLA_NARS_FLAG`), not merely a "less-than" token. The primary analysis keeps
those reported values, because dropping them would preferentially delete the
low-phosphorus end of the gradient where a threshold would sit. Sensitivity S4
drops every flagged value instead. 2017 `ND` (non-detect) records have no
numeric TP result and are therefore absent through the positivity rule, not
through a separate decision.

## 4. Primary sample

The **2017 probability sample** is the primary sample: it is the most recent
survey, it was analysed in one laboratory cycle, and it avoids counting the
same lake twice in the pooled table. The **2007 target sample** is a
replication in a different survey year, reported with the same models. It is
not independent: 216 lakes match one-to-one in the path-dependence sample
(§7), because the 2017 revisit panel re-sampled 2007 sites. A pooled
2017+2007 analysis is sensitivity S2.

## 5. Primary test A — is there a sharp change point?

All models use x = log10(TP in µg/L) and y = log10(CHLA in µg/L), fitted by
ordinary least squares:

| model | form | parameters (BIC) |
|---|---|---|
| **M0** | y = a + b·x | 2 |
| **M1** | y = a + b₁·x + d·(x − τ)₊ (continuous two-segment; slope b₁ below τ, b₁+d above) | 4 (3 coefficients + τ) |
| **M2** | y = a + b·x + c·x² (smooth alternative) | 3 |

τ is found by a grid search over 240 equally spaced values between the 10th and
90th percentiles of x (the central 80% of the observed phosphorus range), with
M1 refitted by least squares at every τ. M1 is continuous at τ by construction.
Reported: τ in log10(µg/L) and µg/L, the two slopes, the slope change d, AIC-free
BIC = n·ln(RSS/n) + k·ln(n), and ΔBIC defined so that **positive favours the
two-segment model**:

- ΔBIC(M1 vs M0) = BIC(M0) − BIC(M1)
- ΔBIC(M1 vs M2) = BIC(M2) − BIC(M1)

**Bootstrap likelihood-ratio test (M1 vs M0).** Parametric residual bootstrap
under M0: 2,000 draws of y* = M0 fit + N(0, σ̂²), M0 and the full M1 grid
refitted on each draw, LR* = n·(ln RSS₀* − ln RSS₁*). The p-value is
(1 + #{LR* ≥ LR_obs}) / (2,001). The naive χ²₁ and χ²₂ p-values are reported as
a check but are not decisive, because τ is not identified under M0.

**τ confidence interval.** 10,000 paired (x, y) bootstrap resamples, M1
refitted on each, percentile 2.5–97.5% interval for τ and for d. The fraction of
resamples landing on either grid edge is reported; if it is material, the
threshold is not located, only bounded.

**Pre-specified decision rule.**

- **SUPPORTED** only if all three hold:
  1. ΔBIC(M1 vs M0) > 6;
  2. bootstrap LRT p < 0.05;
  3. ΔBIC(M1 vs M2) > 0 (M1 beats the smooth quadratic).
- **NOT SUPPORTED** if any of the three fails. This is a failure to confirm the
  threshold signature, not proof that the response is gradual.
- **EVIDENCE AGAINST a sharp switch** if additionally ΔBIC(M1 vs M0) ≤ 0 and
  ΔBIC(M1 vs M2) ≤ 0, i.e. the two-segment model is worse than *both* smooth
  alternatives.

The 2007 replication is judged by exactly the same rule. There is no pooling of
the two surveys in the decisive test.

## 6. Primary test B — conditional bimodality

Residuals are taken from M2 fitted to the primary sample. The statistic is
`lib/stats.dip_like`, the depth of the deepest interior valley of a Gaussian
kernel density, expressed as a fraction of the lower adjacent peak. **It is a
descriptive valley-depth statistic, not Hartigan's dip test** and is reported as
such.

Null calibration: 2,000 parametric draws y* = M2 fit + N(0, σ̂²), M2 refitted,
`dip_like` recomputed. p = (1 + #{depth* ≥ depth_obs}) / 2,001. The null median
and 95th percentile are reported.

**Pre-specified decision rule.** P2 is SUPPORTED if the Benjamini–Hochberg
adjusted p < 0.05 **and** the observed valley depth is at least 0.10 (a 10%
density fall between peaks, fixed as the smallest depth worth calling a valley).
Otherwise P2 is not supported.

A secondary descriptive repeats the statistic on log CHLA itself against a
Gaussian null with the observed mean and SD (equivalent to a lognormal null in
raw units). The two p-values form the **only multiplicity family** in this
test and are BH-adjusted together; every other p-value in this file is a
single pre-specified test or a robustness check and is not corrected again.

## 7. Test C — the "not reversed" clause

The 2007 and 2017 surveys overlap: the 2017 revisit panel re-sampled 2007
probability lakes. Matching is on the NHD `COMID` (2007 `COM_ID`), requiring
exactly one site in each survey per COMID (1:1). Area agreement is reported as
a matching check. If fewer than 50 lakes match, P3 is marked **UNTESTED** and
nothing is substituted for it.

For matched lakes the pre-specified model is

  log10 CHLA₂₀₁₇ = b₀ + b₁·log10 TP₂₀₁₇ + b₂·log10 CHLA₂₀₀₇ + ε

and the path-dependence prediction is **b₂ > 0** with current phosphorus held
fixed. A model adding log10 TP₂₀₀₇ is reported as a check. Inference: OLS
one-sided t-test, a 10,000-draw paired bootstrap 95% CI for b₂, and a
10,000-label permutation p (labels of 2007 CHLA shuffled across lakes). All
three seeds are fixed.

**Pre-specified decision rule.** Path dependence is SUPPORTED if the one-sided
t p < 0.025 and the bootstrap 95% CI for b₂ excludes zero. Otherwise it is not
supported. The clause is UNTESTED if the match count is below 50. A positive b₂
is consistent with hysteresis but does not establish it: persistent lake-level
differences (morphometry, colour, nitrogen) also produce it, and two snapshots
cannot show a lake crossing a boundary twice.

## 8. Depth stratification and alternative definitions

**Depth.** Internal loading is expected to matter more in shallow, polymictic
lakes. M1 vs M0 vs M2 are refitted separately for 2017 lakes with maximum depth
< 3 m (shallow) and ≥ 3 m (deep), with the same criteria. The 3 m cut is fixed
here, before fitting. 2007 is repeated with `DEPTHMAX` for replication. Lakes
missing depth are excluded from this stratification only.

**Alternatives (each reported, verdict compared with the primary):**

- **S1** 2007 target sample instead of 2017 (the replication itself);
- **S2** pooled 2007 + 2017 rows (the 254 matched lakes then contribute twice);
- **S3** 2017 including hand-picked (`HAND`) sites;
- **S4** drop every below-reporting-limit value in both surveys;
- **S5** τ searched over the central 60% rather than 80% of x;
- **S6** 2007 excluding `REF_Lake` reference sites.

The verdict "moves" if the three-part rule in §5 changes state in any of these.

## 9. What this design cannot see

A single visit per lake measures a snapshot, not a trajectory. P1 and P2 test a
cross-sectional signature that a threshold mechanism implies, but the same
signature can be produced without bistability — for example by a smooth but
convex relation, by a change in the limiting nutrient along the gradient, or by
the landscape distribution of lakes. P3 uses two snapshots ten years apart; it
can detect persistence of lake identity but cannot observe a lake returning to
a previous phosphorus concentration, and it cannot separate hysteresis from
fixed lake-level differences. The survey's target population excludes very
small and non-probability water bodies, and each lake is sampled once in the
index season, so within-lake seasonal variation is invisible. No design here
manipulates phosphorus; all statements are associations.

## 10. Reporting conventions

- Two-sided tests unless a direction is stated (P3 is one-sided positive).
- Percentile bootstrap CIs; the number of resamples is stated with each.
- Seed 20261030 for every random draw; Python `numpy` default generator.
- Every number in `README.md` appears in `results.json` under a named key.
- Observational data: "associated with" is the strongest causal verb used.

Included code and results

• tests/05-lake-switch/data/lake_table.csv
• tests/05-lake-switch/results.json
• tests/05-lake-switch/run.py