Tree-Based Scan Statistics (TreeScan)

Most pharmacoepidemiology asks a pre-specified question: does drug X raise the rate of outcome Y? Tree-based scan statistics turn that around. You do not name the outcome in advance.

You hand the method an entire hierarchical classification of outcomes - the ICD diagnosis tree, or the MedDRA System Organ Class -> High Level Group Term -> Preferred Term tree for adverse events - and let it scan all of the outcomes at once, asking at every level of granularity "is the count of events here, in the exposed group, higher than expected?"

The method's job is signal generation: it nominates the specific code (a single MedDRA Preferred Term) or the grouped branch (a whole System Organ Class) where events cluster, while honestly accounting for the fact that it just looked in thousands of places.

The core idea: cuts on a tree

A "cut" is a node of the tree taken together with every leaf beneath it. A leaf cut is one fine-grained code (e.g., the Preferred Term myocarditis); an internal-node cut is the whole branch (e.g., the cardiac System Organ Class, which sums myocarditis + arrhythmia + ...). Scanning the tree means evaluating every cut as a candidate signal.

This is what lets TreeScan catch a real effect whether it is concentrated in one code (sharp leaf signal) or smeared across a related family of codes that no single code makes significant on its own (a branch signal that only appears once you sum the siblings). A flat, one-code-at-a-time screen sees neither the grouping nor shares strength across related outcomes.

The statistic and the multiplicity fix

Each cut gets a log-likelihood-ratio (LLR) comparing the observed count under the cut to its expected count, with the alternative being "more events here than expected." The test statistic for the whole analysis is the maximum LLR over all cuts - the most surprising place on the tree. The hard part is that you tested thousands of overlapping, correlated cuts, so a raw per-cut p-value is meaningless.

TreeScan controls this with a Monte Carlo permutation null: it repeatedly re-generates the data under the no-signal hypothesis (e.g., randomly redistributing the events across the tree in proportion to the expected counts, or permuting exposed/unexposed labels), recomputes the maximum LLR each time, and builds the null distribution of that maximum.

A cut's signal is "significant" only if its observed LLR beats almost all of those simulated maxima. Because the threshold is set on the maximum over the whole tree, the family-wise error is controlled across the entire scan in one shot - no Bonferroni explosion, and the correlation among nested cuts is handled automatically.

Interpreting the output

From the worked example: six cuts evaluated over four MedDRA Preferred Terms under two System Organ Classes. The myocarditis leaf cut has the highest LLR (LLR = 4.07, observed 8 vs expected 3), exceeding the permutation-derived critical boundary of approximately 3.7 (p ≈ 0.017 by Monte Carlo). The cardiac System Organ Class branch cut (observed 10, expected 6) has LLR ≈ 2.10, below the boundary. No GI cut is flagged.

Formal interpretation: The myocarditis Preferred Term cut is the nominated signal — the single location in the outcome tree where the observed-to-expected excess is large enough that it would arise by chance in fewer than 5% of permutation replicates under the null of no excess anywhere on the tree.

The multiplicity correction is achieved through the permutation null on the maximum LLR: because the test statistic is the maximum over all cuts, the family-wise error rate across the entire tree is controlled at alpha without Bonferroni correction and without assuming independence among the overlapping, nested cuts.

The branch cut (cardiac SOC) was not flagged separately because the myocarditis leaf already captures most of the excess; the branch signal does not add independent information in this example.

Practical interpretation: The myocarditis finding is a generated hypothesis — a signal for follow-up investigation — not a confirmed causal relationship. Report the LLR, the permutation p-value, the observed and expected counts, and the tie fraction to the nearest-ancestor branch.

Next steps are: verify the myocarditis case definition PPV (a low-PPV code makes a noise signal), examine confounding (are the exposed patients differently characterized than the reference?), and run a pre-specified confirmatory study (a cohort or SCCS analysis on the myocarditis outcome alone) in an independent data source or time period.

TreeScan produces no effect-size estimate; do not extract a relative risk from the scan output.

Pros, cons, and trade-offs

(specific and comparative).

