All your base R belong to us: small data and expansive analytics in base R

Chad Heilig — R Users Group, 2019-10-24 — 39 slides (38 shown; 1 omitted)

This rendering was constructed from an original PowerPoint presentation. Corrections of original text are indicated as “[corrected: X]”, redactions as “[redacted: X]”, and names withheld as “[generalized: X]”. Handouts and other auxiliary sources mentioned in the presentations are not available. The content of each presentation was believed to be correct and valid at the time of the presentation. This content, including references and links within the slides, has not been checked or updated when preparing this site. This rendering was facilitated by Anthropic Claude agentic tools for automated text extraction through computer vision and for assistance with compiling the completed sets.

Setting: All your base R belong to us: small data and expansive analytics in base R, 2019-10-24. CDC’s R Users Group.

This presentation set up expansive analysis of 2 small datasets from studies in western Kenya. In a cluster-randomized trial, we allocated 6 population-level interventions to 24 geographic areas and compared findings before and after an intervention period. The exact method from these small numbers worked from 262,440,000 possible assignment permutations. In an observational diagnostic study of 32 children with laboratory evidence of tuberculosis (TB), we drew 10,000 pseudosamples of 729 specimen combinations to approximate the multivariate performance of 21,870,000 pseudosample combinations. I performed the entire analysis in base R using careful subscripting. I could perhaps have used other packages, or these days Python, but at the time it was intellectually satisfying to work with efficiency, transparency, and a lack of library dependencies while exercising granular control.

Slide 1All your base R belong to us: small data and expansive analytics in base R

Blue title banner with the HHS and CDC logos at right; a redacted region follows the bullet after the presenter's name.

Chad Heilig, PhD • [redacted: email address]
Associate Director for Data Science
Machine Intelligence and Data Science Team
R Users Group, 2019-10-24

Description

Text in the picture

Center for Surveillance, Epidemiology, and Laboratory Services

Abbreviations: CDC = Centers for Disease Control and Prevention; HHS = Department of Health and Human Services.

Slide 2All your base are belong to us

Pixel-art video game screen: a pale figure with teal hair, a partly mechanical face and a purple cloak on orange, above a caption bar.

https://en.wikipedia.org/wiki/All_your_base_are_belong_to_us

Description

Text in the picture

CATS: ALL YOUR BASE R BELONG
TO US.

Presenter notes

https://en.wikipedia.org/wiki/All_your_base_are_belong_to_us
https://knowyourmeme.com/memes/all-your-base-are-belong-to-us
http://www.allyourbasearebelongtous.com/history/

Slide 3Small data, expansive analyses, base R

Small-data, expansive analyses
- Cluster-randomized trial, 24 geographic areas, western Kenya
  - Compare 6 population interventions before and after intervention period
  - Exact method: 262,440,000 possible assignment permutations
- Observational diagnostic study, 32 children with TB, western Kenya
  - Compare diagnostic yield of 6 specimen types, 3 test types
  - Multivariate performance: 21,870,000 pseudosample/combinations
Base R (no additional packages required*)
- Careful object subscripting
  - For analysis (mostly arrays) and for presentation (base graphics)
- Why base R? Efficiency, ease** (transparent, no dependencies)

* Limited data.table use not required. ** “Ease” is subject to interpretation.

Abbreviations: TB = tuberculosis.

Slide 4Cluster-randomized trial

24 geographic areas in western Kenya
Compare methods for boosting tuberculosis case notifications

Slide 5Cluster-randomized trial

Study design and purpose
- 24 geographic areas in western Kenya, 6 population interventions
- Compare methods for boosting tuberculosis case notifications
  - Performance measure: number/rate of TB case notifications
  - Contrast before and after intervention period, between intervention types
Analytic design
- Exact method: enumerate reference distributions for performance measures; 262,440,000 possible assignment permutations
- Sample* 1,000,000 assignments to verify exact method
- Attribute case notifications to (re)allocated conditions

* without replacement

Abbreviations: TB = tuberculosis.

Slide 6Interventions to boost case notification

Incomplete factorial design
- No facility-based
  - Standard (S)
  - Contact investigation (C)
  - Mobile unit (M)
- Health facility-based
  - Health-facility based (H)
  - H + Contact investigation (HC)
  - HC + Mobile unit (HCM)

3 singular interventions
- Contact investigation (C)
- Mobile unit (M)
- Health-facility based (H)
A complete 2^3 factorial design would have 8 conditions:
- ∅, C, M, H, ~~CM~~, HC, ~~HM~~, HCM [CM and HM struck through]

