One day is not your diet: usual salt intake from two 24-hour recalls

#In plain words

The question. About nine in ten US adults reportedly eat more than the advised 2,300 mg of sodium a day, about a teaspoon of salt. But diets change from day to day: is that right?

What we did. A national survey asked 4,863 adults what they ate the day before, most of them twice. Each person's usual intake, their long-run average, is hidden; each day is that plus the day's ups and downs. Two days per person let the model separate how much people differ from how much one person's days differ, and set the daily part aside.

What we found. Judged by one day, 74 in 100 adults are over the line; judged by their usual diet, about 86 in 100 (85 to 88): many who look under it just had a light day. So the CDC's nine in ten is close; one day would have said three in four.

How sure are we. It is virtually certain (>99%) that more adults are usually over the line than a single day suggests. We have medium confidence in the 86 in 100 itself: two days per person is the minimum, recalls miss salt added at the table, and our ranges are probably a little too narrow.

The animation is saved as E108_usual_salt_intake.gif; part G builds it and shows six of its frames as a still for print.

import os

for v_ in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v_] = "1"  # nutpie runs chains as threads: a BLAS pool per chain would fight over the cores

import logging
import sys
import time
import warnings

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
from matplotlib.colors import LinearSegmentedColormap, to_rgb
from matplotlib.patches import Patch
from scipy import stats

from pymc_challenges import data
from pymc_challenges import public as pub

RANDOM_SEED = 108
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")
logging.getLogger("pymc").setLevel(logging.WARNING)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=UserWarning, module="arviz")
sys.modules.pop("jax", None)  # arviz leaves a lazy jax stub that makes matplotlib import all of JAX
pd.set_option("display.width", 160)
T_START, TIMES = time.time(), {}

HI, LO, INK, MUTED, LIGHT = pub.HI, pub.LO, pub.INK, pub.MUTED, pub.LIGHT
LINE = 2300.0                                   # mg sodium a day: the adult CDRR / Dietary Guidelines limit
STEM = "E108_usual_salt_intake"
print(f"PyMC {pm.__version__}, ArviZ {az.__version__}")

PyMC 6.3.2, ArviZ 1.3.0

#The setting

"Nine in ten Americans eat too much salt" is one of the most repeated numbers in public health. It comes from the National Health and Nutrition Examination Survey (NHANES), in which trained interviewers walk each participant through everything they ate and drank in the previous 24 hours, food by food, and a nutrient database turns the list into milligrams. The limit is 2,300 mg of sodium a day for adults, the intake above which the National Academies (2019) judged that cutting back lowers the risk of chronic disease.

A 24-hour recall measures one day, and one day is a poor guide to a diet. Pizza night and a day of salad are both real days of the same person. If we count everyone whose recall day was over 2,300 mg, we count the person whose usual intake is 2,000 mg but who had pizza, and miss the person whose usual intake is 2,600 mg but who ate lightly. Statisticians in nutrition have worried about this for a long time. Carriquiry (1999) showed why the share of people past a cut-off, computed from single days, is biased, and the remedy has a name: estimate the usual intake distribution. The Iowa State (ISU) method of Nusser, Carriquiry, Dodd & Fuller (1996) transforms intakes towards normality, splits their variance into a between-person and a within-person part using people with a second recall, and keeps only the first. The NCI method (Tooze et al. 2006; Kipnis et al. 2009) adds covariates and handles foods eaten only on some days with a two-part model; Dodd et al. (2006) review the family. NHANES collects a second recall on most participants for exactly this reason.

This is a measurement-error model, and in PyMC it is a hierarchical model. Person $i$ has a hidden usual level $\mu_i$; each recall day $j$ is a noisy reading,

$$y_{ij} = \mu_i + \text{(day effects)} + \varepsilon_{ij},\qquad \varepsilon_{ij}\sim N(0,\sigma_w^2),\qquad \mu_i \sim N(x_i^\top\beta,\ \sigma_b^2),$$

on a scale $y = g(\text{sodium})$ where the noise is close to normal. The between-person spread $\sigma_b$ is what we want; the within-person spread $\sigma_w$ is what one day adds on top, and only people with two days can tell the two apart. The share of the population over the line is then a property of the distribution of $\mu_i$, back on the milligram scale.

The textbook expectation is that one-day shares exaggerate: the noisy distribution is wider than the usual one, so it has more people at both ends. Look at where the line sits, though. If most people's usual intake is above the line, a wider distribution pushes more of them below it, and the one-day share comes out too low. Which way the bias goes depends on where the line cuts the distribution, and the CDC's own "nine in ten" is not a one-day number: Jackson, King, Zhao & Cogswell (2016, MMWR 64:1393-1397) estimated usual intakes from the two NHANES 2009-2012 recalls with the NCI method and found 89% of adults over 2,300 mg (98% of men, 80% of women). Our question is whether the 2017-2018 data, modelled from scratch, agree, and how far a single day would have taken us from that answer.

part	question	tool
A	What do two days of eating look like?	the recalls, survey weights, day 1 against day 2, the one-day share
B	Can a model take the day-to-day noise away?	a toy deconvolution with a known truth
C	A usual-intake model for sodium	a Box-Cox scale, prior predictive check, a funnel and its fix, a prior swap, predictive checks of one day and of two
D	What changes, and for whom?	the back-transform, three lines, a log-scale refit, sex and income
E	What does the second day buy?	person effects integrated out; refits with one day and with 300 second days; income-specific variances
F	Alcohol: zeros that mean "not today"	a correlated two-part model with a never-drinker group; a frequency tree
G	Explaining it to everyone	the HOP GIF, its still, flip rates, an icon-array ladder, a word strip, and the statistician's plot beside them

#A · Two days of eating

#The recalls

The derived file holds one row per adult and reliable recall day: total sodium, potassium, alcohol and energy for the day, the day of the week, age, sex, the family income-to-poverty ratio, race and ethnicity, the alcohol questionnaire, and the dietary sample weights.

data.describe("e108_nhanes_recalls")
raw = data.load("e108_nhanes_recalls")
print(f"{len(raw):,} recall days from {raw.seqn.nunique():,} adults; "
      f"{(raw.groupby('seqn').day.count() == 2).sum():,} with both days")
print("days with zero sodium:", int((raw.sodium_mg <= 0).sum()))

e108_nhanes_recalls: One row per adult and reliable recall day (9,103 rows: 4,863 adults, 4,240 with day 2): seqn, day (1 or 2), dow (1 = Sunday ... 7 = Saturday), sodium_mg, potassium_mg, alcohol_g, kcal, age (years, top-coded 80), sex, pir (family income-to-poverty ratio, top-coded 5, NaN for 572 adults), race_eth (RIDRETH3 code: 1 Mexican American, 2 other Hispanic, 3 NH White, 4 NH Black, 6 NH Asian, 7 other/multi), alq_ever (ALQ111 ever had a drink: 1 yes, 0 never, NaN unknown), alq_freq (ALQ121 past-year drinking frequency code 0-10: 0 none in the past year, also set for lifetime abstainers, 1 every day ... 10 once or twice a year, NaN refused/unknown), alq_drinks (ALQ130 usual drinks on a drinking day, NaN refused/unknown/not asked), wt_day1 / wt_day2 (dietary day-1 and two-day sample weights), psu, stratum (masked design variables).
  source: NHANES 2017-2018 (cycle J), CDC / National Center for Health Statistics: total nutrient intakes from the day-1 (in person) and day-2 (telephone) 24-hour dietary recalls (DR1TOT_J, DR2TOT_J) joined to demographics (DEMO_J) and the alcohol questionnaire (ALQ_J). US federal work, public domain. Adults aged 19+ with a reliable day-1 recall (DR1DRSTZ == 1), plus their day-2 recall when it is reliable too. DERIVED by tools/build_e108_nhanes.py; the url is the first of the four .xpt files, not this file - keep the cached file in data/ (md5 6eb9b657fc13deb94a137be75dad9536).
  url:    https://wwwn.cdc.gov/Nchs/Data/Nhanes/Public/2017/DataFiles/DR1TOT_J.xpt
9,103 recall days from 4,863 adults; 4,240 with both days
days with zero sodium: 1

One adult reports 0 mg of sodium for a whole (reliable) day; we drop that single day, which leaves that person with their second day only. We group people for the variance parts of the model into eight cells of age band by sex, and code income into bands of the income-to-poverty ratio (below the poverty line, 1-2, 2-3, 3-4, 4 or more times it, and unknown for the 572 adults without one). The weekend (Friday to Sunday) gets its own term: NHANES schedules weekend recalls on purpose, and eating differs on those days.

AGE_EDGES, AGE_LAB = [18, 30, 50, 70, 81], ["19-30", "31-50", "51-70", "71+"]
PIR_EDGES, PIR_LAB = [0, 1, 2, 3, 4, 5.01], ["<1", "1-2", "2-3", "3-4", "4+", "unknown"]
RACE = {1: "Mexican Am.", 2: "Other Hisp.", 3: "NH White", 4: "NH Black", 6: "NH Asian", 7: "Other/multi"}
RACES = list(RACE.values())
CELLS = [f"{a} {s}" for a in AGE_LAB for s in ("female", "male")]

