Summary

Key papers (10)

1. Age-Related Differences in Nasopharyngeal Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) Levels in Patients With Mild to Moderate Coronavirus Disease 2019 (COVID-19). Heald-Sargent T, Muller WJ, Zheng X, Rippe J, Patel AB, Kociolek LK; JAMA Pediatrics 174(9):902-903, 2020.
   The origin of the claim under test. In 145 patients tested within a week of symptom onset it reported median Ct 6.5 (IQR 4.8-12.0) in children under 5 years versus 11.1 (6.3-15.7) in children 5-17 and 11.0 (6.9-17.5) in adults on the Abbott m2000 platform, and read the roughly 4.5-cycle gap as a 10- to 100-fold viral excess in the youngest. It supplies the effect size that any replication must be powered to exclude.
2. Nasopharyngeal SARS-CoV-2 viral loads in young children do not differ significantly from those in older children and adults. Madera S, Crawford E, Langelier C, Tran NK, Thornborrow E, Miller S, DeRisi JL; Scientific Reports 11:3044, 2021.
   The largest early US contradiction of the Heald-Sargent finding, using more than five thousand RT-PCR-confirmed cases across multiple centres and reporting no discernible difference in viral nucleic acid between young children and adults. It establishes that the discrepancy is not a matter of sample size alone.
3. Comparison of Symptoms and RNA Levels in Children and Adults With SARS-CoV-2 Infection in the Community Setting. Chung E, Chow EJ, Wilcox NC, et al.; JAMA Pediatrics 175(10):e212025, 2021.
   The most directly comparable published estimate and the model for how the contrast should be reported: an adjusted mean Ct difference of -0.7 cycles (95% CI -2.2 to 0.9) for symptomatic children versus symptomatic adults and -0.6 (95% CI -4.0 to 2.8) for asymptomatic comparisons, against a much larger symptomatic-versus-asymptomatic contrast of about -3 cycles in both age groups. It shows that clinical presentation, not age, structures the Ct distribution.
4. Estimating infectiousness throughout SARS-CoV-2 infection course. Jones TC, Biele G, Muhlemann B, et al.; Science 373(6551):eabi5273, 2021.
   The best-powered study in the field, with 25,381 German cases including 6,110 from test centres attended by presymptomatic, asymptomatic and mildly symptomatic subjects. It places the youngest subjects at most 0.5 log10 units (under about 1.7 cycles) below older subjects and explicitly attributes part of the gap to smaller swab sizes, framing specimen quality as an age-differential measurement problem rather than a nuisance.
5. Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) Viral Load in the Upper Respiratory Tract of Children and Adults With Early Acute Coronavirus Disease 2019 (COVID-19). Baggio S, L'Huillier AG, Yerly S, et al.; Clinical Infectious Diseases 73(1):148-150, 2021.
   Demonstrates the standard confounding control in this literature -- restriction to the first five days after symptom onset -- in 53 children and 352 adults in Geneva, and finds no significant difference. It is the template for handling time since infection when that variable is recorded.
6. Comparison of Upper Respiratory Viral Load Distributions in Asymptomatic and Symptomatic Children Diagnosed with SARS-CoV-2 Infection in Pediatric Hospital Testing Programs. Kociolek LK, Muller WJ, Yee R, et al.; Journal of Clinical Microbiology 59(1):e02593-20, 2020.
   The reference for how much the testing programme itself -- asymptomatic screening versus symptomatic diagnostic testing -- shifts the viral-load distribution within paediatric hospitals. It is the direct justification for treating testing stream as a confounder rather than a design detail.
7. College of American Pathologists (CAP) Microbiology Committee Perspective: Caution Must Be Used in Interpreting the Cycle Threshold (Ct) Value. Rhoads D, Peaper DR, She RC, Nolte FS, Wojewoda CM, Anderson NW, Pritt BS; Clinical Infectious Diseases 72(10):e685-e686, 2021.
   The authoritative statement that Ct is not a quantitative viral load: it depends on assay, target gene, specimen adequacy, collection technique and transport. It is the basis for treating Ct as a proxy measured with error, and for the concern that the error may differ systematically between toddlers and adults.
8. Collider bias undermines our understanding of COVID-19 disease risk and severity. Griffith GJ, Morris TT, Tudball MJ, et al.; Nature Communications 11:5749, 2020.
   The canonical COVID-19 treatment of selection bias in samples defined by having been tested. Since both age and viral load plausibly influence the probability of being tested and of testing positive, conditioning on a positive test can create an age-Ct association in the sample that does not exist in the source population.
9. Bias due to misclassification in the estimation of relative risk. Copeland KT, Checkoway H, McMichael AJ, Holbrook RH; American Journal of Epidemiology 105(5):488-495, 1977.
   Supplies the classical result that non-differential misclassification of a dichotomous exposure biases the estimate toward the null. This is the tool needed when the analytic data have been perturbed for de-identification, because it bounds the direction in which such perturbation can move an estimate.
10. Estimating epidemiologic dynamics from cross-sectional viral load distributions. Hay JA, Kennedy-Shaffer L, Kanjilal S, Lennon NJ, Gabriel SB, Lipsitch M, Mina MJ; Science 373(6552):eabh0635, 2021.
    Shows that the population distribution of Ct values shifts systematically with the epidemic growth rate. It establishes calendar time as a confounder of any Ct comparison spanning a period of changing incidence, and therefore as a variable that must be adjusted for in age contrasts drawn from surveillance or laboratory records.