• vs signal-detection (disproportionality: PRR, ROR, the Bayesian BCPNN/MGPS): Disproportionality works on spontaneous-report databases (FAERS, VigiBase) and asks, pair by pair, whether a drug-event combination is reported more than expected given the report mix. It has no population denominator, no hierarchy, and treats every Preferred Term as an independent test. TreeScan runs on cohort or self-controlled data with real person-time/denominators, exploits the MedDRA/ICD hierarchy to borrow strength across related codes, and delivers one honest family-wise-controlled answer instead of thousands of separately-corrected pair tests. Prefer disproportionality when all you have is a spontaneous-report dump with no denominator; prefer TreeScan when you have an enumerated cohort (claims/EHR) where expected counts are estimable.
• vs a pre-specified single-outcome analysis (a Cox model or sequential MaxSPRT on outcome Y): A targeted analysis is more powerful for the outcome you named and is the right tool for confirmation. TreeScan trades some per-outcome power for breadth - it will surface the outcome you would not have thought to pre-specify. The cost is that a TreeScan hit is a lead, not a verdict: it must be re-tested in an independent, pre-specified, confounding-controlled study. Do not report a scan hit as a confirmed causal effect.
• vs naive multiple testing (scan every code, Bonferroni-correct): Bonferroni over thousands of correlated, nested codes is brutally conservative and ignores the hierarchy; it cannot express a branch-level signal at all. TreeScan's permutation-on-the-maximum is both less conservative (it respects the correlation) and able to test groupings. The price is computation (thousands of permutations x a full tree scan each).

When NOT to use - and when it is actively misleading

• Do not treat a signal as a confirmed effect. TreeScan is a generator. A flagged cut is hypothesis- generating; reporting it as an established harm conflates screening with confirmation and invites over-reaction to noise. Every hit needs an independent, pre-specified, confounder-adjusted follow-up.
• Garbage expected counts -> garbage signals. The whole method rests on credible expected counts (from a comparator cohort, historical rates, or the self-controlled comparison window). If the expected counts are confounded - the exposed and reference groups differ in age, comorbidity, or surveillance intensity - TreeScan will faithfully flag that confounding as a "signal." The unconditional Poisson model especially assumes the expected counts are correct; it does not adjust for confounding on its own. Use a self-controlled design, matching, or stratified expected counts, and pair scanning with negative-control outcomes / empirical calibration to gauge residual bias.
• It does not estimate effect size or do confounding control. TreeScan answers "where is there an excess?", not "how big, and is it causal?". Pulling a relative risk off a scan hit and acting on it - without a designed comparative study - is a misuse.
• Outcome misclassification rides along. The leaves are code-based outcome definitions; if a Preferred Term or ICD code has poor positive predictive value, its "signal" may be a coding artifact, not a real cluster.

Data-source operational depth

• Claims (FFS): The natural substrate - an enumerated cohort with denominators (person-time) and an ICD diagnosis tree. Build the exposed cohort, define the post-exposure risk window, count events at each ICD leaf, and get expected counts from a comparator cohort (active comparator or matched unexposed) or from the same people's comparison window (self-controlled). Watch Medicare Advantage / capitated person-time (incomplete encounter capture deflates counts) and the usual claims outcome-validity caveats - a noisy leaf code makes a noisy signal.
• EHR: Richer outcome detail (labs, vitals, problem lists) can sharpen leaf definitions and let you build a finer outcome tree, but encounter capture is leakier (care outside the system is invisible), so expected counts and risk windows are more fragile. Best when the network is reasonably closed.
• Registry: Adjudicated outcomes make clean leaves, but registries usually track a narrow outcome set, which defeats the point of scanning a broad tree; use when the registry itself is the surveillance target.
• Linked claims-EHR-registry / distributed networks (Sentinel): The ideal substrate - claims for denominators and exposure timing, EHR/registry for outcome refinement, and a common data model so the same tree and the same scan run across many sites. The self-controlled tree-temporal scan is the standard Sentinel tool here because it sidesteps between-person confounding by construction.

Worked surveillance example

A new vaccine is monitored across a claims network. Events in the post-vaccination risk window are counted at four MedDRA Preferred Terms grouped under two System Organ Classes; expected counts come from the matched comparison person-time. The scan evaluates six cuts (four leaves + two branches), scores each with the conditional-Poisson LLR, takes the maximum, and benchmarks it against the permutation null.

The worked_example below carries the exact, hand-checkable arithmetic - the myocarditis leaf wins with LLR 4.07 and a permutation p of about 0.017, nominating it as the candidate signal to confirm in a designed study.

Decision diagram