d = raw[raw.sodium_mg > 0].copy()
per = d.sort_values("day").groupby("seqn").first().reset_index()
per = per[["seqn", "age", "sex", "pir", "race_eth", "alq_freq", "wt_day1"]]
per["age_g"] = pd.cut(per.age, AGE_EDGES, labels=AGE_LAB).astype(str)
pb = pd.cut(per.pir, PIR_EDGES, labels=PIR_LAB[:5], right=False).astype(object)
per["pir_g"] = pb.where(per.pir.notna(), "unknown").astype(str)
per["cell"] = per.age_g + " " + per.sex
per["race"] = per.race_eth.map(RACE)
d["p"] = pd.Index(per.seqn).get_indexer(d.seqn)
d["weekend"] = d.dow.isin([1, 6, 7]).astype(float)       # Sunday, Friday, Saturday
d = d.reset_index(drop=True)
CI = pd.Index(CELLS).get_indexer(per.cell)
PI = pd.Index(PIR_LAB).get_indexer(per.pir_g)
RI = pd.Index(RACES).get_indexer(per.race)
assert (d.p >= 0).all() and (CI >= 0).all() and (PI >= 0).all() and (RI >= 0).all()
W = per.wt_day1.to_numpy() / per.wt_day1.sum()           # day-1 dietary weights, one per adult
N = len(per)

d1 = d[d.day == 1]
x1, p1 = d1.sodium_mg.to_numpy(float), d1.p.to_numpy()
w1 = W[p1]
ONE_DAY = np.average(x1 > LINE, weights=w1)
print(f"{N:,} adults, {len(d):,} sodium days; weekend share of day-1 recalls {d1.weekend.mean():.0%}")
print(f"one-day share over {LINE:,.0f} mg: weighted {ONE_DAY:.3f}, unweighted {(x1 > LINE).mean():.3f}")
wkd1 = d1.weekend.to_numpy() == 1
w_dow = np.where(wkd1, (3 / 7) / wkd1.mean(), (4 / 7) / (1 - wkd1.mean()))   # weekday mix only
print(f"weekend share of day-1 recalls after weighting {np.average(wkd1, weights=w1):.0%}; "
      f"share over re-weighted for day of week only {np.average(x1 > LINE, weights=w_dow):.3f}")
for lab, m_ in (("weekdays", ~wkd1), ("weekends", wkd1)):
    print(f"  {lab}: unweighted {(x1[m_] > LINE).mean():.3f}, weighted {np.average(x1[m_] > LINE, weights=w1[m_]):.3f}")
print(pd.Series(PIR_LAB, name="pir band").to_frame().assign(
    adults=[int((PI == j).sum()) for j in range(6)],
    one_day_over=[round(np.average(x1[PI[p1] == j] > LINE, weights=w1[PI[p1] == j]), 3) for j in range(6)]
).to_string(index=False))


def wquant(x, w, q):
    """Weighted quantiles (q may be an array)."""
    o = np.argsort(x)
    c = np.cumsum(w[o]) / w[o].sum()
    return x[o][np.minimum(np.searchsorted(c, q), len(x) - 1)]


print("weighted one-day percentiles 5/50/95 (mg):", wquant(x1, w1, [0.05, 0.5, 0.95]).round())

4,863 adults, 9,102 sodium days; weekend share of day-1 recalls 63%
one-day share over 2,300 mg: weighted 0.740, unweighted 0.708
weekend share of day-1 recalls after weighting 42%; share over re-weighted for day of week only 0.702
  weekdays: unweighted 0.691, weighted 0.732
  weekends: unweighted 0.718, weighted 0.751
pir band  adults  one_day_over
      <1     810         0.668
     1-2    1196         0.695
     2-3     709         0.755
     3-4     468         0.779
      4+    1108         0.783
 unknown     572         0.696
weighted one-day percentiles 5/50/95 (mg): [1192. 3204. 7041.]

Weighted to the US adult population, 74% of single recall days are over 2,300 mg (2,300 mg of sodium is about a teaspoon of salt). The weights move the share up three points from the 71% unweighted, and not because of the weekend. NHANES schedules weekend recalls on purpose (63% of first recalls, 42% after weighting), and weekend days read higher, so undoing that excess alone would lower the share slightly, to 70%. The weights raise it on weekdays and weekends alike (69% to 73%, 72% to 75%), because they undo the survey's oversampling of groups that report less sodium. We meet this again when the model reproduces the weighted share. The weighted median day is about 3,200 mg, well above the line, and the 5th to 95th percentiles run from about 1,200 to 7,000 mg. Across income bands the one-day share climbs from 67% below the poverty line to 78% at four times it; we come back to that gap in part D.

#Day 1 against day 2

The people with two recalls show how far one day is from the next. On the log scale, two days of the same person should agree to the extent that people differ more than days do; their correlation is the share of the one-day variance that belongs to the person.

both = d.pivot_table(index="p", columns="day", values="sodium_mg").dropna()
l1, l2 = np.log(both[1].to_numpy()), np.log(both[2].to_numpy())
r = np.corrcoef(l1, l2)[0, 1]
within = np.var(l2 - l1 - (l2 - l1).mean()) / 2
between = np.cov(l1, l2)[0, 1]
print(f"{len(both):,} adults with two days: log-scale correlation {r:.2f}")
print(f"method of moments on the log scale: within {within:.3f}, between {between:.3f}, "
      f"ratio {within / between:.1f}; mean log day 2 minus day 1 {np.mean(l2 - l1):+.3f}")
o1, o2 = both[1].to_numpy() > LINE, both[2].to_numpy() > LINE
print(f"over on both days {np.mean(o1 & o2):.0%}, on one day only {np.mean(o1 != o2):.0%}, "
      f"on neither {np.mean(~o1 & ~o2):.0%}  (unweighted)")

fig, axs = plt.subplots(1, 2, figsize=(12, 5.2), width_ratios=[1.1, 1])
ax = axs[0]
ax.scatter(both[1], both[2], s=5, alpha=0.25, color=INK, lw=0)
for f in (ax.axvline, ax.axhline):
    f(LINE, color=LO, lw=1.2)
ax.set_xscale("log"); ax.set_yscale("log")
ax.set_xlim(300, 30000); ax.set_ylim(300, 30000)
ax.set_xlabel("sodium on day 1 (mg, log scale)"); ax.set_ylabel("sodium on day 2 (mg, log scale)")
ax.set_title(f"Two days of one person agree only loosely\n(log correlation {r:.2f})", loc="left")
for (x, y, t) in [(20000, 20000, f"over both days\n{np.mean(o1 & o2):.0%}"),
                  (500, 500, f"under both days\n{np.mean(~o1 & ~o2):.0%}"),
                  (20000, 500, f"over on day 1 only\n{np.mean(o1 & ~o2):.0%}"),
                  (500, 20000, f"over on day 2 only\n{np.mean(~o1 & o2):.0%}")]:
    ax.text(x, y, t, ha="center", va="center", fontsize=10.5, color=INK,
            bbox=dict(fc="white", ec=LIGHT, alpha=0.9))
ax = axs[1]
DOW = ["Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"]
for day, col, off in ((1, INK, -0.18), (2, MUTED, 0.18)):
    g = d[d.day == day].groupby("dow").sodium_mg
    med, n = g.median(), g.size()
    ax.bar(np.arange(7) + off, med.reindex(range(1, 8)).to_numpy(), width=0.36, color=col,
           label=f"day {day} ({'in person' if day == 1 else 'by phone'})")
ax.axhline(LINE, color=LO, lw=1.2)
ax.text(6.6, LINE + 60, "2,300 mg", color=LO, ha="right", va="bottom", fontsize=10,
        bbox=dict(fc="white", ec="none", pad=1))
ax.set_xticks(range(7), DOW)
ax.set_ylabel("median sodium (mg)")
ax.set_ylim(0, 4300)
ax.legend(loc="upper left", fontsize=10, frameon=False, ncols=2)
ax.set_title("The phone day reads lower;\nweekends read a little higher", loc="left");

4,239 adults with two days: log-scale correlation 0.39
method of moments on the log scale: within 0.203, between 0.128, ratio 1.6; mean log day 2 minus day 1 -0.097
over on both days 53%, on one day only 31%, on neither 17%  (unweighted)

The cloud is wide: the log correlation between two days of the same person is 0.39. In the language of the model, about 40% of the one-day variance on this scale is the person and 60% is the day (method-of-moments within / between ratio 1.6). Almost a third of adults (31%) are over the line on one of their two days and under it on the other. The second recall, taken by telephone, also reads about 9% lower on average, a known feature of NHANES (people report less on the phone, or the days differ); the model gets a day-2 term rather than treating it as noise.