1. Idea 1Selected

   Testing indication, not age: a replication attempt of the reported viral-load excess in children under 5 years in 2020 paediatric hospital RT-PCR records

   Hypothesis. Among SARS-CoV-2 RT-PCR-positive specimens collected at a single paediatric hospital laboratory, children under 5 years have a mean cycle threshold at least 3 cycles lower (i.e. viral load at least about 8-fold higher) than adults aged 18 years or more, and this difference persists after adjustment for testing stream (drive-through versus on-site collection, order-set versus individually ordered test, patient class), sex and calendar day. The hypothesis is refuted if the adjusted difference in mean Ct for the under-5 versus 18-and-over contrast has a 95% confidence interval that excludes -3 cycles.

   Dataset. COVID-19 PCR testing at the Children's Hospital of Philadelphia, 2020 chop_covid_testing

   Test plan, novelty and risk

   Test plan. Dataset: chop_covid_testing (15,524 RT-PCR tests, pandemic days 4-107, 2020). Sample filters, applied in this order and each reported with the number of records lost: keep result == 'positive' (865 records; drop the 14,358 negatives, whose ct_result is the fixed non-detection ceiling of 45, and the 301 invalid results); drop records with missing ct_result (5) and with ct_result >= 45 (3), which are non-quantifiable; keep one record per subject_id, the earliest by pan_day and then by test_id, to remove within-person correlation from repeat testing. Expected analytic sample about 780-790 specimens from as many individuals. Exposure: age group from the `age` field, primary contrast 0-4 years versus 18 years and over, with 5-11 and 12-17 years carried as intermediate categories and a per-10-year linear trend as a secondary parameterisation. Outcome: ct_result in cycles, treated as continuous (lower Ct = higher viral load; roughly 3.3 cycles per log10 RNA copies). Covariates for confounding by testing indication and epidemic phase: drive_thru_ind (0/1), orderset (0/1), patient_class collapsed to inpatient / emergency / outpatient-or-observation / not-recorded, gender, and calendar time as a restricted set of pandemic-period indicators (roughly monthly blocks of pan_day). Estimators: (1) crude difference in mean Ct by age group with Welch (unequal-variance) 95% confidence intervals; (2) ordinary least squares of ct_result on age group and covariates with heteroscedasticity-robust (HC3) standard errors, reporting the adjusted mean difference and 95% CI for each age group against the 18-and-over reference; (3) direct standardisation of age-group mean Ct to the joint distribution of testing stream and calendar block in the whole positive sample, with nonparametric bootstrap percentile 95% CIs (2,000 resamples, fixed seed); (4) a dichotomous high-viral-load outcome (Ct < 25) analysed with a log-binomial or Poisson-with-robust-variance model to give a risk ratio and risk difference with 95% CIs; (5) median regression as a distribution-robust check. Stratified analysis: repeat the primary contrast within the drive-through stream and within the on-site stream separately, and report the two stratum-specific estimates alongside the crude one, so that any Simpson-type reversal is visible. Bias analysis: because the released records were de-identified by permutation, apply the classical result for non-differential misclassification of a dichotomous exposure to bound the direction of the attenuation, and report the factor by which the observed contrast would have to be attenuated to be compatible with a -4.5 cycle true difference. Refutation criteria: the hypothesis is refuted if the adjusted under-5 versus adult mean difference and its 95% CI lie entirely above -3 cycles; it is supported if the CI lies entirely below -3; it is inconclusive if the CI spans -3.

   Novelty. The literature contains one small paediatric-hospital study reporting a 4.5-cycle deficit in children under 5, and several larger community or multi-centre studies reporting no material difference, but no one has attempted the replication in a second paediatric hospital laboratory while making the confounding structure explicit. Crucially, no published analysis reports the crude and testing-stream-stratified contrasts side by side, even though in a children's hospital the age variable is almost perfectly confounded with reason for testing: children arrive through emergency and inpatient streams, adults are largely staff and caregivers swabbed at a drive-through site. The contribution is a precise interval rather than another 'no significant difference' statement, plus an explicit statement of the assumptions under which the contrast could be read causally.

   Expected outcome. The crude difference is expected to be small and in the opposite direction to the hypothesis, on the order of +1 cycle (children under 5 having slightly lower viral load than adults), with a 95% CI roughly one and a quarter cycles wide on each side, and the setting-stratified estimates are expected to differ from the crude estimate and possibly to change sign, because drive-through collection and order-set testing are strongly associated with both age and Ct. The most likely conclusion is that the hypothesis is refuted: the data will exclude a 3-cycle deficit in the under-5s, and by extension the 10- to 100-fold excess originally claimed, while being compatible with a difference of at most one to two cycles in either direction.

   Main risk. Days since symptom onset and symptom status are not recorded, so the age contrast remains confounded by where in the infection course each group is sampled; testing stream is only a proxy for that. The released data were de-identified by permutation, which for a dichotomous exposure biases toward the null and so cannot manufacture a positive point estimate where a large negative one is true, but does mean the interval should be read as a bound rather than a precise adjusted estimate. Finally, the sample is everyone tested at one hospital, so conditioning on a positive test can induce collider bias if both age and viral load affect the probability of being tested.