flowchart TD
  Root["Root (all events)<br/>observed 15 / expected 15"]
  Root --> Cardiac["Cardiac SOC (cut)<br/>obs 10 / exp 6<br/>LLR 2.17"]
  Root --> GI["GI SOC (cut)<br/>obs 5 / exp 9<br/>LLR 0"]
  Cardiac --> Myo["myocarditis (leaf cut)<br/>obs 8 / exp 3<br/>LLR 4.07 - max"]
  Cardiac --> Arr["arrhythmia (leaf cut)<br/>obs 2 / exp 3<br/>LLR 0"]
  GI --> Nau["nausea (leaf cut)<br/>obs 3 / exp 5<br/>LLR 0"]
  GI --> Vom["vomiting (leaf cut)<br/>obs 2 / exp 4<br/>LLR 0"]

flowchart LR
  A[Outcome hierarchy<br/>ICD / MedDRA tree] --> B[Count observed events<br/>at every leaf]
  B --> C[Expected counts<br/>comparator or self-control window]
  C --> D[Score every cut<br/>node + descendants by LLR]
  D --> E[Test statistic =<br/>maximum LLR over all cuts]
  E --> F{Beats the Monte Carlo<br/>permutation null?}
  F -- Yes --> G[Candidate signal<br/>hypothesis-generating]
  F -- No --> H[No signal at this cut]
  G --> I[Confirm in a pre-specified<br/>confounder-controlled study]

Worked example

Scenario

A new vaccine is being watched in a claims network. In the weeks after vaccination we count adverse events at four specific MedDRA codes (the leaves), grouped under two broad categories (the branches): a cardiac category holding myocarditis and arrhythmia, and a GI category holding nausea and vomiting. From a matched comparison group we already know how many of each event we would expect if the vaccine were harmless. We want the scan to tell us which code, or which whole category, has a real excess - and to do it honestly even though we are checking six places at once.

Dataset

The counts an analyst feeds the scan - observed events in the vaccinated group and expected events from the matched comparison, for each leaf and each branch (a branch is just the sum of its leaves).

FIG. 1 — DESIGN TIMELINE

Steps

Result

The scan nominates the myocarditis leaf as the candidate signal (maximum LLR = 4.07, permutation p = 0.017); the cardiac branch is elevated (LLR 2.17) only because it contains myocarditis, and the GI side shows nothing. This is a hypothesis to confirm in a pre-specified, confounder-controlled study - not a confirmed harm.

Trade-offs

Runnable example

import math, random

LEAVES   = ["myocarditis", "arrhythmia", "nausea", "vomiting"]
OBSERVED = {"myocarditis": 8, "arrhythmia": 2, "nausea": 3, "vomiting": 2}   # sum = 15
EXPECTED = {"myocarditis": 3, "arrhythmia": 3, "nausea": 5, "vomiting": 4}   # conditioned to sum = 15
CUTS = {                                  # node -> leaves beneath it (leaf cuts + branch cuts)
    "myocarditis": ["myocarditis"],
    "arrhythmia":  ["arrhythmia"],
    "nausea":      ["nausea"],
    "vomiting":    ["vomiting"],
    "cardiac_SOC": ["myocarditis", "arrhythmia"],
    "gi_SOC":      ["nausea", "vomiting"],
}

def cut_llr(c, mu, C):
    # conditional-Poisson LLR for one cut; 0 unless the cut is in excess (c > mu)
    if c <= mu or c >= C:
        return 0.0
    return c * math.log(c / mu) + (C - c) * math.log((C - c) / (C - mu))

def scan(obs):
    C = sum(obs.values())
    best_node, best_llr = None, 0.0
    for node, members in CUTS.items():
        c  = sum(obs[m] for m in members)
        mu = sum(EXPECTED[m] for m in members)
        s  = cut_llr(c, mu, C)
        if s > best_llr:
            best_llr, best_node = s, node
    return best_llr, best_node

obs_llr, node = scan(OBSERVED)          # -> ~4.07 at "myocarditis"

C      = sum(OBSERVED.values())
weights = [EXPECTED[l] for l in LEAVES]
random.seed(20230601)
NSIM, hits = 999, 0
for _ in range(NSIM):
    draw = random.choices(LEAVES, weights=weights, k=C)   # redistribute C events under the null
    sim  = {l: 0 for l in LEAVES}
    for l in draw:
        sim[l] += 1
    if scan(sim)[0] >= obs_llr - 1e-9:
        hits += 1
p = (hits + 1) / (NSIM + 1)
print(f"most likely cut = {node}  LLR = {obs_llr:.2f}  permutation p = {p:.3f}")