Abbreviations: HC = health-facility enhancement plus contact investigation (trial condition H + C); HCM = health-facility enhancement plus contact investigation plus mobile unit (trial condition H + C + M); HM = health-facility enhancement plus mobile unit (trial condition H + M).

Slide 7Geographic study units

24 units, aimed for average of 25,000 persons (2009 census)
Rural units
- 6 in the Siaya DSS* area, Siaya County (labeled S1–S6)
- 6 in the Kombewa DSS area, Kisumu County (W1–W6)
Urban units
- 6 in the western portion of Kisumu Town (L1–L6)
- 6 in the eastern portion of Kisumu Town (M1–M6)

* Demographic Surveillance System

Abbreviations: DSS = demographic surveillance system.

Slide 8Random allocation to geographic units

Rural units, Siaya and Kombewa DSS
- Assign {S, C, M} to a set of 6 and {H, HC, HCM} to the other 6
  - Health-facility conditions with probability 0.5; 2 possibilities
- For 6 units assigned {S, C, M}, permute (S, S, C, C, M, M); 90*
- For other 6 units, permute (H, H, HC, HC, HCM, HCM); 90
- 2 × 90 × 90 = 16,200 (constrained) permutations
Urban units, Kisumu Town west and east—same scheme
- 16,200 (constrained) permutations
Overall, 262,440,000 possible assignment allocations

* C(6, 2) × C(4, 2) × C(2, 2) = 15 × 6 = 90; C(n, m) = n! / [m! (n − m)!]

Presenter notes

In the 6 units with enhanced health-facility conditions, there were 90 ways to assign 3 conditions in pairs: C(6, 2), or 15, ways to assign 1 pair among the 6 units; C(4, 2), or 6, ways to assign another pair among the remaining 4 units, and C(2, 2), or 1, way to assign the third pair among the remaining 2 units, giving C(6, 2) × C(4, 2) × C(2, 2) = 15 × 6 = 90. Similarly, in the remaining 6 rural units, there were 90 ways to assign 3 pairs of interventions that did not involve enhanced case-finding at health facilities. These numbers 2 × 90 × 90 yielded 16,200 possible assignments, each with equal probability. C(6, 2) is notation for the number of way to choose 2 GSUs a set of 6 GSUs, also known as a binomial coefficient.
In factorial notation, C(n, m) = n! / [m! (n − m)!].

Abbreviations: DSS = demographic surveillance system; GSUs = geographic study units; HC = health-facility enhancement plus contact investigation (trial condition H + C); HCM = health-facility enhancement plus contact investigation plus mobile unit (trial condition H + C + M).

Slide 9Random allocation in R

set.seed(19631212) # ISO 8601 date of Kenyan independence
rand_hf <- c(sample(1:2), sample(3:4)) # health facility vs not
# [1] 1 2 3 4 # randomly left with identity permutation!
rand_un <- apply(outer(c("L", "M", "S", "W"), 1:6, paste0), 1, sample)
# 6 rows by 4 columns, random row order within each column
# rand_un[1,] # "L2" "M5" "S1" "W1"
# rand_un[,1] # "L2" "L1" "L3" "L4" "L6" "L5"
# data.frame with 24 rows, columns "unit" (random order) and "condition" (fixed)
rand_all <- data.frame(unit=c(rand_un[, rand_hf]), condition=rep(rep(c(
"0 standard", "1 contact investigation", "2 mobile unit",
"3 enhanced health facility", "4 enhanced HF, CI, MU", "5 enhanced HF, CI"),
each=2), times=2), stringsAsFactors=FALSE)

Abbreviations: CI = contact investigation; HF = health facility; ISO = International Organization for Standardization; MU = mobile unit.

Slide 10Effects of interest