2. Idea 2

   Does the testing stream explain the age gradient in SARS-CoV-2 test positivity at a paediatric hospital?

   Hypothesis. The crude increase in RT-PCR test positivity with age observed across a paediatric hospital's 2020 testing programme is at least halved once testing stream (drive-through, order-set, patient class) and calendar period are adjusted for; that is, the adjusted odds ratio for positivity comparing adults aged 18 and over with children under 5 is less than half the crude odds ratio on the log scale. The hypothesis is refuted if the adjusted log odds ratio retains more than half of the crude log odds ratio.

   Dataset. COVID-19 PCR testing at the Children's Hospital of Philadelphia, 2020 chop_covid_testing

   Test plan, novelty and risk

   Test plan. Dataset: chop_covid_testing. Sample: all 15,223 tests with a valid result after excluding the 301 invalid results; report the count lost. Outcome: positive versus negative. Exposure: age group (0-4, 5-11, 12-17, 18-49, 50 and over). Estimator: logistic regression with cluster-robust standard errors by subject_id, first crude and then adjusted for drive_thru_ind, orderset, patient_class, gender and monthly blocks of pan_day; report crude and adjusted odds ratios with 95% CIs, plus risk differences from marginal standardisation of the fitted adjusted model with bootstrap CIs. Refutation: the adjusted log odds ratio for 18-and-over versus 0-4 remains above half the crude value.

   Novelty. Age gradients in test positivity are routinely reported from laboratory extracts and routinely read as gradients in infection risk, but positivity in a hospital testing programme is jointly determined by who is tested and why. Decomposing how much of the gradient survives adjustment for testing stream is a direct demonstration of indication-driven bias in a familiar surveillance statistic.

   Expected outcome. A crude gradient from roughly 4 percent positivity in children under 5 to roughly 9 percent in adults aged 50 and over, attenuating substantially but probably not by half after adjustment, so the hypothesis is likely to be refuted while still showing that a material part of the gradient is indication-driven.

   Main risk. Testing stream is measured coarsely, so residual confounding by indication is certain, and the analysis cannot separate genuine age differences in infection risk from age differences in exposure and in the threshold for being tested. The estimand is also not a population prevalence, which limits its interest.

3. Idea 3

   Confounding by the sick-quitter effect: how much of the crude association between smoking cessation and 10-year mortality is explained by baseline age and illness?

   Hypothesis. In the NHEFS cohort, the crude 10-year all-cause mortality risk in those who quit smoking exceeds that in continuing smokers by at least 5 percentage points, and adjustment for baseline age, sex, education, smoking history and prevalent disease removes at least three quarters of that crude risk difference. The hypothesis is refuted if the standardised risk difference retains more than a quarter of the crude risk difference.

   Dataset. NHANES I Epidemiologic Follow-up Study smoking cessation cohort nhefs_cohort

   Test plan, novelty and risk

   Test plan. Dataset: nhefs_cohort (1,629 smokers). Exposure: qsmk. Outcome: death by 1982. Estimator: crude risk difference and risk ratio with 95% CIs from a binomial model; then standardisation (parametric g-formula) over baseline age, sex, race, education, smokeintensity, smokeyrs, exercise, active and the disease-history indicators, with nonparametric bootstrap 95% CIs (fixed seed); inverse-probability-of-treatment weighted estimates as a cross-check; E-value for the standardised risk ratio. Refutation: the standardised risk difference retains more than a quarter of the crude value.

   Novelty. The sick-quitter phenomenon is folklore in tobacco epidemiology but is usually asserted rather than quantified as the fraction of a crude association removed by adjustment. Reporting the crude, standardised and IP-weighted risk differences side by side with an E-value gives an explicit accounting.

   Expected outcome. A crude excess mortality risk in quitters of roughly 5 to 10 percentage points, driven almost entirely by quitters being older and sicker, with the standardised risk difference close to zero and a wide interval.

   Main risk. The exposure-outcome contrast is heavily studied with this exact dataset and the direction of confounding is well known, so the study risks being a demonstration rather than a discovery. With about 300 deaths the standardised estimate will be imprecise, and reverse causation, illness causing cessation, is indistinguishable from confounding without information on the timing of quitting.

Ranking rationale

Idea 1 is ranked first because it takes a specific, published, quantified claim (a roughly 4.5-cycle viral-load deficit in children under 5) and subjects it to a genuine replication attempt in independent data of the same type, in a sample large enough that the resulting confidence interval can exclude the claimed effect rather than merely fail to detect it; it also has an internal methodological payoff, since testing stream is almost collinear with age in a children's hospital and the crude and stratified contrasts can be shown side by side. Idea 2 uses the same dataset but asks a weaker question, because a positivity gradient in a hospital convenience sample has no clean target parameter and residual confounding by indication cannot be bounded. Idea 3 is competently answerable and would produce a clean quantitative accounting, but the NHEFS smoking-cessation contrast is the single most analysed exposure-outcome pair in the causal-inference teaching literature, so the direction and rough magnitude of the confounding are already documented and the study would be a demonstration rather than a test.

Data provenance

Result (RESULT_JSON)