leaves   <- c("myocarditis", "arrhythmia", "nausea", "vomiting")
observed <- c(myocarditis = 8, arrhythmia = 2, nausea = 3, vomiting = 2)   # sum = 15
expected <- c(myocarditis = 3, arrhythmia = 3, nausea = 5, vomiting = 4)   # conditioned to sum = 15
cuts <- list(myocarditis = "myocarditis", arrhythmia = "arrhythmia",
             nausea = "nausea", vomiting = "vomiting",
             cardiac_SOC = c("myocarditis", "arrhythmia"),
             gi_SOC      = c("nausea", "vomiting"))

cut_llr <- function(c, mu, C) {
  if (c <= mu || c >= C) return(0)
  c * log(c / mu) + (C - c) * log((C - c) / (C - mu))
}
scan <- function(obs) {
  C <- sum(obs); best_llr <- 0; best_node <- NA_character_
  for (node in names(cuts)) {
    m  <- cuts[[node]]
    s  <- cut_llr(sum(obs[m]), sum(expected[m]), C)
    if (s > best_llr) { best_llr <- s; best_node <- node }
  }
  list(llr = best_llr, node = best_node)
}
obs <- scan(observed)                       # -> ~4.07 at "myocarditis"

C       <- sum(observed)
weights <- expected[leaves] / sum(expected)
set.seed(20230601); NSIM <- 999L; hits <- 0L
for (i in seq_len(NSIM)) {
  draw <- sample(leaves, C, replace = TRUE, prob = weights)
  sim  <- setNames(as.integer(table(factor(draw, levels = leaves))), leaves)
  if (scan(sim)$llr >= obs$llr - 1e-9) hits <- hits + 1L
}
p <- (hits + 1L) / (NSIM + 1L)
cat(sprintf("most likely cut = %s  LLR = %.2f  permutation p = %.3f\n", obs$node, obs$llr, p))

%let C    = 15;     /* total observed = total expected (conditioned)  */
%let nsim = 999;

/* 1) Score the six candidate cuts (leaf cuts + branch-node cuts). */
data cuts;
  length node $12;
  input node $ c mu;
  if c > mu and c < &C
    then llr = c*log(c/mu) + (&C - c)*log((&C - c)/(&C - mu));
    else llr = 0;
  datalines;
myocarditis 8 3
arrhythmia  2 3
nausea      3 5
vomiting    2 4
cardiac_SOC 10 6
gi_SOC      5 9
;
run;

proc sort data=cuts; by descending llr; run;
data _null_;
  set cuts;
  if _n_ = 1 then do;
    call symputx('best_node', node);
    call symputx('best_llr', put(llr, 8.4));
  end;
run;

/* 2) Monte Carlo permutation null: redistribute C events across the 4 leaves
      in proportion to expected counts (3,3,5,4)/15, recompute the maximum LLR. */
data perm;
  call streaminit(20230601);
  array p[4] _temporary_ (0.2 0.2 0.3333333 0.2666667);
  array mu[6] _temporary_ (3 3 5 4 6 9);
  do sim = 1 to &nsim;
    array cnt[4] c1-c4;
    do k = 1 to 4; cnt[k] = 0; end;
    do e = 1 to &C;
      u = rand('TABLE', of p[*]);     /* category 1..4 */
      cnt[u] + 1;
    end;
    array cc[6] cc1-cc6;
    cc1 = c1; cc2 = c2; cc3 = c3; cc4 = c4;   /* leaf cuts   */
    cc5 = c1 + c2; cc6 = c3 + c4;             /* branch cuts */
    maxllr = 0;
    do j = 1 to 6;
      if cc[j] > mu[j] and cc[j] < &C
        then L = cc[j]*log(cc[j]/mu[j]) + (&C - cc[j])*log((&C - cc[j])/(&C - mu[j]));
        else L = 0;
      if L > maxllr then maxllr = L;
    end;
    hit = (maxllr >= &best_llr - 1e-9);
    output;
  end;
  keep sim maxllr hit;
run;

proc sql noprint;
  select (sum(hit) + 1) / (&nsim + 1) into :pval trimmed from perm;
quit;
%put NOTE: most likely cut = &best_node  observed LLR = &best_llr  permutation p = &pval;

Citations

1. [1]Kulldorff M, Fang Z, Walsh SJ. A tree-based scan statistic for database disease surveillance. Biometrics. 2003;59(2):323-331.

1. [2]Russo MW, Wang SV. An open-source implementation of tree-based scan statistics. Pharmacoepidemiology and Drug Safety. 2024;33(3):e5765.

1. [3]Yih WK, Kulldorff M, Dashevsky I, Maro JC. Using the self-controlled tree-temporal scan statistic to assess the safety of live attenuated herpes zoster vaccine. American Journal of Epidemiology. 2019;188(7):1383-1388.