The counts in the corners already hint at the answer. 53% of people are over on both days and only 17% under on both: the bulk sits above the line, and most of the people who cross it are people from above it having a light day.

#B · A toy deconvolution

Before the real model, a world where we know the truth. We simulate 500 people whose usual sodium intakes have a median of 3,300 mg, like real adults, and give each two noisy days with twice as much within-person as between-person variance on the log scale. Then we fit the model

$$\log x_{ij} \sim N(\mu_i, \sigma_w^2), \qquad \mu_i = a + \sigma_b u_i,\quad u_i\sim N(0,1),$$

and compute each person's usual intake on the milligram scale. That last step is easy to get wrong. Usual intake is the long-run mean of a person's days, and the mean of a log-normal day is $E[x \mid \mu_i] = \exp(\mu_i + \sigma_w^2/2)$, not $\exp(\mu_i)$, which is the person's median day.

t0 = time.time()
rs = np.random.default_rng(1)
n_toy, sb_true, sw_true = 500, 0.30, 0.42
mu_true = np.log(3300) - sw_true**2 / 2 + sb_true * rs.standard_normal(n_toy)   # log median day
usual_true = np.exp(mu_true + sw_true**2 / 2)                                     # mean day, mg
toy = np.exp(mu_true[:, None] + sw_true * rs.standard_normal((n_toy, 2)))
with pm.Model(coords={"person": np.arange(n_toy), "day": [1, 2]}) as toy_model:
    a = pm.Normal("a", np.log(3000), 1)
    s_b = pm.HalfNormal("s_b", 1)
    s_w = pm.HalfNormal("s_w", 1)
    u = pm.Normal("u", 0, 1, dims="person")
    pm.Normal("y", (a + s_b * u)[:, None], s_w, observed=np.log(toy), dims=("person", "day"))
    toy_idata = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
post_toy = az.extract(toy_idata, var_names=["a", "s_b", "s_w", "u"], num_samples=1000, random_seed=1)
mu_toy = (post_toy["a"].values + post_toy["s_b"].values * post_toy["u"].values).T      # (S, person)
usual_toy = np.exp(mu_toy + post_toy["s_w"].values[:, None] ** 2 / 2)
print(f"divergences {int(toy_idata.sample_stats['diverging'].sum())}, "
      f"max r_hat {max(float(az.rhat(toy_idata, var_names=[v]).to_dataset()[v].max()) for v in ['a', 's_b', 's_w', 'u']):.3f}")
print(az.summary(toy_idata, var_names=["s_b", "s_w"], round_to=3)[["mean", "sd", "eti89_lb", "eti89_ub", "r_hat"]])
print(f"truth: s_b {sb_true}, s_w {sw_true}; realised in this sample: sd of usual levels "
      f"{mu_true.std():.3f}, sd of day noise {np.std(np.log(toy) - mu_true[:, None]):.3f}")
for lab, f in [("over 2,300 mg", lambda v: v > LINE), ("under 1,500 mg", lambda v: v < 1500),
               ("over 5,000 mg", lambda v: v > 5000)]:
    est = f(usual_toy).mean(1)
    print(f"{lab:15s} true usual {f(usual_true).mean():.3f}   day 1 {f(toy[:, 0]).mean():.3f}   "
          f"model {np.median(est):.3f} ({np.quantile(est, 0.055):.3f}-{np.quantile(est, 0.945):.3f})")
TIMES["toy"] = time.time() - t0

divergences 0, max r_hat 1.008
      mean     sd  eti89_lb  eti89_ub  r_hat
s_b  0.281  0.024     0.243     0.319  1.003
s_w  0.449  0.014     0.427     0.471  1.001
truth: s_b 0.3, s_w 0.42; realised in this sample: sd of usual levels 0.274, sd of day noise 0.442
over 2,300 mg   true usual 0.908   day 1 0.664   model 0.886 (0.850-0.926)
under 1,500 mg  true usual 0.004   day 1 0.080   model 0.002 (0.000-0.010)
over 5,000 mg   true usual 0.058   day 1 0.136   model 0.060 (0.034-0.084)

The model recovers the split. Its between-person sd covers the truth; its within-person sd, 0.45, is above the nominal 0.42 but matches the 0.44 of noise this sample actually drew, which is what two days can see. With them come the shares. In this toy world 91% of people usually eat over 2,300 mg; a single day says 66%, and the model, from the same two days, says about 89% (85-93%). The tails go the other way: one day finds 8% of people under 1,500 mg where the truth is under 1%, and 14% over 5,000 mg where the truth is 6%; the model gets both. Same data, same noise, opposite biases, all depending on which side of the middle the line falls. The picture makes it plain.

fig, ax = plt.subplots(figsize=(11, 4.6))
bins = np.geomspace(500, 20000, 60)
ax.hist(toy[:, 0], bins=bins, color=LIGHT, label="one day (simulated recall)")
ax.hist(usual_true, bins=bins, histtype="step", lw=2, color=INK, label="true usual intake")
for k in range(25):
    h, _ = np.histogram(usual_toy[k], bins=bins)
    ax.stairs(h, bins, color=HI, alpha=0.12, lw=1.2)
ax.plot([], [], color=HI, lw=1.5, label="model's usual intake (25 posterior draws)")
for v, t in ((1500, "1,500"), (LINE, "2,300"), (5000, "5,000")):
    ax.axvline(v, color=LO, lw=1.1, ls="-" if v == LINE else ":")
    ax.text(v * 1.03, ax.get_ylim()[1] * 0.92, t + " mg", color=LO, fontsize=10)
ax.set_xscale("log")
ax.set_xlabel("sodium (mg a day, log scale)"); ax.set_ylabel("people")
ax.legend(loc="upper left", fontsize=10, frameon=False)
ax.set_title("Removing day-to-day noise squeezes the distribution towards its middle, which is over the line",
             loc="left");

#C · A usual-intake model for sodium

#A scale where the noise is normal

The log is the usual first choice for intakes, but on the log scale real sodium days are skewed to the left: a few very light days stretch the lower tail. Nusser et al. (1996) transform intakes with a power before modelling, and so do we. A Box-Cox transform $g_\lambda(x) = (x^\lambda - 1)/\lambda$ includes the log at $\lambda = 0$; we take $\lambda$ from its maximum-likelihood value on day 1, then standardise.

x_all = d.sodium_mg.to_numpy(float)
LAM = float(stats.boxcox_normmax(x1))
zbc = (x_all**LAM - 1) / LAM
M0, S0 = zbc[d.day.to_numpy() == 1].mean(), zbc[d.day.to_numpy() == 1].std()
Y = (zbc - M0) / S0
for lab, v in (("log", np.log(x1)), (f"Box-Cox {LAM:.2f}", (x1**LAM - 1) / LAM)):
    print(f"{lab:13s} skewness {stats.skew(v):+.2f}, excess kurtosis {stats.kurtosis(v):.2f}")


def ginv(zs):
    """Standardised Box-Cox scale -> mg."""
    return np.maximum(LAM * (zs * S0 + M0) + 1, 1e-12) ** (1 / LAM)

log           skewness -1.10, excess kurtosis 6.42
Box-Cox 0.30  skewness +0.09, excess kurtosis 1.14