{
  "sample": {
    "rows_loaded": 15524,
    "result_negative": 14358,
    "result_invalid": 301,
    "result_positive": 865,
    "dropped_ct_missing": 5,
    "dropped_ct_at_ceiling_45": 3,
    "positive_quantifiable_tests": 857,
    "dropped_repeat_positive_tests": 71,
    "analytic_n": 786
  },
  "n_analysed": 786,
  "n_individuals": 786,
  "pan_day_min": 10,
  "pan_day_max": 107,
  "overall_mean_ct": 28.984796437659035,
  "overall_sd_ct": 7.2913358947784115,
  "overall_median_ct": 28.365,
  "descriptives_by_age_group": {
    "0-4": {
      "n": 186,
      "mean_ct": 29.188064516129035,
      "sd_ct": 8.048353299169191,
      "median_ct": 28.99,
      "p25_ct": 21.595,
      "p75_ct": 36.95,
      "pct_drive_thru": 0.3602150537634409,
      "pct_orderset": 0.6827956989247311,
      "pct_emergency_or_inpatient": 0.41935483870967744,
      "pct_male": 0.3709677419354839,
      "median_pan_day": 32,
      "n_hi_vl": 73,
      "risk_hi_vl": 0.3924731182795699
    },
    "5-11": {
      "n": 117,
      "mean_ct": 28.92119658119658,
      "sd_ct": 7.602818121159843,
      "median_ct": 28.19,
      "p25_ct": 22.18,
      "p75_ct": 36.2,
      "pct_drive_thru": 0.49572649572649574,
      "pct_orderset": 0.8290598290598291,
      "pct_emergency_or_inpatient": 0.3504273504273504,
      "pct_male": 0.5128205128205128,
      "median_pan_day": 53,
      "n_hi_vl": 44,
      "risk_hi_vl": 0.37606837606837606
    },
    "12-17": {
      "n": 168,
      "mean_ct": 29.358333333333334,
      "sd_ct": 6.909778861307767,
      "median_ct": 28.57,
      "p25_ct": 23.51,
      "p75_ct": 36.095,
      "pct_drive_thru": 0.4523809523809524,
      "pct_orderset": 0.8452380952380952,
      "pct_emergency_or_inpatient": 0.2916666666666667,
      "pct_male": 0.5297619047619048,
      "median_pan_day": 84.5,
      "n_hi_vl": 57,
      "risk_hi_vl": 0.3392857142857143
    },
    "18+": {
      "n": 315,
      "mean_ct": 28.689174603174603,
      "sd_ct": 6.91443868827219,
      "median_ct": 27.96,
      "p25_ct": 22.744999999999997,
      "p75_ct": 35.065,
      "pct_drive_thru": 0.7968253968253968,
      "pct_orderset": 0.8507936507936508,
      "pct_emergency_or_inpatient": 0.08888888888888889,
      "pct_male": 0.49523809523809526,
      "median_pan_day": 71,
      "n_hi_vl": 120,
      "risk_hi_vl": 0.380952380952381
    }
  },
  "crude_mean_ct_difference_vs_18plus": {
    "0-4": {
      "diff": 0.4988899129544322,
      "lcl": -0.8919687673791592,
      "ucl": 1.8897485932880236,
      "se": 0.7071308833828485
    },
    "5-11": {
      "diff": 0.23202197802197588,
      "lcl": -1.3530744420513237,
      "ucl": 1.8171183980952752,
      "se": 0.8036276920737634
    },
    "12-17": {
      "diff": 0.6691587301587347,
      "lcl": -0.6295771185357264,
      "ucl": 1.967894578853196,
      "se": 0.6602823342946575
    }
  },
  "adjusted_mean_ct_difference_vs_18plus_hc3": {
    "0-4": {
      "diff": -0.723154072664272,
      "lcl": -2.327816122348293,
      "ucl": 0.8815079770197491,
      "se": 0.8187201715650853
    },
    "5-11": {
      "diff": 0.06255102415590841,
      "lcl": -1.6605094407749872,
      "ucl": 1.785611489086804,
      "se": 0.8791286363025937
    },
    "12-17": {
      "diff": 0.7544298105256146,
      "lcl": -0.6384289409010113,
      "ucl": 2.147288561952241,
      "se": 0.7106552785731362
    }
  },
  "standardised_mean_ct_difference_vs_18plus_bootstrap": {
    "0-4": {
      "diff": -0.7231540726642685,
      "lcl": -2.2860773544572948,
      "ucl": 0.8043939292408588
    },
    "5-11": {
      "diff": 0.06255102415590841,
      "lcl": -1.7068562935186748,
      "ucl": 1.6454540614319817
    },
    "12-17": {
      "diff": 0.7544298105256182,
      "lcl": -0.5879355687201818,
      "ucl": 2.1062232909419794
    }
  },
  "adjusted_age_trend_per_10_years": {
    "diff_per_10y": 0.04969227199744353,
    "lcl": -0.3314207832601781,
    "ucl": 0.4308053272550651
  },
  "median_regression_difference_vs_18plus": {
    "0-4": {
      "diff": 0.06500099998226005,
      "lcl": -2.57461669981893,
      "ucl": 2.7046186997834503
    },
    "5-11": {
      "diff": -0.30000000011887096,
      "lcl": -3.198579120519065,
      "ucl": 2.5985791202813227
    },
    "12-17": {
      "diff": 0.5360094762145309,
      "lcl": -1.8910031514170311,
      "ucl": 2.963022103846093
    }
  },
  "high_viral_load_ct_lt_25_risk_ratio_vs_18plus": {
    "0-4": {
      "rr": 1.2246902208206631,
      "lcl": 0.9309446906756556,
      "ucl": 1.6111227143743636
    },
    "5-11": {
      "rr": 1.041027766121899,
      "lcl": 0.7712503449060434,
      "ucl": 1.4051712482135437
    },
    "12-17": {
      "rr": 0.8697152728174202,
      "lcl": 0.6674614426315144,
      "ucl": 1.1332559567631355
    }
  },
  "high_viral_load_ct_lt_25_risk_difference_vs_18plus": {
    "0-4": {
      "rd": 0.07900865055029616,
      "lcl": -0.026674555575315752,
      "ucl": 0.18822945518491585
    },
    "5-11": {
      "rd": 0.014329323038268171,
      "lcl": -0.09401361421521336,
      "ucl": 0.13230651811216815
    },
    "12-17": {
      "rd": -0.048459520367670206,
      "lcl": -0.14062457738971176,
      "ucl": 0.041621085693551454
    }
  },
  "stratified_by_testing_stream_0_4_vs_18plus": {
    "drive_thru": {
      "diff": -0.7839442231075715,
      "lcl": -2.8555335113529425,
      "ucl": 1.2876450651377995,
      "se": 1.0433874950548527,
      "n_0_4": 67,
      "n_18plus": 251,
      "mean_ct_0_4": 27.52,
      "mean_ct_18plus": 28.30394422310757
    },
    "on_site": {
      "diff": -0.07277310924370184,
      "lcl": -2.397294571398226,
      "ucl": 2.2517483529108224,
      "se": 1.1757821223100653,
      "n_0_4": 119,
      "n_18plus": 64,
      "mean_ct_0_4": 30.1272268907563,
      "mean_ct_18plus": 30.200000000000003
    }
  },
  "interaction_age0_4_x_drive_thru": {
    "diff": -1.4828106919184378,
    "lcl": -4.221608253872865,
    "ucl": 1.255986870035989
  },
  "sensitivity_all_positive_tests_cluster_robust": {
    "n_tests": 857,
    "n_subjects": 786,
    "diff": -0.1993512518837406,
    "lcl": -1.6376602095456312,
    "ucl": 1.23895770577815
  },
  "sensitivity_exclude_age_over_100": {
    "n": 785,
    "dropped": 1,
    "diff": -0.678600107562275,
    "lcl": -2.2827690551257738,
    "ucl": 0.925568840001224
  },
  "age_variable_internal_validity_check": {
    "positivity_by_age_group": {
      "0-4": {
        "n_tests": 5754,
        "n_positive": 220,
        "positivity": 0.03823427181091415
      },
      "5-11": {
        "n_tests": 2833,
        "n_positive": 128,
        "positivity": 0.04518178609248147
      },
      "12-17": {
        "n_tests": 2425,
        "n_positive": 175,
        "positivity": 0.07216494845360824
      },
      "18+": {
        "n_tests": 4211,
        "n_positive": 342,
        "positivity": 0.08121586321538826
      }
    },
    "adjusted_rr_positive_0_4_vs_18plus": 0.4931218787204185,
    "lcl": 0.4095220339925216,
    "ucl": 0.5937877991619755
  },
  "bias_analysis": {
    "attenuation_needed_for_true_minus4p5_at_point": 0.1607009050365049,
    "attenuation_needed_for_true_minus4p5_at_lcl": 0.517292471632954,
    "nondifferential_attenuation_could_reconcile_lcl": true,
    "bias_required_to_reach_minus3_cycles": 2.276845927335728,
    "u_ct_effect_needed_if_prevalence_gap_0p5": 4.553691854671456,
    "u_prevalence_gap_needed_if_ct_effect_5cycles": 0.4553691854671456
  },
  "fold_change_translation": {
    "adjusted_fold_higher_viral_load_0_4_vs_18plus": 1.6507871034748185,
    "lcl": 0.5427997731688575,
    "ucl": 5.020448046781037,
    "claimed_fold_in_literature": 22.62741699796952,
    "log10_copies_difference_point": 0.2191375977770521,
    "log10_copies_difference_lcl": -0.267123629399924,
    "log10_copies_difference_ucl": 0.7053988249540283
  },
  "hypothesis_threshold_cycles": -3,
  "claimed_literature_difference_cycles": -4.5,
  "hypothesis_refuted": true,
  "hypothesis_supported": false,
  "n_bootstrap": 2000,
  "seed": 20200401
}

