""" Shrinking clownfish: a scaling-physiology and survival-odds model ================================================================= Science Journaling Club, September 2026. Companion calculation for our field note on Versteeg, M. A., et al. (2025), "Individual clown anemonefish shrink to survive heat stress and social conflict", Science Advances 11(21), eadt7079, doi:10.1126/sciadv.adt7079. THIS IS NOT THE PAPER'S ANALYSIS. Versteeg and colleagues measured 134 wild Amphiprion percula in Kimbe Bay with callipers, six times over five lunar months, and fitted Cox proportional-hazards models to the survival data. Everything below is the club's own arithmetic, built from published scaling relationships and from numbers we typed in by hand. We have no fish, no anemone, no oxygen probe and no thermistor. We have a laptop. WHAT THE MODEL DOES ------------------- PART A - body geometry. A clownfish is idealised as a triaxial ellipsoid: standard length SL = 0.82 x TL, body depth D = 0.45 x SL, body width W = 0.18 x SL. Volume is the exact ellipsoid volume; surface area is the Knud Thomsen approximation (p = 1.6075, error < 1.1% for these axis ratios). Mass is volume times a tissue density of 1.05 g/cm^3. Because the shape is held fixed, surface-to-volume ratio is exactly proportional to 1/TL, and mass is exactly proportional to TL^3. Every ratio we report is therefore scale-free: the calibration constants cancel. PART B - the oxygen squeeze. Oxygen SUPPLY capacity is taken as S(L, T) = k_s * M(L)^d_g * C_O2(T) where d_g is the gill-surface-area scaling exponent and C_O2(T) is the air-equilibrium oxygen concentration of seawater at salinity 35, from the Garcia & Gordon (1992) refit of the Benson & Krause data. Oxygen DEMAND is taken as Dm(L, T) = k_d * M(L)^0.75 * Q10^((T - T0)/10) with Q10 = 2.0, a standard value for teleost resting metabolism. The supply-to-demand RATIO is normalised to 1.0 at the reference state (TL = 67.50 mm, T = 28.5 degC), so R = (L/L0)^(3*(d_g - 0.75)) * [C_O2(T)/C_O2(T0)] / Q10^((T-T0)/10). Baseline d_g = 0.70. That is the gill-limitation-friendly end of the published spread and it is chosen deliberately: it is the assumption under which shrinking helps at all. We sweep d_g from 0.60 to 0.90 and report the sign flip at d_g = 0.75, where the geometric benefit of shrinking vanishes. Two currencies are reported, because they give answers that differ by more than an order of magnitude: (i) the RATIO currency: how much warming a given shrink offsets, if what limits the fish is supply per unit demand; (ii) the ABSOLUTE currency: how much warming a given shrink offsets, if what limits the fish is the total oxygen or food bill, M^0.75. PART C - survival odds. Simulated, not observed. We build a population of breeding pairs, let each fish shrink or not with the rank-specific frequencies the paper reports, give shrinkers a shrink magnitude, and generate survival from logit(p) = b0 + b1 * s + b2 * tandem where s is percent of total length lost and tandem = 1 when BOTH members of a pair shrank. b1 is solved by bisection so that the population-level shrank-versus-did-not odds ratio is exactly 1.78, matching the paper's reported 78% survival benefit. Outcomes are then drawn once and a logistic regression is fitted back by iteratively reweighted least squares, in plain numpy. We recover the odds ratio, the survival-probability curve against shrink magnitude, and a Wald test of tandem against solo shrinking among shrinkers only. A power replicate then asks what a study of 134 fish with the paper's observed 8% mortality could have detected with a plain binary endpoint. ASSUMPTIONS, STATED PLAINLY --------------------------- 1. Shape is conserved under shrinking. A fish that loses 1.4% of its length is assumed geometrically similar to its former self. Nobody has shown this. 2. Gill surface area scales as a single power of mass within an individual. Gills are plastic; this ignores acclimation entirely. 3. Q10 = 2.0 and is constant across 26-34 degC. Near a thermal limit it is not. 4. Oxygen supply is treated as proportional to ambient concentration. Real uptake depends on the partial-pressure gradient across the lamellae, ventilation rate, haemoglobin affinity and boundary-layer thickness, none of which appear here. 5. Demand carries no activity, digestion, gonad or stress component. 6. The survival model is a single binary endpoint. The paper used time-to- event data with repeated measures, which carries more information. 7. Shrinking is treated as a cause of survival. In an observational study it could equally be a marker of a fish that was in good enough condition to afford it. Our simulation assumes causation because we built it in. 8. Salinity fixed at 35, pressure at 1 atm, air saturation assumed. LIMITS ------ This model cannot tell you why clownfish shrink. It can tell you how large the purely geometric benefit is under stated assumptions, and the answer is small. It reproduces an odds ratio we typed in rather than one we measured. Nothing here is evidence about fish. It is arithmetic about a claim. Run: python shrinking-clownfish.py """ import math import numpy as np SEED = 20250521 # the paper's publication date rng = np.random.default_rng(SEED) OUT = [] def say(line=""): print(line) OUT.append(line) def rule(title): say() say("=" * 72) say(title) say("=" * 72) # ---------------------------------------------------------------------- # PART A. Body geometry # ---------------------------------------------------------------------- SL_FRAC = 0.82 # standard length as a fraction of total length DEPTH_FRAC = 0.45 # body depth as a fraction of standard length WIDTH_FRAC = 0.18 # body width as a fraction of standard length RHO_TISSUE = 1.05e-3 # g / mm^3 THOMSEN_P = 1.6075 TL_MIN, TL_MAX = 27.00, 83.50 # observed range, mm TL TL_R1, TL_R2 = 67.50, 52.60 # mean initial TL, rank 1 and rank 2 TL_REF = TL_R1 def axes_mm(tl): """Semi-axes of the idealised body ellipsoid, in mm.""" sl = SL_FRAC * tl return sl / 2.0, (DEPTH_FRAC * sl) / 2.0, (WIDTH_FRAC * sl) / 2.0 def volume_mm3(tl): a, b, c = axes_mm(tl) return (4.0 / 3.0) * math.pi * a * b * c def surface_mm2(tl): a, b, c = axes_mm(tl) p = THOMSEN_P term = (a ** p * b ** p + a ** p * c ** p + b ** p * c ** p) / 3.0 return 4.0 * math.pi * term ** (1.0 / p) def mass_g(tl): return volume_mm3(tl) * RHO_TISSUE def sa_vol(tl): return surface_mm2(tl) / volume_mm3(tl) rule("PART A. Body geometry across the observed length range") say("Idealised body: triaxial ellipsoid, SL = 0.82*TL, depth = 0.45*SL,") say("width = 0.18*SL, tissue density 1.05 g/cm^3. Shape held constant,") say("so SA/V is exactly proportional to 1/TL and mass to TL^3.") say() say(f"{'TL (mm)':>9} {'mass (g)':>10} {'area (mm2)':>12} {'SA/V (1/mm)':>13} {'rel. SA/V':>10}") say("-" * 58) geom_rows = [] for tl in [27.0, 35.0, 45.0, 52.6, 60.0, 67.5, 75.0, 83.5]: r = sa_vol(tl) / sa_vol(TL_REF) geom_rows.append((tl, mass_g(tl), surface_mm2(tl), sa_vol(tl), r)) say(f"{tl:9.2f} {mass_g(tl):10.3f} {surface_mm2(tl):12.1f} {sa_vol(tl):13.5f} {r:10.3f}") say() say(f"Mass range over the observed lengths: {mass_g(TL_MIN):.3f} g at " f"{TL_MIN:.2f} mm to {mass_g(TL_MAX):.2f} g at {TL_MAX:.2f} mm, " f"a factor of {mass_g(TL_MAX)/mass_g(TL_MIN):.1f}.") say(f"SA/V at the rank-1 mean ({TL_R1:.2f} mm): {sa_vol(TL_R1):.5f} /mm") say(f"SA/V at the rank-2 mean ({TL_R2:.2f} mm): {sa_vol(TL_R2):.5f} /mm " f"({100*(sa_vol(TL_R2)/sa_vol(TL_R1) - 1):.1f}% higher)") SHRINKS = [0.5, 1.0, 1.4, 2.0, 2.5, 3.0, 4.0, 5.0] say() say("Effect of a shrink on the rank-1 mean fish:") say(f"{'shrink %TL':>11} {'new TL':>8} {'new mass':>9} {'mass lost %':>12} {'SA/V gain %':>12}") say("-" * 56) for s in SHRINKS: tl2 = TL_REF * (1 - s / 100.0) dm = 100 * (1 - mass_g(tl2) / mass_g(TL_REF)) dsv = 100 * (sa_vol(tl2) / sa_vol(TL_REF) - 1) say(f"{s:11.1f} {tl2:8.2f} {mass_g(tl2):9.3f} {dm:12.2f} {dsv:12.2f}") # ---------------------------------------------------------------------- # PART B. The oxygen squeeze # ---------------------------------------------------------------------- # Garcia & Gordon (1992) combined fit to Benson & Krause, umol/kg, S in PSU. GG_A = [5.80871, 3.20291, 4.17887, 5.10006, -9.86643e-2, 3.80369] GG_B = [-7.01577e-3, -7.70028e-3, -1.13864e-2, -9.51519e-3] GG_C0 = -2.75915e-7 SALINITY = 35.0 def o2_sat_umol_kg(t_c, s=SALINITY): """Air-equilibrium O2 concentration of seawater, umol/kg.""" ts = math.log((298.15 - t_c) / (273.15 + t_c)) lnc = sum(GG_A[i] * ts ** i for i in range(6)) lnc += s * sum(GG_B[i] * ts ** i for i in range(4)) lnc += GG_C0 * s * s return math.exp(lnc) def o2_sat_mg_l(t_c, s=SALINITY): """Same, converted to mg/L using a seawater density of 1023 kg/m^3.""" return o2_sat_umol_kg(t_c, s) * 31.9988e-3 * 1.023 T_REF = 28.5 # club's pre-heatwave reference, degC T_HEAT = 32.5 # reference + the 4 degC anomaly the paper reports Q10 = 2.0 D_GILL = 0.70 # baseline gill-area scaling exponent D_METAB = 0.75 # metabolic scaling exponent def supply_rel(tl, t_c, d_g=D_GILL): return (mass_g(tl) / mass_g(TL_REF)) ** d_g * (o2_sat_umol_kg(t_c) / o2_sat_umol_kg(T_REF)) def demand_rel(tl, t_c): return (mass_g(tl) / mass_g(TL_REF)) ** D_METAB * Q10 ** ((t_c - T_REF) / 10.0) def ratio_rel(tl, t_c, d_g=D_GILL): """Supply-to-demand ratio, normalised to 1.0 at (TL_REF, T_REF).""" return supply_rel(tl, t_c, d_g) / demand_rel(tl, t_c) rule("PART B1. Oxygen solubility falls while demand rises") say(f"Salinity {SALINITY:.0f}, air saturation, 1 atm. Reference state: " f"TL {TL_REF:.2f} mm at {T_REF:.1f} degC.") say() say(f"{'T (degC)':>9} {'O2 (umol/kg)':>13} {'O2 (mg/L)':>11} {'supply rel':>11} " f"{'demand rel':>11} {'S/D ratio':>10}") say("-" * 70) temps = [26.0, 27.0, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.0] squeeze_rows = [] for t in temps: sup = o2_sat_umol_kg(t) / o2_sat_umol_kg(T_REF) dem = Q10 ** ((t - T_REF) / 10.0) squeeze_rows.append((t, o2_sat_umol_kg(t), o2_sat_mg_l(t), sup, dem, sup / dem)) say(f"{t:9.1f} {o2_sat_umol_kg(t):13.2f} {o2_sat_mg_l(t):11.3f} " f"{sup:11.4f} {dem:11.4f} {sup/dem:10.4f}") o2_ref = o2_sat_umol_kg(T_REF) o2_heat = o2_sat_umol_kg(T_HEAT) sup_drop = 100 * (1 - o2_heat / o2_ref) dem_rise = 100 * (Q10 ** ((T_HEAT - T_REF) / 10.0) - 1) squeeze = ratio_rel(TL_REF, T_HEAT) say() say(f"Over the {T_HEAT - T_REF:.1f} degC anomaly, from {T_REF:.1f} to {T_HEAT:.1f} degC:") say(f" dissolved oxygen falls {sup_drop:.2f}% " f"({o2_ref:.2f} to {o2_heat:.2f} umol/kg, {o2_sat_mg_l(T_REF):.3f} to " f"{o2_sat_mg_l(T_HEAT):.3f} mg/L)") say(f" metabolic demand rises {dem_rise:.2f}% at Q10 = {Q10:.1f}") say(f" supply-to-demand ratio for a fish that does not change size: " f"{squeeze:.4f}") say(f" that is a squeeze of {100*(1-squeeze):.2f}% on the oxygen budget.") # local log-slope of the ratio with temperature, used for offset arithmetic def ratio_temp_factor(t_c): return (o2_sat_umol_kg(t_c) / o2_sat_umol_kg(T_REF)) / Q10 ** ((t_c - T_REF) / 10.0) def offset_temp_ratio(shrink_pct, d_g=D_GILL, t0=T_REF): """degC of warming that a given shrink offsets, in the RATIO currency.""" gain = (1 - shrink_pct / 100.0) ** (3.0 * (d_g - D_METAB)) target = 1.0 / gain # ratio_temp_factor falls monotonically with temperature, so the signed # bracket lets d_g above 0.75 return a negative offset: shrinking there # costs the fish warming rather than buying it any. lo, hi = -15.0, 15.0 for _ in range(300): mid = 0.5 * (lo + hi) f = ratio_temp_factor(t0 + mid) / ratio_temp_factor(t0) if f > target: lo = mid else: hi = mid return 0.5 * (lo + hi) def offset_temp_absolute(shrink_pct, t0=T_REF): """degC of warming that a given shrink offsets, in the ABSOLUTE-demand currency: the shrink cuts M^0.75, warming raises Q10^(dT/10).""" cut = (1 - shrink_pct / 100.0) ** (3.0 * D_METAB) return 10.0 * math.log(1.0 / cut) / math.log(Q10) rule("PART B2. How much warming does shrinking actually offset?") say("Two currencies. RATIO asks how much warming is cancelled if the fish is") say("limited by supply per unit demand. ABSOLUTE asks the same question if the") say("fish is limited by its total oxygen and food bill. Same model, same") say("shrink, answers an order of magnitude apart.") say() say(f"{'shrink %TL':>11} {'ratio gain %':>13} {'offset degC':>12} " f"{'demand cut %':>13} {'offset degC':>12}") say(f"{'':>11} {'(d_g=0.70)':>13} {'(ratio)':>12} {'(M^0.75)':>13} {'(absolute)':>12}") say("-" * 65) offset_rows = [] for s in SHRINKS: gain = 100 * ((1 - s / 100.0) ** (3.0 * (D_GILL - D_METAB)) - 1) dt_r = offset_temp_ratio(s) cut = 100 * (1 - (1 - s / 100.0) ** (3.0 * D_METAB)) dt_a = offset_temp_absolute(s) offset_rows.append((s, gain, dt_r, cut, dt_a)) say(f"{s:11.1f} {gain:13.3f} {dt_r:12.4f} {cut:13.2f} {dt_a:12.4f}") S_TYPICAL = 1.4 dt_r14 = offset_temp_ratio(S_TYPICAL) dt_a14 = offset_temp_absolute(S_TYPICAL) say() say(f"A typical {S_TYPICAL:.1f}% shrink offsets {dt_r14:.4f} degC in the ratio") say(f"currency and {dt_a14:.4f} degC in the absolute currency. The absolute") say(f"answer is {dt_a14/dt_r14:.1f} times larger.") say(f"Against the {T_HEAT - T_REF:.1f} degC anomaly, the ratio currency recovers " f"{100*dt_r14/(T_HEAT-T_REF):.2f}% of the heat;") say(f"the absolute currency recovers {100*dt_a14/(T_HEAT-T_REF):.2f}%.") # how much shrinking would it take to fully cancel 4 degC, ratio currency? need = 1.0 / squeeze frac = need ** (1.0 / (3.0 * (D_GILL - D_METAB))) say() say(f"To cancel the full {T_HEAT-T_REF:.1f} degC in the RATIO currency, a fish would have to") say(f"shrink to {100*frac:.2f}% of its length: {TL_REF:.2f} mm down to " f"{TL_REF*frac:.2f} mm.") say(f"That is a {100*(1-frac):.2f}% reduction, which is not a fish.") frac_abs = (1.0 / (Q10 ** ((T_HEAT - T_REF) / 10.0))) ** (1.0 / (3.0 * D_METAB)) say(f"To cancel it in the ABSOLUTE currency: shrink to {100*frac_abs:.2f}% of length, " f"{TL_REF*frac_abs:.2f} mm.") say(f"That is a {100*(1-frac_abs):.2f}% reduction. Still far beyond anything observed,") say("but within the same order of magnitude as a real animal.") rule("PART B3. Sensitivity to the gill-area exponent") say("d_g is the exponent in gill surface area proportional to mass^d_g. The") say("whole geometric argument lives or dies on whether it is below 0.75.") say() say(f"{'d_g':>6} {'ratio exponent':>16} {'gain at 1.4%':>14} {'offset degC':>12} {'verdict':>22}") say("-" * 74) dg_rows = [] for d_g in [0.60, 0.65, 0.67, 0.70, 0.75, 0.80, 0.85, 0.90]: expo = 3.0 * (d_g - D_METAB) gain = 100 * ((1 - S_TYPICAL / 100.0) ** expo - 1) dt = offset_temp_ratio(S_TYPICAL, d_g) if d_g < D_METAB: verdict = "shrinking helps" elif d_g == D_METAB: verdict = "shrinking is neutral" else: verdict = "shrinking hurts" dg_rows.append((d_g, expo, gain, dt, verdict)) say(f"{d_g:6.2f} {expo:16.3f} {gain:14.3f} {dt:12.4f} {verdict:>22}") say() say("At d_g = 0.75 the two exponents cancel and body size drops out of the") say("ratio entirely. Above it, the sign flips and a smaller fish is worse off.") say("Published intraspecific gill exponents straddle 0.75. This is the single") say("softest joint in the whole calculation.") rule("PART B4. The length that restores balance, temperature by temperature") def restoring_length(t_c, d_g=D_GILL): """TL at which the supply/demand ratio returns to 1.0 at temperature t_c.""" expo = 3.0 * (d_g - D_METAB) if abs(expo) < 1e-12: return float("nan") need = 1.0 / ratio_temp_factor(t_c) return TL_REF * need ** (1.0 / expo) say(f"{'T (degC)':>9} {'S/D at 67.50':>13} {'restoring TL':>13} {'shrink needed %':>16}") say("-" * 55) restore_rows = [] for t in [29.0, 29.5, 30.0, 30.5, 31.0, 31.5, 32.0, 32.5]: r = ratio_rel(TL_REF, t) lr = restoring_length(t) restore_rows.append((t, r, lr, 100 * (1 - lr / TL_REF))) say(f"{t:9.1f} {r:13.4f} {lr:13.4f} {100*(1-lr/TL_REF):16.2f}") say() say("Read the last column as a reductio. Half a degree of warming already asks") say("for a shrink of tens of percent under these assumptions. The fish shrank") say("by one or two. Whatever shrinking is doing, the gill-geometry ratio on its") say("own is not a sufficient explanation.") # ---------------------------------------------------------------------- # PART C. Survival odds # ---------------------------------------------------------------------- N_PAIRS = 3000 P_SHRINK_R1 = 0.71 # paper: 71% of rank 1 shrank at least once P_SHRINK_R2 = 0.79 # paper: 79% of rank 2 shrank at least once MEAN_SHRINK = 1.4 # percent TL, club assumption SHRINK_SD_LOG = 0.55 P0_NONSHRINK = 0.72 # baseline survival in the simulated heatwave TANDEM_OR = 1.55 # true tandem bonus built into the generator TARGET_OR = 1.78 # the paper's reported 78% survival benefit def logistic(z): return 1.0 / (1.0 + np.exp(-z)) def norm_sf(z): """Upper-tail normal probability, no scipy.""" return 0.5 * math.erfc(abs(z) / math.sqrt(2.0)) def irls_logit(X, y, iters=60): """Logistic regression by iteratively reweighted least squares.""" n, k = X.shape beta = np.zeros(k) for _ in range(iters): eta = X @ beta mu = logistic(eta) w = np.clip(mu * (1 - mu), 1e-10, None) z = eta + (y - mu) / w XtW = X.T * w beta_new = np.linalg.solve(XtW @ X, XtW @ z) if np.max(np.abs(beta_new - beta)) < 1e-11: beta = beta_new break beta = beta_new eta = X @ beta mu = logistic(eta) w = np.clip(mu * (1 - mu), 1e-10, None) cov = np.linalg.inv((X.T * w) @ X) return beta, np.sqrt(np.diag(cov)) rule("PART C1. Building the simulated population") rank = np.tile([1, 2], N_PAIRS) shrank = np.where( rank == 1, rng.random(2 * N_PAIRS) < P_SHRINK_R1, rng.random(2 * N_PAIRS) < P_SHRINK_R2, ) mag = np.zeros(2 * N_PAIRS) draws = rng.lognormal(mean=math.log(MEAN_SHRINK) - 0.5 * SHRINK_SD_LOG ** 2, sigma=SHRINK_SD_LOG, size=2 * N_PAIRS) mag[shrank] = np.clip(draws[shrank], 0.2, 8.0) pair_both = shrank.reshape(N_PAIRS, 2).all(axis=1) tandem = np.repeat(pair_both, 2).astype(float) say(f"pairs: {N_PAIRS}, individuals: {2*N_PAIRS}") say(f"shrank at least once: {shrank.sum()} ({100*shrank.mean():.2f}%)") say(f" rank 1: {shrank[rank==1].sum()} ({100*shrank[rank==1].mean():.2f}%), " f"target {100*P_SHRINK_R1:.0f}%") say(f" rank 2: {shrank[rank==2].sum()} ({100*shrank[rank==2].mean():.2f}%), " f"target {100*P_SHRINK_R2:.0f}%") say(f"pairs where both shrank (tandem): {pair_both.sum()} " f"({100*pair_both.mean():.2f}% of pairs)") say(f"shrink magnitude among shrinkers: mean {mag[shrank].mean():.3f}% TL, " f"median {np.median(mag[shrank]):.3f}%, max {mag[shrank].max():.3f}%") say(f" in mm on a {TL_REF:.2f} mm fish: mean " f"{TL_REF*mag[shrank].mean()/100:.3f} mm, max " f"{TL_REF*mag[shrank].max()/100:.3f} mm") b0 = math.log(P0_NONSHRINK / (1 - P0_NONSHRINK)) b2 = math.log(TANDEM_OR) def population_or(b1): """Odds ratio, shrank versus did not, averaged over the built covariates.""" p = logistic(b0 + b1 * mag + b2 * tandem) ps = p[shrank].mean() pn = p[~shrank].mean() return (ps / (1 - ps)) / (pn / (1 - pn)) lo, hi = 0.0, 3.0 for _ in range(300): mid = 0.5 * (lo + hi) if population_or(mid) < TARGET_OR: lo = mid else: hi = mid b1 = 0.5 * (lo + hi) say() say(f"generator coefficients (logit scale):") say(f" b0 intercept = {b0:+.6f} (non-shrinker survival " f"{P0_NONSHRINK:.2f})") say(f" b1 per 1% TL shrink = {b1:+.6f} (solved so the population odds") say(f" ratio is exactly {TARGET_OR:.2f})") say(f" b2 tandem pair = {b2:+.6f} (odds ratio {TANDEM_OR:.2f})") say(f" population OR achieved = {population_or(b1):.6f}") p_true = logistic(b0 + b1 * mag + b2 * tandem) survived = (rng.random(2 * N_PAIRS) < p_true).astype(float) say() say(f"drawn outcomes: {int(survived.sum())} survived of {2*N_PAIRS} " f"({100*survived.mean():.2f}%)") rule("PART C2. Fitted logistic regression") X_bin = np.column_stack([np.ones(2 * N_PAIRS), shrank.astype(float)]) beta_bin, se_bin = irls_logit(X_bin, survived) or_bin = math.exp(beta_bin[1]) or_lo = math.exp(beta_bin[1] - 1.96 * se_bin[1]) or_hi = math.exp(beta_bin[1] + 1.96 * se_bin[1]) z_bin = beta_bin[1] / se_bin[1] say("Model 1: survival ~ shrank (binary). The headline comparison.") say(f" intercept {beta_bin[0]:+.5f} (SE {se_bin[0]:.5f})") say(f" shrank {beta_bin[1]:+.5f} (SE {se_bin[1]:.5f})") say(f" ODDS RATIO {or_bin:.4f} 95% CI {or_lo:.4f} to {or_hi:.4f}") say(f" Wald z {z_bin:.3f}, p = {2*norm_sf(z_bin):.3e}") say(f" paper's reported effect: 78% increase, i.e. OR {TARGET_OR:.2f}") say(f" recovery error: {100*(or_bin/TARGET_OR - 1):+.2f}%") X_dose = np.column_stack([np.ones(2 * N_PAIRS), mag, tandem]) beta_dose, se_dose = irls_logit(X_dose, survived) say() say("Model 2: survival ~ shrink magnitude + tandem. The dose-response.") say(f" intercept {beta_dose[0]:+.5f} (SE {se_dose[0]:.5f}) " f"-> survival {100*logistic(beta_dose[0]):.2f}% at zero shrink") say(f" per 1% TL {beta_dose[1]:+.5f} (SE {se_dose[1]:.5f}) " f"-> OR {math.exp(beta_dose[1]):.4f} per percentage point") say(f" tandem {beta_dose[2]:+.5f} (SE {se_dose[2]:.5f}) " f"-> OR {math.exp(beta_dose[2]):.4f}") say(f" generator truth: b1 {b1:+.5f}, b2 {b2:+.5f}") say() say("Fitted survival probability against shrink magnitude:") say(f"{'shrink %TL':>11} {'solo':>9} {'tandem':>9} {'gap (pp)':>10}") say("-" * 42) curve_rows = [] for s in [0.0, 0.5, 1.0, 1.4, 2.0, 2.5, 3.0, 4.0, 5.0]: solo = logistic(beta_dose[0] + beta_dose[1] * s) tand = logistic(beta_dose[0] + beta_dose[1] * s + beta_dose[2]) curve_rows.append((s, solo, tand)) say(f"{s:11.1f} {100*solo:8.2f}% {100*tand:8.2f}% {100*(tand-solo):10.2f}") p_solo_14 = logistic(beta_dose[0] + beta_dose[1] * 1.4) p_tand_14 = logistic(beta_dose[0] + beta_dose[1] * 1.4 + beta_dose[2]) p_zero = logistic(beta_dose[0]) say() say(f"HEADLINE NUMBERS at the club's default settings:") say(f" non-shrinker survival {100*p_zero:.2f}%") say(f" 1.4% solo shrinker survival {100*p_solo_14:.2f}%") say(f" 1.4% tandem shrinker {100*p_tand_14:.2f}%") say(f" fitted binary odds ratio {or_bin:.4f}") rule("PART C3. Does tandem shrinking beat solo shrinking?") sh = shrank X_t = np.column_stack([np.ones(sh.sum()), mag[sh], tandem[sh]]) beta_t, se_t = irls_logit(X_t, survived[sh]) or_t = math.exp(beta_t[2]) z_t = beta_t[2] / se_t[2] p_t = 2 * norm_sf(z_t) say("Shrinkers only, so the comparison is tandem against solo with shrink") say("magnitude held in the model.") say(f" n = {int(sh.sum())} shrinkers, of whom {int(tandem[sh].sum())} were in " f"tandem pairs") say(f" tandem coefficient {beta_t[2]:+.5f} (SE {se_t[2]:.5f})") say(f" ODDS RATIO {or_t:.4f} 95% CI " f"{math.exp(beta_t[2]-1.96*se_t[2]):.4f} to " f"{math.exp(beta_t[2]+1.96*se_t[2]):.4f}") say(f" Wald z {z_t:.3f}, p = {p_t:.3e}") say(f" verdict: tandem {'beats' if (z_t > 0 and p_t < 0.05) else 'does not beat'} " f"solo at the 5% level") raw_none = survived[~sh].mean() raw_solo = survived[sh & (tandem == 0)].mean() raw_tand = survived[sh & (tandem == 1)].mean() say() say("Raw survival rates, no model:") say(f" did not shrink {100*raw_none:.2f}% (n = {int((~sh).sum())})") say(f" shrank, solo {100*raw_solo:.2f}% (n = {int((sh & (tandem==0)).sum())})") say(f" shrank, tandem {100*raw_tand:.2f}% (n = {int((sh & (tandem==1)).sum())})") say(f" tandem minus solo: {100*(raw_tand - raw_solo):+.2f} percentage points") rule("PART C4. Could 134 fish have seen this with a binary endpoint?") N_STUDY = 134 N_REPS = 4000 P0_FIELD = 0.87 # tuned so overall mortality lands near the paper's 8% b0f = math.log(P0_FIELD / (1 - P0_FIELD)) b1f = math.log(TARGET_OR) rng2 = np.random.default_rng(SEED + 1) detect = 0 ors = [] deaths = [] for _ in range(N_REPS): rk = np.tile([1, 2], N_STUDY // 2) sh2 = np.where(rk == 1, rng2.random(N_STUDY) < P_SHRINK_R1, rng2.random(N_STUDY) < P_SHRINK_R2) p2 = logistic(b0f + b1f * sh2.astype(float)) y2 = (rng2.random(N_STUDY) < p2).astype(float) deaths.append(N_STUDY - y2.sum()) if y2.sum() in (0, N_STUDY) or sh2.sum() in (0, N_STUDY): continue a = ((sh2 == 1) & (y2 == 1)).sum() b = ((sh2 == 1) & (y2 == 0)).sum() c = ((sh2 == 0) & (y2 == 1)).sum() d = ((sh2 == 0) & (y2 == 0)).sum() if min(a, b, c, d) == 0: continue lor = math.log((a * d) / (b * c)) se = math.sqrt(1 / a + 1 / b + 1 / c + 1 / d) ors.append(math.exp(lor)) if lor - 1.96 * se > 0: detect += 1 say(f"{N_REPS} simulated studies of {N_STUDY} fish, non-shrinker survival " f"{P0_FIELD:.2f}, true OR {TARGET_OR:.2f}.") say(f" mean deaths per study: {np.mean(deaths):.2f} " f"({100*np.mean(deaths)/N_STUDY:.2f}% mortality; the paper observed 11, 8.21%)") say(f" usable 2x2 tables: {len(ors)} of {N_REPS}") say(f" median estimated OR: {np.median(ors):.3f}") say(f" studies whose 95% CI excluded 1: {detect} " f"({100*detect/N_REPS:.2f}% power)") say() say("With a plain survived-or-not endpoint, a study this size would almost") say("never reach significance. The paper did not use a binary endpoint. It used") say("repeated measures and time-to-event models, which carry far more") say("information per fish. Read this number as a floor on the design's power,") say("not as a criticism of it.") rule("SUMMARY OF QUOTED FIGURES") say(f"A1 mass at TL 27.00 mm {mass_g(27.0):.3f} g") say(f"A2 mass at TL 83.50 mm {mass_g(83.5):.2f} g") say(f"A3 SA/V at TL 67.50 mm {sa_vol(67.5):.5f} /mm") say(f"A4 SA/V gain from a 1.4% shrink " f"{100*(sa_vol(67.5*0.986)/sa_vol(67.5)-1):.3f}%") say(f"A5 mass lost in a 1.4% shrink " f"{100*(1-mass_g(67.5*0.986)/mass_g(67.5)):.2f}%") say(f"B1 O2 at 28.5 degC {o2_sat_umol_kg(28.5):.2f} umol/kg " f"({o2_sat_mg_l(28.5):.3f} mg/L)") say(f"B2 O2 at 32.5 degC {o2_sat_umol_kg(32.5):.2f} umol/kg " f"({o2_sat_mg_l(32.5):.3f} mg/L)") say(f"B3 solubility drop over 4 degC {sup_drop:.2f}%") say(f"B4 demand rise over 4 degC (Q10 2.0) {dem_rise:.2f}%") say(f"B5 supply/demand ratio at 32.5 degC {squeeze:.4f}") say(f"B6 offset of a 1.4% shrink, ratio {dt_r14:.4f} degC") say(f"B7 offset of a 1.4% shrink, absolute {dt_a14:.4f} degC") say(f"B8 absolute demand cut, 1.4% shrink " f"{100*(1-(0.986)**(3*D_METAB)):.2f}%") say(f"B9 shrink needed to cancel 4 degC, ratio {100*(1-frac):.2f}% " f"(to {TL_REF*frac:.2f} mm)") say(f"B10 gill exponent where sign flips {D_METAB:.2f}") say(f"C1 fitted binary odds ratio {or_bin:.4f} " f"[{or_lo:.4f}, {or_hi:.4f}]") say(f"C2 non-shrinker survival, fitted {100*p_zero:.2f}%") say(f"C3 1.4% solo shrinker survival {100*p_solo_14:.2f}%") say(f"C4 1.4% tandem shrinker survival {100*p_tand_14:.2f}%") say(f"C5 tandem odds ratio among shrinkers {or_t:.4f}, p = {p_t:.3e}") say(f"C6 raw tandem minus solo {100*(raw_tand-raw_solo):+.2f} pp") say(f"C7 power of a 134-fish binary study {100*detect/N_REPS:.2f}%") say() say(f"seed = {SEED}; numpy {np.__version__}") with open(__file__.replace(".py", "-output.txt"), "w", encoding="utf-8") as fh: fh.write("\n".join(OUT) + "\n")