The log leaves a skewness of -1.1 and an excess kurtosis of 6.4; $\lambda \approx 0.30$ (between a cube root and a square root) removes the skew (+0.09) and most of the excess kurtosis (1.1 is left, so the tails are still a little heavier than a normal's; we meet them again in the predictive checks and in the first exercise). On this standardised scale one-day values have mean 0 and sd 1, which makes priors easy to state.

#The model

For person $i$ in variance cell $c(i)$ (age band by sex), recall day $j$:

$$y_{ij} \sim N\big(\mu_i + \gamma_2[j=2] + \gamma_w(\text{weekend}_{ij} - 3/7),\ \sigma_{w,c(i)}^2\big),$$ $$\mu_i = a + \beta_{c(i)} + \beta_{\text{income}(i)} + \beta_{\text{race}(i)} + \sigma_{b,c(i)}\,u_i, \qquad u_i \sim N(0, 1).$$

The group effects are zero-sum normals. The day-2 term absorbs the phone-day drop, and the weekend term is centred at 3/7, so $\mu_i$ is the level of an average day of a week. Each cell has its own within-person and between-person sd, and both sets are partially pooled on the log scale: $\log\sigma_{w,c}\sim N(\ell_w, \tau_w^2)$ and $\log\sigma_{b,c}\sim N(\ell_b, \tau_b^2)$. The priors on $\ell_w$ and $\ell_b$ centre the within / between variance ratio at $(0.75/0.6)^2 \approx 1.6$. That centre is not independent of the data: it is the method-of-moments ratio we computed on the log scale in part A. The prior is wide enough (the prior predictive check below gives the ratio a range from about a sixth to 18) that the data should decide, and after the fit we refit with the prior centred at a ratio of 1 to check that they do.

A usual intake in milligrams is then the mean over day noise and over a typical week (four weekdays, three weekend days):

$$T_i = \tfrac47\,E\big[g^{-1}(\mu_i - \tfrac37\gamma_w + \sigma_w\varepsilon)\big] + \tfrac37\,E\big[g^{-1}(\mu_i + \tfrac47\gamma_w + \sigma_w\varepsilon)\big],\qquad \varepsilon\sim N(0,1),$$

computed with 24-node Gauss-Hermite quadrature. The in-person day 1 is the reference: we do not add the phone-day drop to anyone's usual intake.

P_ARR, DAY2 = d.p.to_numpy(), (d.day.to_numpy() == 2).astype(float)
WK = d.weekend.to_numpy() - 3 / 7
COORDS = {"person": per.seqn.to_numpy(), "cell": CELLS, "pir": PIR_LAB, "race": RACES}


def sodium_model(centred_between=False, y=Y, sw_centre=0.75):
    """The sodium model on scale y; sw_centre/0.6 squared is the prior centre of the variance ratio."""
    with pm.Model(coords=COORDS) as m:
        a = pm.Normal("a", 0, 1)
        b_cell = pm.ZeroSumNormal("b_cell", 0.5, dims="cell")
        b_pir = pm.ZeroSumNormal("b_pir", 0.5, dims="pir")
        b_race = pm.ZeroSumNormal("b_race", 0.5, dims="race")
        g_day2 = pm.Normal("g_day2", 0, 0.5)
        g_wk = pm.Normal("g_wk", 0, 0.5)
        ls_w0 = pm.Normal("ls_w0", np.log(sw_centre), 0.5)
        ls_b0 = pm.Normal("ls_b0", np.log(0.6), 0.5)
        t_w = pm.HalfNormal("t_w", 0.3)
        t_b = pm.HalfNormal("t_b", 0.3)
        ls_w = pm.Normal("ls_w", ls_w0, t_w, dims="cell")                     # centred
        if centred_between:
            ls_b = pm.Normal("ls_b", ls_b0, t_b, dims="cell")
        else:                                                                # non-centred
            ls_b = pm.Deterministic("ls_b", ls_b0 + t_b * pm.Normal("e_b", 0, 1, dims="cell"), dims="cell")
        s_w = pm.Deterministic("s_w", pt.exp(ls_w), dims="cell")
        s_b = pm.Deterministic("s_b", pt.exp(ls_b), dims="cell")
        u = pm.Normal("u", 0, 1, dims="person")
        mu = a + b_cell[CI] + b_pir[PI] + b_race[RI] + s_b[CI] * u
        pm.Normal("y", mu[P_ARR] + g_day2 * DAY2 + g_wk * WK, s_w[CI][P_ARR], observed=y)
    return m

#Prior predictive check

What do these priors say about one-day intakes, before seeing any? We draw 300 sets of parameters and recall days, back-transform to milligrams, and look at the median day and the share over the line in each.

t0 = time.time()
m_sod = sodium_model()
with m_sod:
    prior = pm.sample_prior_predictive(300, random_seed=RANDOM_SEED, var_names=["y", "s_w", "s_b"])
yp = prior.prior_predictive["y"].values.reshape(300, -1)[:, d.day.to_numpy() == 1]
xp = ginv(yp)
sw_p, sb_p = prior.prior["s_w"].values.reshape(300, -1), prior.prior["s_b"].values.reshape(300, -1)
ratio_p = (sw_p**2).mean(1) / (sb_p**2).mean(1)
q = lambda v: np.round(np.quantile(v, [0.055, 0.5, 0.945]), 2)
print("prior median day (mg):      ", q(np.median(xp, 1)).round(-1))
print("prior share over 2,300 mg:  ", q((xp > LINE).mean(1)))
print("prior 99th percentile (mg): ", q(np.quantile(xp, 0.99, axis=1)).round(-2))
print("prior within/between ratio: ", q(ratio_p))
print(f"share of prior days at the Box-Cox floor (x = 0): {(LAM * (yp * S0 + M0) + 1 <= 1e-12).mean():.4f}")

fig, axs = plt.subplots(1, 2, figsize=(11, 3.8))
axs[0].hist((xp > LINE).mean(1), bins=30, color=LIGHT, edgecolor="white")
axs[0].axvline(ONE_DAY, color=LO, lw=1.5)
axs[0].text(ONE_DAY, axs[0].get_ylim()[1] * 0.9, " observed\n one-day share", color=LO, fontsize=10)
axs[0].set_xlabel("share of simulated days over 2,300 mg"); axs[0].set_yticks([])
axs[0].set_title("Prior: share of days over the line", loc="left")
axs[1].hist(np.log10(ratio_p), bins=30, color=LIGHT, edgecolor="white")
axs[1].set_xticks(np.log10([1 / 16, 1 / 4, 1, 4, 16, 64]), ["1/16", "1/4", "1", "4", "16", "64"])
axs[1].set_xlabel("within / between variance ratio (prior)"); axs[1].set_yticks([])
axs[1].set_title("Prior: day noise against person spread", loc="left");
TIMES["prior"] = time.time() - t0

prior median day (mg):       [1130. 3140. 6810.]
prior share over 2,300 mg:   [0.22 0.65 0.94]
prior 99th percentile (mg):  [ 6000. 14300. 32100.]
prior within/between ratio:  [ 0.17  1.86 17.8 ]
share of prior days at the Box-Cox floor (x = 0): 0.0020

The priors are wide where we know little and sensible where we know a lot: the median day stays within the range of real diets (1,100-6,800 mg), the share over the line runs from about a fifth to nearly everyone around the observed one, and the variance ratio spans about a sixth to 18. Only 0.2% of simulated days hit the transform's floor at zero milligrams.

#Failure: a funnel in the between-person scales

A first version wrote both scale hierarchies the plain way, $\log\sigma_{b,c} \sim N(\ell_b,\tau_b)$ with $\log\sigma_{b,c}$ sampled directly. Each cell has hundreds of people, so this "centred" form usually samples well (E25 and E94 met the same pattern). Here it does not.

t0 = time.time()
with sodium_model(centred_between=True):
    idata_c = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
div_c = idata_c.sample_stats["diverging"].values
tb_c = idata_c.posterior["t_b"].values
print(f"centred between-person scales: {int(div_c.sum())} divergences (nutpie)")
print(f"t_b at divergent draws: median {np.median(tb_c[div_c]) if div_c.any() else float('nan'):.3f}; "
      f"all draws: median {np.median(tb_c):.3f}, 5% quantile {np.quantile(tb_c, 0.05):.3f}")
print(az.summary(idata_c, var_names=["t_b", "t_w", "ls_b0"], round_to=3)[["mean", "sd", "ess_bulk", "r_hat"]])
TIMES["fit_centred"] = time.time() - t0
del idata_c

centred between-person scales: 22 divergences (nutpie)
t_b at divergent draws: median 0.010; all draws: median 0.070, 5% quantile 0.018
        mean     sd  ess_bulk  r_hat
t_b    0.077  0.046   513.560  1.010
t_w    0.156  0.051  3559.640  1.000
ls_b0 -0.612  0.041   857.258  1.004

The divergences sit at small $\tau_b$. That is the funnel of Neal (2003) in its classic form: the eight cells barely differ in their between-person spread, so $\tau_b$ is pushed towards zero, where the eight $\log\sigma_{b,c}$ must squeeze together and the sampler's step size no longer fits. The within-person scales are different. Each is informed directly by the gaps between two days of hundreds of people, and $\tau_w$ is clearly above zero (young men are much noisier than older women), so they are well identified and centred works for them. Making both hierarchies non-centred fails the other way: in a trial run not repeated here (the same seed and nutpie's default 4 x 1,000 draws; swap ls_w for ls_w0 + t_w * e_w to reproduce it), the within-person level $\ell_w$ mixed so slowly that its $\hat R$ reached 1.9. The fix is a mixed parameterisation, centred where the data are strong and non-centred where the group-level sd is near zero, which is the model above.

#The fit

t0 = time.time()
with m_sod:
    idata = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
TIMES["fit_sodium"] = time.time() - t0
print(f"sampling took {TIMES['fit_sodium']:.0f} s with {idata.posterior.attrs['tuning_steps']} tuning steps; "
      f"divergences {int(idata.sample_stats['diverging'].sum())}")
print("tree depth counts:", np.bincount(idata.sample_stats["depth"].values.ravel()))
SCAL = ["a", "g_day2", "g_wk", "ls_w0", "ls_b0", "t_w", "t_b"]
print(az.summary(idata, var_names=SCAL, round_to=3)[["mean", "sd", "eti89_lb", "eti89_ub", "ess_bulk", "r_hat"]])
for v in ["s_w", "s_b", "u", "e_b", "b_pir", "b_cell", "b_race"]:
    print(f"{v:7s} max r_hat {float(az.rhat(idata, var_names=[v]).to_dataset()[v].max()):.3f}, "
          f"min bulk ESS {float(az.ess(idata, var_names=[v]).to_dataset()[v].min()):.0f}")