"""Are RT-PCR cycle thresholds lower in children under 5 than in adults?
Real data: SARS-CoV-2 RT-PCR records, Children's Hospital of Philadelphia, 2020
(medicaldata::covid_testing via Rdatasets), loaded from disk - never simulated."""
import json
import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.regression.quantile_regression import QuantReg
from dataset import load

SEED, N_BOOT, HI_CUT = 20200401, 2000, 25.0
CLAIMED, REFUTE = -4.5, -3.0
AGE_BINS = [-0.001, 4.999, 11.999, 17.999, 1000.0]
LABS, REF = ["0-4", "5-11", "12-17", "18+"], "18+"
SETLEV = ["not recorded", "emergency", "inpatient", "outpatient/observation", "not applicable"]
BLKLEV = ["d04-33", "d34-63", "d64-93", "d94-107"]
Z = stats.norm.ppf(0.975)


def collapse(x):
    if pd.isna(x):
        return "not recorded"
    return {"inpatient": "inpatient", "emergency": "emergency",
            "not applicable": "not applicable"}.get(x, "outpatient/observation")


def prep(f):
    f = f.copy()
    f["ageg"] = pd.cut(f["age"], bins=AGE_BINS, labels=LABS).astype(str)
    f["setting"] = f["patient_class"].map(collapse)
    f["block"] = pd.cut(f["pan_day"], bins=[3, 33, 63, 93, 108], labels=BLKLEV).astype(str)
    return f


def design(f, age=True, interact=False):
    p = [pd.Series(1.0, index=f.index, name="const")]
    if age:
        p += [(f["ageg"] == a).astype(float).rename("age_" + a) for a in LABS if a != REF]
    p.append(f["drive_thru_ind"].astype(float).rename("drive_thru"))
    p.append(f["orderset"].astype(float).rename("orderset"))
    p += [(f["setting"] == s).astype(float).rename("set_" + s) for s in SETLEV[1:]]
    p.append((f["gender"] == "male").astype(float).rename("male"))
    p += [(f["block"] == b).astype(float).rename("blk_" + b) for b in BLKLEV[1:]]
    if interact:
        p.append(((f["ageg"] == "0-4") * f["drive_thru_ind"]).astype(float).rename("age04_x_dt"))
    X = pd.concat(p, axis=1)
    return X.values.astype(float), list(X.columns)


def ols(X, y, hc3=False):
    A = np.linalg.pinv(X.T @ X)
    b = A @ (X.T @ y)
    if not hc3:
        return b, None
    r = y - X @ b
    h = np.clip(np.einsum("ij,jk,ik->i", X, A, X), 0, 0.999)
    u = r / (1.0 - h)
    V = A @ (X.T @ (X * (u ** 2)[:, None])) @ A
    return b, np.sqrt(np.diag(V))


