--- title: 'Statistical Analysis and Exploratory Data Analysis' source: 'https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis' language: en updated: '2026-10-07T04:31:12+00:00' license: 'CC BY-SA 4.0' lessons: 13 outcomes: - 'Carries out exploratory data analysis in a structured order and writes down what it looks for at each step' - 'Presents distribution and relationship summaries with a measure suited to the scale and shows what the summary hides' - 'Applies a statistical test with its assumptions and reports the result together with its margin of variation' --- # Statistical Analysis and Exploratory Data Analysis Writing next to every reported difference how much that difference can move: the hypothesis-generating job of exploration, the shape a summary hides in a univariate distribution, the pattern the correlation coefficient cannot see in a bivariate relationship, the spurious relationship produced by the number of comparisons in a correlation matrix, the data problems caught at first glance, the number moved by the decision between mean and median, what five numbers hide in a distribution comparison, what a confidence interval says and does not say, permutation-based hypothesis testing and assumption checking, the decision count that comes out of an experiment set, the interpretive limit of a regression coefficient, the retention rate moved by the cohort definition, and counting the findings produced from noise. ## Exploratory Analysis Counting what is looked for at each step of the sequence run to get to know the dataset, and which question each step produces: the distinction between the hypothesis-generating job of exploration and confirmation, the distribution shape hidden by a measure of central tendency in univariate analysis, what a scatter plot and a cross tabulation show beyond the correlation coefficient in bivariate analysis, the spurious relationships produced by the number of comparisons in a correlation matrix and the role of dimensionality reduction, and the heaping, truncation and duplication problems caught in the first visual inspection. ### Purpose of Exploratory Analysis Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/purpose-of-exploratory-analysis The two numbers left by the Data Analysis with Python course are taken up again: the east mean is 24.674 m³, the north mean 20.210 m³, a difference of 4.464 m³. The dataset is not a population but a sample — 1,329 of the 1,400 meters in the network entered the records, 3,199 of 3,987 possible readings were observed, and coverage is 80.2%. When the region labels are randomly redistributed 2,000 times, no run reaches 4.463 m³; the most extreme permutation gives 1.754 m³, so the margin of variation of the observed difference is below 1 in 2,000. When the same procedure is applied to the rows of the north region with odd and even subscriber IDs, the observed difference of 0.232 m³ is exceeded in 1,314 of 2,000 runs. A scan of ten checks produces 8 questions and makes 0 decisions. The Data Analysis with Python course built a toolkit over twenty-two lessons and wrote down every default of it. Two numbers were left at the close: the mean consumption of the east region is 24.674 m³, that of the north 20.210 m³. The report wrote the difference as 4.464 m³ and stopped there. Two questions went unasked: do these numbers summarize all the meters in the network, or only the part that could be read, and is the 4.464 m³ difference a regional difference, or a fluctuation that chance can produce at the same size? This lesson asks those two questions and sets the unit of measure of the course. From here on, every number has a second number next to it: the **margin of variation**, how much that value could move with nothing changing. The rule of the course follows from this: a difference written without its margin of variation is an observation, not a finding. **A number whose margin is not written counts as unmeasured.** - **EA1.** The dataset is **fictional** and carried over from K03–K04: a long-format reading table of 3,199 rows, 1,329 subscribers, three periods, five regions, and tariff tiers of 10/25/40 m³. The seed is written in the code; the numbers are the same from run to run. - **EA2. Exploratory data analysis** **generates** hypotheses; it does not confirm them. Its output is not a decision but a list of questions to be answered; confirmation is a separate job. - **EA3.** This dataset is not a **population** but a **sample**. Every observed summary is a summary of the rows at hand, not of the population, and would have come out differently in another sample. - **EA4. The margin of variation** is measured in this course with **permutation**: group labels are randomly redistributed, the same number is recomputed, and the number of runs in which the observed value is exceeded is counted. No distribution table or critical value is used. - **EA5.** The resolution of the margin is **limited by the run count**. For a difference that is exceeded in none of 2,000 runs, what can be said is not "its margin is zero" but "its margin is below 1 in 2,000". #### Population or Sample ```python # exploration.py — MODEL: the fictional measurement network of K03 and K04 is rebuilt # from the same seed. The dataset is fictional; no real subscriber, meter or region appears. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF PERIODS = ["2026-01", "2026-02", "2026-03"] TARIFFS = [10, 25, 40] REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t: int): """Linear congruential generator: seed written in the code, same output every run.""" x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value() -> float: s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u: float, w: list[float]) -> int: t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 _subscribers, _dropped = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] k = {"subscriber_id": 10001 + i, "region": b[0], "base": b[2], "household_size": pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]), "satisfaction": pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) + 1} (_dropped if r() < 0.046 else _subscribers).append(k) _valid = {a["subscriber_id"] for a in _subscribers} READINGS = [] for k in _subscribers + _dropped: r = rng(SEED + 7000 + k["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor(k["base"] * math.exp( (r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if k["subscriber_id"] not in _valid: r() continue READINGS.append({"subscriber_id": k["subscriber_id"], "region": k["region"], "period": PERIODS[d], "consumption_m3": t, "household_size": k["household_size"], "satisfaction": k["satisfaction"], "tier": 1 if t < TARIFFS[0] else 2 if t < TARIFFS[1] else 3 if t < TARIFFS[2] else 4}) r() def mean(v: list[float]) -> float: return sum(v) / len(v) def region_values(name: str) -> list[float]: return [o["consumption_m3"] for o in READINGS if o["region"] == name] east, north = region_values("east"), region_values("north") POSSIBLE = len(_valid) * len(PERIODS) print(f"fictional dataset (seed {SEED}): {len(READINGS)} rows, {len(_valid)} subscribers") print(f"east {len(east):4d} rows, mean {mean(east):.3f} m3") print(f"north {len(north):4d} rows, mean {mean(north):.3f} m3") print(f"observed difference {mean(east) - mean(north):.3f} m3") print(f"meters in network {RAW}, subscribers in records {len(_valid)}, " f"never in records {RAW - len(_valid)}") print(f"possible readings {POSSIBLE}, observed readings {len(READINGS)}, " f"missing {POSSIBLE - len(READINGS)}, " f"coverage {100 * len(READINGS) / POSSIBLE:.1f} percent") ``` ``` fictional dataset (seed 20260218): 3199 rows, 1329 subscribers east 557 rows, mean 24.674 m3 north 921 rows, mean 20.210 m3 observed difference 4.463 m3 meters in network 1400, subscribers in records 1329, never in records 71 possible readings 3987, observed readings 3199, missing 788, coverage 80.2 percent ``` The last two lines close the question. The network has 1,400 meters, and the number of subscribers that entered the records is 1,329; 71 meters appear in no row. Had the subscribers in the records been read over three periods, there would be 3,987 rows; the observed rows number 3,199 and coverage stays at 80.2%. This dataset is not the whole network but the part of it that **could be read** — that is, a sample, and not a randomly drawn one either, because what selects which meter could be read is not a lottery but the field itself. The consequence is direct: 24.674 m³ is not the mean consumption of the east region but the mean of the **557 rows that could be read** in the east region. 557 rows drawn from the same network in another week would have given another number. The question itself changes here: the question "is the number correct" is replaced by "how much does this number move". The close of the previous course wrote the difference as 4.464 m³, while the running block writes 4.463 m³: the close subtracted two rounded numbers, the block the exact values. The 0.001 m³ between them changes nothing here, but it is the smallest example of the fact that a reported number depends on how it was computed. #### The Margin of Variation of a Difference One way to ask whether a difference is real is to build a world in which the difference **does not exist** and count how often a difference of the same size appears there. If the region label says nothing about consumption, randomly redistributing the labels should not change the result. ```python # margin.py — margin of variation: region labels are randomly redistributed, the same # difference is recomputed, and the runs that exceed the observed value are counted. RUNS = 2000 TOLERANCE = 1e-9 # one difference computed two ways shifts in the last digit; equality is not tested def permutation_margin(a: list[float], b: list[float], seed: int): """The values of both groups are pooled. In each run len(a) items of the pool are picked at random as the first group and the difference is measured again.""" r = rng(seed) pool, na, nb = a + b, len(a), len(b) total, observed = sum(pool), mean(a) - mean(b) exceeded, extreme = 0, 0.0 for _ in range(RUNS): h, s = pool[:], 0.0 for i in range(na): j = i + int(r() * (na + nb - i)) h[i], h[j] = h[j], h[i] s += h[i] diff = s / na - (total - s) / nb if abs(diff) >= abs(observed) - TOLERANCE: exceeded += 1 extreme = max(extreme, abs(diff)) return observed, exceeded, extreme def north_half(odd: bool) -> list[float]: return [o["consumption_m3"] for o in READINGS if o["region"] == "north" and (o["subscriber_id"] % 2 == 1) == odd] COMPARISONS = [("east / north", east, north, SEED + 501), ("north: odd / even subscriber id", north_half(True), north_half(False), SEED + 502)] WIDTHS = [33, 7, 7, 12, 12, 10] y = lambda h: "".join(str(v).rjust(WIDTHS[i]) if i else str(v).ljust(WIDTHS[0]) for i, v in enumerate(h)) print(y(["comparison", "n1", "n2", "observed", "extreme", "exceeded"])) counts = {} for name, a, b, seed in COMPARISONS: observed, exceeded, extreme = permutation_margin(a, b, seed) counts[name] = exceeded print(y([name, len(a), len(b), f"{observed:.3f}", f"{extreme:.3f}", exceeded])) print(f"\n{RUNS} runs; exceeded = runs that produced a difference at least as large as the observed one") for name, exceeded in counts.items(): print(f"{name:<33}margin of variation " + (f"below 1/{RUNS}" if exceeded == 0 else f"{exceeded / RUNS:.3f}, exceeded in " f"{round(100 * exceeded / RUNS)} of 100 runs")) ``` ``` comparison n1 n2 observed extreme exceeded east / north 557 921 4.463 1.754 0 north: odd / even subscriber id 452 469 0.232 1.647 1314 2000 runs; exceeded = runs that produced a difference at least as large as the observed one east / north margin of variation below 1/2000 north: odd / even subscriber id margin of variation 0.657, exceeded in 66 of 100 runs ``` Two rows went through the same procedure and gave opposite answers. None of the 2,000 permutations reached the 4.463 m³ difference between east and north; the most extreme difference is 1.754 m³, less than half of the observed one. The region label says something about consumption. The second row is a check. Whether a subscriber ID is odd or even cannot affect consumption; this breakdown is noise by definition, and the tool says so: the observed difference of 0.232 m³ is exceeded in 1,314 of 2,000 runs, that is, in sixty-six of every hundred runs. Writing such a difference into the report as "odd-numbered subscribers in the north consume more" would mean counting as a finding something that permutation produces two thirds of the time. The resolution limit shows in the first row. A run count of zero exceeded does not mean "the margin is zero"; with 2,000 runs the smallest measurable margin is 1 in 2,000, and the observed difference is below it. Had the run count been raised tenfold, the limit would have dropped tenfold. The margin is itself a measurement and is written together with its own resolution. What the procedure does not say is also written. Permutation answers only the question "how often does a difference of this size appear when the labels are permuted"; it does not say the **cause** of the difference. The tariff base of the east may be higher, household size may differ, or the meters that could be read in the east may systematically be the large subscribers. A small margin does not take the place of an explanation. #### Exploration Does Not Decide, It Produces Questions Measuring the margin of a difference is a step of exploration, not its end. The job of exploration is to walk through the dataset and work out which questions are worth asking; the scan below writes this as a list of measures and turns every measure that crosses its threshold into a question. ```python # question.py — the output of exploration is a list of questions, not a decision. Each # check computes a measure; one that crosses its threshold yields a question, not a decision. from collections import Counter def median(v: list[float]) -> float: s = sorted(v) n = len(s) return s[n // 2] if n % 2 else (s[n // 2 - 1] + s[n // 2]) / 2 T, S = [o["consumption_m3"] for o in READINGS], sorted(o["consumption_m3"] for o in READINGS) region_count = Counter(o["region"] for o in READINGS) period_count = Counter(o["period"] for o in READINGS) tier_count = Counter(o["tier"] for o in READINGS) subscriber_count = Counter(o["subscriber_id"] for o in READINGS) satisfaction_count = Counter(o["satisfaction"] for o in READINGS) share = lambda k: 100 * k / len(READINGS) CHECKS = [ ("consumption", "share of rows written as zero", share(sum(1 for t in T if t == 0)), 1.0, "upper", "is zero a measurement or an unread meter"), ("coverage", "observed / possible readings", 100 * len(READINGS) / POSSIBLE, 95.0, "lower", "do the read differ from the unread"), ("subscriber_id", "share of subscribers with 3 periods", 100 * sum(1 for v in subscriber_count.values() if v == 3) / len(subscriber_count), 90.0, "lower", "are the missing periods randomly distributed"), ("region", "smallest / largest region", 100 * min(region_count.values()) / max(region_count.values()), 60.0, "lower", "how much does the margin of variation of the small region grow"), ("region", "share of rows, empty region label", share(sum(1 for o in READINGS if not o["region"])), 0.0, "upper", "which region do the unlabeled rows belong to"), ("period", "row range between periods", 100 * (max(period_count.values()) - min(period_count.values())) / max(period_count.values()), 10.0, "upper", "is the period difference behavior or the calendar"), ("consumption", "mean - median range", 100 * (mean(T) - median(T)) / median(T), 2.0, "upper", "is the distribution skewed, which summary is reported"), ("consumption", "largest value / third quartile", max(S) / S[int(0.75 * len(S))], 2.5, "upper", "how much do extreme values move the mean"), ("tier", "share of the rarest tier", share(min(tier_count.values())), 5.0, "lower", "how much does the ratio of the rare tier vary"), ("satisfaction", "share of most common satisfaction level", share(max(satisfaction_count.values())), 30.0, "upper", "is a mean legitimate on an ordinal scale"), ] E = [15, 42, 8, 11] z = lambda h: "".join(str(v).ljust(E[i]) if i < 2 else str(v).rjust(E[i]) for i, v in enumerate(h[:4])) + " " + str(h[4]) print(z(["field", "measure", "value", "threshold", "question produced"])) questions = 0 for field, measure, value, threshold, direction, text in CHECKS: fired = value > threshold if direction == "upper" else value < threshold questions += fired print(z([field, measure, f"{value:.1f}", f"{threshold:.1f}", text if fired else "-"])) print(f"\n{len(CHECKS)} checks run, {questions} questions produced, 0 decisions made") ``` ``` field measure value threshold question produced consumption share of rows written as zero 3.7 1.0 is zero a measurement or an unread meter coverage observed / possible readings 80.2 95.0 do the read differ from the unread subscriber_id share of subscribers with 3 periods 66.5 90.0 are the missing periods randomly distributed region smallest / largest region 51.7 60.0 how much does the margin of variation of the small region grow region share of rows, empty region label 0.0 0.0 - period row range between periods 33.5 10.0 is the period difference behavior or the calendar consumption mean - median range 3.9 2.0 is the distribution skewed, which summary is reported consumption largest value / third quartile 2.1 2.5 - tier share of the rarest tier 1.8 5.0 how much does the ratio of the rare tier vary satisfaction share of most common satisfaction level 33.4 30.0 is a mean legitimate on an ordinal scale 10 checks run, 8 questions produced, 0 decisions made ``` Ten checks ran and eight questions came out. There is not a single decision in the table: no row says "these rows should be dropped" or "this region consumes more". That is the output of exploration, and it is exactly where exploration parts from confirmation. Confirmation takes one question and produces one answer for it; exploration walks through the dataset and produces a list of questions to be answered. The two rows that do not cross the threshold carry information as well. There are no rows with an empty region label, because the region comes from the subscriber record and was matched in the Data Collection and Preparation course; the ratio of the largest value to the third quartile is 2.1, that is, there is no extreme break at the upper end. A check that passes is also a measurement result. The thresholds are themselves decisions. Had 75% been written instead of the 95% coverage threshold, the second row would not have produced a question; had 5% been chosen instead of 2% for the mean–median range, the seventh row would have gone silent. The number eight in the table is not a property of the dataset but a number produced jointly by the dataset and the threshold list. That is why the thresholds are written in the code and visible in the table. The cost of the eight questions is paid in the next topic. Every breakdown looked at during exploration is a debt that enlarges the margin of the test: looking at ten measures and testing the most striking one is not the same as looking at a single measure and testing it. This lesson only writes that the debt was opened. #### Summary - The dataset is a sample, not a population: 1,329 of the 1,400 meters entered the records, 3,199 of 3,987 possible readings were observed, and coverage is 80.2%. - The observed difference between east and north is 4.463 m³; when the region labels are randomly redistributed 2,000 times, the most extreme difference stays at 1.754 m³ and no run reaches the observed one. - When the same procedure is applied to the rows of the north with odd and even subscriber IDs, the observed difference of 0.232 m³ is exceeded in 1,314 of 2,000 runs; the tool reports noise as noise. - The resolution of the margin is the run count: for a difference that is never exceeded in 2,000 runs, what can be said is not "its margin is zero" but "its margin is below 1 in 2,000". - The output of exploration is a question, not a decision: the scan of ten checks produced 8 questions and made 0 decisions, and this number is the joint product of the dataset and the threshold list. #### Next Step This lesson looked at the difference between two region means and wrote its margin next to it. The mean itself was not questioned: when "the mean of 557 rows is 24.674 m³" is written, nothing has been said about how those rows are distributed, and distributions that give the same mean can be very different from each other. The next lesson covers univariate analysis: two distributions that give the same mean are built, the five-number summary and bin counts are printed, and it is counted how the range between the mean and the median gives skewness away and how much the outliers that were not removed move the summary. ### Univariate Analysis Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/univariate-examination The 557 readings of the east region are placed beside a symmetric model distribution with the same mean (24.674 m³) and the same standard deviation (9.198), and the two reported numbers are shown to say nothing about shape: the five-number summaries diverge as 0.00 / 19.23 / 24.12 / 29.85 / 51.90 versus 0.66 / 18.37 / 24.61 / 31.43 / 47.79, the interquartile range comes out 10.62 against 13.06, and the range between the mean and the median is 0.554 against 0.064. Bin counts give what the five numbers cannot see: the observation has 20 values in the first bin and 1 in the second, and if a bin width of 10 is chosen this break is buried in a single number, 21. The 30 values outside the fences are not removed, their effect is measured: removal would move the east mean to 25.154 and the difference between east and north from 4.463 to 4.622 m³. The previous lesson looked at the difference between two region means and wrote its margin of variation next to it: the 4.463 m³ difference appeared in none of 2,000 permutations. The mean itself was never questioned. The 24.674 m³ written for the east region is 557 readings compressed into a single number, and it says nothing about how those 557 numbers are arranged. This lesson covers how a column is examined on its own, and its unit of measure is this: **the shape that the summary hides**. Two distributions that give the same mean and the same standard deviation are built, both are reported with the same two numbers, and it is shown that the difference between them never enters the report. - **EA6.** Univariate analysis asks a column three questions: where the center is, how scattered the values are, and what the shape of the distribution is. Unless all three are written together, the column has not been described. - **EA7. The five-number summary** is the minimum, the first quartile, the median, the third quartile and the maximum. The interquartile range is the distance between Q3 and Q1 and was defined in the Data Collection and Preparation course. - **EA8. Skewness** is measured in this course **by the range between the mean and the median**. If the range is positive, the long tail is on the right; the mean follows the tail, the median does not. - **EA9. Bin counts** give the shape of the distribution as numbers: the limits of each bin and the number of values falling into it are printed. **The bin width is a decision** and it changes the shape. - **EA10.** Outliers are **not removed** in this course. Removal was a cleaning decision of the Data Collection and Preparation course and its unit of measure was the row count; here the same values stay in the set, and only their effect on the summary is counted. #### Center, Dispersion and Shape ```python # distribution.py — MODEL: the fictional measurement network of K03 and K04 is rebuilt # from the same seed. The dataset is fictional; no real subscriber, meter or region appears. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF PERIODS = ["2026-01", "2026-02", "2026-03"] TARIFFS = [10, 25, 40] REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t: int): """Linear congruential generator: seed written in the code, same output every run.""" x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value() -> float: s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u: float, w: list[float]) -> int: t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 _subscribers, _dropped = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] k = {"subscriber_id": 10001 + i, "region": b[0], "base": b[2], "household_size": pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]), "satisfaction": pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) + 1} (_dropped if r() < 0.046 else _subscribers).append(k) _valid = {a["subscriber_id"] for a in _subscribers} READINGS = [] for k in _subscribers + _dropped: r = rng(SEED + 7000 + k["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor(k["base"] * math.exp( (r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if k["subscriber_id"] not in _valid: r() continue READINGS.append({"subscriber_id": k["subscriber_id"], "region": k["region"], "period": PERIODS[d], "consumption_m3": t, "household_size": k["household_size"], "satisfaction": k["satisfaction"], "tier": 1 if t < TARIFFS[0] else 2 if t < TARIFFS[1] else 3 if t < TARIFFS[2] else 4}) r() def mean(v: list[float]) -> float: return sum(v) / len(v) def std_dev(v: list[float]) -> float: """Standard deviation; denominator n-1, because what is at hand is a sample.""" m = mean(v) return math.sqrt(sum((x - m) ** 2 for x in v) / (len(v) - 1)) def median(v: list[float]) -> float: s = sorted(v) n = len(s) return s[n // 2] if n % 2 else (s[n // 2 - 1] + s[n // 2]) / 2 def quartile(v: list[float], p: float) -> float: """Linear-interpolation quartile; same as the definition in K03.""" s = sorted(v) i = p * (len(s) - 1) a = int(i) return s[a] if a + 1 >= len(s) else s[a] + (s[a + 1] - s[a]) * (i - a) def region_values(name: str) -> list[float]: return [o["consumption_m3"] for o in READINGS if o["region"] == name] east, north = region_values("east"), region_values("north") print(f"fictional dataset (seed {SEED}): {len(READINGS)} rows, {len(_valid)} subscribers") print(f"east: {len(east)} rows, mean {mean(east):.3f} m3, " f"median {median(east):.2f} m3, standard deviation {std_dev(east):.3f}") print(f"mean - median = {mean(east) - median(east):.3f} m3, " f"{100 * (mean(east) - median(east)) / std_dev(east):.1f} percent of the standard deviation") ``` ``` fictional dataset (seed 20260218): 3199 rows, 1329 subscribers east: 557 rows, mean 24.674 m3, median 24.12 m3, standard deviation 9.198 mean - median = 0.554 m3, 6.0 percent of the standard deviation ``` A column cannot be described with a single number. Central tendency gives two separate numbers: the mean is 24.674 m³ and the median 24.12 m³. **Dispersion** is the third number: the standard deviation is 9.198, that is, the readings wander around the center by more than nine cubic meters on average. These three numbers do not yet tell the shape, but they give the first sign of it. The sign is the parting of the mean and the median. The median looks at the ranking and does not care how far away a value is; the mean carries the distance of every value. The 0.554 m³ range between the two says the long tail is on the right: a few large readings have pulled the mean up and left the median in place. The range is 6.0% of the standard deviation and can be reported directly as the measure of skewness. #### Same Mean, Different Distribution ```python # equivalent.py — a second distribution with the same mean and the same standard deviation. # It is fictional: the sum of three uniform draws is symmetric, then it is fitted to the mean # and standard deviation of the east. The aim is comparison, not simulation. def symmetric(n: int, target_mean: float, target_std: float, seed: int) -> list[float]: r = rng(seed) raw = [r() + r() + r() - 1.5 for _ in range(n)] m, s = mean(raw), std_dev(raw) v = [round(target_mean + (x - m) * target_std / s, 2) for x in raw] shift = target_mean - mean(v) # rounding shifts the mean a little return [round(x + shift, 2) for x in v] MODEL = symmetric(len(east), mean(east), std_dev(east), SEED + 601) def five_numbers(v: list[float]) -> list[float]: return [min(v), quartile(v, 0.25), median(v), quartile(v, 0.75), max(v)] HEADERS = ["minimum", "Q1", "median", "Q3", "maximum", "IQR", "mean", "skewness"] WIDTHS = [17] + [10] * 8 y = lambda h: "".join(str(v).ljust(WIDTHS[i]) if not i else str(v).rjust(WIDTHS[i]) for i, v in enumerate(h)) print(y(["distribution"] + HEADERS)) for name, v in [("east (observed)", east), ("model (symmetric)", MODEL)]: b = five_numbers(v) print(y([name] + [f"{x:.2f}" for x in b] + [f"{b[3] - b[1]:.2f}", f"{mean(v):.3f}", f"{mean(v) - median(v):.3f}"])) print(f"\nboth distributions have a mean of {mean(east):.3f} and a standard deviation of {std_dev(east):.3f}") print(f"mean - median range: observed {mean(east) - median(east):.3f}, " f"model {mean(MODEL) - median(MODEL):.3f}") ``` ``` distribution minimum Q1 median Q3 maximum IQR mean skewness east (observed) 0.00 19.23 24.12 29.85 51.90 10.62 24.674 0.554 model (symmetric) 0.66 18.37 24.61 31.43 47.79 13.06 24.674 0.064 both distributions have a mean of 24.674 and a standard deviation of 9.198 mean - median range: observed 0.554, model 0.064 ``` Two rows are reported with the same two numbers. The mean is the same to the third decimal, and so is the standard deviation. The sentence "Mean consumption in the east region is 24.674 m³, standard deviation 9.198" is true for both and does not tell them apart. The five-number summary does. The minimum of the observation is 0.00, that of the model 0.66. The third quartile of the observation is 29.85, that of the model 31.43. The real difference is in the interquartile range: 10.62 in the observation, 13.06 in the model. Of two distributions that give the same standard deviation, one has a markedly narrower middle half. This means that what enlarges the standard deviation of the observation is not the dispersion in the middle but a few values at the ends. The last column reduces the measure to a single number. Skewness is 0.554 in the observation and 0.064 in the model — more than eight times. In a symmetric distribution the mean and the median are expected to coincide, and the model does this; the remaining 0.064 is the residue of the draw. In the observation the range is no longer a residue, it is the shape of the distribution itself. #### Shape as Numbers ```python # bins.py — the shape of the distribution is printed as numbers: how many values in each bin. # No chart is drawn; the bin width is a decision and it changes the shape. def count_bins(v: list[float], width: int, top: int = 55) -> list[int]: k = [0] * (top // width) for x in v: k[min(int(x // width), len(k) - 1)] += 1 return k WIDTH = 5 a5, m5 = count_bins(east, WIDTH), count_bins(MODEL, WIDTH) print(f"{'bin':>12}{'east':>8}{'model':>8} east bar") for i, (x, y2) in enumerate(zip(a5, m5)): print(f"[{i * WIDTH:2d},{(i + 1) * WIDTH:2d})".rjust(12) + f"{x:8d}{y2:8d} " + "#" * round(x / 4)) a10 = count_bins(east, 2 * WIDTH) print(f"\nwidth {WIDTH}: first two bins {a5[0]} and {a5[1]}; " f"width {2 * WIDTH}: first bin {a10[0]}") print(f"readings written as zero {sum(1 for x in east if x == 0)}, " f"real readings in the [{WIDTH},{2 * WIDTH}) bin {a5[1]}") ``` ``` bin east model east bar [ 0, 5) 20 7 ##### [ 5,10) 1 32 [10,15) 35 51 ######### [15,20) 111 83 ############################ [20,25) 135 112 ################################## [25,30) 119 98 ############################## [30,35) 71 92 ################## [35,40) 35 57 ######### [40,45) 16 23 #### [45,50) 11 2 ### [50,55) 3 0 # width 5: first two bins 20 and 1; width 10: first bin 21 readings written as zero 20, real readings in the [5,10) bin 1 ``` Bin counts show what the five numbers cannot see. There are 20 values in the first bin and 1 value in the second. This is not a tail, it is a break: between the 20 records heaped at zero and the region of real readings sits an almost empty bin. The count says this because the 20 records heaped at zero are rows where zero was written since the meter could not be read; they are not a measured consumption but a record format. The model has no such break: 7 and 32 are consecutive and run smoothly. The bin width is a decision, and changing this decision changes the shape. When width 10 is chosen, the first bin holds 21 values; 20 and 1 melt into a single number, and the break disappears from the table. The same data, the same procedure, two different shapes. This is why the bin width is written inside the report, not in a footnote. #### Outliers Are Not Removed, They Are Measured ```python # outlier.py — values outside the fences are NOT removed here. In the Data Collection and # Preparation course removal was a cleaning decision; here the values stay, and the effect is measured. def fences(v: list[float]) -> tuple[float, float]: q1, q3 = quartile(v, 0.25), quartile(v, 0.75) iqr = q3 - q1 return q1 - 1.5 * iqr, q3 + 1.5 * iqr WIDTHS = [9, 6, 12, 12, 10, 10, 10, 11, 8] row = lambda h: "".join(str(v).ljust(WIDTHS[i]) if not i else str(v).rjust(WIDTHS[i]) for i, v in enumerate(h)) print(row(["dataset", "n", "lower fence", "upper fence", "below", "above", "mean", "kept mean", "diff"])) kept = {} for name, v in [("east", east), ("north", north), ("model", MODEL)]: lower, upper = fences(v) kept[name] = [x for x in v if lower <= x <= upper] print(row([name, len(v), f"{lower:.2f}", f"{upper:.2f}", sum(1 for x in v if x < lower), sum(1 for x in v if x > upper), f"{mean(v):.3f}", f"{mean(kept[name]):.3f}", f"{mean(kept[name]) - mean(v):+.3f}"])) print(f"\neast / north difference: with all rows {mean(east) - mean(north):.3f} m3, " f"if the outside-fence values were removed {mean(kept['east']) - mean(kept['north']):.3f} m3") print(f"east standard deviation {std_dev(east):.3f} -> {std_dev(kept['east']):.3f}, " f"median {median(east):.2f} -> {median(kept['east']):.2f}") ``` ``` dataset n lower fence upper fence below above mean kept mean diff east 557 3.30 45.78 20 10 24.674 25.154 +0.480 north 921 2.41 37.96 32 18 20.210 20.532 +0.321 model 557 -1.22 51.02 0 0 24.674 24.674 +0.000 east / north difference: with all rows 4.463 m3, if the outside-fence values were removed 4.622 m3 east standard deviation 9.198 -> 7.416, median 24.12 -> 24.22 ``` The values outside the fences are not deleted here; the table only counts what would happen if they were. In the east there are 30 rows outside the fences, and 20 of them are outside the **lower** fence. Removal does not lower the mean, it **raises** it from 24.674 to 25.154, because most of the removed rows are readings written as zero. This is where counting an outlier as a "large value" misleads. The same decision moves the mean in the north by 0.321 m³, and since the two movements do not go the same way, the regional difference changes too: 4.622 m³ instead of 4.463 m³. The difference measured in the previous lesson has shifted by 0.159 m³ with no new data, through a single decision about the fence. This is a variation separate from the margin that comes from the sample, and added to it. The model row makes clear what the measure is. In a symmetric distribution there is **no** value outside the fences, and removal leaves the mean exactly where it is. The number of values outside the fences is a consequence not of the data being "faulty" but of the shape of the distribution; the number 30 measures the tails of the east, not its quality. The last row shows why the median is called resistant. When thirty rows are taken out, the standard deviation falls from 9.198 to 7.416, that is, by 19.4%; the median moves from 24.12 to 24.22, that is, by 0.10. The measure of dispersion depends on the ends, the measure of ranking does not. #### Summary - A column is described with three measures: for the east, a mean of 24.674 m³, a median of 24.12 m³ and a standard deviation of 9.198; the 0.554 m³ range between them is the first sign of shape. - A symmetric model distribution that gives the same mean and the same standard deviation diverges in the five-number summary: an interquartile range of 10.62 against 13.06, skewness of 0.554 against 0.064. - Bin counts give what the five numbers cannot see: the observation has 20 values in the first bin and 1 in the second, and this break comes from the readings written as zero. - The bin width is a decision: when width 10 is chosen, the first bin becomes 21 and the break is not visible in the table. - The 30 rows outside the fences are not removed; removal would have moved the east mean to 25.154 and the difference between east and north from 4.463 to 4.622 m³, and lowered the standard deviation by 19.4%. #### Next Step This lesson looked at a single column and wrote its shape as numbers. Most of the questions in the dataset cannot be answered with a single column, though: do consumption and household size rise together, is there a link between tariff tier and region, how does the satisfaction level relate to consumption. Examining two columns together brings a new summary — the correlation coefficient — and that coefficient, just like the mean, compresses everything into a single number. The next lesson covers bivariate analysis: two different relationships that give the same coefficient are built, a cross tabulation shows what the coefficient cannot see, and it is written why rank correlation is the right measure on an ordinal scale. ### Bivariate Analysis Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/bivariate-examination The consumption column is paired with four columns and the margin of variation is printed next to each coefficient: the coefficients +0.010 with household size and −0.010 with satisfaction are exceeded in 1,189 and 1,154 of 2,000 permutations, the −0.111 that comes out with period order is exceeded in none, and the +0.883 with tariff tier is not a discovery but a definition. Three fictional relationships that give the same +0.600 coefficient are placed side by side: the slice means of the linear one climb smoothly from 6.6 to 24.9, those of the curvilinear one rise to 39.8 and return to 37.1, those of the leverage one stay constant around 20, and when a single extreme point is removed the coefficient drops to −0.057. The cross tabulation of region and tariff tier counts that the fourth tier never appears in the west and that the east has three times its dataset share. When the coding changes on an ordinal scale, the coefficient moves by 79% in absolute value on satisfaction, while rank correlation stays at −0.0102. The previous lesson looked at a single column and wrote its shape with five numbers and bin counts. Most questions ask for two columns, though: does consumption rise together with household size, does it fall as the periods go on, is it related to the satisfaction level. When two columns are examined together a new summary comes out — the correlation coefficient — and it too compresses a pile of numbers into a single number. The unit of measure of this lesson is **the pattern the coefficient cannot see**. Three separate relationships that give the same coefficient are built, all three are reported with the same number, and it is shown that the difference between them never shows in the coefficient. - **EA11.** Bivariate analysis considers two columns together. The **correlation coefficient** is the ratio of the joint variation to the product of the two separate variations and measures only **linear** joint variation. - **EA12.** The margin is written next to the coefficient: one of the columns is randomly rearranged and it is counted in how many runs a coefficient of the same size appears. - **EA13.** The coefficient is not a substitute for the **scatter plot**. The scatter is not drawn here; it is printed with **slice means**: the driver column is divided into five equal parts and the mean of the second column is written for each part. - **EA14. Cross tabulation** is the joint count of two categorical columns. On a nominal scale the correlation coefficient is not defined; the relationship is read from cell counts. - **EA15. Rank correlation** looks only at the ranking. The numbers given to levels on an ordinal scale are a coding decision — the scale definitions were set up in the Introduction to Data Analytics course and are not repeated here. The coefficient is affected by this decision; rank correlation is not. #### The Coefficient and Its Margin ```python # exploration.py — MODEL: the fictional measurement network of K03 and K04 is rebuilt # from the same seed. The dataset is fictional; no real subscriber, meter or region appears. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF PERIODS = ["2026-01", "2026-02", "2026-03"] TARIFFS = [10, 25, 40] REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t: int): """Linear congruential generator: seed written in the code, same output every run.""" x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value() -> float: s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u: float, w: list[float]) -> int: t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 _subscribers, _dropped = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] k = {"subscriber_id": 10001 + i, "region": b[0], "base": b[2], "household_size": pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]), "satisfaction": pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) + 1} (_dropped if r() < 0.046 else _subscribers).append(k) _valid = {a["subscriber_id"] for a in _subscribers} READINGS = [] for k in _subscribers + _dropped: r = rng(SEED + 7000 + k["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor(k["base"] * math.exp( (r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if k["subscriber_id"] not in _valid: r() continue READINGS.append({"subscriber_id": k["subscriber_id"], "region": k["region"], "period": PERIODS[d], "consumption_m3": t, "household_size": k["household_size"], "satisfaction": k["satisfaction"], "tier": 1 if t < TARIFFS[0] else 2 if t < TARIFFS[1] else 3 if t < TARIFFS[2] else 4}) r() def mean(v: list[float]) -> float: return sum(v) / len(v) def corr(x: list[float], y: list[float]) -> float: """Correlation coefficient: the ratio of the joint variation to the product of the two separate variations.""" mx, my = mean(x), mean(y) sx = math.sqrt(sum((a - mx) ** 2 for a in x)) sy = math.sqrt(sum((b - my) ** 2 for b in y)) return sum((a - mx) * (b - my) for a, b in zip(x, y)) / (sx * sy) def ranks(v: list[float]) -> list[float]: """Equal values take the average rank; this is what an ordinal scale requires.""" positions: dict[float, list[int]] = {} for k, i in enumerate(sorted(range(len(v)), key=lambda j: v[j])): positions.setdefault(v[i], []).append(k + 1) return [mean(positions[x]) for x in v] def rank_corr(x: list[float], y: list[float]) -> float: return corr(ranks(x), ranks(y)) RUNS = 2000 TOLERANCE = 1e-9 # the same coefficient computed two ways shifts in the last digit def corr_margin(x: list[float], y: list[float], seed: int) -> tuple[int, float]: """The y column is randomly permuted; the coefficient is re-measured in every run.""" r, observed = rng(seed), corr(x, y) mx, my = mean(x), mean(y) divisor = math.sqrt(sum((a - mx) ** 2 for a in x) * sum((b - my) ** 2 for b in y)) dx, dy = [a - mx for a in x], [b - my for b in y] exceeded, extreme = 0, 0.0 for _ in range(RUNS): h = dy[:] for i in range(len(h) - 1, 0, -1): j = int(r() * (i + 1)) h[i], h[j] = h[j], h[i] c = sum(a * b for a, b in zip(dx, h)) / divisor exceeded += abs(c) >= abs(observed) - TOLERANCE extreme = max(extreme, abs(c)) return exceeded, extreme CONSUMPTION = [o["consumption_m3"] for o in READINGS] PAIRS = [("household_size", [o["household_size"] for o in READINGS]), ("satisfaction", [o["satisfaction"] for o in READINGS]), ("period order", [PERIODS.index(o["period"]) for o in READINGS]), ("tier", [o["tier"] for o in READINGS])] def align(w: list[int], left: int = 1): """The first `left` columns are left-aligned, the rest right-aligned.""" return lambda h: "".join(str(v).ljust(w[i]) if i < left else str(v).rjust(w[i]) for i, v in enumerate(h)) y = align([30, 7, 11, 12, 10, 10]) print(y(["pair", "n", "coef.", "rank coef.", "extreme", "exceeded"])) for i, (name, v) in enumerate(PAIRS): exceeded, extreme = corr_margin(CONSUMPTION, v, SEED + 710 + i) print(y([f"consumption / {name}", len(v), f"{corr(CONSUMPTION, v):+.3f}", f"{rank_corr(CONSUMPTION, v):+.3f}", f"{extreme:.3f}", exceeded])) print(f"\n{RUNS} runs; extreme = the largest coefficient seen in the permuted sets") ``` ``` pair n coef. rank coef. extreme exceeded consumption / household_size 3199 +0.010 +0.014 0.069 1189 consumption / satisfaction 3199 -0.010 -0.010 0.058 1154 consumption / period order 3199 -0.111 -0.113 0.061 0 consumption / tier 3199 +0.883 +0.813 0.065 0 2000 runs; extreme = the largest coefficient seen in the permuted sets ``` If only the coefficient column is looked at, the first two rows can be read as "a weak but existing relationship": +0.010 and −0.010 are not zero. The margin column closes this. The permuted sets produce a coefficient of this size in 1,189 and in 1,154 of 2,000 runs. Both rows are noise. The third row goes the other way. The −0.111 between period order and consumption is small too, but the largest coefficient seen in the permuted sets is 0.061 and none of the 2,000 runs reaches the observed one. What speaks is not the size of the coefficient but where it stands against the permutation distribution. The fourth row shows that the margin alone is not enough either. The coefficient between tariff tier and consumption is +0.883 and its margin is zero, but the tier is computed from consumption: the limits of 10, 25 and 40 m³ determine it. The coefficient here measures a definition, not a discovery. #### Same Coefficient, Three Different Relationships ```python # same_coefficient.py — FICTIONAL: three different relationships that give the same coefficient. # The free parameter (noise weight, distance of the extreme point) is fitted to the target by bisection. N, TARGET = 557, 0.600 def driver(seed: int): r = rng(seed) return ([round(4 + 16 * r(), 2) for _ in range(N)], [r() + r() + r() - 1.5 for _ in range(N)]) def tune(f, upper: float) -> float: """f(p) gives the coefficient and varies monotonically with p; the target is approached by bisection.""" lower, increasing = 0.0, f(upper) > f(0.0) for _ in range(50): mid = (lower + upper) / 2 if (f(mid) < TARGET) == increasing: lower = mid else: upper = mid return (lower + upper) / 2 def mix(x, base, noise): w = tune(lambda p: corr(x, [t + p * g for t, g in zip(base, noise)]), 200.0) return [round(t + w * g, 2) for t, g in zip(base, noise)] x1, g1 = driver(SEED + 801) y1 = mix(x1, [1.4 * a for a in x1], g1) x2, g2 = driver(SEED + 802) y2 = mix(x2, [-0.16 * (a - 14.0) ** 2 + 40 for a in x2], g2) # peak at 14 m3 r3 = rng(SEED + 803) x3 = [round(4 + 16 * r3(), 2) for _ in range(N - 1)] y3 = [round(20 + 9 * (r3() + r3() + r3() - 1.5), 2) for _ in range(N - 1)] far = tune(lambda k: corr(x3 + [20 + k], y3 + [20 + 3 * k]), 400.0) x3, y3 = x3 + [round(20 + far, 2)], y3 + [round(20 + 3 * far, 2)] def slice_mean(x: list[float], y: list[float], d: int = 5) -> list[float]: s = sorted(zip(x, y)) return [mean([b for _, b in s[i * len(s) // d:(i + 1) * len(s) // d]]) for i in range(d)] rel = align([13, 11, 12]) print(rel(["relationship", "coef.", "rank coef."]) + " slice means (y)") for name, x, y in [("linear", x1, y1), ("curvilinear", x2, y2), ("leverage", x3, y3)]: print(rel([name, f"{corr(x, y):+.3f}", f"{rank_corr(x, y):+.3f}"]) + " " + " ".join(f"{v:5.1f}" for v in slice_mean(x, y))) print(f"\nin the leverage relationship the extreme point is ({x3[-1]}, {y3[-1]}); the remaining {N - 1} points have an " f"x range of 4-20") print(f"when that single point is removed the coefficient becomes {corr(x3[:-1], y3[:-1]):+.3f}") ``` ``` relationship coef. rank coef. slice means (y) linear +0.600 +0.594 6.6 12.6 15.4 21.4 24.9 curvilinear +0.600 +0.559 27.8 34.9 38.8 39.8 37.1 leverage +0.600 -0.052 20.1 20.0 19.9 19.4 21.7 in the leverage relationship the extreme point is (104.4, 273.2); the remaining 556 points have an x range of 4-20 when that single point is removed the coefficient becomes -0.057 ``` The coefficients of the three relationships are the same to the third decimal, and the sentence "correlation coefficient +0.600" does not tell the three apart. What does is the last column: when the driver column is divided into five slices and the mean of the second column is written for each slice, three different patterns come out. In the linear relationship the slice means climb smoothly from 6.6 to 24.9. In the curvilinear one they rise to 39.8 from 27.8 and return to 37.1: the second column increases up to a point, then decreases. The coefficient cannot see this turn, because it measures linear joint variation and the two sides of the turn partly cancel each other. Rank correlation cannot see it either; the difference between 0.559 and 0.594 does not count as a warning on its own. The third row is the harshest. The slice means hold along the line 20.1 — 20.0 — 19.9 — 19.4 and jump to 21.7 in the last slice: not a trend, but the trace of a single point. That point carries the coefficient — while 556 points sit between 4 and 20, one is at (104.4, 273.2) — and when that point is removed the coefficient falls to −0.057. Rank correlation says this from the start, because ranking does not care about the distance of a point. The summary of the measure is this: the coefficient reduces a two-dimensional shape to a single number and does not say what it reduced; five slice means tell all three relationships apart. #### Categorical and Ordinal Columns ```python # categorical.py — a cross tabulation counts two categorical columns together; region is # a nominal scale and the correlation coefficient is not defined for this pair. from collections import Counter cell = Counter((o["region"], o["tier"]) for o in READINGS) REGIONS = sorted({o["region"] for o in READINGS}) TIERS = [1, 2, 3, 4] c = align([9, 6] + [6] * 4 + [7] * 4) print(c(["region", "n"] + [f"t{k}" for k in TIERS] + [f"%t{k}" for k in TIERS])) for b in REGIONS + ["dataset"]: count = [sum(cell[(x, k)] for x in (REGIONS if b == "dataset" else [b])) for k in TIERS] print(c([b, sum(count)] + count + [f"{100 * s / sum(count):.1f}" for s in count])) # On an ordinal scale the numbers given to levels are a decision: when the coding changes # the coefficient changes; rank correlation does not, since it looks only at the ranking. COLUMNS = [("tier", [o["tier"] for o in READINGS], [("equal interval", {1: 1, 2: 2, 3: 3, 4: 4}), ("doubling", {1: 1, 2: 2, 3: 4, 4: 8}), ("tier midpoint", {1: 5.0, 2: 17.5, 3: 32.5, 4: 45.0})]), ("satisfaction", [o["satisfaction"] for o in READINGS], [("equal interval", {d: d for d in range(1, 6)}), ("doubling", {1: 1, 2: 2, 3: 4, 4: 8, 5: 16})])] o_ = align([13, 23, 12, 12], 2) print("\n" + o_(["column", "coding", "coef.", "rank coef."])) for col_name, col, codings in COLUMNS: coefs = [corr(CONSUMPTION, [m[d] for d in col]) for _, m in codings] for (name, m), value in zip(codings, coefs): print(o_([col_name, name, f"{value:+.4f}", f"{rank_corr(CONSUMPTION, [m[d] for d in col]):+.4f}"])) shift = 100 * (max(map(abs, coefs)) / min(map(abs, coefs)) - 1) print(f"{col_name}: coefficient moved {shift:.0f} percent in absolute value, " f"rank correlation did not move") ``` ``` region n t1 t2 t3 t4 %t1 %t2 %t3 %t4 central 600 30 352 198 20 5.0 58.7 33.0 3.3 east 557 21 281 225 30 3.8 50.4 40.4 5.4 north 921 42 661 210 8 4.6 71.8 22.8 0.9 south 645 59 538 48 0 9.1 83.4 7.4 0.0 west 476 88 377 11 0 18.5 79.2 2.3 0.0 dataset 3199 240 2209 692 58 7.5 69.1 21.6 1.8 column coding coef. rank coef. tier equal interval +0.8832 +0.8125 tier doubling +0.8417 +0.8125 tier tier midpoint +0.8800 +0.8125 tier: coefficient moved 5 percent in absolute value, rank correlation did not move satisfaction equal interval -0.0103 -0.0102 satisfaction doubling -0.0184 -0.0102 satisfaction: coefficient moved 79 percent in absolute value, rank correlation did not move ``` Region is a nominal scale: there is no order between north and east, and the correlation coefficient cannot be computed. The relationship is read from cell counts. In the west 79.2% of the 476 rows are in the second tier, and there are **no** rows in the fourth tier; in the east 5.4% of the 557 rows are in the fourth tier. In the whole dataset this share is 1.8%, so the share of the east is close to three times as large. The table's own margin shows in the cell counts. There are 8 rows in the fourth tier of the north and 11 in the third of the west; in a cell of 8 rows a single row moves the ratio by about one eighth. Unless the percentage is written together with how many observations it rests on, the table cannot be read. The lower table measures the decision of the ordinal scale. The tariff tier can be numbered 1, 2, 3, 4; a doubling sequence and the midpoints of the tiers give the same ranking too. The coefficient nevertheless changes as 0.8832 — 0.8417 — 0.8800. On satisfaction the effect is larger: 79% in absolute value. Rank correlation stays in place on both columns. The reason is in the definition of the scale. On an ordinal scale the **distance** between levels is not measured, only the order is known. Since the coefficient multiplies the level numbers, it counts the range between the first and second tier as equal to the range between the third and fourth; the scale does not support this. Rank correlation is independent of the coding. #### Summary - The coefficient is read with its margin: the coefficients +0.010 and −0.010 are exceeded in 1,189 and 1,154 of 2,000 permutations, and −0.111 in none; size alone does not tell them apart. - A coefficient whose margin is zero may not be a finding either: +0.883 measures a definition, since the tier is computed from consumption. - Three relationships that give the same +0.600 coefficient are told apart by slice means: a rise from 6.6 to 24.9, a curve that rises to 39.8 and returns to 37.1, and a flat line around 20. - In the leverage relationship a single point carries the coefficient; when that point is removed the coefficient falls to −0.057. - On an ordinal scale the coding moves the coefficient — by 79% in absolute value on satisfaction — and does not move rank correlation. #### Next Step This lesson looked at a single pair and wrote its margin next to the coefficient. In a real dataset the number of pairs does not stay at two: four numeric columns produce six pairs, six columns produce fifteen pairs. Looking at each pair separately is asking a new question each time, and as questions multiply, the chance that coincidence produces a high coefficient in at least one of them grows. The next lesson builds the correlation matrix, counts how many pairs give a high coefficient in columns produced entirely independently, and computes how much of the variation how many components carry. ### Multivariate Analysis Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/multivariate-examination A correlation matrix of six columns carries 15 pairs, and the threshold comes not from a distribution table but from permutation: in 2,000 permutations the 95th percentile of the coefficient distribution comes out as 0.0354, and 5 of the 15 pairs exceed this threshold — but three of the five exceed it by definition and the other two because of the subscriber number column. In columns generated completely independently of one another, when the number of pairs rises from 6 to 190 the number of pairs exceeding the threshold rises from 0 to 12, and the expected value is 9.5; the highest coefficient, 0.049, is clearly above the threshold. Dimensionality reduction is reduced to a single question: in the measurement network columns the first two components carry 49.3% of the variation, in the dataset generated from two hidden drivers 82.1%, and in independent columns only 35.0%. The previous lesson examined two columns together and wrote its margin of variation next to the coefficient of a single pair. In a real analysis, though, pairs are not asked one by one; a table is taken and the coefficients of all numeric columns are computed at the same time. Six columns produce 15 pairs, twenty columns 190 pairs. Each pair is a separate question, and each question gives chance one more shot. The unit of measure of this lesson is **the spurious relationship produced by the number of comparisons**. In columns generated completely independently of one another, the number of pairs exceeding the threshold is counted and it is shown how this number grows with the number of columns. The second measure is the gain of dimensionality reduction, and it is reduced to a single question: how many components carry how much of the variation. - **EA16.** Multivariate analysis looks at more than two columns at once. A **correlation matrix** with $k$ columns carries $k(k-1)/2$ pairs; the number of pairs grows quadratically with the number of columns. - **EA17.** Every cell in the matrix is a separate question. A threshold that is defensible for a single pair automatically inflates the **number of pairs exceeding the threshold** when it is applied to all pairs at once. - **EA18.** Here the threshold comes not from a distribution table but from **permutation**: a column is permuted within itself 2,000 times, and the 95th percentile of the resulting coefficient distribution is taken as the threshold. - **EA19. Dimensionality reduction** is covered in this lesson through a single question: **how many components carry how much of the variation**. Choice of method, rotation and interpretation of the components are not written. - **EA20.** Component shares are read from the **eigenvalues of the correlation matrix**; the sum of the eigenvalues equals the number of columns, and when written as shares the total is 100. If the columns do not repeat one another, the shares are distributed close to equal. #### The Matrix and the Threshold ```python # exploration.py — MODEL: the fictional measurement network of K03 and K04 is rebuilt # from the same seed. The dataset is fictional; no real subscriber, meter or region appears. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF PERIODS = ["2026-01", "2026-02", "2026-03"] TARIFFS = [10, 25, 40] REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t: int): """Linear congruential generator: seed written in the code, same output every run.""" x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value() -> float: s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u: float, w: list[float]) -> int: t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 _subscribers, _dropped = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] k = {"subscriber_id": 10001 + i, "region": b[0], "base": b[2], "household_size": pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]), "satisfaction": pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) + 1} (_dropped if r() < 0.046 else _subscribers).append(k) _valid = {a["subscriber_id"] for a in _subscribers} READINGS = [] for k in _subscribers + _dropped: r = rng(SEED + 7000 + k["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor(k["base"] * math.exp( (r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if k["subscriber_id"] not in _valid: r() continue READINGS.append({"subscriber_id": k["subscriber_id"], "region": k["region"], "period": PERIODS[d], "consumption_m3": t, "household_size": k["household_size"], "satisfaction": k["satisfaction"], "tier": 1 if t < TARIFFS[0] else 2 if t < TARIFFS[1] else 3 if t < TARIFFS[2] else 4}) r() def mean(v: list[float]) -> float: return sum(v) / len(v) def corr(x: list[float], y: list[float]) -> float: mx, my = mean(x), mean(y) sx = math.sqrt(sum((a - mx) ** 2 for a in x)) sy = math.sqrt(sum((b - my) ** 2 for b in y)) return sum((a - mx) * (b - my) for a, b in zip(x, y)) / (sx * sy) def align(w: list[int], left: int = 1): """The first `left` columns are left-aligned, the rest right-aligned.""" return lambda h: "".join(str(v).ljust(w[i]) if i < left else str(v).rjust(w[i]) for i, v in enumerate(h)) COLUMNS = [("consumption", [o["consumption_m3"] for o in READINGS]), ("tier", [float(o["tier"]) for o in READINGS]), ("household_size", [float(o["household_size"]) for o in READINGS]), ("satisfaction", [float(o["satisfaction"]) for o in READINGS]), ("period", [float(PERIODS.index(o["period"])) for o in READINGS]), ("subscriber_id", [float(o["subscriber_id"]) for o in READINGS])] K = len(COLUMNS) m = align([15] + [11] * K) print(m(["column"] + [a[:9] for a, _ in COLUMNS])) for i, (name, x) in enumerate(COLUMNS): print(m([name] + [f"{corr(x, y):+.3f}" for _, y in COLUMNS[:i + 1]])) RUNS = 2000 def threshold_95(x: list[float], seed: int) -> float: """The column is permuted against itself; the 95th percentile of the |coefficient| distribution is the threshold.""" r, n = rng(seed), len(x) mx = mean(x) d = [a - mx for a in x] total_sq = sum(a * a for a in d) values = [] for _ in range(RUNS): h = d[:] for i in range(n - 1, 0, -1): j = int(r() * (i + 1)) h[i], h[j] = h[j], h[i] values.append(abs(sum(a * b for a, b in zip(d, h)) / total_sq)) return sorted(values)[int(0.95 * RUNS)] THRESHOLD = threshold_95(COLUMNS[0][1], SEED + 901) pairs = [(COLUMNS[i][0], COLUMNS[j][0], corr(COLUMNS[i][1], COLUMNS[j][1])) for i in range(K) for j in range(i + 1, K)] exceeded = [c for c in pairs if abs(c[2]) >= THRESHOLD] print(f"\n{K} columns, {len(pairs)} pairs; in {RUNS} permutations the 95th percentile of the " f"|coefficient| distribution is {THRESHOLD:.4f}") print(f"pairs exceeding the threshold: {len(exceeded)} / {len(pairs)}") for a, b, v in sorted(exceeded, key=lambda c: -abs(c[2])): print(f" {a} / {b}: {v:+.3f}") ``` ``` column consumpti tier household satisfact period subscribe consumption +1.000 tier +0.883 +1.000 household_size +0.010 +0.007 +1.000 satisfaction -0.010 -0.015 -0.004 +1.000 period -0.111 -0.116 +0.000 -0.006 +1.000 subscriber_id +0.031 +0.041 +0.046 -0.001 -0.005 +1.000 6 columns, 15 pairs; in 2000 permutations the 95th percentile of the |coefficient| distribution is 0.0354 pairs exceeding the threshold: 5 / 15 consumption / tier: +0.883 tier / period: -0.116 consumption / period: -0.111 household_size / subscriber_id: +0.046 tier / subscriber_id: +0.041 ``` Six columns produced 15 pairs and the threshold came out as 0.0354. The threshold was not read from a table: the consumption column was permuted within itself 2,000 times, the coefficient was recomputed in every run, and the 95th percentile of the resulting distribution was taken as the threshold. The meaning of this threshold is narrow and is exactly this — when there is no relationship at all, in a **single pair**, a coefficient of this size is seen in five of a hundred runs. Five pairs exceed the threshold, but not all five are of the same kind. The +0.883 between consumption and tier is not a discovery but a definition; tier is computed from consumption with the limits 10, 25 and 40 m³. The −0.116 between tier and period is the same as the −0.111 between consumption and period, because tier is a rewriting of consumption — two lines count the same fact twice. What remains are the two pairs that include the subscriber number: the subscriber number is a record number and can be the cause of nothing. Five of the 15 pairs exceeded the threshold and the number of independent findings that will enter the report is at most one. The matrix itself does not produce a finding; it enlarges the number of questions to be asked, and the larger it grows, the easier it becomes to find a cell that exceeds the threshold. #### The Number of Pairs Produces Findings ```python # multiplicity.py — FICTIONAL: columns with no relation to one another are generated and the # same threshold is applied to all pairs. What is sought is the number of pairs that exceed it. def independent(k: int, n: int, seed: int) -> list[list[float]]: r = rng(seed) return [[r() + r() + r() for _ in range(n)] for _ in range(k)] INDEPENDENT = independent(20, len(READINGS), SEED + 902) _m = [mean(s) for s in INDEPENDENT] _d = [[a - t for a in s] for s, t in zip(INDEPENDENT, _m)] _n = [math.sqrt(sum(a * a for a in s)) for s in _d] pairwise = lambda i, j: sum(a * b for a, b in zip(_d[i], _d[j])) / (_n[i] * _n[j]) g = align([8, 8, 12, 12, 12]) print(g(["columns", "pairs", "exceeded", "expected", "largest"]) + " largest pair") for k in (4, 6, 8, 12, 20): combos = [(i, j) for i in range(k) for j in range(i + 1, k)] values = {(i, j): abs(pairwise(i, j)) for i, j in combos} best = max(values, key=values.get) print(g([k, len(combos), sum(1 for v in values.values() if v >= THRESHOLD), f"{0.05 * len(combos):.1f}", f"{values[best]:.3f}"]) + f" column{best[0] + 1} / column{best[1] + 1}") print(f"\nthreshold {THRESHOLD:.4f}; none of the columns is derived from another") ``` ``` columns pairs exceeded expected largest largest pair 4 6 0 0.3 0.021 column1 / column2 6 15 1 0.8 0.037 column2 / column6 8 28 2 1.4 0.048 column1 / column8 12 66 5 3.3 0.049 column6 / column12 20 190 12 9.5 0.049 column6 / column12 threshold 0.0354; none of the columns is derived from another ``` None of the columns in this table is derived from another; all of them come from separate draws and there is no relationship among them. Even so, 0 pairs in four columns, 1 in six columns, 2 in eight columns and **12** in twenty columns exceed the threshold. Where the number comes from is written in the definition of the threshold. The threshold is the value exceeded in five of a hundred runs for a single pair; when the same threshold is applied to 190 pairs, about 190 × 0.05 = 9.5 pairs are expected to exceed it. The observed 12 is near the expected 9.5. So the pairs exceeding the threshold are not a defect but a direct consequence of the definition of the threshold itself. The danger lies in pulling a single row out of this table. The sentence "A relationship of 0.049 was found between the sixth and the twelfth column, the threshold is 0.0354" contains correct numbers and is entirely made up: that pair was selected **because it is the largest** among 190 pairs. Next to a coefficient, not only its margin but also **how many pairs it was selected from** must be written. "One of a single pair" and "the largest of 190 pairs" are not the same sentence. One way to reduce the number of comparisons is not to raise the threshold but to lower the **number of columns**. This is the subject of the next section. #### How Many Components ```python # component.py — the single question of dimensionality reduction: how many components carry how # much of the variation. Shares are read from the eigenvalues of the correlation matrix. def eigenvalues(A: list[list[float]]) -> list[float]: """Eigenvalues of a symmetric matrix: rotated until the largest off-diagonal entry drops to zero.""" n, M = len(A), [s[:] for s in A] for _ in range(200): p, q, largest = 0, 1, 0.0 for i in range(n): for j in range(i + 1, n): if abs(M[i][j]) > largest: p, q, largest = i, j, abs(M[i][j]) if largest < 1e-12: # the rotation stops when the off-diagonal entries fall below this threshold break th = (M[q][q] - M[p][p]) / (2 * M[p][q]) t = (1.0 if th >= 0 else -1.0) / (abs(th) + math.sqrt(th * th + 1)) c = 1 / math.sqrt(t * t + 1) s = t * c for k in range(n): kp, kq = M[k][p], M[k][q] M[k][p], M[k][q] = c * kp - s * kq, s * kp + c * kq for k in range(n): pk, qk = M[p][k], M[q][k] M[p][k], M[q][k] = c * pk - s * qk, s * pk + c * qk return sorted((M[i][i] for i in range(n)), reverse=True) def shared_driver(seed: int) -> list[list[float]]: """FICTIONAL: six columns are generated from two hidden drivers and noise.""" r = rng(seed) f1 = [r() + r() + r() - 1.5 for _ in range(len(READINGS))] f2 = [r() + r() + r() - 1.5 for _ in range(len(READINGS))] WEIGHTS = [(1.0, 0.0), (0.9, 0.2), (0.8, 0.3), (0.1, 1.0), (0.2, 0.9), (0.3, 0.8)] return [[a * x + b * y + 0.55 * (r() + r() + r() - 1.5) for x, y in zip(f1, f2)] for a, b in WEIGHTS] DATASETS = [("measurement network", [x for _, x in COLUMNS]), ("shared driver", shared_driver(SEED + 903)), ("independent", independent(6, len(READINGS), SEED + 904))] b = align([20] + [9] * 6) print(b(["dataset"] + [f"c{i + 1}" for i in range(6)]) + " first two components") for name, S in DATASETS: eig = eigenvalues([[corr(a, c) for c in S] for a in S]) total = sum(eig) print(b([name] + [f"{100 * v / total:.1f}" for v in eig]) + f" {100 * (eig[0] + eig[1]) / total:.1f} percent") print(f"\nshares sum to 100 in every row; the sum of eigenvalues is the number of columns ({len(COLUMNS)})") ``` ``` dataset c1 c2 c3 c4 c5 c6 first two components measurement network 31.9 17.4 16.7 16.2 15.9 1.9 49.3 percent shared driver 55.2 26.9 4.9 4.7 4.3 4.0 82.1 percent independent 17.8 17.2 16.7 16.4 16.3 15.5 35.0 percent shares sum to 100 in every row; the sum of eigenvalues is the number of columns (6) ``` The three rows show when dimensionality reduction does work. In the dataset generated from a shared driver, six columns come from two hidden variables and the first two components carry 82.1% of the variation; writing two numbers instead of six columns means leaving 17.9% of the variation behind. In the independent columns the shares are almost equal, from 17.8 to 15.5, and the first two components carry only 35.0%: there is nothing to reduce, because the columns never repeat one another. The measurement network row lies between the two. The first component carries 31.9%, while the next four are squeezed between 17.4 and 15.9. The real information is in the last column: the share of the sixth component has dropped to 1.9%. This drop is a sign of repetition, and what the repetition is can be known from the first table — consumption and tier write the same thing twice. The share table does not say which columns repeat one another, but it gives **whether there is repetition** at a glance. The two sections meet here. Six columns produce 15 pairs; if the columns are reduced to two components the number of pairs falls to 1. This is the measurable gain of dimensionality reduction: it lowers the number of cells to look at and thus also lowers the chance that coincidence exceeds the threshold. Its cost is that the components no longer have a name like "consumption" or "household size". Choice of method and interpretation of the components are outside this lesson. #### Summary - A matrix of six columns carries 15 pairs; the threshold comes not from a distribution table but from permutation, and as the 95th percentile of the 2,000-run distribution it comes out as 0.0354. - Of the five pairs exceeding the threshold, three exceed it by definition (tier is computed from consumption) and two come from the subscriber number column; the number of independent findings is at most one. - In columns with no relation to one another, when the number of pairs rises from 6 to 190 the number of pairs exceeding the threshold rises from 0 to 12, and the expected value is 9.5. - Next to a coefficient, its margin is written together with **how many pairs it was selected from**; "one of a single pair" and "the largest of 190 pairs" are not the same claim. - Component shares tell the gain of dimensionality reduction in advance: the first two components carry 82.1% in the shared-driver dataset, 49.3% in the measurement network and 35.0% in the independent columns. #### Next Step The four lessons of this topic approached the dataset always through computation: a mean was taken, a distribution was counted, a coefficient was computed, a matrix was built. Their common assumption was that the values in the dataset are measured values. Yet during exploration the first things caught are often not found by computation but seen while looking at the table: the repetition of the same value with suspicious frequency, records piling up at an upper limit, field combinations that should never occur together, one period being markedly sparser than the others. The next lesson takes up this first inspection and counts which summary each caught problem moves and by how much. ### Visual Inspection of Data Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/visual-inspection-of-data A six-rule scan finds four problems in the raw record layer and each one distorts a different summary: the value 0.00 repeats 119 times and holds the dataset mean at 19.693, whereas it would be 20.454 if the rows that do not count as readings were left out; 20 records have piled up at the 45.00 m³ ceiling and the mean of the largest 1% appears as 44.663 instead of 46.597; 213 rows carry measured consumption in a zero-person household and consumption per person comes out as 11.204 or 10.617 depending on the decision; the row range between periods is 33.5% and the first-period mean rises from 20.644 to 20.868 when balanced subscribers are taken. No row is deleted. At the close of the topic the debt opened by exploration is counted: five breakdown fields produce 23 single and 209 pair cells, 232 separate subsets in total, and 33 of them have fewer than 30 rows. The first four lessons of this topic always approached the dataset through computation: a mean was taken, a distribution was counted, a coefficient was computed, a matrix was built. Their common assumption was that the values in the record were measured values. Yet during exploration the first things caught are often not found by computation but seen while looking at the table: the repetition of the same value with suspicious frequency, records piling up at a limit, field combinations that should never occur together, one period being markedly sparser than the others. The unit of measure of this lesson is **the number of problems the inspection catches and the summary each one distorts**. Four problems are found, a single question is asked for each — which number does it move and by how much — and no row is deleted. - **EA21.** Visual inspection here is **not drawing** but a scan turned into numbers: how many times each value is repeated, where records pile up, which field combination cannot exist, which slice is sparse. - **EA22. Heaping** is the repetition of the same value more than expected; its measure is the row share of the most frequent value. The same number appearing hundreds of times in a continuous quantity is not a measurement but a recording convention. - **EA23. Truncation** is the piling up of records at an upper or lower limit; its measure is the number of records at that limit. Truncation **cannot be repaired**: the truncated value is not in the record. The only things that can be counted are how many records were truncated and which summary was distorted by how much. - **EA24. An impossible combination** is two individually valid fields being invalid together; its measure is the number of rows that break the rule. - **EA25.** In this lesson no row is **deleted**. Cleaning was the job of the Data Collection and Preparation course and the measure there was the number of rows; here, for each problem two computations are printed side by side and the difference between them is written. #### Six Rules, Four Findings ```python # exploration.py — MODEL: the fictional measurement network of K03 and K04 is rebuilt # from the same seed. The dataset is fictional; no real subscriber, meter or region appears. import math from collections import Counter SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF PERIODS = ["2026-01", "2026-02", "2026-03"] TARIFFS = [10, 25, 40] REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t: int): """Linear congruential generator: seed written in the code, same output every run.""" x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value() -> float: s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u: float, w: list[float]) -> int: t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 _subscribers, _dropped = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] k = {"subscriber_id": 10001 + i, "region": b[0], "base": b[2], "household_size": pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]), "satisfaction": pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) + 1} (_dropped if r() < 0.046 else _subscribers).append(k) _valid = {a["subscriber_id"] for a in _subscribers} READINGS = [] for k in _subscribers + _dropped: r = rng(SEED + 7000 + k["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor(k["base"] * math.exp( (r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if k["subscriber_id"] not in _valid: r() continue READINGS.append({"subscriber_id": k["subscriber_id"], "region": k["region"], "period": PERIODS[d], "consumption_m3": t, "household_size": k["household_size"], "satisfaction": k["satisfaction"], "tier": 1 if t < TARIFFS[0] else 2 if t < TARIFFS[1] else 3 if t < TARIFFS[2] else 4}) r() def mean(v: list[float]) -> float: return sum(v) / len(v) def align(w: list[int], left: int = 1): """The first `left` columns are left-aligned, the rest right-aligned.""" return lambda h: "".join(str(v).ljust(w[i]) if i < left else str(v).rjust(w[i]) for i, v in enumerate(h)) # The raw record layer is a MODEL: the field device cannot read above 45.00 m3 and writes # the ceiling. This is the table that reaches the analyst; READINGS exists only for comparison. CEILING = 45.0 RECORDS = [dict(o, consumption_m3=min(o["consumption_m3"], CEILING)) for o in READINGS] T = [k["consumption_m3"] for k in RECORDS] top = Counter(T).most_common(1)[0] period_counts = Counter(k["period"] for k in RECORDS).values() duplicates = Counter((k["subscriber_id"], k["period"]) for k in RECORDS).values() RULES = [ ("heaping", "row share of the most frequent value (percent)", 100 * top[1] / len(T), 1.0, f"value {top[0]:.2f} repeated {top[1]} times"), ("truncation", "records at the largest value", float(Counter(T)[max(T)]), 5.0, f"{Counter(T)[max(T)]} records piled up at the {max(T):.2f} m3 upper limit"), ("impossible combination", "rows with consumption > 0 when household is 0", float(sum(1 for k in RECORDS if k["household_size"] == 0 and k["consumption_m3"] > 0)), 1.0, "measured consumption in a zero-person household"), ("period gap", "row range between periods", 100 * (max(period_counts) - min(period_counts)) / max(period_counts), 10.0, f"sparsest period {min(period_counts)} rows, densest {max(period_counts)} rows"), ("duplicate record", "rows of the same subscriber-period pair", float(max(duplicates)), 1.5, "-"), ("undefined value", "rows with negative consumption", float(sum(1 for t in T if t < 0)), 1.0, "-"), ] d = align([24, 46, 9, 10], 2) print(d(["problem", "measured", "value", "threshold"]) + " finding") problems = 0 for name, measure, value, threshold, text in RULES: hit = value >= threshold problems += hit print(d([name, measure, f"{value:.1f}", f"{threshold:.1f}"]) + " " + (text if hit else "-")) print(f"\n{len(RULES)} rules run, {problems} problems found; none deletes a row") ``` ``` problem measured value threshold finding heaping row share of the most frequent value (percent) 3.7 1.0 value 0.00 repeated 119 times truncation records at the largest value 20.0 5.0 20 records piled up at the 45.00 m3 upper limit impossible combination rows with consumption > 0 when household is 0 213.0 1.0 measured consumption in a zero-person household period gap row range between periods 33.5 10.0 sparsest period 838 rows, densest 1260 rows duplicate record rows of the same subscriber-period pair 1.0 1.5 - undefined value rows with negative consumption 0.0 1.0 - 6 rules run, 4 problems found; none deletes a row ``` The record layer is a model and what it models is clear: the field device cannot read above 45.00 m³ and writes the ceiling. This is the table that reaches the analyst; the untruncated values are kept only for comparison. The four findings are of four different kinds. The value 0.00 appears in 119 of 3,199 rows, that is 3.7%; the repetition of a number at this frequency in a continuous quantity does not come from measurement but from a recording convention that writes zero for a meter that could not be read. At the largest value 20 records have piled up: the exact same largest value appearing twenty times in a quantity with two decimals is the signature of the ceiling. 213 rows carry measured consumption in a zero-person household; the `household_size` field is valid on its own, the `consumption_m3` field is valid on its own, their combination is invalid. The row range between periods is 33.5%. The two rules that did not exceed the threshold give information too. The same subscriber and period pair does not appear more than twice in any row and there is no row that writes negative consumption. These are problems that were searched for and not found, and the fact that they were searched for is written down; when it is not written, the reader assumes the inspection never asked that question. #### The Number Each Problem Moves ```python # effect.py — the same question for each problem: which summary does it move and by how much. No # row is deleted; two computations are printed side by side and the difference is written. def per_person(records, household_zero_is_one: bool) -> float: v = [k["consumption_m3"] / (max(k["household_size"], 1) if household_zero_is_one else k["household_size"]) for k in records if household_zero_is_one or k["household_size"] > 0] return mean(v) def period_mean(records, period: str, subscribers=None) -> float: return mean([k["consumption_m3"] for k in records if k["period"] == period and (subscribers is None or k["subscriber_id"] in subscribers)]) complete = {a for a, c in Counter(k["subscriber_id"] for k in RECORDS).items() if c == len(PERIODS)} TOP = round(0.01 * len(T)) EFFECT = [ ("heaping", top[1], "dataset mean (m3)", mean(T), mean([t for t in T if t > 0]), "if rows where zero is written are not counted as readings"), ("truncation", Counter(T)[max(T)], "mean of the largest 1% (m3)", mean(sorted(T)[-TOP:]), mean(sorted(o["consumption_m3"] for o in READINGS)[-TOP:]), "had there been no ceiling (known only in the model)"), ("impossible combination", sum(1 for k in RECORDS if k["household_size"] == 0), "consumption per person (m3)", per_person(RECORDS, True), per_person(RECORDS, False), "household-0 rows left out"), ("period gap", len(RECORDS) - 3 * len(complete), "first-period mean (m3)", period_mean(RECORDS, PERIODS[0]), period_mean(RECORDS, PERIODS[0], complete), "only subscribers with all three periods"), ] e = align([24, 7]) print(e(["problem", "rows"]) + " " + "distorted summary".ljust(30) + "recorded".rjust(10) + "problem-free".rjust(13) + "diff".rjust(9)) for name, n, summary, recorded, clean, _ in EFFECT: print(e([name, n]) + " " + summary.ljust(30) + f"{recorded:10.3f}{clean:13.3f}{clean - recorded:+9.3f}") print() for name, _, _, recorded, clean, explanation in EFFECT: print(f"{name:<24}{explanation}; {100 * (clean - recorded) / recorded:+.1f} percent") ``` ``` problem rows distorted summary recorded problem-free diff heaping 119 dataset mean (m3) 19.693 20.454 +0.761 truncation 20 mean of the largest 1% (m3) 44.663 46.597 +1.934 impossible combination 222 consumption per person (m3) 11.204 10.617 -0.586 period gap 685 first-period mean (m3) 20.644 20.868 +0.224 heaping if rows where zero is written are not counted as readings; +3.9 percent truncation had there been no ceiling (known only in the model); +4.3 percent impossible combination household-0 rows left out; -5.2 percent period gap only subscribers with all three periods; +1.1 percent ``` Four problems distort four different summaries and none of them substitutes for another. The zero heaping holds the dataset mean down: instead of 19.693, it would be 20.454 if the rows that do not count as readings were left out, that is 3.9%. This number is inside every region mean that is reported too, because the zeros are distributed over all regions. The number that truncation moves is not the mean but the upper end. The dataset mean is almost unaffected by the ceiling, but the mean of the largest 1% stays at 44.663; without the ceiling it would be 46.597. This second number is known only in the model. In the field a truncated value is not in the record and cannot be brought back; the only thing that can be done is to write the sentence "20 of the 32 highest readings rest on the ceiling" next to the number. In a truncated column the largest value, every summary above the quartiles and every decision concerning the upper end rest on a clipped number, not a measured one. The impossible combination is different: here there are two defensible decisions and both produce a number. If a zero-person household is counted as 1, consumption per person is 11.204 m³; if those 222 rows are left out it comes out as 10.617 m³. The difference is 5.2%, and the data does not say which is right; what says it is how many people the household actually has, and that is not in the record either. The number that the period gap moves is the smallest but its direction is instructive. The first-period mean is 20.644 m³ in the record; when only the subscribers with all three periods are taken it becomes 20.868 m³. The 0.224 m³ between them is not a behavior difference but a **composition** difference: the subscribers who drop out early consumed below the average in the first period. A sentence tying the decline between periods to behavior cannot be written without taking this difference into account. #### The Debt Exploration Opens ```python # debt.py — the number of breakdowns opened in the exploration stage: each cell is a subset that # can be looked at, and each look is a debt that enlarges the margin of the test. FIELDS = [("region", 5), ("period", 3), ("tier", 4), ("satisfaction", 5), ("household_size", 6)] single = sum(k for _, k in FIELDS) field_pairs = [(i, j) for i in range(len(FIELDS)) for j in range(i + 1, len(FIELDS))] double = sum(FIELDS[i][1] * FIELDS[j][1] for i, j in field_pairs) small, empty = 0, 0 for i, j in field_pairs: a, b = FIELDS[i][0], FIELDS[j][0] count = Counter((k[a], k[b]) for k in RECORDS) empty += FIELDS[i][1] * FIELDS[j][1] - len(count) small += sum(1 for v in count.values() if v < 30) print(f"{len(FIELDS)} breakdown fields, {len(field_pairs)} field pairs") print(f"single breakdown cells : {single}") print(f"double breakdown cells : {double}") print(f"total viewable subsets: {single + double}") print(f"{small} of them have fewer than 30 rows, {empty} are empty") print(f"the dataset is 1 table; exploration can extract {single + double} separate summaries") ``` ``` 5 breakdown fields, 10 field pairs single breakdown cells : 23 double breakdown cells : 209 total viewable subsets: 232 33 of them have fewer than 30 rows, 2 are empty the dataset is 1 table; exploration can extract 232 separate summaries ``` Five breakdown fields open 23 cells on their own, and their pair combinations 209 cells more. In total 232 subsets, all coming out of the same table, and for each of them a mean, a rate or a difference can be computed. 33 of them have fewer than 30 rows, and in a small cell the margin grows: the mean of a subset of twenty rows moves noticeably when a single record comes or goes. This number is not a defect but the nature of exploration. The exploration stage exists precisely to walk through these cells, look at which one draws attention and produce questions from there. But every cell walked is a debt added to the margin of the test to be done later — because the question to be tested was selected **because it was the one that drew the most attention** among the 232 cells. #### Summary - The six-rule scan found four problems: heaping (the value 0.00 in 3.7% of the rows), truncation (20 records at 45.00 m³), impossible combination (213 rows) and period gap (33.5%). - Heaping holds the dataset mean at 19.693; if the zeros are not counted as readings it becomes 20.454, so a range of 3.9% is inside every region mean. - Truncation hardly moves the mean at all but shows the mean of the largest 1% as 44.663 instead of 46.597 and **cannot be repaired**; the only thing that can be counted is how many records rest on the ceiling. - The impossible combination leaves two defensible decisions: consumption per person is 11.204 or 10.617 m³, a difference of 5.2%. - The period gap carries the first-period mean from 20.644 to 20.868; this is not a behavior but a composition difference. #### Next Step In five lessons this topic did one thing: it got to know the dataset and produced questions. It was seen that the 4.463 m³ difference between regions could not be produced by permutation, it was measured that the east distribution is skewed and that the 30 rows outside the fences that were not dropped moved the mean by 0.480 m³, it was shown that three different relationships giving the same coefficient separate by their slice means, it was counted that 12 of the 190 pairs in twenty independent columns exceed the threshold, and it was written that four record problems distorted four different summaries. The common gap is this: **none of these questions was tested**. Exploration says that a difference seems to exist, it does not say whether that difference is worth reporting; it says a summary is skewed, it does not decide which summary will be reported. Moreover, exploration itself left a debt: 232 breakdown cells were walked and the question to be tested was selected from among them. The next topic, Applied Statistics, begins by paying this debt — first it measures the choice of the summary to be reported and the number that choice moves, then it builds the difference's own interval and test. ## Applied Statistics Testing with a permutation test whether a difference is real or noise and reporting the result with its margin of variation: the decision between mean and median moving the reported number, what box and density summaries squeeze into five numbers and what they hide, what a confidence interval says and does not say, test choice and assumption checking, the number of decisions an experiment set yields through the garden of forking paths, the interpretive limit of a regression coefficient and how much of the effect it explains, the retention rate moved by the cohort definition, selection bias, and counting the difference found by trying many things. ### Choosing Descriptive Summaries Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/choosing-descriptive-summaries Three defensible summaries come out of the same 557 readings: the mean is 24.674 m³, the median 24.120 m³ and the 10% trimmed mean 24.581 m³, with a largest difference of 0.554 m³ among them. Bootstrap measures the margin of variation of all three and the widths come out as 1.552, 1.762 and 1.458 m³: the difference the choice moves is a third of the margins, and all three intervals contain the 25 m³ tariff step. When the same computation drops from the 557-reading group to the 142-reading group, the choice difference rises to 1.679 m³ and the margin width to 2.917 m³, and this time the decision cannot be read independently of the margin. The exploration phase did not produce a decision. Over five lessons the distributions were opened, the breakdowns were tried and the relationships were scanned, and what remained is a list of questions; none of the questions on the list has been tested yet. Every breakdown that was looked at during exploration also left a debt to testing, and that debt will be counted at the end of the topic. The first step of testing is not, as one might think, choosing a test. Before it comes a step that is often skipped: with which summary the number to be tested will be written. When the "mean consumption" of the east region's readings is asked for, there are at least three defensible answers, and the three give three different numbers from the same dataset. The measure of this lesson is how many cubic meters this choice moves the number in the report, and how much that movement amounts to next to the number's own margin. - **AS1. A descriptive summary is a decision.** For the same dataset the mean, the median and the trimmed mean give three different numbers; none of them is an approximation of another, and the three answer three different questions. - **AS2. Every summary has a margin of variation**, and the margin is measured with **bootstrap**: the dataset is drawn from itself many times, at the same size and with replacement, and the summary is recomputed at every draw. The draw count is written down; the margin depends on that count too. - **AS3.** A summary with a narrow margin is not the right summary by itself. Narrowness does not change **what** the summary answers; it only says how much that answer moves from sample to sample. - **AS4.** The difference the choice moves is read next to the margin. An argument over a difference smaller than the margin is an argument over noise. - **AS5.** The smaller the group, the larger the margin. Next to every reported summary is written **how many observations it rests on**; the margin of a summary without an `n` cannot be read. #### Three Summaries, Three Numbers ```python # dataset.py — MODEL: the measurement network of K03-K04 is regenerated with the same seed. # A long-format reading table of 3,199 rows, 1,329 subscribers, three periods, five regions. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF PERIODS = ["2026-01", "2026-02", "2026-03"] REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t): x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value(): s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u, w): t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 _kept, _dropped = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]) pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) (_dropped if r() < 0.046 else _kept).append( {"subscriber_id": 10001 + i, "region": b[0], "base": b[2]}) VALID = {a["subscriber_id"] for a in _kept} READINGS = [] for a in _kept + _dropped: r = rng(SEED + 7000 + a["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor( a["base"] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if a["subscriber_id"] not in VALID: r() continue READINGS.append({"subscriber_id": a["subscriber_id"], "region": a["region"], "period": PERIODS[d], "consumption_m3": t}) r() def mean(v): return sum(v) / len(v) def median(v): s = sorted(v) n = len(s) return s[n // 2] if n % 2 else (s[n // 2 - 1] + s[n // 2]) / 2 def trimmed(v, p=0.10): s = sorted(v) k = int(len(s) * p) return mean(s[k:len(s) - k]) east = [o["consumption_m3"] for o in READINGS if o["region"] == "east"] U = [mean(east), median(east), trimmed(east)] print(f"fictional dataset (seed {SEED}): {len(READINGS)} rows, {len(VALID)} subscribers") print(f"east region: {len(east)} readings, {sum(1 for x in east if x == 0.0)} zero readings") print(f"mean {U[0]:.3f} | median {U[1]:.3f} | trimmed(10%) {U[2]:.3f}") print(f"largest difference among the three {max(U) - min(U):.3f} m3") ``` ``` fictional dataset (seed 20260218): 3199 rows, 1329 subscribers east region: 557 readings, 20 zero readings mean 24.674 | median 24.120 | trimmed(10%) 24.581 largest difference among the three 0.554 m3 ``` The dataset is fictional and is regenerated with the same seed from the measurement network of K03–K04; the 3,199 rows and 1,329 subscribers are the same numbers the previous course left behind. The east region's mean of 24.674 m³ also comes from there, and it is the number asked about at the opening of this course. Now two more numbers stand next to it. The median is 24.120 m³: the value that stands in the middle when the readings are sorted, that is, the number below which half of the dataset stays. The 10% trimmed mean is 24.581 m³: the mean of the 447 values left after the smallest 55 and the largest 55 readings are removed. All three are correct, all three come from the same 557 readings, and the largest difference among them is 0.554 m³. That the mean is larger than the median is a piece of information: the distribution is **skewed** to the right, that is, a few large readings pull the mean up. The twenty zero readings, on the other hand, work in the opposite direction and press the mean down. These two effects partly cancel each other inside the same number, and someone who looks only at the mean cannot see either of them. #### Measuring the Margin ```python # margin.py — the margin of variation is measured with bootstrap: the dataset is drawn from itself # 1,000 times at the same size and with replacement, the summary is recomputed at every draw, # and the middle 95% of the 1,000 values is taken. DRAWS = 1000 def percentile(s, p): i = p * (len(s) - 1) lower = int(i) upper = min(lower + 1, len(s) - 1) return s[lower] + (s[upper] - s[lower]) * (i - lower) def margin(v, summary, seed, draws=DRAWS): rand = rng(seed) n = len(v) d = sorted(summary([v[int(rand() * n)] for _ in range(n)]) for _ in range(draws)) return percentile(d, 0.025), percentile(d, 0.975) E = [12, 10, 10, 10, 10] y = lambda h: "".join(str(v).rjust(E[k]) if k else str(v).ljust(E[0]) for k, v in enumerate(h)) SUMMARIES = (("mean", mean), ("median", median), ("trimmed", trimmed)) print(y(["summary", "number", "lower", "upper", "width"])) for name, f in SUMMARIES: a, b = margin(east, f, SEED + 101) print(y([name, round(f(east), 3), round(a, 3), round(b, 3), round(b - a, 3)])) print(f"\n{len(east)} observations, {DRAWS} draws, all three on the same draws") print("the 25 m3 tariff step stays inside all three intervals") ``` ``` summary number lower upper width mean 24.674 23.889 25.441 1.552 median 24.12 23.248 25.01 1.762 trimmed 24.581 23.876 25.334 1.458 557 observations, 1000 draws, all three on the same draws the 25 m3 tariff step stays inside all three intervals ``` Bootstrap answers a single question: if another set of readings of the same size had come from the same network, what would this summary have produced. Since there is no other dataset at hand, the dataset itself is put in the place of the population, 557 values are drawn from it with replacement, the summary is recomputed, and this is repeated 1,000 times. The middle 95% of the resulting 1,000 values is the summary's **margin of variation**. The three summaries were computed on the same draws; the difference between them comes from the summary, not from the luck of the draw. The widths lie between 1.458 and 1.762 m³. The difference the choice moved, on the other hand, was 0.554 m³, that is, a third of even the narrowest margin. The argument over which of the three summaries to write is, in this dataset, an argument that stays inside the margin. In the report this amounts to a single sentence: is the east region's consumption below the 25 m³ tariff step. All three numbers are below it, but all three intervals contain 25. The claim of staying below the step cannot be defended without writing the margin, whichever summary is chosen. In a report that does write the margin, the claim softens by itself and takes its correct form: the three observed numbers are all below the step, but this data cannot rule out that another set of readings from the same network would exceed the step. #### Which Summary Answers What ```python # choice.py — what the summary choice writes into the report: total volume and trimmed rows. k = int(len(east) * 0.10) print(f"mean x n = {mean(east) * len(east):.2f} m3 reported volume") print(f"median x n = {median(east) * len(east):.2f} m3 same product, different number") print(f"difference between the two volumes {(mean(east) - median(east)) * len(east):.2f} m3") print(f"trimmed mean leaves {2 * k} readings out; " f"{sum(1 for x in sorted(east)[:k] if x == 0.0)} of them are zero readings") nonzero = [x for x in east if x > 0.0] print(f"if zero readings are not counted at all the mean becomes {mean(nonzero):.3f} m3") ``` ``` mean x n = 13743.34 m3 reported volume median x n = 13434.84 m3 same product, different number difference between the two volumes 308.50 m3 trimmed mean leaves 110 readings out; 20 of them are zero readings if zero readings are not counted at all the mean becomes 25.593 m3 ``` That the three summaries answer three different questions is visible here. The mean carries the **total** of the dataset: multiplied by 557 it gives 13,743.34 m³, and this is the volume that passed through the network itself. The median gives 13,434.84 m³ in the same product, and this number is the total of nothing; the 308.50 m³ in between is the measure of the median not answering the volume question. When the volume to be billed is asked, the median is the wrong summary. The trimmed mean answers the third question: a center that resists extreme readings. It does so by leaving 110 readings out, and 20 of the ones it leaves out are zero readings. This is the summary with the narrowest margin, but the source of the narrowness is not the dataset itself, it is that 20% of the dataset does not enter the computation. A narrow margin here is not an advantage, it is a sign of a narrowed coverage. The last line shows what the coverage decision does by itself. If the zero readings are not counted at all, the mean becomes 25.593 m³ and exceeds the tariff step. What the zero readings are was dealt with in the cleaning lesson of M26/K03; no elimination is done here, it is only recorded that the coverage decision moved the summary by 0.919 m³ and that this movement is larger than the choice difference. #### When the Group Shrinks ```python # small.py — the same three summaries and the same margin, recomputed on a small breakdown. sub = [o["consumption_m3"] for o in READINGS if o["region"] == "east" and o["period"] == "2026-03"] for name, v in (("east / all", east), ("east / 2026-03", sub)): s = [f(v) for _, f in SUMMARIES] g = [b - a for a, b in (margin(v, f, SEED + 101) for _, f in SUMMARIES)] print(f"{name}: n={len(v)}, choice difference {max(s) - min(s):.3f} m3, " f"margin width between {min(g):.3f} and {max(g):.3f} m3") print(f" mean {s[0]:.3f} | median {s[1]:.3f} | trimmed {s[2]:.3f}") ``` ``` east / all: n=557, choice difference 0.554 m3, margin width between 1.458 and 1.762 m3 mean 24.674 | median 24.120 | trimmed 24.581 east / 2026-03: n=142, choice difference 1.679 m3, margin width between 2.527 and 2.917 m3 mean 22.394 | median 20.715 | trimmed 22.112 ``` The same region, a single period: 142 readings. The choice difference rose from 0.554 m³ to 1.679 m³, that is, to three times. The widths of the margins grew too and settled in a band of 2.5–2.9 m³. The two increases are read together: in the small breakdown the summaries both move away from each other and the variation of each one grows. The source of the divergence is visible. In this period the mean is 22.394 m³ and the median is 20.715 m³, and the range of 1.679 m³ in between comes from a few large readings pulling the mean in the single-period breakdown. The trimmed mean, at 22.112 m³, stays next to the mean, because the readings doing the pulling are more than the trimmed 10%. For the report the result is this: in the 557-reading group the summary choice is a detail that stays inside the margin, whereas in the 142-reading group the choice difference exceeds half of the margin and the decision can no longer be read independently of the margin. In both, what must be written is the same — next to the number, the `n` and the margin. #### The Margin Is Also a Number ```python # draws.py — the margin comes from a run, not from a calculation; when the run count changes, the # bounds of the margin move too. In this course the upper limit is 1,000 draws. for c in (200, 500, 1000): a, b = margin(east, mean, SEED + 101, c) print(f"{c:5d} draws: {a:.3f} - {b:.3f}, width {b - a:.3f} m3") ``` ``` 200 draws: 23.955 - 25.406, width 1.451 m3 500 draws: 23.949 - 25.426, width 1.477 m3 1000 draws: 23.889 - 25.441, width 1.552 m3 ``` The margin comes not from a formula but from a run, and when the run count changes the bounds of the margin move too. With two hundred draws the width is 1.451 m³, with a thousand draws it is 1.552 m³; the 0.101 m³ in between is the margin's own uncertainty. The ends of the interval are the part that moves the most, because the ends rest on the outermost twenty-fifth values of the 1,000, and in that region there are few draws. In this course the draw count does not exceed 1,000. The limit is not a convenience but a record: a reported margin is not complete without writing how many draws it came from. The same rule will hold for the permutation count in later lessons, and there the run count will directly determine the resolution of the margin. #### Summary - The same 557 readings give three defensible summaries: the mean is 24.674 m³, the median 24.120 m³, the 10% trimmed mean 24.581 m³; the largest difference the choice moves is 0.554 m³. - Bootstrap puts the dataset in the place of the population, draws from it 1,000 times with replacement and measures the margin of each summary: the widths come out as 1.552, 1.762 and 1.458 m³, that is, three times the choice difference. - The summary with the narrowest margin is the trimmed mean, but the narrowness comes from 110 readings not entering the computation; a narrow margin does not mean the right summary, it means the coverage narrowed. - The three summaries answer three different questions: the mean carries the total and its product with 557 is 13,743.34 m³, while the median's same product, 13,434.84 m³, corresponds to no volume. - When the breakdown drops from 557 readings to 142, the choice difference rises to 1.679 m³ and the margin width to 2.917 m³; this is why the number of observations it rests on is written next to every summary. #### Next Step A single number was chosen and its margin was written, but a single number does not describe a distribution. The real question in a report is often the comparison of two groups, and by habit this comparison is made over five numbers: the minimum, the first quartile, the median, the third quartile, the maximum. The next lesson shows that these five numbers do not determine a distribution — two datasets whose five numbers are identical to each other are built, the bin counts are printed side by side, and it comes out in numbers that one has a single peak and the other has two. The notch interval of the median is again produced with bootstrap. ### Distribution Comparison Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/distribution-comparison The 921 readings of the north region come down to five numbers, and a second dataset that gives the same five numbers is built: both box summaries are 0.00 | 15.74 | 19.66 | 24.63 | 45.47, the interquartile range is 8.89 m³ and the fences are 2.405 and 37.965 m³. Bin counts separate the two — in north a single peak in the 16–20 m³ bin, in the twin two peaks in the 12–16 and 24–28 m³ bins with the median left in the trough. The notch of the same median comes out as 0.920 m³ and 3.146 m³ with bootstrap, and the number of observations outside the fences varies between 50 and 35. The previous lesson chose a single number and wrote its margin. Most of the questions in a report ask not for a single number but for a comparison: is the consumption distribution of two regions different from each other. By habit this comparison is made over five numbers — the minimum, the first quartile, the median, the third quartile, the maximum — and the five numbers are put side by side as a box summary. In this lesson the box summary and the **density display** are not drawn; both are built as data structures and printed as numbers. Drawing, axis and design decisions are the job of M26/K06. The question here comes earlier: do five numbers determine a distribution. The measure is how many peaks two datasets whose five numbers are identical to each other give in their bin counts. - **AS6. The box summary is five numbers and two fences.** The interquartile range is the difference between the third quartile and the first quartile; the fences stand outside by 1.5 times this range. In M26/K03 these five numbers were an **elimination** decision; here there is no elimination, the numbers are used to describe the shape of the distribution. - **AS7. Five numbers do not determine a distribution.** There are several datasets that give the same five numbers, and the difference among them is not reflected in the box summary at all. - **AS8. A bin count is the numeric counterpart of density.** The value axis is divided into equal-width **bin**s and the observations falling into each bin are counted; the number of peaks is read from these counts. - **AS9. The notch is the margin of the median** and is produced with bootstrap: the median is recomputed in 1,000 draws, and the middle 95% of the resulting values is the notch. - **AS10.** The width of the notch depends on the number of observations **and** on the density around the median. The same median gives a much wider notch in a dataset where the neighborhood of the median is sparse. - **AS11.** Whether two notches intersect is a display cue, not a test; what the intersection of intervals means is the subject of the next lesson. #### The Box Summary Is a Data Structure ```python # dataset.py — MODEL: the measurement network of K03-K04 is regenerated with the same seed. # A long-format reading table of 3,199 rows, 1,329 subscribers, three periods, five regions. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF PERIODS = ["2026-01", "2026-02", "2026-03"] REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t): x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value(): s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u, w): t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 _kept, _dropped = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]) pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) (_dropped if r() < 0.046 else _kept).append( {"subscriber_id": 10001 + i, "region": b[0], "base": b[2]}) VALID = {a["subscriber_id"] for a in _kept} READINGS = [] for a in _kept + _dropped: r = rng(SEED + 7000 + a["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor( a["base"] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if a["subscriber_id"] not in VALID: r() continue READINGS.append({"subscriber_id": a["subscriber_id"], "region": a["region"], "period": PERIODS[d], "consumption_m3": t}) r() def percentile2(v, p): s = sorted(v) i = p * (len(s) - 1) a = int(i) b = min(a + 1, len(s) - 1) return s[a] + (s[b] - s[a]) * (i - a) def five_numbers(v): return [min(v), percentile2(v, 0.25), percentile2(v, 0.50), percentile2(v, 0.75), max(v)] def box(v): b = five_numbers(v) iqr = b[3] - b[1] return b, iqr, b[1] - 1.5 * iqr, b[3] + 1.5 * iqr north = sorted(o["consumption_m3"] for o in READINGS if o["region"] == "north") b, iqr, fence_low, fence_high = box(north) print(f"north: {len(north)} readings") print("five-number summary: " + " | ".join(f"{x:.2f}" for x in b)) print(f"interquartile range {iqr:.2f} m3, fences {fence_low:.3f} and {fence_high:.3f} m3") print(f"{sum(1 for x in north if x < fence_low or x > fence_high)} readings outside the fences") ``` ``` north: 921 readings five-number summary: 0.00 | 15.74 | 19.66 | 24.63 | 45.47 interquartile range 8.89 m3, fences 2.405 and 37.965 m3 50 readings outside the fences ``` The box summary here is not a drawing but a record of seven numbers: five numbers, the interquartile range and two fences. The 921 readings of the north region came down to these seven. The median is 19.66 m³, the interquartile range is 8.89 m³, the fences are between 2.405 and 37.965 m³, and 50 readings stay outside these limits. That the smallest value of the dataset is zero suggests that most of the readings below the lower fence are zero readings. M26/K03 had dealt with these readings as a cleaning decision; here none of them is eliminated, only how many of them stay outside the fences is counted. The fence is not an elimination threshold but a measure of shape. Seven numbers take little space in a report and speed up the comparison. Their price is measured in the next section: these seven numbers say much less than one would think about how the 921 readings are arranged. #### The Same Five Numbers, Another Dataset ```python # twin.py — SETUP: the sorted values of north are pushed, inside each quarter slice, toward the end # of the slice. Since the 0th, 230th, 460th, 690th and 920th values in the sorted array stay in # place, the five-number summary is preserved exactly, the inner distribution changes. def push(part, end, p=4): a, z = part[0], part[-1] u = [(x - a) / (z - a) for x in part] return [z - (z - a) * (1 - t) ** p if end == "upper" else a + (z - a) * t ** p for t in u] C = [0, 230, 460, 690, len(north) - 1] twin = (push(north[C[0]:C[1] + 1], "upper")[:-1] + push(north[C[1]:C[2] + 1], "lower")[:-1] + push(north[C[2]:C[3] + 1], "upper")[:-1] + push(north[C[3]:], "lower")) b2, iqr2, low2, high2 = box(twin) print(f"twin dataset: {len(twin)} values (fictional, derived from the sorted values of north)") print("five-number summary: " + " | ".join(f"{x:.2f}" for x in b2)) print(f"interquartile range {iqr2:.2f} m3, fences {low2:.3f} and {high2:.3f} m3") print(f"{sum(1 for x in twin if x < low2 or x > high2)} values outside the fences") print(f"are the two box summaries the same: {[round(x, 2) for x in b] == [round(x, 2) for x in b2]}") ``` ``` twin dataset: 921 values (fictional, derived from the sorted values of north) five-number summary: 0.00 | 15.74 | 19.66 | 24.63 | 45.47 interquartile range 8.89 m3, fences 2.405 and 37.965 m3 35 values outside the fences are the two box summaries the same: True ``` The twin dataset is fictional and how it was built is in the code: the sorted array is divided into four quarter slices, and the values inside each slice are pushed toward one end of the slice. The push does not break the ordering and the values that coincide with the slice ends stay in place, so the five numbers are preserved exactly. Because the push direction of the lower and the upper half was chosen opposite, the values pile up on the two sides of the median. The two box summaries are the same: the same median, the same quartiles, the same interquartile range, the same fences. A report that puts the box summaries side by side cannot report any difference between the two datasets. The only thing it can report is in the last line, and that is not among the five numbers either: the number of observations outside the fences has dropped from 50 to 35, that is, the same fences leave a different number of values outside in the two datasets. This is not a defect of the box summary but its scope. The five numbers report **five points of the ordering**; they say nothing about how the values are distributed between those points, and they were not designed to. The information the box summary carries is independent of the position in the ordering of all but four of the 921 readings: every arrangement inside a slice gives the same seven numbers. It must be written once more that the twin dataset is a fictional setup with no counterpart in the real network. Its purpose is not to imitate a region but to show the size of the difference the box summary cannot distinguish. The same effect appears in real datasets in a more moderate form, and it usually arises from two different subscriber groups being collected into a single row. #### Bin Counts ```python # bins.py — the density display as numbers: the value axis is divided into 4 m3 bins, # the observations in each bin are counted. The text bar is one mark per ten observations. WIDTH, N = 4.0, 12 def counts(v): c = [0] * N for x in v: c[min(int(x / WIDTH), N - 1)] += 1 return c def peaks(c): return [i for i in range(1, N - 1) if c[i] > c[i - 1] and c[i] > c[i + 1] and c[i] >= 0.3 * max(c)] cn, ct = counts(north), counts(twin) for i in range(N): print(f"{i * WIDTH:5.0f}-{(i + 1) * WIDTH:<3.0f} {cn[i]:4d} {'#' * round(cn[i] / 10):<34}" f" {ct[i]:4d} {'#' * round(ct[i] / 10)}") print(f"peak count: north {len(peaks(cn))}, twin {len(peaks(ct))}") print(f"the peaks of the twin are in the {[int(i * WIDTH) for i in peaks(ct)]} m3 bins, " f"while the median {b2[2]:.2f} m3 is in the trough") ``` ``` 0-4 32 ### 32 ### 4-8 0 0 8-12 41 #### 0 12-16 175 ################## 318 ################################ 16-20 234 ####################### 114 ########### 20-24 182 ################## 110 ########### 24-28 130 ############# 329 ################################# 28-32 66 ####### 11 # 32-36 37 #### 4 36-40 16 ## 0 40-44 5 2 44-48 3 1 peak count: north 1, twin 2 the peaks of the twin are in the [12, 24] m3 bins, while the median 19.66 m3 is in the trough ``` The two datasets that the box summary could not separate were separated with twelve pairs of numbers. In north the counts rise to a single peak in the 16–20 m³ bin and fall off to both sides. In the twin there are two peaks: 318 values in the 12–16 m³ bin, 329 values in the 24–28 m³ bin, and between them stands a trough dropping to 114 and 110. The place of the trough matters: the median is 19.66 m³ and falls exactly into this trough. In the twin dataset the values closest to the median are the sparsest ones. The sentence "the typical subscriber consumes 19.66 m³" is defensible for north, but for the twin dataset it is true and still misleading — there is almost nobody around that value. The bin width is a decision, and this decision itself moves the number: if 10 m³ had been chosen instead of 4 m³, the two peaks would fall into a single bin and the difference would again be invisible. The visual effect and the choice of bin width are the subject of M26/K06; the only thing recorded here is that the chosen width is written under the output. #### Notch Interval ```python # notch.py — the margin of the median: drawn 1,000 times with replacement, the median is # recomputed at every draw, and the middle 95% of the resulting values is the notch interval. DRAWS = 1000 def percentile(s, p): i = p * (len(s) - 1) a = int(i) z = min(a + 1, len(s) - 1) return s[a] + (s[z] - s[a]) * (i - a) def notch(v, seed, draws=DRAWS): rand = rng(seed) n = len(v) d = sorted(percentile2([v[int(rand() * n)] for _ in range(n)], 0.50) for _ in range(draws)) return percentile(d, 0.025), percentile(d, 0.975) central = [o["consumption_m3"] for o in READINGS if o["region"] == "central"] for name, v in (("north", north), ("twin", twin), ("central", central)): a, z = notch(v, SEED + 211) print(f"{name:7s} n={len(v):4d} median {percentile2(v, 0.50):6.3f} " f"notch {a:.3f} - {z:.3f} width {z - a:.3f}") na, nz = notch(north, SEED + 211) ca, cz = notch(central, SEED + 211) print(f"do the notches of north and central intersect: {na < cz and ca < nz}") ``` ``` north n= 921 median 19.660 notch 19.230 - 20.150 width 0.920 twin n= 921 median 19.660 notch 18.203 - 21.349 width 3.146 central n= 600 median 21.535 notch 20.795 - 22.330 width 1.535 do the notches of north and central intersect: False ``` The notch is the margin of variation of the median and is produced with the very same bootstrap as in the previous lesson. The medians of north and the twin dataset are the same number and so are their observation counts, yet their notches are 0.920 and 3.146 m³ wide. The three and a half times in between come not from the number of observations but from the density around the median: a median standing in a trough goes much farther from draw to draw. This is the numeric expression of what the box summary lacks. While two boxes that put the five numbers side by side cannot be told apart, the margin of the same median differs by more than three times. A report that writes the notch next to the median reports, even without drawing anything, that the two distributions are not the same. The notch here is the sixth and seventh number added on top of the five, and its cost is a single run of 1,000 draws. The last line is a real comparison: the notches of the north and central regions do not intersect. Whether this means "the medians of the two regions are different" is not decided here. The interpretation built on intersection and non-intersection is the first question of the next lesson. #### Summary - The box summary is a record, not a drawing: five numbers, the interquartile range and two fences; the 921 readings of north come down to the values 0.00 | 15.74 | 19.66 | 24.63 | 45.47 and an interquartile range of 8.89 m³. - A second dataset that gives the same five numbers can be built; the two box summaries are exactly the same and the difference shows only in the number of observations outside the fences, 50 against 35. - Bin counts separate the two datasets: in north a single peak at 16–20 m³, in the twin two peaks in the 12–16 and 24–28 m³ bins, and the median is inside the trough with counts of 114 and 110. - The notch of the median is produced with bootstrap and comes out 0.920 m³ and 3.146 m³ wide for the same median; the width depends on the number of observations and on the density around the median. - The notches of north and central do not intersect; this is a display cue and by itself not a test result. #### Next Step In this lesson intervals were produced twice and widths were written twice, but what an interval says was never defined. The confidence interval was used as a tool in M15/K02; the next lesson takes the interval itself as its subject. How the width depends on the number of observations is measured — what happens to the width when the number of observations is quadrupled is counted — and what the interval does not say is written out plainly. In the end the intersection of two intervals is dealt with: intersection does not mean no difference, and the proof of that is computing the difference's own interval. ### Interpreting Confidence Intervals Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/interpreting-confidence-intervals The east region's mean of 24.674 m³ settles into a 23.894–25.431 m³ interval with bootstrap and the width comes out as 1.537 m³, but only 37 of the 557 readings are inside this interval: the interval belongs to the summary, not to the observations. When the observation count rises from 35 to 139 and from 139 to 557, the width falls from 6.601 to 3.009 and then to 1.459 m³. In 113 of 120 samples of forty observations the interval contains the population mean. The two periods of north give two intervals that overlap by 0.148 m³, while the interval of the difference itself, 0.281–2.612 m³, leaves zero outside. In the previous two lessons intervals were produced twice and widths were written twice, but what an interval says was not defined. The confidence interval is not new in this curriculum: it was used as a **tool** in the Measurement topic of the User Experience and Behavior Design course, and the sample size and the multiple comparison correction were computed there. Those computations are not repeated here. In this lesson the interval itself is the subject. Here the interval is produced not with a formula but with **bootstrap**. This choice is not a convenience: the formula route asks for an assumption about the distribution, and that assumption is disputable for reading data, while the draw route asks only for the dataset itself. The measure is gathered under two headings: what the width depends on and what the interval does not say. The most expensive item of the second heading is left for the end — the overlap of two intervals does not mean that two means are not different, and the proof of this is given with a number. - **AS12. A confidence interval is not a measurement but the output of a procedure.** The dataset is drawn from 1,000 times at the same size and with replacement, the summary is recomputed at every draw, and the middle 95% of the resulting values is written as the interval. - **AS13. The interval belongs to the summary, not to the observations.** It says nothing about where the individual readings stand, and the vast majority of the readings are outside the interval. - **AS14. The width depends on the observation count**, and when the observation count is quadrupled the width falls to about half. The price of halving the margin of error is four times the observations. - **AS15. The 95% belongs not to a single interval but to the procedure.** About 95% of the intervals built with the same procedure contain the population value; this is what can be said for the single interval at hand. - **AS16.** The interval assumes that the dataset represents the population. It neither measures nor corrects selection bias; in a biased dataset a narrow interval is a narrow interval around the wrong number. - **AS17.** The overlap of two intervals does not say "no difference". The difference is a separate number. - **AS18.** The place where the decision is read is the interval of the difference itself: it is checked whether zero is inside. #### How the Interval Is Produced ```python # dataset.py — MODEL: the measurement network of K03-K04 is regenerated with the same seed. # A long-format reading table of 3,199 rows, 1,329 subscribers, three periods, five regions. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF PERIODS = ["2026-01", "2026-02", "2026-03"] REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t): x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value(): s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u, w): t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 _kept, _dropped = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]) pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) (_dropped if r() < 0.046 else _kept).append( {"subscriber_id": 10001 + i, "region": b[0], "base": b[2]}) VALID = {a["subscriber_id"] for a in _kept} READINGS = [] for a in _kept + _dropped: r = rng(SEED + 7000 + a["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor( a["base"] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if a["subscriber_id"] not in VALID: r() continue READINGS.append({"subscriber_id": a["subscriber_id"], "region": a["region"], "period": PERIODS[d], "consumption_m3": t}) r() def mean(v): return sum(v) / len(v) DRAWS = 1000 def percentile(s, p): i = p * (len(s) - 1) a = int(i) z = min(a + 1, len(s) - 1) return s[a] + (s[z] - s[a]) * (i - a) def interval(v, seed, draws=DRAWS): rand = rng(seed) n = len(v) d = sorted(mean([v[int(rand() * n)] for _ in range(n)]) for _ in range(draws)) return percentile(d, 0.025), percentile(d, 0.975) east = [o["consumption_m3"] for o in READINGS if o["region"] == "east"] a, z = interval(east, SEED + 301) print(f"east mean {mean(east):.3f} m3, n={len(east)}, {DRAWS} draws") print(f"95% confidence interval {a:.3f} - {z:.3f} m3, width {z - a:.3f} m3") inside = sum(1 for x in east if a <= x <= z) print(f"{inside} of the 557 readings ({100 * inside / len(east):.2f} percent) are inside this interval") ``` ``` east mean 24.674 m3, n=557, 1000 draws 95% confidence interval 23.894 - 25.431 m3, width 1.537 m3 37 of the 557 readings (6.64 percent) are inside this interval ``` The procedure is a single sentence: since there is no population at hand, the dataset is put in the place of the population, 1,000 draws of the same size are made from it with replacement, the mean is recomputed at every draw, and the middle 95% of the resulting 1,000 means is written as the interval. For the east region the result is 23.894–25.431 m³, and the margin of error, that is, the half-width of the interval, is 0.769 m³. The last line turns the most commonly confused side of the interval into a number. The interval contains only 37 of the 557 readings, that is, 6.64%. The interval reports not where the readings are but where the **mean** may be. The sentence "95% of the subscribers consume between 23.9 and 25.4 m³" cannot be derived from this output and is false in the dataset. In report language the distinction is kept like this: the interval is always written under the name of a **summary** — the interval of the mean, the notch of the median, the interval of the ratio. When the name is not written, the reader taking it for an interval of the observations is an expected misunderstanding. The only assumption of the procedure also stands in the open: the dataset is accepted as representing the population. Bootstrap puts the dataset at hand in the place of the population that is not at hand, and the price of this substitution grows in small datasets. A value that is not in the dataset cannot appear in any draw, so the interval sees only as much as the dataset can see. #### What the Width Depends On ```python # width.py — the readings of the same region are first shuffled, then three nested # subsets are taken: 35, 139, 557. At every step the observation count is roughly quadrupled. rand = rng(SEED + 304) shuffled = east[:] for i in range(len(shuffled) - 1, 0, -1): j = int(rand() * (i + 1)) shuffled[i], shuffled[j] = shuffled[j], shuffled[i] previous = None for n in (35, 139, 557): sub = shuffled[:n] a2, z2 = interval(sub, SEED + 305) g = z2 - a2 extra = "" if previous is None else f", {g / previous:.3f} times the previous width" print(f"n={n:3d} mean {mean(sub):6.3f} interval {a2:6.3f} - {z2:6.3f} " f"width {g:.3f}{extra}") previous = g ``` ``` n= 35 mean 24.242 interval 20.825 - 27.426 width 6.601 n=139 mean 24.115 interval 22.507 - 25.515 width 3.009, 0.456 times the previous width n=557 mean 24.674 interval 23.985 - 25.445 width 1.459, 0.485 times the previous width ``` The three subsets are nested: the 35 readings are inside the 139, and the 139 inside the 557. The observation count is roughly quadrupled at every step, and the width falls to roughly half at every step — from 6.601 to 3.009 and then to 1.459 m³. The ratios come out as 0.456 and 0.485, that is, around halving. The counterpart of this in planning is direct: a measurement that wants to halve the margin of error must collect four times the observations. In a breakdown of thirty-five readings the mean has an interval 6.6 m³ wide, and in that breakdown nothing can be said about the 25 m³ tariff step; the observation count needed to be able to answer the same question is sixteen times as large. The price of cutting a breakdown finer is also visible here: every split lowers the observation count, and the falling observation count grows the width in inverse proportion to the square root. It should be noticed that the means move too: 24.242, 24.115 and 24.674. All three numbers are correct and all three come from the same region. The difference among them is exactly what the interval measures. The second determinant of the width is the dataset's own **dispersion**. At the same observation count, a breakdown whose values stand close to each other gives a narrow interval, and a widely dispersed breakdown gives a wide one. The analyst has only one of these two levers: the observation count is the subject of the measurement plan, while dispersion is a property of the thing measured. This is why a narrow interval is not an achievement of analysis but mostly the result of a collection decision. #### The 95% Is the 95% of What ```python # coverage.py — the 557 readings are accepted as the population, samples of 40 observations # are drawn from it and a short interval is built for every sample. What is counted is # how many of the intervals contain the population mean. POPULATION = mean(east) rand = rng(SEED + 308) TRIALS, SAMPLE, SHORT = 120, 40, 200 covering = 0 for _ in range(TRIALS): sample = [east[int(rand() * len(east))] for _ in range(SAMPLE)] d = sorted(mean([sample[int(rand() * SAMPLE)] for _ in range(SAMPLE)]) for _ in range(SHORT)) if percentile(d, 0.025) <= POPULATION <= percentile(d, 0.975): covering += 1 print(f"population mean {POPULATION:.3f} m3 accepted; {TRIALS} samplings of " f"{SAMPLE} observations, interval built with {SHORT} draws") print(f"{covering} of the intervals ({100 * covering / TRIALS:.1f} percent) contain the " f"population mean") ``` ``` population mean 24.674 m3 accepted; 120 samplings of 40 observations, interval built with 200 draws 113 of the intervals (94.2 percent) contain the population mean ``` The number 95 does not belong to the single interval at hand; it belongs to the **procedure**. To see this, the 557 readings of the east region were taken as the population, 120 separate samples of 40 observations were drawn from it, and an interval was built for each. Of the 120 intervals built, 113 contain the population mean, that is, 94.2%. This is the correct form of the sentence: about 95% of the intervals built with this procedure contain the population value. For the single interval at hand it cannot be said that "the population mean is inside with 95% probability", because the population mean is a number and it is either inside or not; what moves is the interval, not the population value. That the resulting number is 94.2 instead of 95 is also a piece of information and not a deviation to be corrected. In samples of forty observations the bootstrap builds the interval a little narrow; this ratio, measured with 120 trials, has a margin of its own of about two points. If the run count were enlarged the ratio would approach 95, but the speed of approach also obeys the very rule this lesson measures. What the interval does not say also stays outside this run. If the dataset does not represent the population — if the network has only the meters that could be read — the interval neither sees nor corrects this deficiency; a narrow interval is only narrow around the wrong number. Selection bias is the job of the last lesson of this topic. #### The Overlap of Two Intervals ```python # overlap.py — the two periods of north: first two separate intervals, then the interval # of the difference itself. The difference is recomputed at every draw by drawing # separately from the two groups. north2 = [o["consumption_m3"] for o in READINGS if o["region"] == "north" and o["period"] == "2026-02"] north3 = [o["consumption_m3"] for o in READINGS if o["region"] == "north" and o["period"] == "2026-03"] a2, z2 = interval(north2, SEED + 306) a3, z3 = interval(north3, SEED + 306) print(f"north 2026-02: n={len(north2)} mean {mean(north2):.3f} interval {a2:.3f} - {z2:.3f}") print(f"north 2026-03: n={len(north3)} mean {mean(north3):.3f} interval {a3:.3f} - {z3:.3f}") print(f"the two intervals overlap for {min(z2, z3) - max(a2, a3):.3f} m3") rand = rng(SEED + 307) n2, n3 = len(north2), len(north3) f = sorted(mean([north2[int(rand() * n2)] for _ in range(n2)]) - mean([north3[int(rand() * n3)] for _ in range(n3)]) for _ in range(DRAWS)) fa, fz = percentile(f, 0.025), percentile(f, 0.975) print(f"interval of the difference itself: {mean(north2) - mean(north3):.3f} m3 " f"({fa:.3f} - {fz:.3f}), is zero inside: {fa <= 0 <= fz}") ``` ``` north 2026-02: n=319 mean 20.324 interval 19.491 - 21.193 north 2026-03: n=239 mean 18.879 interval 18.068 - 19.639 the two intervals overlap for 0.148 m3 interval of the difference itself: 1.446 m3 (0.281 - 2.612), is zero inside: False ``` The two periods of the north region were put side by side. The two intervals overlap by 0.148 m³, and a reader who looks at the overlap draws the conclusion "there is no difference between the two periods". This conclusion is wrong, and its wrongness is not an interpretive preference but something that can be shown by computation. When the interval of the difference itself is computed, the observed difference of 1.446 m³ settles into the 0.281–2.612 m³ interval, and zero stays outside this interval. That is, the same data says "uncertain" when looking at the two separate intervals and "away from zero" when looking at the interval of the difference. Of the two, the correct one is the second, because the question asked is what the difference is, and the difference is a separate number. The reason in short is this: two separate intervals put two separate uncertainties side by side **at their full length**, whereas the uncertainty of the difference is not as large as the sum of the two. Non-overlap always shows that the difference is away from zero, but overlap shows nothing. If a decision is to be made, the number to look at is single: the interval of the difference and the place of zero in that interval. The rule in practice follows from this. If the intervals of two groups are printed side by side in a report, that is not the place where the comparison is made; the interval of the difference is computed separately and written separately. A two-interval display is useful for introducing each group on its own, not for comparison. When this distinction is not written, the same output produces two different decisions, and which one is given depends on which line the reader looks at. #### Summary - The confidence interval was produced not with a formula but with bootstrap: the east region's mean of 24.674 m³ settled into the 23.894–25.431 m³ interval in 1,000 draws, that is, into a margin of error of 0.769 m³. - The interval belongs to the summary, not to the observations; only 37 of the 557 readings, 6.64%, are inside this interval. - When the observation count rose from 35 to 139 and from 139 to 557, the width fell from 6.601 to 3.009 and then to 1.459 m³; the price of halving the margin of error is four times the observations. - The 95% belongs not to a single interval but to the procedure: in 113 of 120 samples of forty observations, that is, in 94.2%, the interval contained the population mean. - The two periods of north gave two intervals that overlap by 0.148 m³, but the interval of the difference itself came out as 0.281–2.612 m³ and left zero outside; overlap does not mean no difference. #### Next Step The interval of the difference left zero outside, and this is the output that stands closest to a decision. Still, it is not a test by itself: the question not asked is "how often a difference of this size would appear when there is no difference at all". The next lesson answers this question with a **permutation test**, defines the null hypothesis, the p-value and the significance level through permutation, and tests the same two groups by counting them both as paired and as independent. The p-value is printed twice, the effect size stands next to it, and what the p-value is not is shown with numbers. ### Applying Hypothesis Testing Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/applying-hypothesis-testing Among the 1,101 subscribers read in both periods the observed difference is 0.957 m³, and this single difference gives two separate p-values: in the independent assumption 24 of 2,000 permutations exceed the difference and p comes out as 0.01249, in the paired assumption 5 exceed it and p becomes 0.00300. The assumption check shows with a correlation of 0.213 that the pairing is real; when the pairing is broken p returns to 0.01099. The effect sizes, 0.112 and 0.089, stand next to p. In a quarter of the same dataset, 275 pairs, the difference is larger at 1.132 m³, yet p comes out thirty-five times larger at 0.10695. The previous lesson computed the interval of the difference itself and saw that zero stayed outside the interval. That stood close to a decision but it was not a test. The question a test asks is different: when there is no difference at all, **in how many runs out of how many** does a difference of the observed size appear. In this course the answer to that question is not read from a table; it is counted by running. The first lesson of the exploratory analysis topic randomly redistributed the region labels and counted in how many runs the difference between east and north appeared; here the same method is named a test, and the main measure is this: **the choice of test is an assumption, and when the assumption changes the p-value moves.** The same two groups, the same observed difference, two p-values. - **AS19. The null hypothesis is not a sentence but a procedure.** The assumption "there is no difference between the two groups" means that the labels are interchangeable among the observations; the **permutation test** does exactly this and recomputes the difference at every shuffle. - **AS20. The p-value** is how many of the shuffled datasets produced a difference as large as the observed one or larger. The observed run also takes part in the count, so the numerator and the denominator are each increased by one. - **AS21. The significance level** is a threshold written **before** the test. A threshold chosen after the result is seen is not a threshold but a result. - **AS22.** The choice of test is an assumption, and the assumption can be checked. Whether the pairing is real is measured with the correlation of the paired readings. - **AS23.** The same observed difference gives two p-values under two assumptions. The difference does not change, the p changes. - **AS24. The p-value is not the size of the effect.** The effect size, the difference divided by the standard deviation, is written next to p as a separate number. - **AS25. The p-value is not the probability of the null hypothesis either.** Moreover its resolution depends on the run count: the smallest p-value that can be written with 2,000 permutations is 0.00050. #### The Null Hypothesis Is a Procedure ```python # dataset.py — MODEL: the fictional measurement network of K03-K04 is regenerated with the same seed. # A long-format reading table of 3,199 rows, 1,329 subscribers, three periods, five regions. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF PERIODS = ["2026-01", "2026-02", "2026-03"] REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t): x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value(): s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u, w): t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 _kept, _dropped = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]) pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) (_dropped if r() < 0.046 else _kept).append( {"subscriber_id": 10001 + i, "region": b[0], "base": b[2]}) VALID = {a["subscriber_id"] for a in _kept} READINGS = [] for a in _kept + _dropped: r = rng(SEED + 7000 + a["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor( a["base"] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if a["subscriber_id"] not in VALID: r() continue READINGS.append({"subscriber_id": a["subscriber_id"], "region": a["region"], "period": PERIODS[d], "consumption_m3": t}) r() def mean(v): return sum(v) / len(v) PERMUTATIONS = 2000 by_subscriber = {} for o in READINGS: by_subscriber.setdefault(o["subscriber_id"], {})[o["period"]] = o["consumption_m3"] PAIRS = [(v["2026-01"], v["2026-02"]) for v in by_subscriber.values() if "2026-01" in v and "2026-02" in v] first = [x for x, _ in PAIRS] second = [y for _, y in PAIRS] diff = [x - y for x, y in PAIRS] OBSERVED = mean(diff) print(f"subscribers read in both periods: {len(PAIRS)}") print(f"2026-01 mean {mean(first):.3f} m3, 2026-02 mean {mean(second):.3f} m3") print(f"observed difference {OBSERVED:.3f} m3") def shuffle(v, draw): h = v[:] for i in range(len(h) - 1, 0, -1): j = int(draw() * (i + 1)) h[i], h[j] = h[j], h[i] return h draw = rng(SEED + 401) everything, n1 = first + second, len(first) exceeded = 0 for _ in range(PERMUTATIONS): h = shuffle(everything, draw) if abs(mean(h[:n1]) - mean(h[n1:])) >= abs(OBSERVED): exceeded += 1 p_independent = (exceeded + 1) / (PERMUTATIONS + 1) print(f"independent assumption: {exceeded} of {PERMUTATIONS} permutations exceeded the observed difference, " f"p = {p_independent:.5f}") ``` ``` subscribers read in both periods: 1101 2026-01 mean 20.638 m3, 2026-02 mean 19.681 m3 observed difference 0.957 m3 independent assumption: 24 of 2000 permutations exceeded the observed difference, p = 0.01249 ``` The question is a report question: did consumption fall from the first period to the second. There are 1,101 subscribers who could be read in both periods, the means are 20.638 and 19.681 m³, and the observed difference is 0.957 m³. The dataset is fictional and regenerated with the same seed from the measurement network of K03–K04; that only 1,101 of the 1,329 subscribers could be read in both periods is the dataset's own deficiency, not a selection. The null hypothesis is set up here not as a sentence but as a procedure. Saying "the period label has no effect" means accepting that which period each of the 2,202 readings belongs to is random. The code does exactly this: it puts all the readings in a single bag, shuffles it, counts the first 1,101 as the first period, and recomputes the difference. This is repeated 2,000 times. The p-value is a **count** from this run: 24 of the permutations produced a difference as large as the observed 0.957 m³ or larger. One is added to the numerator and the denominator, because the observed arrangement is itself one of the possible arrangements; p comes out as 0.01249. If a significance level of 5% was written before the test, this result is below the threshold. That the significance level is written **before** is the load-bearing part of this sentence. When the threshold is chosen after the result, it stops being a decision rule and turns into a label that repeats the result itself; in that case the test has eliminated nothing. The number of the threshold is also a convention, not a measurement: in this dataset the choice between 5% and 1% does more work than which of the two assumptions below is looked at. #### The Pairing Assumption ```python # paired.py — the same dataset, another assumption: every row is the two readings of one subscriber. # This time the null hypothesis means "within every pair the sign of the difference is random". draw = rng(SEED + 403) exceeded2 = 0 for _ in range(PERMUTATIONS): t = sum(f if draw() < 0.5 else -f for f in diff) / len(diff) if abs(t) >= abs(OBSERVED): exceeded2 += 1 p_paired = (exceeded2 + 1) / (PERMUTATIONS + 1) def std_dev(v): m = mean(v) return math.sqrt(sum((x - m) ** 2 for x in v) / (len(v) - 1)) pooled = math.sqrt((std_dev(first) ** 2 + std_dev(second) ** 2) / 2) ma, mb = mean(first), mean(second) cov = sum((x - ma) * (y - mb) for x, y in PAIRS) / (len(PAIRS) - 1) print(f"paired assumption: {exceeded2} permutations exceeded the observed difference, p = {p_paired:.5f}") print(f"assumption check: correlation of the paired readings {cov / (std_dev(first) * std_dev(second)):.3f}") print(f"the difference is {OBSERVED:.3f} m3 in both; effect size " f"{OBSERVED / pooled:.3f} (independent) and {OBSERVED / std_dev(diff):.3f} (paired)") print(f"p ratio {p_independent / p_paired:.2f} times") ``` ``` paired assumption: 5 permutations exceeded the observed difference, p = 0.00300 assumption check: correlation of the paired readings 0.213 the difference is 0.957 m3 in both; effect size 0.112 (independent) and 0.089 (paired) p ratio 4.17 times ``` The same 2,202 readings were counted with a second assumption. This time every row is the two readings of one subscriber, and the null hypothesis takes the form "within every pair the sign of the difference is random". The shuffle changes accordingly: the labels swap places not between the groups but within every pair. The result is p 0.00300. The observed difference is 0.957 m³ under both assumptions and has not moved by a single cubic meter. The only number that moved is the p-value, and it shrank 4.17 times. The whole difference between them comes from the assumption taken while choosing the test. Both assumptions stay below the 5% threshold, but at the 1% threshold the decisions part ways. The assumption check is in the second line. For the pairing to bring a gain, the two readings of the same subscriber must be related to each other; the measured correlation is 0.213. This means that the subscribers' consumption level is partly preserved across periods, and it explains why the paired assumption gave a smaller p: the dispersion between subscribers is subtracted from the computation. The effect size stands next to p on the same line: 0.112 and 0.089. Both are the difference divided by one standard deviation, and both are small. What is being talked about is a decline of 4.6% of mean consumption; the smallness of the p-value does not enlarge this decline. #### Pairing Does Not Come from the Procedure ```python # fake.py — pairing is not a procedure but a fact: if the second-period readings # are matched with other subscribers the difference stays the same, the gain disappears. draw = rng(SEED + 404) fake_second = shuffle(second, draw) fake_diff = [x - y for x, y in zip(first, fake_second)] draw = rng(SEED + 405) exceeded3 = 0 for _ in range(PERMUTATIONS): t = sum(f if draw() < 0.5 else -f for f in fake_diff) / len(fake_diff) if abs(t) >= abs(mean(fake_diff)): exceeded3 += 1 print(f"fake pairing: same difference {mean(fake_diff):.3f} m3, " f"p = {(exceeded3 + 1) / (PERMUTATIONS + 1):.5f}") print(f"standard deviation of the differences: with the real pairing {std_dev(diff):.3f}, " f"with the fake {std_dev(fake_diff):.3f} m3") ``` ``` fake pairing: same difference 0.957 m3, p = 0.01099 standard deviation of the differences: with the real pairing 10.753, with the fake 12.275 m3 ``` It can be thought that the small p-value of the paired test was produced by the test itself. This run tests that: the second-period readings are shuffled and matched with other subscribers, and then exactly the same paired test is applied. The observed difference does not change, because the mean of the two groups has stayed the same: again 0.957 m³. But p rises to 0.01099 and returns to the same region as the 0.01249 that the independent count gave. The reason is in the last line: the standard deviation of the differences has risen from 10.753 to 12.275 m³, that is, the dispersion the real pairing had lowered has come back. The result gives the rule about the choice of test. The paired test is not a more powerful method; **if the pairing is real** it is more powerful. Inventing a pairing that does not exist does not shrink the p-value, it only moves the computation to the wrong place. This is why the assumption check is written not as an appendix to the test but as its precondition. The check itself is not expensive either: here a single correlation number and a single standard deviation comparison were enough, and both are too small to be counted next to the cost of the test. #### What the p-Value Is Not ```python # scale.py — the same effect, a different observation count. The subsets are chosen # nested and the same paired test is applied to each. draw = rng(SEED + 406) order = shuffle(list(range(len(PAIRS))), draw) for divisor in (4, 2, 1): sub = [diff[i] for i in sorted(order[:len(PAIRS) // divisor])] r2 = rng(SEED + 403) over4 = 0 for _ in range(PERMUTATIONS): t = sum(f if r2() < 0.5 else -f for f in sub) / len(sub) if abs(t) >= abs(mean(sub)): over4 += 1 print(f"n={len(sub):4d} difference {mean(sub):6.3f} m3 effect {mean(sub) / std_dev(sub):.3f} " f"p = {(over4 + 1) / (PERMUTATIONS + 1):.5f}") print(f"smallest p that can be written with {PERMUTATIONS} permutations = {1 / (PERMUTATIONS + 1):.5f}") ``` ``` n= 275 difference 1.132 m3 effect 0.099 p = 0.10695 n= 550 difference 1.192 m3 effect 0.110 p = 0.01149 n=1101 difference 0.957 m3 effect 0.089 p = 0.00300 smallest p that can be written with 2000 permutations = 0.00050 ``` The three lines show directly what the p-value is not. In the subset of a quarter the observed difference is 1.132 m³ — **larger** than the 0.957 m³ of the full dataset — and the effect size, 0.099, is also above 0.089. Even so the p-value is 0.10695, that is, thirty-five times the full dataset's 0.00300 and above the 5% threshold. The same three lines are also read the other way: as p shrank, the effect did not grow. The line that gives the smallest p is the line that gives the smallest effect. The p-value measures not the size of the difference but how far the difference, together with the observation count, is from chance. This is why p is not reported alone; the observation count and the effect size are written next to it. The second wrong reading closes here too: p is not the probability that the null hypothesis is true. What it measures is how often a difference this large is produced when the null hypothesis is true — the number 0.00300 does not mean "the chance of no difference is three in a thousand", it means "when there is no difference, a result like this came up in 5 of 2,000 permutations". The last line writes the limit of the number itself: the smallest p that can be produced with 2,000 permutations is 0.00050. In a test where no shuffle could exceed the observed difference, what can be written is not "p is zero" but "p is smaller than 0.00050". The run count sets the resolution of the margin. #### Summary - The null hypothesis in the permutation test is a procedure: the period labels are randomly redistributed 2,000 times, the difference is recomputed each time, and the runs that exceed the observed difference are counted. - In the same 1,101 pairs the observed difference is 0.957 m³ and does not change; the p-value comes out as 0.01249 in the independent assumption and 0.00300 in the paired assumption, that is, the assumption moves p by 4.17 times. - The pairing assumption is checked: the correlation of the paired readings is 0.213. When the pairing is broken the standard deviation of the differences rises from 10.753 to 12.275 m³ and p returns to the independent assumption with 0.01099. - The p-value is not the effect size: in the subset of 275 pairs the difference is 1.132 m³ and the effect is larger at 0.099, yet p is 0.10695, thirty-five times that of the full dataset. - The p-value is not the probability of the null hypothesis either, and its resolution depends on the run count; the smallest p that can be written with 2,000 permutations is 0.00050. #### Next Step In this lesson a single question was asked and a single p-value was written: is there a difference between the two periods. Once the question is chosen, the test itself is orderly work and, together with the assumption check, gives a defensible result. In a real analysis, however, the question is not single. Which metric will be counted from the same experiment dataset, which breakdown will be opened, which observation window will be taken and what will be done about an unbalanced assignment are chosen one by one; every choice means a new test, and the number of questions enlarges the margin. The next lesson counts how many different decisions come out of the same experiment dataset and measures how much of this margin pre-registration closes. ### A/B Test Analysis Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/ab-test-analysis Counting the four analytic degrees of freedom in front of the analyst once the dataset of a finished-design experiment is in hand: the choice of metric, the breakdown, the observation window and the handling of the unbalanced assignment produce 108 combinations. The single pre-registered decision gives a 0.123 m³ difference, a margin of variation of [-0.698, 0.952] and a p-value of 0.786; of the 108 combinations walked on the same dataset, 4 fall below the 0.05 threshold, 9 below 0.10, 45 point toward an increase and 63 toward a decrease, and the difference in the mean metric wanders between -2.113 and 1.119 m³. The true difference between the groups is exactly zero. Pre-registration brings the number of combinations from 108 down to 1. The previous lesson ran a single test: two groups, one p-value, the effect size next to it. The answer given was correct but it belonged to a single question. Once the dataset of a finished-design experiment is in hand, the question is not single. Which metric will be reported, will the breakdown be looked at, where will the observation window close, what will be done if the assignment stayed unbalanced. Every choice is defensible separately and every choice produces a different number. The number of this lesson is **the number of decisions that come out of the same experiment dataset**. The design of the experiment is not set up here: how many subscribers were needed, why the experiment was not stopped early and how many comparisons tightened the limit and how were computed in the A/B and Multivariate Testing lesson of the User Experience and Behavior Design course. This lesson starts after the design is finished: the assignment is made, the readings are collected, the dataset is in hand. - **AS26.** The dataset is inherited from K03 and is **fictional**: a reading table of 3,199 rows, 1,329 subscribers, three periods, five regions. The seed is written in the code. - **AS27.** The experiment is **fictional** too: period 1 is the baseline reading, the reminder notification was sent after period 1, and the assignment is unbalanced 70/30. The assignment is drawn independently of the readings and changes no reading; **the true difference between the groups is exactly zero.** - **AS28.** The test is produced by **permutation** (500 permutations), the interval by **bootstrap** (1,000 draws). No distribution table, critical value or ready-made test function is used. - **AS29.** The "significant" threshold is 0.05. This is not a measurement but a decision, and its share is written at the end of the lesson. - **AS30.** The four analytic degrees of freedom and their option counts: metric 2, breakdown 9, observation window 3, handling of the unbalanced assignment 2. #### The Single Decision Written Before the Experiment A subscriber who enters the experiment dataset is one who has a period 1 reading. The primary decision is written in a single line: metric mean consumption, no breakdown, window periods 2 and 3, assignment left raw. ```python # ab_test.py — MODEL. The fictional dataset of K03 is produced with the same seed, and a fictional # experiment assignment is added on top of it. The assignment is drawn independently of the readings # and changes no reading: the true difference between the groups is exactly zero. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t): x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value(): s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u, w): t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 SUBSCRIBERS, DROPPED = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]) pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) (DROPPED if r() < 0.046 else SUBSCRIBERS).append( {"subscriber_id": 10001 + i, "region": b[0], "base": b[2]}) VALID = {k["subscriber_id"] for k in SUBSCRIBERS} READINGS = {} for k in SUBSCRIBERS + DROPPED: r = rng(SEED + 7000 + k["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor( k["base"] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if k["subscriber_id"] in VALID: READINGS.setdefault(k["subscriber_id"], {})[d + 1] = t r() # The experiment is FICTIONAL: period 1 is the baseline reading, the reminder notification was # sent after period 1 and the assignment is unbalanced 70/30 (two of every three subscribers are B). EXPERIMENT = [] for k in SUBSCRIBERS: d = READINGS.get(k["subscriber_id"], {}) if 1 not in d: continue r = rng(SEED + 31000 + k["subscriber_id"]) EXPERIMENT.append({"region": k["region"], "group": "B" if r() < 0.70 else "A", "baseline_reading": d[1], "d2": d.get(2), "d3": d.get(3)}) PERMUTATIONS, RESAMPLES = 500, 1000 rs = rng(SEED + 55000) def shuffle(items): h = items[:] for i in range(len(h) - 1, 0, -1): j = int(rs() * (i + 1)) h[i], h[j] = h[j], h[i] return h def permutation_test(a, b): """Permutation test: the labels are shuffled, the difference is recomputed.""" na, nb = len(a), len(b) observed = sum(b) / nb - sum(a) / na exceeded = 0 for _ in range(PERMUTATIONS): h = shuffle(a + b) if abs(sum(h[na:]) / nb - sum(h[:na]) / na) >= abs(observed) - 1e-12: exceeded += 1 return observed, exceeded / PERMUTATIONS, na, nb def interval(a, b): """Bootstrap: the difference is computed on a dataset redrawn 1,000 times.""" na, nb, diffs = len(a), len(b), [] for _ in range(RESAMPLES): ya = sum(a[int(rs() * na)] for _ in range(na)) / na yb = sum(b[int(rs() * nb)] for _ in range(nb)) / nb diffs.append(yb - ya) diffs.sort() return diffs[24], diffs[974] def value(x, fields, metric): d = [x[a] for a in fields if x[a] is not None] if not d: return None return sum(d) / len(d) if metric == "mean" else (1.0 if max(d) >= 25 else 0.0) def split_groups(sub, fields, metric): va, vb = [], [] for x in sub: v = value(x, fields, metric) if v is not None: (va if x["group"] == "A" else vb).append(v) return va, vb print(f"fictional dataset (seed {SEED}): {len(SUBSCRIBERS)} subscribers, " f"{sum(len(v) for v in READINGS.values())} readings") print(f"experiment dataset: {len(EXPERIMENT)} subscribers, " f"A {sum(1 for x in EXPERIMENT if x['group'] == 'A')} / " f"B {sum(1 for x in EXPERIMENT if x['group'] == 'B')}") # The single decision written before the experiment: metric mean, no breakdown, window # periods 2-3, assignment left raw. va, vb = split_groups(EXPERIMENT, ("d2", "d3"), "mean") observed, p, na, nb = permutation_test(va, vb) low, high = interval(va, vb) print(f"\nprimary decision mean consumption difference (B - A)") print(f" observed difference : {observed:.3f} m3 ({na} + {nb} subscribers)") print(f" margin of variation : [{low:.3f}, {high:.3f}] m3 (1,000 bootstrap draws)") print(f" p-value : {p:.3f} ({PERMUTATIONS} permutations)") ``` ``` fictional dataset (seed 20260218): 1329 subscribers, 3199 readings experiment dataset: 1260 subscribers, A 402 / B 858 primary decision mean consumption difference (B - A) observed difference : 0.123 m3 (344 + 757 subscribers) margin of variation : [-0.698, 0.952] m3 (1,000 bootstrap draws) p-value : 0.786 (500 permutations) ``` This is the line to be reported: a difference of 0.123 m³, next to it a margin of variation 1.650 m³ wide, and a p-value of 0.786. The interval contains zero; the permutation test produces a difference of the same size in four of every five runs. The number comes from the mean of 344 and 757 subscribers, and the group sizes are not equal. The only information this line conveys is that the difference stays inside the margin of variation. This is where the experiment ends. What comes after is analysis. #### Four Analytic Degrees of Freedom On the same dataset there are four choices, and each of the four has more than one defensible value. **Metric:** mean consumption, or the proportion of subscribers entering the upper tariff tier. **Breakdown:** the whole dataset, one of the five regions, or one of three tariff tiers by baseline reading. **Observation window:** both periods together, period 2 only, period 3 only. **Balance:** the unbalanced assignment is left raw or the B group is trimmed to the size of A. ```python # Four degrees of freedom: the metric, the breakdown, the observation window, how the # unbalanced assignment is handled. All the combinations are walked and it is counted how many # come out "significant". METRICS = ("mean", "upper tier rate") BREAKDOWN = ("whole dataset", "north", "south", "east", "west", "central", "under 10", "10-25", "over 25") WINDOWS = {"periods 2-3": ("d2", "d3"), "period 2 only": ("d2",), "period 3 only": ("d3",)} BALANCE = ("raw", "trimmed") def tier(v): return "under 10" if v < 10 else ("10-25" if v < 25 else "over 25") def selected(brk): if brk == "whole dataset": return EXPERIMENT return [x for x in EXPERIMENT if x["region"] == brk or tier(x["baseline_reading"]) == brk] SWEEP = [] for metric in METRICS: for brk in BREAKDOWN: for name, fields in WINDOWS.items(): for balance in BALANCE: sub = [x for x in selected(brk) if value(x, fields, metric) is not None] if balance == "trimmed": a = [x for x in sub if x["group"] == "A"] sub = a + [x for x in sub if x["group"] == "B"][:len(a)] va, vb = split_groups(sub, fields, metric) if len(va) < 2 or len(vb) < 2: continue SWEEP.append((metric, brk, name, balance) + permutation_test(va, vb)) significant = [s for s in SWEEP if s[5] <= 0.05] positive = sum(1 for s in SWEEP if s[4] > 0) means = [s[4] for s in SWEEP if s[0] == "mean"] print(f"combinations run : {len(SWEEP)} " f"({len(METRICS)} metric x {len(BREAKDOWN)} breakdown x {len(WINDOWS)} window " f"x {len(BALANCE)} balance)") print(f"p <= 0.05 : {len(significant)}") print(f"p <= 0.10 : {sum(1 for s in SWEEP if s[5] <= 0.10)}") print(f"direction : {positive} combinations up, {len(SWEEP) - positive} combinations down") print(f"difference range for the mean metric: {min(means):.3f} to {max(means):.3f} m3") print("\nmetric breakdown window balance difference p n(A) n(B)") for s in sorted(significant, key=lambda s: s[5]): print(f"{s[0]:<18} {s[1]:>8} {s[2]:<15} {s[3]:<10} {s[4]:>7.3f} {s[5]:>6.3f} " f"{s[6]:>5} {s[7]:>5}") ``` ``` combinations run : 108 (2 metric x 9 breakdown x 3 window x 2 balance) p <= 0.05 : 4 p <= 0.10 : 9 direction : 45 combinations up, 63 combinations down difference range for the mean metric: -2.113 to 1.119 m3 metric breakdown window balance difference p n(A) n(B) upper tier rate over 25 period 3 only trimmed -0.195 0.006 77 77 upper tier rate over 25 period 3 only raw -0.160 0.014 77 157 mean west period 2 only trimmed -1.857 0.040 57 57 upper tier rate central period 3 only trimmed -0.189 0.040 53 53 ``` While the true difference is zero, 4 of the 108 combinations fall below the 0.05 threshold and 9 below 0.10. All four are products of the analysis; there is no effect in the dataset that corresponds to them. The real information of the table is not the threshold count but the distribution of two columns. In the mean metric the observed difference wanders between -2.113 and 1.119 m³; the 0.123 m³ given by the primary decision stands in the middle of this range. Of the combinations, 45 point toward an increase and 63 toward a decrease. From the same dataset both the sentences "the notification reduced consumption" and "the notification increased consumption" can be supported with numbers. The data does not speak; which combination is written speaks. What the four significant rows have in common shows in the third and fourth columns. All four come from the narrowest breakdowns: the observation counts are between 53 and 157, whereas the primary decision worked on 1,101 subscribers. Three are in the narrowest observation window, three in the trimmed balance. In a small dataset the margin of variation grows, drawing a large observed difference from inside a large margin becomes easy, and the threshold is crossed more often. Trimming works in the same direction: reducing the B group of 858 subscribers to 402 makes the sample smaller. #### Pre-registration Since the cost of the freedom can be counted, closing it can be counted too. The four axes are written before the experiment one by one, and at every step it is looked at how many combinations remain. ```python # From how many to how many does pre-registration bring the # number of combinations. WRITTEN = (("metric", 0, "mean"), ("window", 2, "periods 2-3"), ("breakdown", 1, "whole dataset"), ("balance", 3, "raw")) remaining = SWEEP print("written before experiment combinations left p <= 0.05 smallest p") print(f"{'none':<27} {len(remaining):>17} {sum(1 for s in remaining if s[5] <= 0.05):>9}" f" {min(s[5] for s in remaining):>10.3f}") written = [] for name, i, choice in WRITTEN: written.append(name) remaining = [s for s in remaining if s[i] == choice] print(f"{'+ ' + name:<27} {len(remaining):>17} {sum(1 for s in remaining if s[5] <= 0.05):>9}" f" {min(s[5] for s in remaining):>10.3f}") print(f"\nclosed combinations: {len(SWEEP)} -> {len(remaining)} " f"({len(SWEEP) - len(remaining)} choices are eliminated before the experiment)") smallest = min(SWEEP, key=lambda s: s[5]) primary = [s for s in SWEEP if s[:4] == ("mean", "whole dataset", "periods 2-3", "raw")][0] print(f"\ntwo legitimate analyses of the same dataset:") print(f" {primary[0]} / {primary[1]} / {primary[2]} / {primary[3]}: " f"difference {primary[4]:.3f}, p {primary[5]:.3f}, {primary[6]} + {primary[7]} subscribers") print(f" {smallest[0]} / {smallest[1]} / {smallest[2]} / {smallest[3]}: " f"difference {smallest[4]:.3f}, p {smallest[5]:.3f}, {smallest[6]} + {smallest[7]} subscribers") print("true difference: 0.000 (the assignment did not change the readings)") ``` ``` written before experiment combinations left p <= 0.05 smallest p none 108 4 0.006 + metric 54 1 0.040 + window 18 0 0.214 + breakdown 2 0 0.758 + balance 1 0 0.772 closed combinations: 108 -> 1 (107 choices are eliminated before the experiment) two legitimate analyses of the same dataset: mean / whole dataset / periods 2-3 / raw: difference 0.123, p 0.772, 344 + 757 subscribers upper tier rate / over 25 / period 3 only / trimmed: difference -0.195, p 0.006, 77 + 77 subscribers true difference: 0.000 (the assignment did not change the readings) ``` The ladder descends in four steps. Writing the metric alone brings 108 combinations to 54 and the number of combinations crossing the threshold from 4 to 1. When the observation window is written too, 18 combinations remain and none crosses the threshold; the smallest p rises from 0.006 to 0.214. When all four are written, the number of combinations is 1 and 107 choices have been eliminated before the experiment. What is closed is not a flaw of method. Every value of the four axes is defensible on its own; the problem is that the value can be chosen **after the result is seen**. A decision written before the experiment reverses this order and makes the choice independent of the result. The last block also shows the margin of the margin itself. The p-value of the primary combination came out as 0.786 in the first block and 0.772 in this one. The two numbers are two runs of the same test; the difference between them is the own margin of variation of 500 permutations. If the permutation count is increased this difference narrows, it does not vanish. If a p-value is not written together with how many permutations produced it, the third digit of that p-value should not be read. The rule of the course grows by one item here: a p-value that does not carry next to it **how many combinations were walked** is not the p-value of a single test. #### Summary - The dataset of a finished-design experiment has four analytic degrees of freedom: metric (2), breakdown (9), observation window (3), handling of the unbalanced assignment (2). The number of combinations is 108. - The single pre-registered decision gives a 0.123 m³ difference, a margin of variation of [-0.698, 0.952] m³ and a p-value of 0.786; the number comes from the mean of 344 and 757 subscribers. - While the true difference between the groups is zero, 4 of the 108 combinations fall below 0.05 and 9 below 0.10; 45 combinations point toward an increase and 63 toward a decrease, and the difference in the mean metric wanders between -2.113 and 1.119 m³. - The four combinations that cross the threshold rest on 53 to 157 observations; in a small breakdown the margin of variation grows and the threshold is crossed more often. - Pre-registration brings the combinations from 108 to 1; 107 choices are eliminated and the smallest p rises from 0.006 to 0.772. - Two runs of the same test give 0.786 and 0.772; this difference is the own margin of 500 permutations. #### Next Step This lesson counted in how many combinations a difference was measured, but what it counted was always the difference between two groups. How consumption changes not with a notification but together with a continuous variable was never asked. The next lesson models the relationship directly: least squares is written inside the lesson, a coefficient is computed, and the margin is printed next to the coefficient with bootstrap. Then a confounder is added to the dataset and the sign of the same coefficient turns around; it is shown side by side that the coefficient of determination does not announce this turning. ### Modeling Relationships with Regression Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/modeling-relationships-with-regression Least squares is written inside the lesson, and what a coefficient says is separated from what it does not say: the univariate coefficient of meter age comes out as -0.605 m³/year, the margin of variation [-0.714, -0.503] and the coefficient of determination 0.098. When the base level of the region is added as a confounder, the same coefficient turns to 1.825, its margin becomes [1.629, 2.003] and the two intervals never overlap. In each of the five regions the within-group slope is positive (between 0.771 and 2.416). The coefficient of determination rises from 0.098 to 0.505 but does not announce the sign; when a completely random column is added it rises again, from 0.505363 to 0.505377. The regional mean of the residuals wanders between -3.020 and 2.481 in the univariate model and drops to between -0.229 and 0.295 when the confounder is added. The previous lesson counted in how many combinations a difference was measured, but what it counted was always the difference between two groups. How consumption changes together with a continuous variable was never asked. When a line is asked about instead of two groups, the number that comes out is not a difference but a **coefficient**: by how many units consumption changes on average when the variable grows by one unit. The number of this lesson is **the margin of variation of the coefficient**. A second number stands next to it and it is the real subject of the lesson: the same coefficient changes sign when a single variable is added to the dataset. Prediction accuracy measures are not computed here; mean absolute error, bias and the band were measured in the Estimation and Planning lesson of the Process, Team and Delivery course. Here the subject is the **interpretation and the limit** of the coefficient. - **AS31.** The dataset is inherited from K03 and is **fictional**. It is reduced to the subscriber level: the value of each subscriber is the mean of that subscriber's periods. Of 1,329 subscribers, 1,260 have at least one reading. - **AS32.** The meter age field is **fictional** too and carries two components together: a component tied to the region and a component tied to the within-region consumption deviation. The aim is to collect two directions in a single variable. - **AS33.** Least squares is **written inside the lesson**: the normal equations are built and solved by elimination with partial pivoting. No ready-made solver is used. - **AS34.** The margin of variation of the coefficient is produced by **drawing with replacement**: the subscriber list is redrawn 1,000 times and the model is rebuilt at every draw. - **AS35.** The confounder is **the base level of the region**: the base consumption of the region determines both the meter age and the consumption. #### A Coefficient and Its Margin Least squares looks for the coefficients that make the sum of the squares of the residuals smallest. This condition gives the normal equations, and the system of equations is solved by elimination. The block below builds the dataset, writes the solver and computes the univariate model. ```python # regression.py — MODEL. The fictional dataset of K03 is produced with the same seed and is # reduced to the subscriber level. The meter age field is FICTIONAL too: it carries a component # tied to the region and a component tied to the within-region consumption deviation together. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t): x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value(): s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u, w): t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 SUBSCRIBERS, DROPPED = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]) satisfied = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) (DROPPED if r() < 0.046 else SUBSCRIBERS).append( {"subscriber_id": 10001 + i, "region": b[0], "base": b[2], "household": household, "satisfaction": satisfied}) VALID = {k["subscriber_id"] for k in SUBSCRIBERS} READINGS = {} for k in SUBSCRIBERS + DROPPED: r = rng(SEED + 7000 + k["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor( k["base"] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if k["subscriber_id"] in VALID: READINGS.setdefault(k["subscriber_id"], []).append(t) r() DATA = [dict(k, consumption=sum(READINGS[k["subscriber_id"]]) / len(READINGS[k["subscriber_id"]])) for k in SUBSCRIBERS if k["subscriber_id"] in READINGS] REGION_MEAN = {} for x in DATA: REGION_MEAN.setdefault(x["region"], []).append(x["consumption"]) REGION_MEAN = {b: sum(v) / len(v) for b, v in REGION_MEAN.items()} for x in DATA: r = rng(SEED + 91000 + x["subscriber_id"]) y = (14.0 - 0.75 * (x["base"] - 20.5) + 0.15 * (x["consumption"] - REGION_MEAN[x["region"]]) + (r() + r() + r() - 1.5) * 2.4) x["age"] = max(1, min(30, int(y + 0.5))) def solve(A, b): """Solves the normal equations by elimination with partial pivoting.""" n = len(b) M = [A[i][:] + [b[i]] for i in range(n)] for i in range(n): p = max(range(i, n), key=lambda k: abs(M[k][i])) M[i], M[p] = M[p], M[i] for k in range(i + 1, n): f = M[k][i] / M[i][i] for j in range(i, n + 1): M[k][j] -= f * M[i][j] x = [0.0] * n for i in range(n - 1, -1, -1): x[i] = (M[i][n] - sum(M[i][j] * x[j] for j in range(i + 1, n))) / M[i][i] return x def ols(rows, fields): """Least squares: the coefficients and the coefficient of determination.""" p = len(fields) + 1 A, b = [[0.0] * p for _ in range(p)], [0.0] * p for s in rows: v = [1.0] + [float(s[a]) for a in fields] for i in range(p): for j in range(p): A[i][j] += v[i] * v[j] b[i] += v[i] * s["consumption"] coefs = solve(A, b) avg = sum(s["consumption"] for s in rows) / len(rows) sst = sum((s["consumption"] - avg) ** 2 for s in rows) sse = sum((s["consumption"] - sum(k * x for k, x in zip(coefs, [1.0] + [float(s[a]) for a in fields]))) ** 2 for s in rows) return coefs, 1 - sse / sst RESAMPLES = 1000 rs = rng(SEED + 77000) def coefficient_interval(rows, fields, index): """Drawing with replacement: the subscriber list is redrawn 1,000 times.""" n, drawn = len(rows), [] for _ in range(RESAMPLES): sample = [rows[int(rs() * n)] for _ in range(n)] drawn.append(ols(sample, fields)[0][index]) drawn.sort() return drawn[24], drawn[974] print(f"fictional dataset (seed {SEED}): {len(SUBSCRIBERS)} subscribers, " f"{sum(len(v) for v in READINGS.values())} readings") print(f"with at least one reading: {len(DATA)} subscribers") coefs, r2 = ols(DATA, ["age"]) low, high = coefficient_interval(DATA, ["age"], 1) print(f"\nconsumption = {coefs[0]:.3f} + ({coefs[1]:.3f}) x meter_age") print(f" meter age coefficient : {coefs[1]:.3f} m3/year") print(f" margin of variation : [{low:.3f}, {high:.3f}] ({RESAMPLES} draws with replacement)") print(f" explained variation : {r2:.4f} ({len(DATA)} observations)") ``` ``` fictional dataset (seed 20260218): 1329 subscribers, 3199 readings with at least one reading: 1260 subscribers consumption = 28.352 + (-0.605) x meter_age meter age coefficient : -0.605 m3/year margin of variation : [-0.714, -0.503] (1000 draws with replacement) explained variation : 0.0984 (1260 observations) ``` This is the line to be reported: when the meter age grows by one year, mean consumption falls by 0.605 m³, the margin of variation is [-0.714, -0.503]. The interval is far from zero and narrow; it rests on 1,260 observations. This number looks as if it satisfies the rule of the course. The third number next to the coefficient is the coefficient of determination, 0.098. The meter age carries about a tenth of the variation in consumption. This low value is often read as "the relationship is weak" and turns into a one-line warning in the report. That reading is incomplete. The coefficient of determination says something about the strength of the relationship; it says nothing about the **direction** of the relationship. #### A Confounder Turns the Coefficient Around The base consumption level of the region determines two things at once: the higher the consumption in a region, the newer the meters, and the higher the consumption. This definition is a confounder. When it enters the model, what the coefficient is can be measured. ```python # Confounder: the base level of the region determines both the meter age and the consumption. for x in DATA: r = rng(SEED + 99000 + x["subscriber_id"]) x["random"] = r() * 10 MODELS = [("age", ["age"]), ("age + base", ["age", "base"]), ("age + base + satisfaction", ["age", "base", "satisfaction"]), ("age + base + random", ["age", "base", "random"])] print("model age coefficient explained variation") for name, fields in MODELS: k, r2 = ols(DATA, fields) print(f"{name:<28} {k[1]:>15.3f} {r2:>19.6f}") low2, high2 = coefficient_interval(DATA, ["age", "base"], 1) k1, _ = ols(DATA, ["age"]) k2, _ = ols(DATA, ["age", "base"]) print(f"\ncoefficient observed margin of variation") print(f"univariate {k1[1]:>8.3f} [{low:.3f}, {high:.3f}]") print(f"confounder added {k2[1]:>8.3f} [{low2:.3f}, {high2:.3f}]") print("\nregion n mean meter age mean consumption within-group slope") for b in sorted(REGION_MEAN, key=lambda b: -REGION_MEAN[b]): sub = [x for x in DATA if x["region"] == b] e = ols(sub, ["age"])[0][1] print(f"{b:<7} {len(sub):>3} {sum(x['age'] for x in sub)/len(sub):>14.2f}" f" {REGION_MEAN[b]:>16.3f} {e:>19.3f}") ``` ``` model age coefficient explained variation age -0.605 0.098357 age + base 1.825 0.505363 age + base + satisfaction 1.827 0.505764 age + base + random 1.825 0.505377 coefficient observed margin of variation univariate -0.605 [-0.714, -0.503] confounder added 1.825 [1.629, 2.003] region n mean meter age mean consumption within-group slope east 219 9.95 24.814 2.416 central 231 12.13 22.483 2.154 north 363 13.60 20.351 1.777 south 258 16.69 16.728 1.550 west 189 18.84 13.935 0.771 ``` The coefficient turned from -0.605 to 1.825. The two margins of variation, [-0.714, -0.503] and [1.629, 2.003], do not overlap at all. Both numbers are correct and both come from the same dataset; the only thing that differs is which variables are in the model. The third table shows the mechanism in the open. Between regions the direction is clear: while the mean age of the meters rises from 9.95 to 18.84, mean consumption falls from 24.814 to 13.935. **Within** the regions, the slope is positive in all five of the five regions; the lowest is 0.771, the highest 2.416. The univariate model collects these two directions in a single number, and because the direction between regions weighs more, the result comes out negative. When the confounder enters the model, the direction between regions is separated out and the direction within regions remains. The first table gives two separate ways in which the coefficient of determination does not speak. The first: while the coefficient of determination rose from 0.098 to 0.505, the coefficient changed sign. The rise of the coefficient of determination does not say which sign is right; it says only that the residuals got smaller. The second: when a completely randomly generated column is added to the model, the coefficient of determination **rises** from 0.505363 to 0.505377. Because least squares can give every new column a coefficient of zero at worst, this coefficient never falls. A coefficient of determination whose variable count is not written cannot be compared. A coefficient is read as "while the other variables in the model are held constant", and the whole of this reading depends on what is inside the model. What the coefficient says is a relationship; it is not causation, and it carries no information about what a variable that is not in the model does. #### Residuals and the Width of the Margin The assumption check is done here. If the model is built correctly, the residuals should gather around zero in every breakdown and their dispersion should not wander together with the prediction level. ```python # Assumption check: the mean of the residuals in a breakdown and their dispersion against the prediction. def residuals(rows, fields): coefs = ols(rows, fields)[0] result = [] for s in rows: v = [1.0] + [float(s[a]) for a in fields] prediction = sum(k * x for k, x in zip(coefs, v)) result.append((s, prediction, s["consumption"] - prediction)) return result R1, R2 = residuals(DATA, ["age"]), residuals(DATA, ["age", "base"]) print("region n residual mean (age) residual mean (age + base)") for b in sorted(REGION_MEAN, key=lambda b: -REGION_MEAN[b]): d1 = [a for s, _, a in R1 if s["region"] == b] d2 = [a for s, _, a in R2 if s["region"] == b] print(f"{b:<7} {len(d1):>3} {sum(d1)/len(d1):>19.3f} {sum(d2)/len(d2):>26.3f}") def std_dev(d): o = sum(d) / len(d) return math.sqrt(sum((x - o) ** 2 for x in d) / len(d)) print("\nprediction band n residual mean residual dispersion") ordered = sorted(R2, key=lambda t: t[1]) step = len(ordered) // 4 for i in range(4): chunk = ordered[i * step:(i + 1) * step if i < 3 else len(ordered)] d = [a for _, _, a in chunk] print(f"{chunk[0][1]:>6.1f} - {chunk[-1][1]:>5.1f} {len(chunk):>5} {sum(d)/len(d):>16.3f}" f" {std_dev(d):>19.3f}") # The same coefficient, with a quarter of the observations r4 = rng(SEED + 88000) QUARTER = [x for x in DATA if r4() < 0.25] low3, high3 = coefficient_interval(QUARTER, ["age", "base"], 1) k3 = ols(QUARTER, ["age", "base"])[0][1] print(f"\n observations coefficient margin of variation width") print(f"{len(DATA):>13} {k2[1]:>11.3f} [{low2:.3f}, {high2:.3f}] {high2 - low2:>8.3f}") print(f"{len(QUARTER):>13} {k3:>11.3f} [{low3:.3f}, {high3:.3f}] {high3 - low3:>8.3f}") ``` ``` region n residual mean (age) residual mean (age + base) east 219 2.481 -0.229 central 231 1.468 0.295 north 363 0.227 0.031 south 258 -1.526 -0.118 west 189 -3.020 0.005 prediction band n residual mean residual dispersion 8.7 - 16.5 315 0.293 3.891 16.5 - 19.7 315 -0.168 3.988 19.7 - 23.3 315 -0.346 4.852 23.3 - 34.3 315 0.221 4.893 observations coefficient margin of variation width 1260 1.825 [1.629, 2.003] 0.374 313 1.665 [1.341, 1.989] 0.648 ``` The first table catches the flaw of the univariate model without looking at the coefficient. The regional mean of the residuals wanders between -3.020 and 2.481: the model systematically underestimates the east region and systematically overestimates the west region. A pattern that remains in a set of residuals by breakdown is the signature of a variable the model left out. When the confounder is added, the same numbers drop to between -0.229 and 0.295. The examination of residuals says which model is missing something before the result is seen. The second table shows that an assumption is not met. The dispersion of the residuals grows together with the prediction level: 3.891 in the lowest band, 4.893 in the highest. The coefficient can still be computed in this case, but a single margin of variation is not valid for the whole prediction range; among high-consumption subscribers the margin is wider. Drawing with replacement carries this irregularity by itself, because it redraws the dataset as it is. The third table gives what the margin depends on. When the same model is built with a quarter of the observations, the coefficient slips from 1.825 to 1.665 and the width of the margin of variation rises from 0.374 to 0.648. The observation count fell to a quarter, the width roughly doubled. Reporting a coefficient computed in a small breakdown without writing its margin is showing the two rows of this table as a single row. #### Summary - The univariate coefficient of meter age is -0.605 m³/year, the margin of variation [-0.714, -0.503] and the coefficient of determination 0.098; the number rests on 1,260 observations. - When the base level of the region is added as a confounder, the same coefficient becomes 1.825 and its margin comes out as [1.629, 2.003]; the two intervals never overlap. - The mechanism is visible in the breakdown: between regions the mean age rises while mean consumption falls, whereas in all five of the five regions the within-group slope is positive (between 0.771 and 2.416). - While the coefficient of determination rose from 0.098 to 0.505, the coefficient changed sign; the coefficient of determination does not announce the direction, and when a completely random column is added it rises again from 0.505363 to 0.505377. - The regional mean of the residuals wanders between -3.020 and 2.481 in the univariate model and drops to between -0.229 and 0.295 when the confounder is added; the dispersion of the residuals grows with the prediction from 3.891 to 4.893. - With a quarter of the observations, the margin of the same coefficient grows in width from 0.374 to 0.648. #### Next Step This lesson showed that a coefficient changes direction by breakdown, but the breakdown was always the region and time never entered the work. Subscribers do not join the network at the same time and are not observed for the same length of time. The next lesson splits the dataset by entry period: the retention rate is computed, then by changing the entry window and the observation window it is counted how many different rates come out of the same dataset. The same change of direction appears here once more, this time not in a coefficient but in a rate: the aggregated rate shows one direction, the rate of every cohort shows the opposite. ### Cohort and Funnel Analysis Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/cohort-and-funnel-analysis Counting how many separate retention rates the cohort definition produces from the same dataset: the rule that counts the unobservable cohort as zero gives 0.619, the rule that excludes it gives 0.867, and outside the three-period window 843 of the 3,199 readings remain. Taking the entry window as a period gives two observable cohorts and a rate range of 0.009; taking it as five weeks gives six cohorts and a range of 0.040; since the narrowest margin of variation is 0.087, the whole range fits inside the margin of a single cohort. The funnel enters only as far as the aggregation reversal: while the aggregated difference is -0.237, the difference in both of the two observable cohorts is +0.083 and +0.028. The previous lesson showed that a coefficient changes direction by breakdown, but the breakdown was always the region and time never entered the work. Subscribers do not join the network at the same time and are not observed for the same length of time. Still being able to read a subscriber one period after the registration opened is a success, and the rate of this success is called **retention**. The number of this lesson is **the retention rate that the cohort definition moves**. A cohort is the set built by entry period, and its definition contains two decisions: how wide the entry window is and where the observation window closes. Both decisions are defensible, and both produce a different rate from the same dataset. The step computation of the conversion funnel is not repeated here; step conversion, drop-off rate and cumulative conversion were measured in the Conversion Funnel lesson of the User Experience course. The funnel enters this lesson only as far as the **aggregation reversal**. - **AS36.** The dataset is inherited from K03 and is **fictional**: 3,199 readings, 1,329 subscribers, three periods, five regions. - **AS37.** The entry week field is **fictional** too: the measurement network did not reach the regions at the same time, so the distribution of the entry period depends on the region. - **AS38.** In this lesson the reading index gives **the period after entry**; the calendar period is the sum of the entry period and the index minus one, and the calendar **closes at three periods**. - **AS39.** Retention is the share of the subscribers of a cohort who still produce a reading $K$ periods after entry. If the period $K$ later falls outside the calendar, the rate **cannot be observed**. - **AS40.** Two rules can be applied to an unobservable subscriber: the **zero-counting** rule counts the subscriber as not retained, the **excluding** rule removes the subscriber from the denominator. - **AS41.** The margin of variation is produced by **drawing with replacement**: the cohort list is redrawn 1,000 times. - **AS42.** The reading method field is **fictional** too and carries two relationships together: remote reading has become widespread in the late periods and has fewer interruptions. #### The Observation Window Cuts the Cohort The calendar closes at three periods. For a subscriber who joined the network in the third period, the answer to the question "can this subscriber still be read one period later" is not in the dataset. This subscriber has not failed to be retained; the subscriber **has not been measured yet**. ```python # cohort.py — MODEL. The fictional dataset of K03 is produced with the same seed. The entry # week field is FICTIONAL too: the measurement network did not reach the regions at the same # time, so the distribution of the entry period depends on the region. The reading index gives # the period after entry; calendar period = entry period + index - 1 and the calendar closes # at three periods. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] ENTRY = {"east": [0.60, 0.28, 0.12], "central": [0.52, 0.30, 0.18], "north": [0.40, 0.34, 0.26], "south": [0.28, 0.36, 0.36], "west": [0.16, 0.32, 0.52]} def rng(t): x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value(): s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u, w): t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 SUBSCRIBERS, DROPPED = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]) pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) (DROPPED if r() < 0.046 else SUBSCRIBERS).append( {"subscriber_id": 10001 + i, "region": b[0], "base": b[2]}) VALID = {k["subscriber_id"] for k in SUBSCRIBERS} READINGS = {} for k in SUBSCRIBERS + DROPPED: r = rng(SEED + 7000 + k["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor( k["base"] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 if k["subscriber_id"] in VALID: READINGS.setdefault(k["subscriber_id"], []).append(t) r() JOINED = [] for k in SUBSCRIBERS: if k["subscriber_id"] not in READINGS: continue r = rng(SEED + 61000 + k["subscriber_id"]) e = pick(r(), ENTRY[k["region"]]) week = e * 13 + 1 + int(r() * 13) JOINED.append({"id": k["subscriber_id"], "region": k["region"], "week": week, "cohort": e + 1, "block": (week - 1) // 5 + 1, "n": len(READINGS[k["subscriber_id"]])}) WINDOW = 3 def observable(x, K): """Does the cohort stay inside the three-period calendar K periods later.""" return x["cohort"] + K <= WINDOW def retained(x, K): return x["n"] >= K + 1 and observable(x, K) def rate(group, K, rule): denominator = group if rule == "zero-counting" else [x for x in group if observable(x, K)] if not denominator: return None, 0 return sum(1 for x in denominator if retained(x, K)) / len(denominator), len(denominator) print(f"fictional dataset (seed {SEED}): {len(SUBSCRIBERS)} subscribers, " f"{sum(len(v) for v in READINGS.values())} readings") print("entry cohort (period): " + ", ".join( f"{k}: {sum(1 for x in JOINED if x['cohort'] == k)}" for k in (1, 2, 3)) + f" ({len(JOINED)} subscribers with at least one reading)") print(f"readings outside the three-period window: " f"{sum(max(0, x['cohort'] + x['n'] - 1 - WINDOW) for x in JOINED)}" f" / {sum(x['n'] for x in JOINED)}") for K in (1, 2): print(f"\nretention ({K} period after entry)") print("cohort n observable zero-counting excluding") for coh in (1, 2, 3): group = [x for x in JOINED if x["cohort"] == coh] o1, n1 = rate(group, K, "zero-counting") o2, n2 = rate(group, K, "excluding") print(f"{coh:>6} {len(group):>6} {n2:>14} {o1:>13.3f} " + (f"{o2:>9.3f}" if o2 is not None else f"{'--':>9}")) o1, _ = rate(JOINED, K, "zero-counting") o2, n2 = rate(JOINED, K, "excluding") print(f"{'total':>6} {len(JOINED):>6} {n2:>14} {o1:>13.3f} {o2:>9.3f}") ``` ``` fictional dataset (seed 20260218): 1329 subscribers, 3199 readings entry cohort (period): 1: 488, 2: 412, 3: 360 (1260 subscribers with at least one reading) readings outside the three-period window: 843 / 3199 retention (1 period after entry) cohort n observable zero-counting excluding 1 488 488 0.871 0.871 2 412 412 0.862 0.862 3 360 0 0.000 -- total 1260 900 0.619 0.867 retention (2 period after entry) cohort n observable zero-counting excluding 1 488 488 0.648 0.648 2 412 0 0.000 -- 3 360 0 0.000 -- total 1260 488 0.251 0.648 ``` The number to enter the report is 0.619 or 0.867 for one-period retention. Both come from the same dataset and the same definition; the only place where they part is what the 360 unobservable subscribers are counted as. The zero-counting rule takes them as not retained and pulls the rate down by 0.248. For two-period retention the difference is even larger: 0.251 against 0.648. Counting an unobservable subscriber as zero is not a measurement but an assumption, and this assumption is not met in the dataset. For the 360 subscribers of cohort 3, the reading of the next period has **never been collected**; had it been collected, it would be expected to give a number close to the rates 0.871 and 0.862 of the first two cohorts. The third row gives the size of this cutting: 843 of the 3,199 readings fall outside the three-period calendar. The cohort definition cuts the dataset, and the part that is cut off is not small. #### The Entry Window Changes the Number of Cohorts The second decision of the cohort is the width of the entry window. If blocks of five weeks are taken instead of periods, the same subscribers distributed over eight cohorts. Next to the rate of every cohort, its margin is printed too. ```python # When the entry window changes, how many separate rates come out of the same dataset. RESAMPLES = 1000 rs = rng(SEED + 44000) def rate_interval(group, K): """Drawing with replacement: the cohort list is redrawn 1,000 times.""" n, drawn = len(group), [] for _ in range(RESAMPLES): drawn.append(sum(1 for _ in range(n) if retained(group[int(rs() * n)], K)) / n) drawn.sort() return drawn[24], drawn[974] def scan(field, label): keys = sorted({x[field] for x in JOINED}) print(f"\nentry window: {label} -> {len(keys)} cohorts " "(excluding rule, 1 period)") print("cohort n retention margin of variation width") widths, rates = [], [] for a in keys: group = [x for x in JOINED if x[field] == a and observable(x, 1)] if len(group) < 10: print(f"{a:>6} {len(group):>6} {'unobservable':>12}") continue o = sum(1 for x in group if retained(x, 1)) / len(group) d, u = rate_interval(group, 1) widths.append(u - d) rates.append(o) print(f"{a:>6} {len(group):>6} {o:>10.3f} [{d:.3f}, {u:.3f}] {u - d:>10.3f}") print(f" rate range {max(rates) - min(rates):.3f}, " f"narrowest margin {min(widths):.3f}, widest margin {max(widths):.3f}") return rates d3 = scan("cohort", "period (13 weeks)") d8 = scan("block", "five weeks") print(f"\nnumber of retention rates from the same dataset: {len(d3)} and {len(d8)}") print(f"range with the period definition {min(d3):.3f} - {max(d3):.3f}, " f"with the five-week definition {min(d8):.3f} - {max(d8):.3f}") ``` ``` entry window: period (13 weeks) -> 3 cohorts (excluding rule, 1 period) cohort n retention margin of variation width 1 488 0.871 [0.840, 0.898] 0.057 2 412 0.862 [0.828, 0.896] 0.068 3 0 unobservable rate range 0.009, narrowest margin 0.057, widest margin 0.068 entry window: five weeks -> 8 cohorts (excluding rule, 1 period) cohort n retention margin of variation width 1 181 0.878 [0.829, 0.923] 0.094 2 185 0.854 [0.805, 0.903] 0.097 3 184 0.891 [0.848, 0.935] 0.087 4 174 0.856 [0.805, 0.908] 0.103 5 155 0.852 [0.800, 0.903] 0.103 6 21 0.857 [0.667, 1.000] 0.333 7 0 unobservable 8 0 unobservable rate range 0.040, narrowest margin 0.087, widest margin 0.333 number of retention rates from the same dataset: 2 and 6 range with the period definition 0.862 - 0.871, with the five-week definition 0.852 - 0.891 ``` When the entry window narrows, the number of reportable rates rises from 2 to 6 and the range between the rates widens from 0.009 to 0.040. A person looking at a cohort matrix may read this widening as a pattern: the third block is 0.891, the fifth block 0.852, and the 4-point drop between them looks like a trend. The column of margins closes this reading. The narrowest margin of variation is 0.087 wide; the whole of the 0.040 range between the six blocks **fits inside the margin of a single cohort**. The pairwise differences between the blocks also stay inside the margins. The thin cohort definition did not produce new information; it divided the same information into more numbers and enlarged the margin of each number. The sixth block is the extreme example of this. The rate resting on twenty-one observable subscribers is 0.857 and its margin is [0.667, 1.000]; its width is 0.333, that is, more than three times the margin of the five cohorts. The seventh and eighth blocks stay entirely outside the observation window and produce no rate. As the cohort definition thins, the number of unobservable cohorts grows too. #### The Aggregated Rate and the Rate of the Cohorts What a cohort is good for shows when it is placed next to an aggregated rate. There are two reading methods in the dataset and the two steps of the funnel are these: entry to the cohort, reading in the next period. ```python # The reading method field is FICTIONAL too and carries two relationships together: remote # reading has become widespread in the late periods and has fewer interruptions. METHOD_BASE = {1: 0.28, 2: 0.55, 3: 0.80} for x in JOINED: r = rng(SEED + 71000 + x["id"]) p = METHOD_BASE[x["cohort"]] + (0.10 if x["n"] >= 2 else -0.10) x["method"] = "remote" if r() < p else "on-site" METHODS = ("on-site", "remote") def difference_interval(a, b, K): na, nb, drawn = len(a), len(b), [] for _ in range(RESAMPLES): oa = sum(1 for _ in range(na) if retained(a[int(rs() * na)], K)) / na ob = sum(1 for _ in range(nb) if retained(b[int(rs() * nb)], K)) / nb drawn.append(ob - oa) drawn.sort() return drawn[24], drawn[974] print("funnel: entry to the cohort -> reading in the next period") print("\ncohort on-site reading remote reading difference margin of variation") for coh in (1, 2, 3): group = [x for x in JOINED if x["cohort"] == coh] a = [x for x in group if x["method"] == "on-site"] b = [x for x in group if x["method"] == "remote"] if not observable(group[0], 1): print(f"{coh:>6} {'--':>5} ({len(a):>3}) {'--':>5} ({len(b):>3})" f" -- unobservable") continue oa = sum(1 for x in a if retained(x, 1)) / len(a) ob = sum(1 for x in b if retained(x, 1)) / len(b) d, u = difference_interval(a, b, 1) print(f"{coh:>6} {oa:.3f} ({len(a):>3}) {ob:.3f} ({len(b):>3})" f" {ob - oa:>+.3f} [{d:+.3f}, {u:+.3f}]") a = [x for x in JOINED if x["method"] == "on-site"] b = [x for x in JOINED if x["method"] == "remote"] oa = sum(1 for x in a if retained(x, 1)) / len(a) ob = sum(1 for x in b if retained(x, 1)) / len(b) d, u = difference_interval(a, b, 1) print(f"{'total':>6} {oa:.3f} ({len(a):>3}) {ob:.3f} ({len(b):>3})" f" {ob - oa:>+.3f} [{d:+.3f}, {u:+.3f}]") print("\ncohort composition of the method (row share)") for y in METHODS: group = [x for x in JOINED if x["method"] == y] print(f" {y:<8} " + " ".join( f"cohort {k}: {sum(1 for x in group if x['cohort'] == k) / len(group):.3f}" for k in (1, 2, 3))) ``` ``` funnel: entry to the cohort -> reading in the next period cohort on-site reading remote reading difference margin of variation 1 0.843 (325) 0.926 (163) +0.083 [+0.028, +0.139] 2 0.844 (154) 0.872 (258) +0.028 [-0.041, +0.099] 3 -- ( 56) -- (304) -- unobservable total 0.755 (535) 0.519 (725) -0.237 [-0.284, -0.184] cohort composition of the method (row share) on-site cohort 1: 0.607 cohort 2: 0.288 cohort 3: 0.105 remote cohort 1: 0.225 cohort 2: 0.356 cohort 3: 0.419 ``` The aggregated row says a single sentence: retention is 0.237 points lower for remotely read meters and its margin is [-0.284, -0.184], away from zero. This number and this margin have been written in the form the rule of the course asks for, and they are wrong. The cohort rows show the opposite. In the first cohort remote reading is 0.083 points **ahead** and its margin is [+0.028, +0.139], not containing zero. In the second cohort it is 0.028 points ahead and its margin contains zero. The third cohort cannot be observed. **No cohort shows** the direction the aggregated rate shows: two say the opposite, one says nothing. This phenomenon is called the **aggregation reversal**, and its mechanism is in the last table. 60.7% of the on-site read meters are in the first cohort; 41.9% of the remotely read ones are in the third cohort. The third cohort is counted as zero because of the observation window, so the aggregated rate of remote reading comes mostly from subscribers who **could not be measured**. What has been measured is not the difference between the two methods but the cohort composition of the two methods. An aggregated rate describes no cohort if there is both a rate difference and a composition difference between the cohorts. An aggregated rate reported without its composition written is not a measurement but a decision to take an average. #### Summary - The rule applied to the unobservable cohort produces two rates from the same dataset: one-period retention is 0.619 with the zero-counting rule and 0.867 with the excluding rule; for two periods it is 0.251 and 0.648. - Outside the three-period observation window 843 of the 3,199 readings remain; the cohort definition cuts the dataset. - When the entry window goes down from a period to five weeks, the number of reportable rates rises from 2 to 6 and the rate range from 0.009 to 0.040; since the narrowest margin of variation is 0.087, the whole range fits inside the margin of a single cohort. - In the sixth block, which rests on twenty-one subscribers, the rate is 0.857, the margin [0.667, 1.000] and the width 0.333; the seventh and eighth blocks are entirely outside the observation window. - Aggregation reversal: while the aggregated difference is -0.237 and its margin is away from zero, the difference in the two observable cohorts is +0.083 and +0.028. 60.7% of those read on-site are in the first cohort, 41.9% of those read remotely are in the third cohort. #### Next Step In three lessons the margin was written next to the number and the three sources that enlarge the margin were counted: the number of combinations, the variable that is not in the model, the cohort definition. In all three the dataset was taken as given. The last lesson of the course removes this assumption: the dataset is made up of meters that could already be read and of subscribers who produced at least one reading, that is, what is measured is not the network itself but the visible part of the network. In the same lesson it is counted how many findings can be produced from noise, and this number is compared with the expected value. ### Misleading Analyses Source: https://academia.sh/en/courses/statistical-analysis-and-exploratory-data-analysis/misleading-analyses Counting the fact that the dataset itself is a choice and that noise can produce findings: the east-north difference is 4.659 m³ in the population of all meters in the network, 4.463 m³ in the sample of readable meters, 4.655 m³ when zero readings are dropped and 4.454 m³ among subscribers who stayed only three periods; the range of the four decisions is 0.206 m³. Of the 1,329 subscribers, 69 appear in no table. Ten scans of twenty breakdowns each run 200 tests on a completely random outcome variable and give 11 findings; the expected number is 10.0 and at least one finding appears in 7 of the ten scans. In the course wrap-up the observed number, margin of variation and whether the decision changes for the thirteen lessons are gathered in a single table. In three lessons the margin was written next to the number and the three sources that enlarge the margin were counted: the number of combinations, the variable that is not in the model, the cohort definition. In all three the dataset was taken as given. This lesson removes that assumption. The dataset is not the network itself but **the visible part of the network**: it is made up of meters that could be read and of subscribers who produced at least one reading. The numbers of this lesson are two. The first is **how much the elimination decisions move the reported difference**, the second is **the number of findings that can be produced from noise**. Both come down to the same thing: the margin of a number comes not only from the sample but also from how many times a choice was made on the way to that number. - **AS43.** The dataset is inherited from K03 and is **fictional**. In this lesson the readings of the eliminated meters are kept too, so the **population** and the **sample** can be put side by side. - **AS44.** The population is the readings of all meters in the network; the sample is the readings of the meters that could enter the record. The elimination in between is part of the production of the dataset and is **drawn independently of the measured value**. - **AS45.** Two more elimination decisions are tested: dropping the readings written as zero and counting only the subscribers who stayed for all three periods. - **AS46.** For the noise sweep the outcome variable is drawn **completely at random**; the breakdowns are taken from the real dataset. There is no real effect in any breakdown. - **AS47.** Twenty breakdowns, ten scans, 200 permutations, threshold 0.05. The threshold is a decision, and the expected number of findings is computed from this threshold. #### The Dataset Itself Is a Choice The anchor of the course is the difference between east and north. This difference changes with which rows are counted. ```python # misleading.py — MODEL. The fictional dataset of K03 is produced with the same seed, but this # time the readings of the eliminated meters are kept too: the population and the sample # are put side by side. import math SEED, RAW, M32 = 20260218, 1400, 0xFFFFFFFF REGION = [("north", 0.28, 21), ("south", 0.22, 17), ("east", 0.18, 26), ("west", 0.14, 14), ("central", 0.18, 23)] def rng(t): x = ((t ^ (t >> 16)) * 2246822507) & M32 x = ((x ^ (x >> 13)) * 3266489909) & M32 s = [(x ^ (x >> 16)) & M32] def next_value(): s[0] = (s[0] * 1664525 + 1013904223) & M32 return s[0] / 4294967296 return next_value def pick(u, w): t = 0.0 for i, x in enumerate(w): t += x if u < t: return i return len(w) - 1 SUBSCRIBERS, DROPPED = [], [] for i in range(RAW): r = rng(SEED + i) b = REGION[pick(r(), [x[1] for x in REGION])] household = 1 + pick(r(), [0.06, 0.24, 0.30, 0.24, 0.11, 0.05]) satisfaction = 1 + pick(r(), [0.08, 0.14, 0.27, 0.34, 0.17]) (DROPPED if r() < 0.046 else SUBSCRIBERS).append( {"subscriber_id": 10001 + i, "region": b[0], "base": b[2], "household": household, "satisfaction": satisfaction}) VALID = {k["subscriber_id"] for k in SUBSCRIBERS} POPULATION, SAMPLE = [], [] for k in SUBSCRIBERS + DROPPED: r = rng(SEED + 7000 + k["subscriber_id"]) for d in range(pick(r(), [0.05, 0.12, 0.21, 0.62])): t = 0.0 if r() < 0.038 else math.floor( k["base"] * math.exp((r() + r() + r() - 1.5) * 0.62) * (1 - d * 0.05) * 100 + 0.5) / 100 row = {"subscriber_id": k["subscriber_id"], "region": k["region"], "period": d + 1, "consumption_m3": t} POPULATION.append(row) if k["subscriber_id"] in VALID: SAMPLE.append(row) r() BY_SUBSCRIBER = {} for s in SAMPLE: BY_SUBSCRIBER.setdefault(s["subscriber_id"], []).append(s["consumption_m3"]) def summarize(rows, label): t, c = {}, {} for s in rows: t[s["region"]] = t.get(s["region"], 0.0) + s["consumption_m3"] c[s["region"]] = c.get(s["region"], 0) + 1 d, k = t["east"] / c["east"], t["north"] / c["north"] print(f"{label:<42} {len(rows):>6} {d:>8.3f} {k:>8.3f} {d - k:>8.3f}") return d - k print(f"{'dataset':<42} {'rows':>6} {'east':>8} {'north':>8} {'gap':>8}") gaps = [summarize(POPULATION, "population (all meters in the network)"), summarize(SAMPLE, "sample (meters that could be read)"), summarize([s for s in SAMPLE if s["consumption_m3"] > 0], "zero readings dropped"), summarize([s for s in SAMPLE if len(BY_SUBSCRIBER[s["subscriber_id"]]) == 3], "only subscribers who stayed three periods")] print(f"\nrange of the gap produced by the four decisions: {max(gaps) - min(gaps):.3f} m3") print(f"dropped subscribers: {len(DROPPED)}, dropped readings: {len(POPULATION) - len(SAMPLE)}") print(f"subscribers who produced no readings: {len(SUBSCRIBERS) - len(BY_SUBSCRIBER)} / {len(SUBSCRIBERS)}") print(f"zero readings: {sum(1 for s in SAMPLE if s['consumption_m3'] == 0)} / {len(SAMPLE)}") print(f"\n{'readings':>12} {'subscribers':>11} {'subscriber mean':>18} {'first-period mean':>22}") for n in (1, 2, 3): group = [v for v in BY_SUBSCRIBER.values() if len(v) == n] print(f"{n:>12} {len(group):>11} {sum(sum(v)/len(v) for v in group)/len(group):>18.3f}" f" {sum(v[0] for v in group)/len(group):>22.3f}") overall = [sum(v) / len(v) for v in BY_SUBSCRIBER.values()] print(f"{'overall':>12} {len(overall):>11} {sum(overall)/len(overall):>18.3f}") ``` ``` dataset rows east north gap population (all meters in the network) 3364 24.784 20.125 4.659 sample (meters that could be read) 3199 24.674 20.210 4.463 zero readings dropped 3080 25.593 20.938 4.655 only subscribers who stayed three periods 2514 24.604 20.150 4.454 range of the gap produced by the four decisions: 0.206 m3 dropped subscribers: 71, dropped readings: 165 subscribers who produced no readings: 69 / 1329 zero readings: 119 / 3199 readings subscribers subscriber mean first-period mean 1 159 20.858 20.858 2 263 19.776 19.818 3 838 19.627 20.895 overall 1260 19.813 ``` The first table separates the two faces of **selection bias**. The difference between the population and the sample goes down from 4.659 to 4.463 m³: 165 readings and 71 subscribers are missing. But this elimination was drawn independently of consumption in the production of the dataset, so the 0.196 m³ in between is not a bias but the sample's own margin. What decides whether an elimination produces bias is not the size of the elimination but **whether it is related to the measured value**. The third row shows a related elimination. When the 119 readings written as zero are counted as "invalid readings" and dropped, the east mean rises from 24.674 to 25.593 m³. The decision is defensible — a zero may not be a consumption but a record of failure to read — but because the zeros are not randomly distributed, it pushes the result in one direction. The range of the four decisions is 0.206 m³, and this is about one twentieth of the anchor difference of the course. The second table gives **survivorship bias**. The mean of the 838 subscribers who stayed for three periods is 19.627 m³, the mean of the whole dataset is 19.813 m³. The difference is small and its source is in the last column: the **first-period** mean of the subscribers who stayed three periods, 20.895 m³, is almost equal to the 20.858 m³ of those with a single period. So the survivors are not subscribers who were different from the start; what creates the difference is that the later-period readings of the long-stayers enter the mean too. The most important number lies outside the table: **69** of the 1,329 subscribers **never produced a reading** and appear in no subscriber table built in this lesson. Every summary computed per subscriber belongs to 1,260 subscribers. The quietest form of survivorship bias is that the denominator itself is invisible. #### Noise Produces Findings The second number is the garden of forking paths. The outcome variable is drawn completely at random, the breakdowns are taken from the real dataset, and a permutation test is run for every breakdown. There is nothing to find. ```python # Producing findings from noise: the outcome variable is drawn completely at random, the # breakdowns are taken from the real dataset. There is no real effect in any breakdown. RECORDS = {k["subscriber_id"]: k for k in SUBSCRIBERS if k["subscriber_id"] in BY_SUBSCRIBER} DATA = [{"region": k["region"], "household": k["household"], "satisfaction": k["satisfaction"], "tier": ("under 10" if BY_SUBSCRIBER[n][0] < 10 else ("10-25" if BY_SUBSCRIBER[n][0] < 25 else "over 25")), "zero": min(BY_SUBSCRIBER[n]) == 0.0, "consumption": sum(BY_SUBSCRIBER[n]) / len(BY_SUBSCRIBER[n])} for n, k in RECORDS.items()] N = len(DATA) BREAKDOWNS = [] for field in ("region", "household", "satisfaction", "tier"): for a in sorted({x[field] for x in DATA}, key=str): BREAKDOWNS.append((f"{field}={a}", [i for i, x in enumerate(DATA) if x[field] == a])) BREAKDOWNS.append(("with zero reading", [i for i, x in enumerate(DATA) if x["zero"]])) PERMUTATIONS, THRESHOLD, SCANS = 200, 0.05, 10 rs = rng(SEED + 33000) PERM = [] for _ in range(PERMUTATIONS): h = list(range(N)) for i in range(N - 1, 0, -1): j = int(rs() * (i + 1)) h[i], h[j] = h[j], h[i] PERM.append(h) def scan(values): """For each breakdown: permutation test of the difference between inside and outside the breakdown.""" total, findings = sum(values), [] for label, idx in BREAKDOWNS: m = len(idx) inside = sum(map(values.__getitem__, idx)) observed = inside / m - (total - inside) / (N - m) exceeded = 0 for h in PERM: p = sum(values[h[i]] for i in idx) if abs(p / m - (total - p) / (N - m)) >= abs(observed) - 1e-12: exceeded += 1 if exceeded / PERMUTATIONS <= THRESHOLD: findings.append((label, exceeded / PERMUTATIONS)) return findings rg = rng(SEED + 22000) print(f"{N} subscribers, {len(BREAKDOWNS)} breakdowns, {PERMUTATIONS} permutations, threshold {THRESHOLD}") print("noise sweep (outcome variable random, no real effect)") total, at_least_one, most = 0, 0, 0 for t in range(SCANS): b = scan([rg() for _ in range(N)]) total += len(b) at_least_one += 1 if b else 0 most = max(most, len(b)) print(f" scan {t + 1:>2}: {len(b)} found") if b: print(" " + ", ".join(f"{label} p={p:.3f}" for label, p in b)) tests = SCANS * len(BREAKDOWNS) print(f"\ntotal tests {tests}, observed findings {total}, expected {tests * THRESHOLD:.1f}") print(f"scans with at least one finding {at_least_one}/{SCANS}, " f"expected {SCANS * (1 - (1 - THRESHOLD) ** len(BREAKDOWNS)):.1f}") print(f"most findings in a single scan: {most}") print(f"same scan on real consumption: {len(scan([x['consumption'] for x in DATA]))} found") ``` ``` 1260 subscribers, 20 breakdowns, 200 permutations, threshold 0.05 noise sweep (outcome variable random, no real effect) scan 1: 1 found region=east p=0.050 scan 2: 0 found scan 3: 1 found with zero reading p=0.040 scan 4: 1 found region=north p=0.050 scan 5: 5 found region=north p=0.030, household=1 p=0.010, household=6 p=0.010, tier=10-25 p=0.000, tier=over 25 p=0.010 scan 6: 0 found scan 7: 1 found tier=over 25 p=0.015 scan 8: 0 found scan 9: 1 found satisfaction=4 p=0.000 scan 10: 1 found region=east p=0.035 total tests 200, observed findings 11, expected 10.0 scans with at least one finding 7/10, expected 6.4 most findings in a single scan: 5 same scan on real consumption: 8 found ``` The two hundred tests gave 11 findings; since the threshold is 0.05, the expected number is 10.0. The observed and the expected are almost the same, because the number of findings comes not from the dataset but from the threshold and the number of tests. The number of scans with at least one finding is 7, expected 6.4. The probability that a scan of twenty breakdowns finds nothing is smaller than one third even when there is no effect at all. The fifth scan deserves to be read on its own: the same random column crosses the threshold in five separate breakdowns, and the p value of one of them is 0.000. These five rows can enter a report as "significant differences were found in the region, household size and tariff tier breakdowns". Behind all five there is nothing. The last line gives the thing that makes the distinction impossible. The same scan finds 8 findings on real consumption and 1.1 per scan on average on noise. Looked at in aggregate, the real structure can be told apart; but a person looking at **a single finding** cannot see which dataset it came from. The only thing that makes the distinction is that the number of breakdowns tried has been written down. A finding whose number of tried breakdowns is not written cannot be told apart from one row of a noise scan. This is also the place to pay the debt the exploration topic left behind. The exploratory analysis had walked through 232 cells in five breakdown fields and had counted that 33 of them had fewer than 30 rows; the question to be tested was chosen from among those 232 cells. If the scan above produces 11 findings from 200 tests on 20 breakdowns, the expected number of findings an exploratory scan of 232 cells would produce at the same threshold is 11.6 — and none of those findings was counted as a test in the exploration stage. The debt is this: when the number of cells walked through in the exploration is not written in the report, a single test done afterwards shows its own margin narrower than it is. Its payment is not a computation but a sentence: it is written how many breakdowns were walked through. Selection bias and p-hacking come together here. The first is choosing which rows enter the dataset, the second is choosing which question is asked of the dataset. Both can be done after the result has been seen, and both stay invisible unless they are written down. #### Summary - The same anchor difference takes four values with four decisions: 4.659 in the population, 4.463 in the sample, 4.655 when zero readings are dropped and 4.454 m³ among subscribers who stayed only three periods; the range is 0.206 m³. - What produces bias is not the size of the elimination but the relationship of the elimination with the measured value: the elimination of 71 subscribers, drawn independently of consumption, moves the difference by 0.196 m³, and dropping the zero readings moves the east mean from 24.674 to 25.593 m³. - Survivorship bias is in two places: the mean of the 838 subscribers who stayed three periods, 19.627 m³, is lower than the 19.813 m³ of the dataset, and 69 of the 1,329 subscribers appear in no subscriber table. - On noise, 10 scans of 20 breakdowns run 200 tests and give 11 findings; since the threshold is 0.05, the expected number is 10.0 and at least one finding appears in 7 of the ten scans. - A single scan can produce five findings and there is no effect behind any of them; a finding whose number of tried breakdowns is not written cannot be told apart from a row of a noise sweep. #### Course Wrap-Up The account of the thirteen lessons stands in a single table, and every number was read from its own lesson. | Lesson | Observed number | Margin of variation | Does the decision change | |---|---|---|---| | Exploration 01 Purpose of Exploratory Analysis | east - north 4.463 m³ | not exceeded in 2,000 runs | no | | Exploration 02 Univariate Analysis | east mean 24.674 m³ | fence-based removal 25.154 m³ | yes, difference 4.622 | | Exploration 03 Bivariate Analysis | coefficient +0.600 | -0.057 when a single point is removed | yes, the relationship disappears | | Exploration 04 Multivariate Analysis | 5 of 15 pairs exceed the threshold | threshold 0.0354, 2,000 runs | yes, independent findings 1 | | Exploration 05 Visual Inspection of Data | dataset mean 19.693 m³ | 20.454 if zeros are not counted | yes, 3.9% range | | Application 01 Choosing Descriptive Summaries | mean 24.674 m³ | margin width 1.552 m³ | no, yes at 142 readings | | Application 02 Distribution Comparison | median 19.66 m³ | notch 0.920 ↔ 3.146 m³ | yes, two peaks are hidden | | Application 03 Interpreting Confidence Intervals | east mean 24.674 m³ | 23.894–25.431 m³ | yes, overlap misleads | | Application 04 Applying Hypothesis Testing | difference 0.957 m³ | p 0.01249 ↔ 0.00300 | no, yes in a quarter of cases | | Application 05 A/B Test Analysis | difference 0.123 m³ | [-0.698, 0.952] m³ | yes, 4 of 108 combinations | | Application 06 Modeling Relationships with Regression | coefficient -0.605 m³/year | [-0.714, -0.503] | yes, coefficient 1.825 | | Application 07 Cohort and Funnel Analysis | retention 0.619 | 0.867 when the rule changes | yes, the direction reverses | | Application 08 Misleading Analyses | 11 findings / 200 tests | expected 10.0 findings | yes, none of them is real | The third column is the rule of the course: **a difference whose margin of variation is not written next to it is not a finding but an observation; a number whose margin is not written is counted as not measured.** The fourth column says what the margin is good for. **In ten of the thirteen lessons the decision changes** — the number stays the same, but with a second defensible decision the sentence that enters the report becomes different. In three it does not change, and these three rows are the real output of the course: a number whose margin is written can also show that the decision **does not change**. The margin is not a warning sign but the basis of the decision. In thirteen lessons the margin was written next to every number. But all these numbers stayed on the analyst's screen as output blocks. The person who will make the decision saw neither the margin nor the number; making the finding **visible and transferable** was not dealt with at all. How the margin itself will stand in a display was not asked either: how an interval is drawn, how the overlap of two intervals is shown to a reader, where in a report the number of breakdowns tried is written. The next course, Visualization, Reporting and Analysis at Scale, takes up these questions.