swb = idata.posterior["s_w"].values.reshape(-1, 8) ** 2, idata.posterior["s_b"].values.reshape(-1, 8) ** 2
ratio = (swb[0][:, CI] @ W) / (swb[1][:, CI] @ W)
print(f"within / between variance ratio (weighted over adults): {q(ratio)}")
print(pd.DataFrame({"s_w": idata.posterior["s_w"].mean(("chain", "draw")).values,
                    "s_b": idata.posterior["s_b"].mean(("chain", "draw")).values}, index=CELLS).round(2).T)

sampling took 10 s with 400 tuning steps; divergences 0
tree depth counts: [   0    0    0    0   29 3971]
         mean     sd  eti89_lb  eti89_ub  ess_bulk  r_hat
a      -0.039  0.017    -0.066    -0.011  2580.992  1.002
g_day2 -0.148  0.017    -0.176    -0.120  4940.756  1.001
g_wk    0.047  0.018     0.019     0.076  5575.766  1.002
ls_w0  -0.289  0.060    -0.387    -0.195  3720.716  1.001
ls_b0  -0.612  0.040    -0.676    -0.549   889.689  1.007
t_w     0.157  0.051     0.096     0.252  3890.099  1.000
t_b     0.078  0.046     0.014     0.161   588.815  1.000
s_w     max r_hat 1.002, min bulk ESS 1426
s_b     max r_hat 1.004, min bulk ESS 797

u       max r_hat 1.010, min bulk ESS 5517
e_b     max r_hat 1.002, min bulk ESS 1528
b_pir   max r_hat 1.003, min bulk ESS 3680
b_cell  max r_hat 1.002, min bulk ESS 3414
b_race  max r_hat 1.003, min bulk ESS 2520
within / between variance ratio (weighted over adults): [1.82 2.   2.22]
     19-30 female  19-30 male  31-50 female  31-50 male  51-70 female  51-70 male  71+ female  71+ male
s_w          0.81        0.93          0.77        0.82          0.67        0.77        0.63      0.66
s_b          0.50        0.56          0.54        0.57          0.51        0.57        0.53      0.56

No divergences under nutpie, $\hat R \le 1.01$ for every parameter printed above and thousands of effective draws for the person effects. The phone day reads lower ($\gamma_2 \approx -0.15$ sd) and weekend days a little higher ($\gamma_w \approx +0.05$). The within-person sds range from about 0.63 to 0.93 across cells, with young men the least predictable from day to day; the between-person sds sit at 0.50-0.57 in every cell. On this scale the day adds about twice as much variance as the person (on the log scale of part A the raw ratio was 1.6).

#Does the prior choose the split?

The prior centre of the variance ratio came from the data, so we refit with the within-person prior moved down to centre the ratio at 1, an even split, and compare the posterior ratio.

t0 = time.time()
with sodium_model(sw_centre=0.6):
    idata_r1 = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
sw1, sb1 = idata_r1.posterior["s_w"].values.reshape(-1, 8) ** 2, idata_r1.posterior["s_b"].values.reshape(-1, 8) ** 2
print(f"prior centred at ratio 1: divergences {int(idata_r1.sample_stats['diverging'].sum())}; "
      f"posterior ratio {q((sw1[:, CI] @ W) / (sb1[:, CI] @ W))} against {q(ratio)} with the main prior")
del idata_r1, sw1, sb1
TIMES["fit_prior1"] = time.time() - t0

prior centred at ratio 1: divergences 0; posterior ratio [1.81 2.   2.21] against [1.82 2.   2.22] with the main prior

The posterior ratio does not move (2.0, 89% interval 1.8-2.2, under either prior), so the data, not the prior's centre, set the split.

#Posterior predictive checks

We replicate every day-1 recall from the posterior, with each person's own $\mu_i$, a fresh day effect and their recorded weekday, and compare quantiles with the data. A usual-intake model lives or dies by its tails, so we look at the 1st and 99th percentiles and at the single largest day.

t0 = time.time()
th = idata.isel(draw=slice(None, None, 4))                       # 1,000 draws for derived quantities
X = lambda n: az.extract(th, var_names=[n]).values
A_, BC_, BP_, BR_ = X("a"), X("b_cell"), X("b_pir"), X("b_race")
SW_, SB_, GWK_, G2_, U_ = X("s_w").T, X("s_b").T, X("g_wk"), X("g_day2"), X("u").T   # (S, cell), (S,), (S, person)
MFIX = (A_ + BC_[CI] + BP_[PI] + BR_[RI]).T                                    # (S, person) fixed part
MU = MFIX + SB_[:, CI] * U_
S = MU.shape[0]
r1 = np.random.default_rng(2)
k_rep = r1.choice(S, 200, replace=False)
wk1 = d1.weekend.to_numpy() - 3 / 7
rep = ginv(MU[k_rep][:, p1] + GWK_[k_rep, None] * wk1 + SW_[k_rep][:, CI[p1]] * r1.standard_normal((200, len(p1))))
rows = []
for qq in [0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99]:
    lo_, hi_ = np.quantile(np.quantile(rep, qq, axis=1), [0.05, 0.95])
    rows.append((f"{qq:.0%}", np.quantile(x1, qq), lo_, hi_))
ppc = pd.DataFrame(rows, columns=["quantile", "observed", "rep 5%", "rep 95%"]).round(0)
print(ppc.to_string(index=False))
print(f"largest day: observed {x1.max():,.0f} mg; replicated maxima 5-95%: "
      f"{np.quantile(rep.max(1), 0.05):,.0f}-{np.quantile(rep.max(1), 0.95):,.0f} mg")
print(f"unweighted share over the line: observed {(x1 > LINE).mean():.3f}, replicated "
      f"{np.quantile((rep > LINE).mean(1), 0.05):.3f}-{np.quantile((rep > LINE).mean(1), 0.95):.3f}")

fig, ax = plt.subplots(figsize=(11, 4.2))
bins = np.geomspace(150, 30000, 70)
for k in range(30):
    h, _ = np.histogram(rep[k], bins=bins)
    ax.stairs(h, bins, color=HI, alpha=0.12)
ax.stairs(np.histogram(x1, bins=bins)[0], bins, color=INK, lw=2, label="observed day 1")
ax.plot([], [], color=HI, label="30 replicated day-1 data sets")
ax.axvline(LINE, color=LO, lw=1.2)
ax.annotate(f"largest real day\n{x1.max():,.0f} mg", (x1.max(), 1), xytext=(x1.max() * 0.45, 150),
            fontsize=10, color=INK, arrowprops=dict(arrowstyle="->", color=INK))
ax.set_xscale("log"); ax.set_xlabel("sodium on day 1 (mg, log scale)"); ax.set_ylabel("adults")
ax.legend(frameon=False, fontsize=10)
ax.set_title("Replicated recall days match the real ones, except for one extreme day", loc="left");
TIMES["ppc"] = time.time() - t0

quantile  observed  rep 5%  rep 95%
      1%     637.0   628.0    735.0
      5%    1176.0  1092.0   1175.0
     25%    2135.0  2063.0   2150.0
     50%    3092.0  3032.0   3126.0
     75%    4309.0  4293.0   4426.0
     95%    6928.0  6768.0   7076.0
     99%    9633.0  9016.0   9727.0
largest day: observed 25,949 mg; replicated maxima 5-95%: 14,005-19,248 mg
unweighted share over the line: observed 0.708, replicated 0.688-0.711

The quantiles from the 1st to the 99th fall inside or at the edge of the replicated ranges (the 5th percentile, 1,176 mg, sits just above its replicated range, so the model puts slightly too many people on very light days), and so does the share over the line. One real day does not: someone reported 25,949 mg, beyond the largest of every replicated data set. One extreme day among 4,862 hardly moves a usual-intake distribution, but it is a reminder that the normal day error is a choice; a Student-t day error is the first exercise.