Area-specific case-notification rate =
1e5 * (annualized # TB cases registered in designated period) / (# persons in 2009 census)
Δ rate = (rate after intervention) − (rate before intervention)
1st-order effects: Δ rate for S, C, M, H, HC, HCM
2nd-order effects, same facility condition:
C – S, M – S, M – C, HC – H, HCM – H, HCM – HC
2nd-order effects, mixed facility conditions:
H – S, HC – S, HCM – S, …, H – M, HC – M, HCM – M

Slide 11Observed effects

Model terms: study unit, intervention assigned, 2009 census,
pre-intervention cases, post-intervention cases
24 observations, 1 for each study unit
Calculate rates and differences of rates, combining unit values in straightforward way
But how can we know whether observed values are more extreme than might be expected by chance?
Exact analysis, randomization-based reference distribution

Slide 12Exact randomization-based analysis (1)

Fisher’s exact test is a randomization-based analysis
- Consider observed result in a contingency table
- Fix observed values; consider all the ways of permuting labels
  - Called a reference distribution
- Compare observed value to reference distribution
  - Is observed value extreme relative to reference distribution?
  - If yes, then “significant”; if no, then “nonsignificant”.
- See the “lady tasting tea” https://en.wikipedia.org/wiki/Lady_tasting_tea
  - Reference distribution of 70 values
  - Observed value exceeded 69 of the values of the reference

Slide 13Exact randomization-based analysis (2)

Recall 262,440,000 possible allocations
- We could compute all of them, but we don’t have to
- Reference distribution is much smaller for most statistics of interest
  - 1st-order effects (S, C, M, H, HC, HCM): 900 values, single ref dist’n
  - 2nd-order effects, same facility condition (C – S, …, HCM – HC): 32,400
  - 2nd-order effects, mixed facility conditions (H – S, …, HCM – M): 202,500
- Build up reference distributions from smaller sets

Slide 14Building up 1st-order effect distribution

i1 <- rep(1:5, 5:1)
i2 <- as.integer(c(2:6, 3:6, 4:6, 5:6, 6))
i15e <- rep(1:15, each=15); i15t <- rep(1:15, times=15)
# raw_diff is 4 × 6 matrix of raw differences in annualized case counts
diff1a <- raw_diff[, i1] + raw_diff[, i2] # 4 × 15
diff1b <- diff1a[c(1,1,2,2), i15e] + diff1a[c(3,4,3,4), i15t] # 4 × 225 -> 900

Slide 15Approximate randomization analysis

Use built-up, full reference distribution of 262,440,000
Subsample 1,000,000 assignments as a check on first method
(1e6 might be overkill, esp considering cardinality of exact distributions)
- Calculate effects of interest over sampled assignments
- Compare to efficiently enumerated distributions

Presenter notes

https://en.wikipedia.org/wiki/Truncatable_prime; https://oeis.org/A024770

Slide 16Approximate randomization analysis in R

set.seed(73939133) # largest right-truncatable prime, oeis.org/A024770
nsim <- 250000
s1 <- sort(sample(65610000, nsim)) # replace=FALSE; also s2, s3, s4
i8100e <- rep(1:8100, each=8100); i8100t <- rep(1:8100, times=8100)
i6 <- c(4:6, 1:3)
lsmw_samp_diff <- rbind(ls_diff[, i8100e[s1]], mw_diff[, i8100t[s1]])
samp_diff <- cbind(lsmw_samp_diff, lwms_samp_diff, mwls_samp_diff, mslw_samp_diff)

Presenter notes

https://en.wikipedia.org/wiki/Truncatable_prime; https://oeis.org/A024770

Slide 17Reference distribution for 1st-order effects

Density plot of the reference distribution for 1st-order effects, exact and sampled, with effects marked along the bottom.

HCM effect:
∆ rate = −26.2/1e5, quantile = 0.918

From the pre-intervention to the post-intervention period, case notification rates went down all around.

Description

Panels: 1
Elements: 2 density curves (legend: exact, thick black line; sampled, thin gray line) depict density over rate difference; 6 rotated orange text labels stand just above the horizontal axis, each naming an effect with a value in parentheses; a yellow leader line runs from the HCM label to the callout text at right of the plot
Horizontal: Rate difference per 100,000; ticks from -200 to 0 in steps of 50; curves extend from about -225 to about 5, slightly beyond the outer ticks
Vertical: Density; ticks from 0.000 to 0.012 in steps of 0.002; the sampled curve's highest spike reaches about 0.013, just above the top tick
Visible pattern: Both curves rise from near zero at the far left, peak between about -85 and -70, and fall back to near zero at the right end, just past 0. The left tail is long and gradual, the right side shorter and steeper. The exact curve is fairly smooth and peaks near 0.0096; the sampled curve tracks it closely but is jagged throughout, with its largest swings near the peak. Four labels cluster between about -120 and -95, on the rising left side, and 2 stand apart on the falling right side, near -40 and -25.

Text in the picture

exact
sampled
M (0.186)
H (0.199)
S (0.281)
HC (0.369)
C (0.854)
HCM (0.918)

Slide 18Reference distribution for 2nd-order effects, same facility condition

Density plot of the reference distribution for 2nd-order effects, same facility condition, exact and sampled, with effects marked.

HCM−H effect:
∆ rate = 91.6/1e5, quantile = 0.897

This reference distribution is necessarily symmetric.

Description

Panels: 1
Elements: 2 density curves (legend: exact, thick black line; sampled, thin gray line) depict density over rate difference; 6 rotated orange text labels stand just above the horizontal axis, each naming an effect with a value in parentheses, 2 of them overlapping; a yellow leader line runs from the HCM-H label to the callout text at right of the plot
Horizontal: Rate difference per 100,000; ticks from -200 to 200 in steps of 100; curves extend from about -225 to about 225, slightly beyond the outer ticks
Vertical: Density; ticks from 0.000 to 0.005 in steps of 0.001; curves peak at about 0.0054 (exact) and 0.0057 (sampled)
Visible pattern: Both curves are bell-shaped and symmetric about 0, rising from near zero at about -225 to a single peak at 0 and falling to near zero at about 225. They nearly coincide in the tails and shoulders; near the peak, from about -50 to 50, the sampled curve runs slightly above the exact one and is slightly ragged. One label on the left near -85 and one near -15; the other 4 bunch on the right between about 20 and 95, where 2 of them overlap.

Text in the picture

exact
sampled
M-C (0.126)
M-S (0.419)
HC-H (0.624)
[illegible]
HCM-H (0.897)

Slide 19Reference distribution for 2nd-order effects, mixed facility conditions

Density plot of the reference distribution for 2nd-order effects, mixed facility conditions, exact and sampled, with effects marked.

HCM−M effect:
∆ rate = 95.0/1e5, quantile = 0.933

This reference distribution is necessarily symmetric.

Description

Panels: 1
Elements: 2 density curves (legend: exact, thick black line; sampled, thin gray line) depict density over rate difference; rotated orange text labels stand just above the horizontal axis, each naming an effect with a value in parentheses, several overlapping near the center; a yellow leader line runs from the HCM-M label to the callout text at right of the plot
Horizontal: Rate difference per 100,000; ticks from -200 to 200 in steps of 100; curves extend from about -220 to about 220, slightly beyond the outer ticks
Vertical: Density; ticks from 0.000 to 0.006 in steps of 0.001; curves peak at about 0.0064 (exact) and 0.0062 (sampled), just above the top tick
Visible pattern: Both curves are bell-shaped and mirror each other about 0, rising from near zero at about -220 to a single peak at 0 and falling to near zero at about 220. The exact curve runs slightly above the sampled one around the peak, from about -75 to 75, and the 2 nearly coincide in the tails. Two labels on the left near -85 and -60, a cluster of overlapping labels sits between about -15 and 25, and 2 labels on the right near 80 and 95.

Text in the picture

exact
sampled
H-C (0.106)
HC-C (0.190)
H-S (0.430)
[illegible]
HC-M (0.656)
HCM-S (0.895)
HCM-M (0.933)

Slide 20Extreme values and the reference dist’n

1st-order effects:
only observed HCM effect is “extreme” relative to reference
2nd-order, same facility:
none is extreme
2nd-order, mixed facility:
only HCM – M is extreme

	Δ rate	F
First-order effects
S	-106.2	0.281
C	-38.1	0.854
M	-121.1	0.186
H	-117.8	0.199
HC	-94.8	0.369
HCM	-26.2	0.918↑
Second-order effects, same facility condition
C - S	68.1	0.826
M - S	-14.9	0.419
M - C	-83.0	0.126
HC - H	22.9	0.624
HCM - H	91.6	0.897
HCM - HC	68.6	0.828
Second-order effects, mixed facility conditions
H - S	-11.6	0.430
HC - S	11.4	0.569
HCM - S	80.0	0.895
H - C	-79.6	0.106
HC - C	-56.7	0.190
HCM - C	11.9	0.573
H - M	3.4	0.520
HC - M	26.3	0.656
HCM - M	95.0	0.933↑

Slide 21Observational study of TB diagnostics

32 children with TB in western Kenya
Compare performance of 6 specimen types, 3 test methods

Abbreviations: TB = tuberculosis.

Slide 22Observational study of diagnostics

Study design and purpose
- 32 children with TB on at least 1 of 8 different specimen types
- Compare Dx performance of 6 specimen types, 3 test methods
  - Performance measure: number of children with lab-confirmed TB
  - Compare # with TB for combinations of specimen and test types
Analytic design
- Resample* 31 children with TB; construct 10,000 pseudosamples
  - Each pseudosample: yield of 729 specimen combinations, 3 tests
  - Preserve multivariate dependencies between specimen types
- Compare combination performance head to head [21,870,000]

* with replacement

Abbreviations: Dx = diagnosis; TB = tuberculosis.

Slide 23Pediatric TB diagnostic results

Dot plot of test results for 32 numbered children across 6 specimen types, rows grouped by HIV status, at the right of the slide.

Specimen types: gastric aspirate (GA), nasopharyngeal aspirate (NPA), induced sputum (IS), string (ST), stool, urine
Test types: nucleic acid amplification (Xpert),
broth culture (MGIT), union of Xpert and MGIT
Results presented in dot plot
- White = negative, black = positive, no dot = no result
- Child in row 32 excluded from further analysis

Exclude blood, lymph node aspirate specimens, fluorescence microscopy results

Description

Panels: 1
Elements: Dot plot of 32 numbered rows by 6 specimen-type columns; each cell holds up to 2 circles, filled black or open white, or is empty; rows in 2 bracketed groups labeled HIV+ (rows 1 to 7) and HIV– (rows 8 to 32); a fraction label at the right end of each row, ranging from 6/6 to 1/1
Horizontal: 6 specimen-type columns, labeled vertically below as GA, NPA, IS, ST, stool, urine, bracketed together under the label Xpert or MGIT; up to 2 dots side by side per column
Vertical: 32 rows numbered 1 to 32 from top to bottom, with a gap between rows 7 and 8 separating the 2 HIV groups
Visible pattern: Within each HIV group the rows run from mostly black at the top to mostly white at the bottom, and the fraction labels at the right descend in step, from 6/6 near the top of each group to 1/6 near its bottom. Black dots are densest in the GA, NPA and IS columns; the ST, stool and urine columns are mostly white below the first several rows of each group. Some cells hold a single dot rather than a pair, and some are empty. The last 2 rows carry the fewest dots, the last only a single black dot, under NPA.

Abbreviations: GA = gastric aspirate; HIV = human immunodeficiency virus; HIV+ = HIV-positive; IS = induced sputum; MGIT = Mycobacteria Growth Indicator Tube (liquid culture system); NPA = nasopharyngeal aspirate; ST = string test (specimen collection method); TB = tuberculosis; Xpert = Xpert MTB/RIF (Cepheid GeneXpert cartridge nucleic acid amplification test for TB and rifampicin resistance).

Slide 24Pediatric TB diagnostic results, expanded

Full-width dot plot of test results for 32 numbered children across 4 test groups and up to 8 specimen types, rows grouped by HIV status.

Description

Panels: 1
Elements: Dot plot of 32 numbered rows by 27 test-and-specimen columns in 4 bracketed column groups; each cell holds up to 2 circles, filled black (positive), open white (negative), or small filled gray (contaminated or indeterminate), or is empty (no result); rows in 2 bracketed groups labeled HIV+ (rows 1 to 7) and HIV– (rows 8 to 32); a fraction label at the right end of each row, ranging from 6/6 to 1/1
Horizontal: 4 column groups, left to right: FM (GA, NPA, IS, ST, stool, urine); Xpert (GA, NPA, IS, ST, stool, urine, FNA); MGIT (GA, NPA, IS, ST, stool, urine, blood, FNA); Xpert or MGIT (GA, NPA, IS, ST, stool, urine); column labels set vertically below, group labels under brackets
Vertical: 32 rows numbered 1 to 32 from top to bottom, with a gap between rows 7 and 8 separating the 2 HIV groups
Visible pattern: The FM group is almost entirely white, with black dots in only a handful of rows, mostly near the top of each HIV group and in the last row. The Xpert and MGIT groups carry many more black dots, concentrated in the upper rows of each HIV group and in the GA, NPA and IS columns, thinning to mostly white toward the bottom. Small gray dots appear almost only in the MGIT group, scattered down many rows and most often in the stool column. The FNA and blood columns are sparse, with dots in only a few rows. The rightmost group, Xpert or MGIT, is the densest in black, and like the others runs from mostly black at the top of each HIV group to mostly white at the bottom; the fraction labels descend from 6/6 to 1/6 within each group.

Presenter notes

Specimen types: gastric aspirate (GA), nasopharyngeal aspirate (NPA), induced sputum (IS), string (ST), stool, urine, blood, lymph node aspirate (FNA)
Test types: fluorescence microscopy (FM), nucleic acid amplification (Xpert), broth culture (MGIT), union of Xpert and MGIT

Abbreviations: FM = fluorescence microscopy; FNA = fine-needle aspirate (of a lymph node); GA = gastric aspirate; HIV = human immunodeficiency virus; HIV+ = HIV-positive; IS = induced sputum; MGIT = Mycobacteria Growth Indicator Tube (liquid culture system); NPA = nasopharyngeal aspirate; ST = string test (specimen collection method); TB = tuberculosis; Xpert = Xpert MTB/RIF (Cepheid GeneXpert cartridge nucleic acid amplification test for TB and rifampicin resistance).

Slide 25Dot plot in base R (1)

Numeric matrix (i.e., a vector with 2-D dim attribute)
- 32 rows: children with any TB, in study ID order
- 53 columns: each test/specimen type
- Values 1 (positive), 2 (negative), 3 (indeterminate), 4:6
Compute x- and y-coordinates, plot character/colors
Create blank canvas
Fill in points (the dots of the dot plot)
Add annotations: 4 kinds of row labels, 3 kinds of column labels

Abbreviations: TB = tuberculosis.

Slide 26Dot plot in base R (2)

# row order by HIV, number of { positive types, positive specimens, positive X/M }
dot_row_order <- order(cbind(hiv, npostype, ntype, npos[, , "xm"])[anypos])
dot_pos_values <- c(dot_pos[dot_row_order, ]) # determines what is plotted
dot_x_values; dot_x_values_rep # constants to achieve desired horizontal spacing
dot_y_values; dot_y_values_rep # constants to achieve desired vertical spacing
# blank plot region
plot(c(1, 53), c(0, 34), bty="n", xlab="", ylab="", xaxt="n", yaxt="n", type="n")
points(dot_x_values_rep, dot_y_values_rep, # x- and y-coordinates
pch=c(21, 21, 21, NA, NA, NA)[dot_pos_values], # plot characters
col=c("black", "black", "white", NA, NA, NA)[dot_pos_values], # fill color
bg=c("black", "white", "gray", NA, NA, NA)[dot_pos_values]) # border color

Abbreviations: HIV = human immunodeficiency virus.

Slide 27Compare diagnostic yields

Example questions
- What’s the most feasible combination of specimens for diagnosing TB in children in low-resource settings?
- Can minimally invasive combinations perform well?
Analytic method
- For a given combination of specimens, it’s easy to count children deemed positive by that combo [examples]
- Just as easy to compare 2 combinations
- How do we separate signal from noise?
  - Resample children with replacement: the 31* with TB, 263 without
  - Calculate yield of various combinations for each pseudosample
  - Summarize across pseudosamples

Description

Panels: 1
Elements: Dot plot of 32 numbered rows by 6 specimen-type columns; each cell holds up to 2 circles, filled black or open white, or is empty; rows in 2 bracketed groups labeled HIV+ (rows 1 to 7) and HIV– (rows 8 to 32); a fraction label at the right end of each row, ranging from 6/6 to 1/1
Horizontal: 6 specimen-type columns, labeled vertically below as GA, NPA, IS, ST, stool, urine, bracketed together under the label Xpert or MGIT; up to 2 dots side by side per column
Vertical: 32 rows numbered 1 to 32 from top to bottom, with a gap between rows 7 and 8 separating the 2 HIV groups
Visible pattern: Within each HIV group the rows run from mostly black at the top to mostly white at the bottom, and the fraction labels at the right descend in step, from 6/6 near the top of each group to 1/6 near its bottom. Black dots are densest in the GA, NPA and IS columns; the ST, stool and urine columns are mostly white below the first several rows of each group. Some cells hold a single dot rather than a pair, and some are empty. The last 2 rows carry the fewest dots, the last only a single black dot, under NPA.

Slide 28Analysis step 1: 729 combinations in 1

Consider a child with {0,1,2} results each for 6 specimen types
- These 12 results determine status on 729 possible combinations
- {0,1,2} × {GA, NPA, IS, string, stool, urine} → 3^6 = 729
  - 0 of each (the null combination),
  - 1 GA specimen and no others, …
  - 2 of each specimen type
- If 1 of a given type, use 1st or 2nd? Use 1st, but randomly permute.
For any collection of children, we know how many would be deemed positive on each of these 729 combinations

Abbreviations: GA = gastric aspirate; IS = induced sputum; NPA = nasopharyngeal aspirate.

Slide 29Analysis step 1 in R: 729 combos

# template for {0,1,2}^6 = 729 combinations by specimen type
# panel6 -- panel of 6 specimen types
# cumu -- cumulative result for each combination
panel6_cumu <- expand.grid(urine=0:2, stool=0:2, strng=0:2, inspt=0:2, npasp=0:2,
gsasp=0:2)[, 6:1]
panel6_cumu <- data.table( # could also use data.frame or similar
pat=apply(panel6_cumu, 1, paste, collapse=""), # length-6 pattern label for combo
panel6_cumu, # result of expand.grid
ntype=rowSums(panel6_cumu[, 1:6] > 0), # number of specimen types in combo
nspec=rowSums(panel6_cumu[, 1:6])) # number of specimens in combo
# pat gsasp npasp inspt strng stool urine ntype nspec
# 2: 000001 0 0 0 0 0 1 1 1
# 3: 000002 0 0 0 0 0 2 1 2
# 4: 000010 0 0 0 0 1 0 1 1
# 5: 000011 0 0 0 0 1 1 2 2

Slide 30Analysis step 2: resample

Imagine a matrix with 32 rows and 53 columns
- Rows 1–31: observed results for 31 children with TB
- Row 32: all 0s, representing other 263 TB-negative children
Sampling scheme
- Sample 32 rows with replacement:
multinomial(294; 1/294, …, 1/294, 263/294)
- For each specimen, favor 1st or 2nd result? binom(2, 1/2)
- Efficiently construct 10,000 pseudosamples of size 294
- This is the full extent of randomness used

Abbreviations: TB = tuberculosis.

Slide 31Analysis step 2 in R: resample

# set.seed(24601) # prison code of Jean Valjean in Hugo's Les Misérables
# 10,000 multinomial pseudosamples: 32 x 1e4 integer matrix
pt_rep_sim <- rmultinom(1e4, size = 294, prob = rep(c(1, 263), c(31, 1))/294)
# row 32 is number in each pseudosample with no positive results
# number positive in each pseudosample
n_anypos_sim <- 294 - pt_rep_sim[32, ] # same as colSums(pt_rep_simx[-32, ])
# equal in distribution to binomial(294, 31/294), mostly 14–50
# randomize order within specimen pairs for each pseudosample
spec_ord_sim <- rbinom(sum(n_anypos_sim) * 6, size=1, prob=0.5)
# it remains only to analyze the pseudosamples; no other randomness needed

Presenter notes

https://en.wikipedia.org/wiki/24,601 “It is the first prison code of the character Jean Valjean in the novel (his only code in the musical) Les Misérables. It was chosen by Victor Hugo when he believed that he was conceived on 24 June 1801 (that is, 24-6-01).” The other code is 9430. https://en.wikipedia.org/wiki/Jean_Valjean

Slide 32Analysis step 3: crunch (pseudo)samples

For observed samples
- Deterministically select 1st or 2nd of all specimen pairs
- Compute 2 matrices, 729 × 31, with each child’s result on each combo
- Report yield or ∆ yield for all specimen combinations of interest
For pseudosamples, use explicit for loop*
- Stochastically select 1st or 2nd of all specimen pairs (spec_ord_sim)
- Compute 10,000 matrices, 729 × n_anypos, with each pseudochild’s result on each combo
- Use distribution of yield or ∆ yield for all specimen combinations of interest, for example, to report deciles or quartiles

* Because there’s too much variability to use an *apply function efficiently

Slide 33Cumulative and incremental yield by specimen type

Chart of point estimates with 2-level intervals: cumulative and incremental counts of TB-positive participants for 6 specimen types and 3 tests.

Description

Panels: 1
Elements: 6 groups of 9 point markers, each marker with a thick black inner interval and a thinner light gray outer interval; marker shape and fill per legend: circles for 1 Xpert (white), 1 MGIT (gray), 1 Xpert/MGIT (black); squares for 2 Xpert (white), 2 MGIT (gray), 2 Xpert/MGIT (black); diamonds for Xpert increment (white), MGIT increment (gray), Xpert/MGIT increment (black). Within each group, circle-square pairs sit side by side at 3 tick positions (white, gray, black) with the diamonds below them; numeric labels beside the black markers; legend in 3 columns across the top
Horizontal: 6 specimen-type categories, each with a 3-tick bracket: gastric aspirate, nasopharyngeal aspirate, induced sputum, string, stool, urine; axis title Specimen type
Vertical: 0 to 30, ticks every 5; axis title Number of participants with TB–positive result
Visible pattern: The circle and square markers are highest for gastric and nasopharyngeal aspirate (black squares about 23 to 24), lower for induced sputum (about 20), lower again for string and stool (about 14 to 15), and lowest for urine, where all markers cluster between 0 and about 4. In stool the gray markers sit well below the white and black ones, at about 6 to 7. The increment diamonds lie along the bottom of every group, between 0 and about 5.

Text in the picture

24.0 21.0 3.0
23.0 20.5 2.5
20.0 17.0 3.0
15.0 12.5 2.5
14.0 12.5 1.5
4.0 0.0
"1 specimen" means either specimen, and "2 specimens" means up to 2.

Abbreviations: MGIT = Mycobacteria Growth Indicator Tube (liquid culture system); TB = tuberculosis; Xpert = Xpert MTB/RIF (Cepheid GeneXpert cartridge nucleic acid amplification test for TB and rifampicin resistance).

Slide 34Average cumulative yield by number of specimens

Chart of point estimates with 2-level intervals: TB-positive participant counts by number of specimens, for 3 tests.

Description

Panels: 1
Elements: 12 triplets of diamonds, white (Xpert), gray (MGIT) and black (Xpert/MGIT) per legend at right, each with a thick black inner interval and a thinner light gray outer interval; numeric labels beside the black diamonds
Horizontal: Number of specimens in combination, 1 to 12
Vertical: 0 to 30, ticks every 5; axis title Number of participants with TB–positive result; the black diamonds and their intervals extend above the top tick to about 31
Visible pattern: All 3 series rise from left to right, steeply over the first few specimens and then flattening. Black stays above gray, and gray above or level with white, throughout: black from about 15 to 31, gray from about 12 to 26, white from about 10 to 26, with white and gray meeting at about 26 at 12 specimens. Intervals are widest at 1 specimen, where the gray outer interval spans roughly 2 to 20, and narrow steadily to the right; at 12 specimens the black diamond shows no interval.

Text in the picture

14.6 19.4 22.6 24.6 26.2 27.4 28.4 29.2 29.8 30.3 30.7 31.0

Slide 35All your base R, redux

Efficiency, ease; transparent, no dependencies

Description

Text in the picture

CATS: ALL YOUR BASE R BELONG
TO US.

Slide 36Expansive analysis, bespoke

Small data → large-sample statistical methods?
3 examples of “expansive analysis”—build up probability model for small-data effects
- Full randomization distribution
- Sampled approximation to full distribution
- Resampling to preserve complex multivariate distributions

Used R-native tools* to implement methods
esp multidimensional arrays, lists
R packages exist to generate
- Permutations
- Pseudosamples
- But none satisfied at time of need
Intellectually satisfying to exploit native features and exercise granular control

* Analogs exist in other environments, of course.

Slide 37End

A redacted region after the presenter's name; HHS and CDC logos at lower right on a blue banner.

Chad Heilig, PhD • [redacted: email address]
The views presented are the author’s own and do not necessarily represent the official position of the Centers for Disease Control and Prevention.

Description

Text in the picture

www.cdc.gov

Abbreviations: CDC = Centers for Disease Control and Prevention; HHS = Department of Health and Human Services.

Slide 38 is omitted.

Slide 39"I love it when you call me Big Data"

Description

Text in the picture

I LOVE IT WHEN
YOU CALL ME BIG DATA

Abbreviations

CDC = Centers for Disease Control and Prevention; CI = contact investigation; DSS = demographic surveillance system; Dx = diagnosis; FM = fluorescence microscopy; FNA = fine-needle aspirate (of a lymph node); GA = gastric aspirate; GSUs = geographic study units; HC = health-facility enhancement plus contact investigation (trial condition H + C); HCM = health-facility enhancement plus contact investigation plus mobile unit (trial condition H + C + M); HF = health facility; HHS = Department of Health and Human Services; HIV = human immunodeficiency virus; HIV+ = HIV-positive; HM = health-facility enhancement plus mobile unit (trial condition H + M); IS = induced sputum; ISO = International Organization for Standardization; MGIT = Mycobacteria Growth Indicator Tube (liquid culture system); MU = mobile unit; NPA = nasopharyngeal aspirate; ST = string test (specimen collection method); TB = tuberculosis; Xpert = Xpert MTB/RIF (Cepheid GeneXpert cartridge nucleic acid amplification test for TB and rifampicin resistance).

References, compiled for this page

Full citations, with current links, for selected references and related works.

Song R, Click ES, McCarthy KD, Heilig CM, Mchembere W, Smith JP, Fajans M, Musau SK, Okeyo E, Okumu A, et al. Sensitive and feasible specimen collection and testing strategies for diagnosing tuberculosis in young children. JAMA Pediatr. 2021;175(5):e206069. DOI 10.1001/jamapediatrics.2020.6069.
Yuen CM, Agaya J, Mchembere W, Okelloh D, Achola M, Opole J, Cowden J, Heilig CM, Borgdorff MW, Cain KP. Optimizing the efficiency of tuberculosis active case-finding in health facilities and communities. Int. J. Tuberc. Lung Dis. 2019;23(7):844–849. DOI 10.5588/ijtld.18.0700.