def logit_fit(X, y, iters=60):
    b = np.zeros(X.shape[1])
    for _ in range(iters):
        p = 1.0 / (1.0 + np.exp(-np.clip(X @ b, -30, 30)))
        w = np.clip(p * (1 - p), 1e-8, None)
        nb = b + np.linalg.pinv(X.T @ (X * w[:, None])) @ (X.T @ (y - p))
        if np.max(np.abs(nb - b)) < 1e-9:
            return nb
        b = nb
    return b


def poisson_robust(X, y, iters=80):
    b = np.zeros(X.shape[1])
    b[0] = np.log(max(y.mean(), 1e-6))
    for _ in range(iters):
        mu = np.exp(np.clip(X @ b, -30, 30))
        nb = b + np.linalg.pinv(X.T @ (X * mu[:, None])) @ (X.T @ (y - mu))
        if np.max(np.abs(nb - b)) < 1e-10:
            b = nb
            break
        b = nb
    mu = np.exp(np.clip(X @ b, -30, 30))
    A = np.linalg.pinv(X.T @ (X * mu[:, None]))
    V = A @ (X.T @ (X * ((y - mu) ** 2)[:, None])) @ A
    return b, np.sqrt(np.diag(V))


def gcomp(X, b, idx, link):
    Xc = X.copy()
    for j in idx.values():
        Xc[:, j] = 0.0
    f = (lambda M: M @ b) if link == "id" else (lambda M: 1 / (1 + np.exp(-np.clip(M @ b, -30, 30))))
    out = {REF: float(np.mean(f(Xc)))}
    for a, j in idx.items():
        Xg = Xc.copy()
        Xg[:, j] = 1.0
        out[a] = float(np.mean(f(Xg)))
    return out


def welch(a, b):
    va, vb = np.var(a, ddof=1) / len(a), np.var(b, ddof=1) / len(b)
    m, se = float(np.mean(a) - np.mean(b)), float(np.sqrt(va + vb))
    dfree = (va + vb) ** 2 / (va ** 2 / (len(a) - 1) + vb ** 2 / (len(b) - 1))
    t = stats.t.ppf(0.975, dfree)
    return {"diff": m, "lcl": m - t * se, "ucl": m + t * se, "se": se}


def fold(diff):
    return float(2.0 ** (-diff))


raw = load()
counts = {"rows_loaded": int(len(raw)),
          "result_negative": int((raw["result"] == "negative").sum()),
          "result_invalid": int((raw["result"] == "invalid").sum()),
          "result_positive": int((raw["result"] == "positive").sum())}
d = raw[raw["result"] == "positive"].copy()
counts["dropped_ct_missing"] = int(d["ct_result"].isna().sum())
d = d[d["ct_result"].notna()]
counts["dropped_ct_at_ceiling_45"] = int((d["ct_result"] >= 45).sum())
d = d[d["ct_result"] < 45]
counts["positive_quantifiable_tests"] = int(len(d))
d = d.sort_values(["pan_day", "subject_id"], kind="mergesort")
counts["dropped_repeat_positive_tests"] = int(len(d) - d["subject_id"].nunique())
d = prep(d.drop_duplicates("subject_id", keep="first").reset_index(drop=True))
n = len(d)
counts["analytic_n"] = int(n)
d["hi"] = (d["ct_result"] < HI_CUT).astype(float)
y = d["ct_result"].values.astype(float)
print("analysed", n, "specimens from", d["subject_id"].nunique(), "individuals")

desc = {}
for a in LABS:
    s = d[d["ageg"] == a]
    desc[a] = {"n": int(len(s)), "mean_ct": float(s["ct_result"].mean()),
               "sd_ct": float(s["ct_result"].std(ddof=1)), "median_ct": float(s["ct_result"].median()),
               "p25_ct": float(np.percentile(s["ct_result"], 25)),
               "p75_ct": float(np.percentile(s["ct_result"], 75)),
               "pct_drive_thru": float(s["drive_thru_ind"].mean()),
               "pct_orderset": float(s["orderset"].mean()),
               "pct_emergency_or_inpatient": float(s["setting"].isin(["emergency", "inpatient"]).mean()),
               "pct_male": float((s["gender"] == "male").mean()),
               "median_pan_day": float(s["pan_day"].median()),
               "n_hi_vl": int(s["hi"].sum()), "risk_hi_vl": float(s["hi"].mean())}

ref_ct = d.loc[d["ageg"] == REF, "ct_result"].values
crude = {a: welch(d.loc[d["ageg"] == a, "ct_result"].values, ref_ct) for a in LABS if a != REF}

X, nm = design(d)
idx = {a: nm.index("age_" + a) for a in LABS if a != REF}
b, se = ols(X, y, hc3=True)
adj = {a: {"diff": float(b[j]), "lcl": float(b[j] - Z * se[j]), "ucl": float(b[j] + Z * se[j]),
           "se": float(se[j])} for a, j in idx.items()}

keep = [nm.index(c) for c in nm if c.startswith(("drive", "orderset", "set_", "male", "blk_"))]
Xt = np.column_stack([np.ones(n), d["age"].values / 10.0, X[:, keep]])
bt, st = ols(Xt, y, hc3=True)
trend = {"diff_per_10y": float(bt[1]), "lcl": float(bt[1] - Z * st[1]), "ucl": float(bt[1] + Z * st[1])}

std_pt = gcomp(X, b, idx, "id")
bh = logit_fit(X, d["hi"].values)
risk_pt = gcomp(X, bh, idx, "logit")