This check has a blind spot, and part E will show why: day-1 values alone say nothing about how their variance splits into person and day, which is exactly what moves the headline. A check of the split has to look at the two days together. For the 4,239 people with two recalls we simulate fresh people with the same covariates and weekdays, two days each, and compare the correlation between the days, the spread of the day-2 minus day-1 difference, and whether that difference is larger for people who eat more (the model assumes a constant day-to-day sd in each cell on this scale; changing within-person spread is what the ISU method's second transformation handles).

t0 = time.time()
cnt2 = d.groupby("p").day.transform("count").to_numpy() == 2
ra, rb = d[cnt2 & (d.day == 1)].sort_values("p"), d[cnt2 & (d.day == 2)].sort_values("p")
pp2 = ra.p.to_numpy()
assert (pp2 == rb.p.to_numpy()).all()
YA, YB = Y[ra.index.to_numpy()], Y[rb.index.to_numpy()]
wka, wkb = ra.weekend.to_numpy() - 3 / 7, rb.weekend.to_numpy() - 3 / 7


def pair_stats(a, b):
    """Correlation of the two days, quantiles of b - a, and mean |b - a| by quintile of (a + b) / 2."""
    df, mid = b - a, (a + b) / 2
    bin_ = np.searchsorted(np.quantile(mid, [0.2, 0.4, 0.6, 0.8]), mid)
    return [np.corrcoef(a, b)[0, 1], *np.quantile(df, [0.05, 0.25, 0.75, 0.95]),
            *[np.abs(df[bin_ == j]).mean() for j in range(5)]]


r3 = np.random.default_rng(12)
reps2 = []
for k in r3.choice(S, 200, replace=False):
    mu_f = MFIX[k, pp2] + SB_[k, CI[pp2]] * r3.standard_normal(len(pp2))       # a fresh person
    sw_f = SW_[k, CI[pp2]]
    reps2.append(pair_stats(mu_f + GWK_[k] * wka + sw_f * r3.standard_normal(len(pp2)),
                            mu_f + G2_[k] + GWK_[k] * wkb + sw_f * r3.standard_normal(len(pp2))))
reps2 = np.array(reps2)
lab2 = (["day 1 - day 2 correlation"] + [f"difference, {v}th pct" for v in (5, 25, 75, 95)]
        + [f"mean |difference|, intake quintile {j + 1}" for j in range(5)])
print(pd.DataFrame({"observed": pair_stats(YA, YB), "rep 5%": np.quantile(reps2, 0.05, 0),
                    "rep median": np.median(reps2, 0), "rep 95%": np.quantile(reps2, 0.95, 0)},
                   index=lab2).round(3).to_string())
y_line = ((LINE**LAM - 1) / LAM - M0) / S0
print(f"2,300 mg sits at the {np.mean((YA + YB) / 2 < y_line):.0%} point of the two-day averages "
      f"(intake quintile 1 is below {np.quantile((YA + YB) / 2, 0.2):.2f}, the line at {y_line:.2f})")
TIMES["ppc_pairs"] = time.time() - t0

                                      observed  rep 5%  rep median  rep 95%
day 1 - day 2 correlation                0.400   0.375       0.404    0.432
difference, 5th pct                     -1.939  -1.997      -1.925   -1.863
difference, 25th pct                    -0.821  -0.929      -0.879   -0.832
difference, 75th pct                     0.491   0.498       0.550    0.598
difference, 95th pct                     1.514   1.510       1.592    1.670
mean |difference|, intake quintile 1     0.840   0.789       0.833    0.866
mean |difference|, intake quintile 2     0.782   0.808       0.844    0.884
mean |difference|, intake quintile 3     0.772   0.822       0.856    0.892
mean |difference|, intake quintile 4     0.816   0.833       0.871    0.910
mean |difference|, intake quintile 5     0.927   0.865       0.908    0.947
2,300 mg sits at the 29% point of the two-day averages (intake quintile 1 is below -0.76, the line at -0.53)

The part that matters most passes. The correlation between a person's two days, which is the between-person share of the variance, is reproduced, and so are the 5th and 95th percentiles of the day-to-day difference. Two things are not. The central half of the real differences is narrower than the model's (the observed 25th and 75th percentiles, -0.82 and +0.49, fall just short of the replicated ranges), with the same outer tails: real day-to-day changes are more often small and occasionally large, the heavier tails the Box-Cox kurtosis already showed. And the real differences are smaller in the middle three intake quintiles (0.77-0.82) than at the ends (0.84 and 0.93), and in those middle three they fall below the model's whole replicated range; the model's own curve rises only gently, from about 0.83 to 0.91. Both point to the Student-t day error of the first exercise, and to a day-to-day spread that varies with intake. The second matters for the headline more than it first seems, because the line is not in the middle of the distribution: 2,300 mg sits near the 14th percentile of usual intake, the 26th of single days and the 29th of the two-day averages, inside the second quintile, where real day-to-day changes are smaller than the model assumes (0.78 against 0.84). A spread that changes with intake shifts a deconvolved share near the line by an amount, and in a direction, that this model cannot tell us; it joins the caveats at the end.

#D · What changes, and for whom?

#The share over the line

The population share over the line is a property of usual intakes across the US adult population, not of the 4,863 people who happened to be sampled. For each posterior draw we give every sampled adult a fresh person effect $u \sim N(0,1)$ around their covariates, compute that new person's usual intake, and average with the survey weights. This is a weighted post-stratification over the sample's covariate mix. For comparison we also use the sampled people's own posterior usual intakes, and the shortcut that back-transforms $\mu_i$ directly, which gives a person's median day instead of their mean.

GH_X, GH_W = np.polynomial.hermite_e.hermegauss(24)
GH_W = GH_W / GH_W.sum()


def usual(mu, sw, gwk):
    """Usual intake (mg): mean over day noise (Gauss-Hermite) and a typical week; arrays (S, person)."""
    out = 0.0
    for wkv, wt in ((-3 / 7, 4 / 7), (4 / 7, 3 / 7)):
        out = out + wt * (ginv((mu + gwk * wkv)[..., None] + sw[..., None] * GH_X) * GH_W).sum(-1)
    return out


def chunked(fn, *arrs, step=100):
    return np.concatenate([fn(*(a[k:k + step] for a in arrs)) for k in range(0, arrs[0].shape[0], step)])


r2 = np.random.default_rng(3)
MU_NEW = MFIX + SB_[:, CI] * r2.standard_normal(MFIX.shape)          # fresh people, same covariates
SWP = SW_[:, CI]
T_SAMPLED = chunked(usual, MU, SWP, np.repeat(GWK_[:, None], N, 1))
T_NEW = chunked(usual, MU_NEW, SWP, np.repeat(GWK_[:, None], N, 1))
T_MEDIAN_DAY = ginv(MU_NEW)
share = lambda T, f=lambda v: v > LINE, m=slice(None): (f(T[:, m]) * W[m]).sum(1) / W[m].sum()
USUAL = share(T_NEW)
P_HIGHER = (USUAL > ONE_DAY).mean()
print(f"one-day share over the line (weighted):          {ONE_DAY:.3f}")
print(f"usual intake, fresh people (population):         {q(USUAL)}   P(usual > one-day) = {P_HIGHER:.3f}")
print(f"usual intake, the sampled people:                {q(share(T_SAMPLED))}")
print(f"'median day' back-transform, fresh people:       {q(share(T_MEDIAN_DAY))}")
print("weighted 5/50/95 percentiles of usual intake (mg):",
      np.median([wquant(t, W, [0.05, 0.5, 0.95]) for t in T_NEW[:200]], axis=0).round(-1))
TIMES["derived"] = time.time() - t0
del idata, th, U_, MU_NEW                   # derived arrays are all we need from here

one-day share over the line (weighted):          0.740
usual intake, fresh people (population):         [0.85 0.86 0.88]   P(usual > one-day) = 1.000
usual intake, the sampled people:                [0.86 0.87 0.88]
'median day' back-transform, fresh people:       [0.8  0.82 0.84]
weighted 5/50/95 percentiles of usual intake (mg): [1890. 3360. 5680.]

About 86% of US adults (89% interval 85-88%) usually eat more than 2,300 mg of sodium, against 74% judged by one day. In every one of the 1,000 posterior draws the usual-intake share is higher than the one-day share. The usual-intake median is about 3,400 mg and the middle 90% of adults sit between about 1,900 and 5,700 mg, a much narrower range than the 1,200-7,000 mg of single days. The sampled people's own usual intakes give almost the same share (87%); the fresh-person version is the population estimate, and slightly wider, because it does not lean on the particular 4,863 people.

The back-transform is not a detail. Using each person's median day, $g^{-1}(\mu_i)$, instead of their mean day gives 82%. A naive analysis that models on a transformed scale and then reports $g^{-1}$ of the fitted level understates everyone's intake, because a skewed day distribution has a mean above its median.

#Three lines, two directions

The direction of the one-day error depends on where the line sits, exactly as in the toy. Here are the lower limit sometimes quoted for people with high blood pressure (1,500 mg), the guideline (2,300 mg), and a heavy-intake mark (5,000 mg).

rows = []
for lab, f in [("under 1,500 mg", lambda v: v < 1500), ("over 2,300 mg", lambda v: v > LINE),
               ("over 5,000 mg", lambda v: v > 5000)]:
    us = share(T_NEW, f)
    rows.append((lab, np.average(f(x1), weights=w1), np.median(us), np.quantile(us, 0.055), np.quantile(us, 0.945)))
print(pd.DataFrame(rows, columns=["share of adults", "one day", "usual", "usual 5.5%", "usual 94.5%"]
                   ).round(3).to_string(index=False))

fig, ax = plt.subplots(figsize=(11, 4.6))
bins = np.geomspace(300, 30000, 70)
ax.hist(x1, bins=bins, weights=w1, color=LIGHT, label="one day (weighted)")
for k in range(30):
    h, _ = np.histogram(T_NEW[k], bins=bins, weights=W)
    ax.stairs(h, bins, color=HI, alpha=0.15, lw=1.2)
ax.plot([], [], color=HI, label="usual intake (30 posterior draws)")
ymax = ax.get_ylim()[1]
for v, t, (od, us), (ha, yy) in zip((1500, LINE, 5000), ("1,500", "2,300", "5,000"),
                                    [(r[1], r[2]) for r in rows], [("right", 0.93), ("right", 0.6), ("left", 0.93)]):
    ax.axvline(v, color=LO, lw=1.1, ls="-" if v == LINE else ":")
    side = "under" if v == 1500 else "over"
    ax.text(v * (0.96 if ha == "right" else 1.04), ymax * yy, f"{t} mg\n{side}: one day {od:.0%},\nusual {us:.0%}",
            color=LO, fontsize=9.5, va="top", ha=ha, bbox=dict(fc="white", ec="none", alpha=0.85))
ax.set_xscale("log"); ax.set_yticks([])
ax.set_xlabel("sodium (mg a day, log scale)")
ax.legend(loc="center right", fontsize=10, frameon=False)
ax.set_title("Real adults: the squeeze lifts the share over 2,300 mg and shrinks both tails", loc="left");

share of adults  one day  usual  usual 5.5%  usual 94.5%
 under 1,500 mg    0.091  0.012       0.008        0.018
  over 2,300 mg    0.740  0.864       0.847        0.880
  over 5,000 mg    0.174  0.111       0.096        0.126

Under 1,500 mg, the one-day share is 9% and the usual share about 1%: almost nobody usually eats that little, though many people have a light day. Over 5,000 mg, one day says 17% and the usual intake about 11%. Over 2,300 mg, the error flips sign. A rule of thumb follows: a single day exaggerates how many people are in a tail and understates how many are on the side of the line where the middle lies.

The tail shares depend more on the model than the headline does. To see how much, we refit the same model on the plain log scale, where the usual intake has the closed form $E[x\mid\mu_i] = \exp(\mu_i + \sigma_w^2/2)$ (with $\mu_i$ and $\sigma_w$ on the log scale).

t0 = time.time()
LX = np.log(x_all)
LM, LS = LX[d.day.to_numpy() == 1].mean(), LX[d.day.to_numpy() == 1].std()
with sodium_model(y=(LX - LM) / LS):
    idata_log = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
tl = idata_log.isel(draw=slice(None, None, 4))
XL = lambda n: az.extract(tl, var_names=[n]).values
mu_l = ((XL("a") + XL("b_cell")[CI] + XL("b_pir")[PI] + XL("b_race")[RI]).T
        + XL("s_b").T[:, CI] * np.random.default_rng(13).standard_normal((S, N)))
sw_l, gw_l = XL("s_w").T[:, CI], XL("g_wk")[:, None]
T_LOG = sum(wt * np.exp(LM + LS * (mu_l + gw_l * wkv) + (LS * sw_l) ** 2 / 2)
            for wkv, wt in ((-3 / 7, 4 / 7), (4 / 7, 3 / 7)))
print(f"log scale: divergences {int(idata_log.sample_stats['diverging'].sum())}, "
      f"max r_hat (s_w, s_b) {max(float(az.rhat(idata_log, var_names=[v]).to_dataset()[v].max()) for v in ['s_w', 's_b']):.3f}")
for lab, f in [("under 1,500 mg", lambda v: v < 1500), ("over 2,300 mg", lambda v: v > LINE),
               ("over 5,000 mg", lambda v: v > 5000)]:
    print(f"{lab:15s} Box-Cox {q(share(T_NEW, f))}   log {q(share(T_LOG, f))}")
del idata_log, tl, mu_l, sw_l, T_LOG
TIMES["fit_log"] = time.time() - t0

log scale: divergences 0, max r_hat (s_w, s_b) 1.004
under 1,500 mg  Box-Cox [0.01 0.01 0.02]   log [0.01 0.01 0.02]
over 2,300 mg   Box-Cox [0.85 0.86 0.88]   log [0.84 0.85 0.87]
over 5,000 mg   Box-Cox [0.1  0.11 0.13]   log [0.13 0.14 0.16]

The headline hardly cares about the scale, but the share over 5,000 mg rises from about 11% to about 14% on the log scale, whose heavier upper tail goes with its left-skewed days. The transformation is a modelling choice that matters most in the tails.

#Sex and income

Jackson et al. (2016) found 98% of men and 80% of women over the line in 2009-2012. And the one-day shares differed by 11 points between the poorest and the richest income bands. Does that gap survive once the noise is gone?

male = (per.sex == "male").to_numpy()
for lab, m in (("men", male), ("women", ~male)):
    mm = m[p1]
    print(f"{lab:6s} one-day {np.average(x1[mm] > LINE, weights=w1[mm]):.3f}   usual {q(share(T_NEW, m=m))}")
BAND = {}
for j, lab in enumerate(PIR_LAB):
    m = PI == j
    mm = m[p1]
    BAND[lab] = (np.average(x1[mm] > LINE, weights=w1[mm]), share(T_NEW, m=m))
    print(f"income {lab:8s} n={m.sum():5d}  one-day {BAND[lab][0]:.3f}   usual {q(BAND[lab][1])}")
gap_one = BAND["4+"][0] - BAND["<1"][0]
gap = BAND["4+"][1] - BAND["<1"][1]
P_GAP = (gap > 0).mean()
print(f"gap 4+ minus <1: one day {gap_one:+.3f}; usual {q(gap)}; P(poorest band lower) = {P_GAP:.3f}")

# can the model reproduce the observed one-day shares? one simulated recall day per fresh person
r4 = np.random.default_rng(14)
wk1 = d1.weekend.to_numpy() - 3 / 7
rep1 = []
for k in range(0, S, 5):
    o_ = ginv(MFIX[k, p1] + SB_[k, CI[p1]] * r4.standard_normal(len(p1)) + GWK_[k] * wk1
              + SW_[k, CI[p1]] * r4.standard_normal(len(p1))) > LINE
    rep1.append([o_.mean(), np.average(o_, weights=w1)]
                + [np.average(o_[PI[p1] == j], weights=w1[PI[p1] == j]) for j in (0, 4)]
                + [o_[PI[p1] == j].mean() for j in (0, 4)])
rep1 = np.array(rep1)
gap_one_uw = (x1[PI[p1] == 4] > LINE).mean() - (x1[PI[p1] == 0] > LINE).mean()
print(f"one-day share, fresh people: unweighted {q(rep1[:, 0])} (observed {(x1 > LINE).mean():.3f}); "
      f"weighted {q(rep1[:, 1])} (observed {ONE_DAY:.3f})")
print(f"one-day gap 4+ minus <1, fresh people: weighted {q(rep1[:, 3] - rep1[:, 2])} (observed {gap_one:+.3f}); "
      f"unweighted {q(rep1[:, 5] - rep1[:, 4])} (observed {gap_one_uw:+.3f})")
print("two-day method of moments on the model scale (unadjusted):",
      {lab: (round(float(np.var((YB - YA)[PI[pp2] == j]) / 2), 2),
             round(float(np.cov(YA[PI[pp2] == j], YB[PI[pp2] == j])[0, 1]), 2)) for j, lab in enumerate(PIR_LAB[:5])},
      "(within, between)")
print("share of adults under 31, by income band:",
      {lab: round(float(np.average(per.age.to_numpy()[PI == j] < 31, weights=W[PI == j])), 2)
       for j, lab in enumerate(PIR_LAB[:5])})

men    one-day 0.831   usual [0.94 0.95 0.96]
women  one-day 0.655   usual [0.75 0.78 0.81]
income <1       n=  810  one-day 0.668   usual [0.79 0.83 0.86]
income 1-2      n= 1196  one-day 0.695   usual [0.82 0.85 0.88]
income 2-3      n=  709  one-day 0.755   usual [0.84 0.88 0.91]
income 3-4      n=  468  one-day 0.779   usual [0.84 0.88 0.92]
income 4+       n= 1108  one-day 0.783   usual [0.85 0.88 0.91]
income unknown  n=  572  one-day 0.696   usual [0.77 0.82 0.87]
gap 4+ minus <1: one day +0.115; usual [0.01 0.06 0.1 ]; P(poorest band lower) = 0.989
one-day share, fresh people: unweighted [0.69 0.7  0.71] (observed 0.708); weighted [0.7  0.71 0.74] (observed 0.740)
one-day gap 4+ minus <1, fresh people: weighted [0.02 0.06 0.11] (observed +0.115); unweighted [0.03 0.06 0.1 ] (observed +0.099)
two-day method of moments on the model scale (unadjusted): {'<1': (0.7, 0.46), '1-2': (0.55, 0.4), '2-3': (0.64, 0.31), '3-4': (0.53, 0.31), '4+': (0.46, 0.36)} (within, between)
share of adults under 31, by income band: {'<1': 0.32, '1-2': 0.26, '2-3': 0.24, '3-4': 0.25, '4+': 0.14}

By sex the picture has the 2009-2012 shape: 95% of men and 78% of women are usually over the line (98% and 80% then), against 83% and 66% on one day. By income, the usual-intake shares are all in the 80s, from about 83% below the poverty line to about 88% at twice it and above. The usual gap between the poorest and the richest band is about 6 points (89% interval 1-10), and it is very likely (98.9%) that the poorest band is usually over less often.

It is tempting to say the gap shrank from 11 points to 6 because noise was removed. The replicated recall days say otherwise. Simulated one-day recalls for fresh people with the same covariates reproduce the unweighted one-day share (70% against 71%) but come out at about 71% weighted (70-74%), against the observed 74%. Their one-day income gap is about 6 points (2-11 weighted, 3-10 unweighted), against an observed 10 points unweighted and 11.5 weighted. So the weights explain only a small part of the shortfall; most of it is the model. This model understates the income gap, and its usual-intake 6 points is best read as a low-leaning estimate whose direction is solid and whose size is not: from it alone we cannot tell whether the usual gap is smaller than the one-day gap at all. The same misfit leaves the weighted one-day share about 3 points low, so the population headline may be a point or two low too.

Two features of the model are candidates. Its income effect is additive, the same on top of every age, sex and race group, with no interactions. And its variance split is shared across income bands, although the poorest band is younger (32% under 31, against 14% in the richest), and the unadjusted two-day moments give it more spread both within people (0.70 against 0.46 on the model's scale) and between them (0.46 against 0.36). Part E lets the split vary by income and checks whether that reproduces the one-day gap. Income bands differ in their age, sex and race mix, so these are differences between groups as they are, not effects of income.

#E · What does the second day buy?

#Integrating the person out

Everything above rests on the second recall. To see that directly we rewrite the model so that the person effects are integrated out, which leaves only the quantities a replicate can inform. Let $s_c$ be the total one-day sd in cell $c$ and $f$ the share of one-day variance that is within-person, so $\sigma_w^2 = f s^2$ and $\sigma_b^2 = (1-f)s^2$. Write $r_{ij}$ for the residual of day $j$ after the covariates and day effects. For a person with one recall,

$$r_{i1} \sim N(0,\ \sigma_b^2 + \sigma_w^2) = N(0,\ s^2),$$

and $f$ does not appear at all. For a person with two recalls, the half-sum and the difference are independent normals,

$$\tfrac12(r_{i1}+r_{i2}) \sim N\big(0,\ \sigma_b^2 + \tfrac12\sigma_w^2\big),\qquad r_{i2}-r_{i1} \sim N(0,\ 2\sigma_w^2),$$

and the difference speaks to $\sigma_w$ alone. With only first days, the likelihood is flat in $f$, and its posterior must be exactly its prior, $f \sim \text{Beta}(4, 3)$. The marginal model has no 4,863 person effects, so it samples quickly, and the usual-intake share still follows from each draw by simulating fresh people.

def marginal_fit(dd, label):
    """Person effects integrated out; dd holds the recall rows to use."""
    cnt = dd.groupby("p").day.transform("count").to_numpy()
    one = dd[cnt == 1]
    two = dd[cnt == 2].sort_values(["p", "day"])
    t1, t2 = two[two.day == 1], two[two.day == 2]
    assert (t1.p.to_numpy() == t2.p.to_numpy()).all()

    def parts(rows):
        return (rows.p.to_numpy(), (rows.day.to_numpy() == 2) * 1.0, rows.weekend.to_numpy() - 3 / 7,
                Y[rows.index.to_numpy()])

    with pm.Model(coords=COORDS) as m:
        a = pm.Normal("a", 0, 1)
        bc = pm.ZeroSumNormal("b_cell", 0.5, dims="cell")
        bp = pm.ZeroSumNormal("b_pir", 0.5, dims="pir")
        br = pm.ZeroSumNormal("b_race", 0.5, dims="race")
        g2, gw = pm.Normal("g_day2", 0, 0.5), pm.Normal("g_wk", 0, 0.5)
        s_tot = pm.LogNormal("s_tot", 0, 0.3, dims="cell")
        f = pm.Beta("f", 4, 3)
        fix = a + bc[CI] + bp[PI] + br[RI]
        mean = lambda pp, dd2, ww: fix[pp] + g2 * dd2 + gw * ww
        po, do, wo, yo = parts(one)
        pm.Normal("y_one", mean(po, do, wo), s_tot[CI[po]], observed=yo)
        if len(t1):
            pa, da, wa, ya = parts(t1)
            pb_, db, wb, yb = parts(t2)
            st = s_tot[CI[pa]]
            pm.Normal("half_sum", (mean(pa, da, wa) + mean(pb_, db, wb)) / 2,
                      st * pt.sqrt(1 - f + f / 2), observed=(ya + yb) / 2)
            pm.Normal("diff", mean(pb_, db, wb) - mean(pa, da, wa), st * pt.sqrt(2 * f), observed=yb - ya)
        idm = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
    rh = max(float(az.rhat(idm, var_names=[v]).to_dataset()[v].max()) for v in ["a", "s_tot", "f", "b_pir"])
    fq = idm.posterior["f"].values.ravel()
    # usual-intake share from fresh people, 1,000 draws
    tk = idm.isel(draw=slice(None, None, 4))
    Z = lambda n: az.extract(tk, var_names=[n]).values
    fx = (Z("a") + Z("b_cell")[CI] + Z("b_pir")[PI] + Z("b_race")[RI]).T
    st, ff, gwk = Z("s_tot").T[:, CI], Z("f"), Z("g_wk")
    mu_new = fx + st * np.sqrt(1 - ff)[:, None] * np.random.default_rng(4).standard_normal(fx.shape)
    Tn = chunked(usual, mu_new, st * np.sqrt(ff)[:, None], np.repeat(gwk[:, None], N, 1))
    sh = share(Tn)
    print(f"{label:22s} {len(t1):5d} people with 2 days  divergences {int(idm.sample_stats['diverging'].sum())}  "
          f"max r_hat {rh:.3f}  f {np.round(np.quantile(fq, [0.055, 0.5, 0.945]), 2)}  "
          f"share over line {np.round(np.quantile(sh, [0.055, 0.5, 0.945]), 2)}")
    return fq, sh


t0 = time.time()
print("prior f 89% interval:", np.round(stats.beta(4, 3).ppf([0.055, 0.5, 0.945]), 2))
r5 = np.random.default_rng(5)
two_ids = d.groupby("p").day.count()
keep300 = r5.choice(two_ids.index[two_ids == 2].to_numpy(), 300, replace=False)
FITS = {"both days": marginal_fit(d, "both days"),
        "day 1 only": marginal_fit(d[d.day == 1], "day 1 only"),
        "300 second days": marginal_fit(d[(d.day == 1) | d.p.isin(keep300)], "300 second days")}
TIMES["part_E"] = time.time() - t0

fig, axs = plt.subplots(1, 2, figsize=(12, 4.2))
grid = np.linspace(0, 1, 300)
axs[0].fill_between(grid, stats.beta(4, 3).pdf(grid), color=LIGHT, label="prior Beta(4, 3)")
for (lab, (fq, _)), col in zip(FITS.items(), (HI, MUTED, LO)):
    axs[0].hist(fq, bins=np.linspace(0, 1, 101), density=True, histtype="step", lw=2, color=col, label=lab)
axs[0].set_xlabel("f: within-person share of one-day variance"); axs[0].set_yticks([])
axs[0].legend(frameon=False, fontsize=10)
axs[0].set_title("With day 1 only, the data say nothing about f", loc="left")
for i, ((lab, (_, sh)), col) in enumerate(zip(FITS.items(), (HI, MUTED, LO))):
    lo_, md_, hi_ = np.quantile(sh, [0.055, 0.5, 0.945])
    axs[1].plot([lo_, hi_], [i, i], color=col, lw=6, solid_capstyle="butt")
    axs[1].plot([md_], [i], "o", color="white", mec=col, ms=8)
    axs[1].text(hi_ + 0.005, i, f"{md_:.0%} ({lo_:.0%}-{hi_:.0%})", va="center", fontsize=10.5)
axs[1].axvline(ONE_DAY, color=INK, ls=":", lw=1.2)
axs[1].text(ONE_DAY + 0.004, -0.75, "one-day share", fontsize=10, ha="left", va="center")
axs[1].set_yticks(range(3), list(FITS)); axs[1].set_ylim(2.5, -1.0)
axs[1].set_xlim(0.7, 1.0); axs[1].set_xlabel("usual-intake share over 2,300 mg (89% interval)")
axs[1].set_title("...so the headline inherits the prior's width", loc="left");