Misleading Analyses

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.

# 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.

# 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.

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.