rng = np.random.default_rng(SEED)
bs = {a: np.empty(N_BOOT) for a in idx}
brd = {a: np.empty(N_BOOT) for a in idx}
for k in range(N_BOOT):
    db = d.iloc[rng.integers(0, n, n)]
    Xb, _ = design(db)
    bb, _ = ols(Xb, db["ct_result"].values.astype(float))
    g = gcomp(Xb, bb, idx, "id")
    gr = gcomp(Xb, logit_fit(Xb, db["hi"].values), idx, "logit")
    for a in idx:
        bs[a][k] = g[a] - g[REF]
        brd[a][k] = gr[a] - gr[REF]
std_res = {a: {"diff": float(std_pt[a] - std_pt[REF]), "lcl": float(np.percentile(bs[a], 2.5)),
               "ucl": float(np.percentile(bs[a], 97.5))} for a in idx}
rd_res = {a: {"rd": float(risk_pt[a] - risk_pt[REF]), "lcl": float(np.percentile(brd[a], 2.5)),
              "ucl": float(np.percentile(brd[a], 97.5))} for a in idx}

bp, sp = poisson_robust(X, d["hi"].values)
rr = {a: {"rr": float(np.exp(bp[j])), "lcl": float(np.exp(bp[j] - Z * sp[j])),
          "ucl": float(np.exp(bp[j] + Z * sp[j]))} for a, j in idx.items()}

qr = QuantReg(y, X).fit(q=0.5)
ci = np.asarray(qr.conf_int())
medreg = {a: {"diff": float(qr.params[j]), "lcl": float(ci[j][0]), "ucl": float(ci[j][1])}
          for a, j in idx.items()}

strata = {}
for nmst, msk in (("drive_thru", d["drive_thru_ind"] == 1), ("on_site", d["drive_thru_ind"] == 0)):
    aa = d.loc[msk & (d["ageg"] == "0-4"), "ct_result"].values
    rr2 = d.loc[msk & (d["ageg"] == REF), "ct_result"].values
    w = welch(aa, rr2)
    w.update({"n_0_4": int(len(aa)), "n_18plus": int(len(rr2)),
              "mean_ct_0_4": float(np.mean(aa)), "mean_ct_18plus": float(np.mean(rr2))})
    strata[nmst] = w
Xi, ni = design(d, interact=True)
bi, si = ols(Xi, y, hc3=True)
ji = ni.index("age04_x_dt")
inter = {"diff": float(bi[ji]), "lcl": float(bi[ji] - Z * si[ji]), "ucl": float(bi[ji] + Z * si[ji])}

dall = prep(raw[(raw["result"] == "positive") & raw["ct_result"].notna() & (raw["ct_result"] < 45)])
Xa, na = design(dall)
ya = dall["ct_result"].values.astype(float)
ba, _ = ols(Xa, ya)
ra = ya - Xa @ ba
A = np.linalg.pinv(Xa.T @ Xa)
meat = np.zeros((Xa.shape[1], Xa.shape[1]))
for _, g in pd.DataFrame(Xa * ra[:, None]).groupby(dall["subject_id"].values):
    s = g.values.sum(axis=0)
    meat += np.outer(s, s)
Vc = A @ meat @ A
ja = na.index("age_0-4")
sens_all = {"n_tests": int(len(dall)), "n_subjects": int(dall["subject_id"].nunique()),
            "diff": float(ba[ja]), "lcl": float(ba[ja] - Z * np.sqrt(Vc[ja, ja])),
            "ucl": float(ba[ja] + Z * np.sqrt(Vc[ja, ja]))}

dp = d[d["age"] <= 100]
Xp, npm = design(dp)
bpp, spp = ols(Xp, dp["ct_result"].values.astype(float), hc3=True)
jp = npm.index("age_0-4")
sens_age = {"n": int(len(dp)), "dropped": int(n - len(dp)), "diff": float(bpp[jp]),
            "lcl": float(bpp[jp] - Z * spp[jp]), "ucl": float(bpp[jp] + Z * spp[jp])}

# internal validity of the de-identified age field: does age still predict positivity?
dv = prep(raw[raw["result"].isin(["positive", "negative"])].copy())
dv["pos"] = (dv["result"] == "positive").astype(float)
pos_by_age = {a: {"n_tests": int((dv["ageg"] == a).sum()),
                  "n_positive": int(dv.loc[dv["ageg"] == a, "pos"].sum()),
                  "positivity": float(dv.loc[dv["ageg"] == a, "pos"].mean())} for a in LABS}
Xv, nv = design(dv)
bv, sv = poisson_robust(Xv, dv["pos"].values)
jv = nv.index("age_0-4")
age_check = {"positivity_by_age_group": pos_by_age,
             "adjusted_rr_positive_0_4_vs_18plus": float(np.exp(bv[jv])),
             "lcl": float(np.exp(bv[jv] - Z * sv[jv])), "ucl": float(np.exp(bv[jv] + Z * sv[jv]))}

a04 = adj["0-4"]
gap = float(a04["diff"] - REFUTE)
bias = {"attenuation_needed_for_true_minus4p5_at_point": float(a04["diff"] / CLAIMED),
        "attenuation_needed_for_true_minus4p5_at_lcl": float(a04["lcl"] / CLAIMED),
        "nondifferential_attenuation_could_reconcile_lcl": bool(0.0 < a04["lcl"] / CLAIMED <= 1.0),
        "bias_required_to_reach_minus3_cycles": gap,
        "u_ct_effect_needed_if_prevalence_gap_0p5": float(gap / 0.5),
        "u_prevalence_gap_needed_if_ct_effect_5cycles": float(gap / 5.0)}

folds = {"adjusted_fold_higher_viral_load_0_4_vs_18plus": fold(a04["diff"]),
         "lcl": fold(a04["ucl"]), "ucl": fold(a04["lcl"]),
         "claimed_fold_in_literature": fold(CLAIMED),
         "log10_copies_difference_point": float(-a04["diff"] / 3.3),
         "log10_copies_difference_lcl": float(-a04["ucl"] / 3.3),
         "log10_copies_difference_ucl": float(-a04["lcl"] / 3.3)}

refuted = bool(a04["lcl"] > REFUTE)
supported = bool(a04["ucl"] <= REFUTE)
print("crude 0-4 vs 18+ =", round(crude["0-4"]["diff"], 3),
      "| adjusted =", round(a04["diff"], 3), (round(a04["lcl"], 3), round(a04["ucl"], 3)))

out = {"sample": counts, "n_analysed": int(n), "n_individuals": int(d["subject_id"].nunique()),
       "pan_day_min": int(d["pan_day"].min()), "pan_day_max": int(d["pan_day"].max()),
       "overall_mean_ct": float(np.mean(y)), "overall_sd_ct": float(np.std(y, ddof=1)),
       "overall_median_ct": float(np.median(y)),
       "descriptives_by_age_group": desc,
       "crude_mean_ct_difference_vs_18plus": crude,
       "adjusted_mean_ct_difference_vs_18plus_hc3": adj,
       "standardised_mean_ct_difference_vs_18plus_bootstrap": std_res,
       "adjusted_age_trend_per_10_years": trend,
       "median_regression_difference_vs_18plus": medreg,
       "high_viral_load_ct_lt_25_risk_ratio_vs_18plus": rr,
       "high_viral_load_ct_lt_25_risk_difference_vs_18plus": rd_res,
       "stratified_by_testing_stream_0_4_vs_18plus": strata,
       "interaction_age0_4_x_drive_thru": inter,
       "sensitivity_all_positive_tests_cluster_robust": sens_all,
       "sensitivity_exclude_age_over_100": sens_age,
       "age_variable_internal_validity_check": age_check,
       "bias_analysis": bias, "fold_change_translation": folds,
       "hypothesis_threshold_cycles": REFUTE, "claimed_literature_difference_cycles": CLAIMED,
       "hypothesis_refuted": refuted, "hypothesis_supported": supported,
       "n_bootstrap": int(N_BOOT), "seed": int(SEED)}
print("RESULT_JSON: " + json.dumps(out))

Results and analysis

Key numbers

Caveats

• Days since symptom onset and symptom status are not recorded. Because Ct is a single snapshot of a rising-then-falling trajectory, any age difference in how quickly people reach testing confounds the contrast, and testing stream is only a proxy for it. This is the largest single threat and it cannot be removed with these variables.
• The sample is everyone tested at one paediatric hospital and the analysis conditions on a positive result. If both age and viral load influence the probability of being tested and of testing positive, conditioning on positivity can induce a collider association between them; no design feature here breaks that path.
• Cycle threshold is not a viral load. It depends on assay, target gene, specimen adequacy and collection technique, and swab quality plausibly differs between a struggling toddler and a cooperating adult, so the measurement error may be differential by age in exactly the direction that matters.
• The released records were de-identified by time-shifting and permutation. The internal-validity check shows the age field retains a strong association with positivity, but the analysis cannot rule out some attenuation of the age-Ct contrast, and the estimate should be read as a bound rather than an exact adjusted difference.
• Precision is limited: 186 under-5 and 315 adult specimens give an interval about 1.6 cycles wide on each side of the point estimate. The study can exclude a 3-cycle deficit but cannot distinguish a difference of one cycle from none, and the high-viral-load risk ratio of 1.22 is compatible with anything between a 7% deficit and a 61% excess.
• Age groups are compared across different calendar periods (median pandemic day 32 for under-5s, 71 for adults, 84.5 for 12-17 year olds). Calendar blocks are adjusted for, but with four blocks the control for epidemic phase is coarse and residual confounding by time remains possible.
• Adults tested at this paediatric hospital are largely staff, caregivers and clients rather than a general adult patient population, so the adult reference group is not representative of adults in Philadelphia and the contrast may not transport to other settings.
• The results are from a single institution over a single 98-day period in 2020, before the emergence of variants that changed viral-load distributions, so they speak to the original claim in its own era and not to later variants.
• All estimates are associations conditional on presenting for testing. A causal reading would additionally require no unmeasured confounding by time since onset and severity given the measured covariates, no age-differential specimen-quality error, and no selection into the tested-and-positive sample depending jointly on age and viral load; none of these can be verified here.

≈ 2988 words · AI-drafted manuscript · not peer-reviewed

A careful, honest replication attempt that reaches a clear negative answer and is faithful to the data it claims to analyse. Using 786 first-positive specimens from deidentified Children's Hospital of Philadelphia RT-PCR records, the adjusted difference in mean cycle threshold between children under 5 and adults is -0.72 cycles (95% CI -2.33 to 0.88), which excludes the pre-set -3 cycle threshold and, comfortably, the -4.5 cycle difference behind the widely cited 10- to 100-fold claim. The paper's most useful contribution is not the null itself but the demonstration that the crude and adjusted contrasts have opposite signs because testing stream is nearly collinear with age in a hospital laboratory extract; that message is well supported by the stratified analysis. Every number I checked against the executed script and the recorded provenance matched, the data are correctly described as real and deidentified, and the reporting follows this journal's conventions on measures of association and restraint in causal language. The remaining problems are presentational and interpretive rather than analytic: the descriptive distribution of the outcome is too thin for a skewed measure, the viral-load translation is untabulated, the misclassification argument is stated more firmly than a permutation of unspecified fields warrants, and selection bias is diagnosed but not quantified alongside the confounding analysis. None of these requires new computation, and all can be handled in revision.