Epidemiologist: MQL recipes

Recipes

SEIR outbreak in a partly vaccinated population
From early growth to R and to the final attack rate
Weekly incidence report with exact Poisson limits
Hospital admissions: delay from cases and short-term forecast
Case-control study: odds ratio and attributable fraction
Diagnostic test: Wilson interval, likelihood ratios, predictive values
Incubation period and quarantine length
Occupational cohort: SMR with exact Poisson test
Vaccine effectiveness from a test-negative study
Direct age standardization of mortality

executed on MentDB

SEIR outbreak in a partly vaccinated population

A measles-like virus (R0 = 12, latent 8 days, infectious 5 days) enters a population with 85 % coverage of a 95 % effective vaccine. The script computes the herd immunity threshold 1 - 1/R0 (91.7 %), starts the immune fraction in R and integrates the SEIR model over one year to see the outbreak that still happens.

Pipeline

01 · ratesbeta = R0 x gamma, sigma, gammamath calc→02 · thresholdHerd immunity 1 - 1/R0, effective coveragemath calc→03 · simulateSEIR system over 365 daysnum ode→04 · plotS, E, I, R against timeplot ode

Script

admin

#1. Measles-like SEIR: R0 = 12, latent 8 days, infectious 5 days;
-> "[R0]" 12; -> "[sigma]" (math calc "1/8"); -> "[gamma]" (math calc "1/5");
-> "[beta]" (math calc "[R0]*[gamma]");
#2. Herd immunity threshold 1 - 1/R0, coverage 85 % with a 95 % effective vaccine;
-> "[pc]" (math calc "1-1/[R0]"); -> "[v]" (math calc "0.85*0.95");
#3. SEIR over one year starting from 10 cases per 100 000;
plot ode (num ode "{\"S\": \"-[beta]*S*I\", \"E\": \"[beta]*S*I-[sigma]*E\", \"I\": \"[sigma]*E-[gamma]*I\", \"R\": \"[gamma]*I\"}" (concat "{\"S\": " (math calc "1-[v]-0.0001") ", \"E\": 0, \"I\": 0.0001, \"R\": " [v] "}") 0 365 "{\"points\": 366}")

Result

mentdb

The chart is shown in the Plot tab of the editor.

Go further

Print the effective reproduction number R0 x S(0) with math calc: above 1 the outbreak takes off
Add an event in the num ode options to get the day when I falls back under its initial value
Search the coverage that keeps Reff below 1 with num solve

executed on MentDB

From early growth to R and to the final attack rate

Fit the exponential growth rate of the first two weeks of daily cases, turn it into a reproduction number with a gamma generation interval (Wallinga-Lipsitch: R = 1 / M(-r), integral computed with the gamma function) and solve the final size equation z = 1 - exp(-R z).

Pipeline

01 · fitExponential growth a exp(r t)stat fit_model→02 · R1 / integral of exp(-r t) g(t), g gamma lawnum integrate · gamma(x)→03 · final sizez = 1 - exp(-R z)num solve→04 · projectGrowth factor over 7 daysbigdec pow

Script

admin

#1. Daily confirmed cases during the first two weeks of the outbreak;
-> "[d]" "[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]";
-> "[c]" "[12, 15, 19, 23, 30, 37, 45, 57, 70, 88, 108, 135, 166, 205]";
#2. Exponential growth rate r and doubling time;
-> "[r]" (json load select "rx_gr" (stat fit_model "a*exp(r*x)" "x" [d] [c] "{\"a\": 10, \"r\": 0.2}") "/parameters/r");
-> "[td]" (math calc "ln(2)/[r]");
#3. Generation interval: gamma law, mean 5.2 days, sd 1.7 days;
-> "[k]" (math calc "(5.2/1.7)^2"); -> "[th]" (math calc "1.7^2/5.2");
#4. Wallinga-Lipsitch: R = 1 / integral of exp - r t times g of t;
-> "[M]" (num integrate (concat "exp(-" [r] "*t)*t^(" [k] "-1)*exp(-t/" [th] ")/(gamma(" [k] ")*" [th] "^" [k] ")") "t" 0 "inf");
-> "[R]" (math calc "1/[M]");
#5. Final size without control: z = 1 - exp - R z, and 7-day growth factor;
-> "[z]" (num solve (concat "1-exp(-" [R] "*x)-x") "x" 0.01 0.999);
-> "[g7]" (bigdec pow (math exp [r]) 7 6);
concat "r = " (math round [r] 4) " per day, doubling time " (math round [td] 2) " days, R = " (math round [R] 3) " (closed form " (math round (math calc "(1+[r]*[th])^[k]") 3) "), expected attack rate " (math round (math calc "100*[z]") 1) " %, growth factor over 7 days x" [g7]

Result

mentdb

r = 0.2149 per day, doubling time 3.22 days, R = 2.874 (closed form 2.874), expected attack rate 93.1 %, growth factor over 7 days x4.50232

Go further

Get the uncertainty of r from the field errors of stat fit_model and propagate it to R with num uncertainty
Compare with an exponential generation interval: R = 1 + r Tg
Draw the fit with plot fit

executed on MentDB

Weekly incidence report with exact Poisson limits

Four weeks of daily notifications show a weekend under-reporting artefact. The script measures the weekly cycle (autocorrelation at lag 7 and spectrum at 1/7 per day on the detrended series), smooths with a 7-day moving average, and reports the 7-day incidence per 100 000 with exact Poisson limits (Garwood) from the regularized incomplete gamma function.

Pipeline

01 · detrendRemove the linear trenddsp detrend→02 · cycleLag-7 autocorrelation and 1/7 componentdsp autocorr · dsp spectrum→03 · smooth7-day moving averagedsp moving_average→04 · rateExact sum and incidence per 100 000bigdec sum · bigdec div · bigdec compare→05 · CIPoisson limits with gamma_incnum solve · gamma_inc(a, x)

Script

admin

#1. Daily notified cases over 4 weeks: low on weekends;
-> "[cases]" "[40, 46, 51, 48, 44, 21, 18, 47, 55, 58, 53, 50, 24, 20, 56, 62, 66, 61, 57, 27, 23, 64, 70, 73, 69, 63, 30, 26]";
#2. 7-day moving average and weekly periodicity of the detrended series;
-> "[ma]" (dsp moving_average [cases] 7);
-> "[dt]" (dsp detrend [cases]);
-> "[r7]" (json load select "rx_ac" (dsp autocorr [dt] 7) "/[7]");
-> "[a7]" (json load select "rx_sp" (dsp spectrum [dt] 1) "/amplitude[4]");
#3. Last 7 days and incidence per 100 000 for a district of 245 000 people;
-> "[k]" (bigdec sum "[64, 70, 73, 69, 63, 30, 26]");
-> "[inc]" (bigdec div (math calc "[k]*100000") 245000 4);
#4. Exact Poisson 95 % limits with the regularized incomplete gamma function;
-> "[lo]" (num solve (concat "gamma_inc(" [k] ", x)-0.025") "x" 1 [k]);
-> "[hi]" (num solve (concat "1-gamma_inc(" [k] "+1, x)-0.025") "x" [k] (math calc "3*[k]"));
concat "weekly autocorrelation r7 = " (math round [r7] 3) ", amplitude of the 1/7 per day component " (math round [a7] 1) " cases" ", last 7-day average " (math round (json load select "rx_ma" [ma] "/[21]") 1) ", " [k] " cases = " [inc] " per 100 000, " (if (> (bigdec compare [inc] 150) 0) {"above"} {"below"}) " the alert threshold of 150 (95 % CI " (math round (math calc "[lo]*100000/245000") 1) " to " (math round (math calc "[hi]*100000/245000") 1) ")"

Result

mentdb

weekly autocorrelation r7 = 0.721, amplitude of the 1/7 per day component 20.3 cases, last 7-day average 56.4, 395 cases = 161.2 per 100 000, above the alert threshold of 150 (95 % CI 145.7 to 177.9)

Go further

Compare two districts with a rate ratio and its exact interval
Look for the outbreak signal with dsp peaks on the moving average
Draw the daily cases and the average with plot chart

executed on MentDB

Hospital admissions: delay from cases and short-term forecast

Admissions follow cases with a delay. The script measures the lag of the maximum cross-correlation between the two centered series, discretizes a gamma delay distribution (mean 6 days, sd 2.5 days) with the regularized incomplete gamma function, and convolves the cases with it (8 % of the cases are hospitalized) to fit today's admissions and forecast the next days from the cases already reported.

Pipeline

01 · centerRemove the meansdsp detrend→02 · lagCross-correlation and its maximumdsp correlate · dsp peaks · stat max→03 · delayP(d < D <= d+1) for the gamma lawgamma_inc(a, x) · math calc→04 · forecastCases convolved with the delaydsp convolve

Script

admin

#1. Daily cases and hospital admissions of a region over 24 days;
-> "[cases]" "[100, 130, 165, 210, 265, 330, 400, 475, 550, 610, 650, 665, 650, 610, 550, 480, 410, 345, 285, 230, 185, 145, 115, 90]";
-> "[adm]" "[0, 0, 0, 2, 4, 5, 7, 14, 15, 20, 27, 30, 37, 40, 44, 45, 48, 49, 45, 43, 39, 33, 31, 26]";
#2. Lag between the two curves: maximum of the cross-correlation of the centered series;
-> "[cc]" (dsp correlate (dsp detrend [adm] "constant") (dsp detrend [cases] "constant"));
-> "[lag]" (- (json load select "rx_cc" (dsp peaks [cc] (concat "{\"height\": " (math calc (concat (stat max [cc]) "*0.999")) "}")) "/indices[0]") 23);
#3. Delay from case to admission: gamma law mean 6 days, sd 2.5 days, discretized with gamma_inc;
-> "[k]" (math calc "(6/2.5)^2"); -> "[th]" (math calc "2.5^2/6");
-> "[pmf]" "["; -> "[sep]" "";
for (-> "[d]" 0) (< [d] 15) (++ "[d]") {
  -> "[pmf]" (concat [pmf] [sep] (math calc "gamma_inc([k], ([d]+1)/[th])-gamma_inc([k], [d]/[th])")); -> "[sep]" ",";
};
-> "[pmf]" (concat [pmf] "]");
#4. Expected admissions = 8 % of the cases convolved with the delay: fit of today, forecast of tomorrow;
json load "rx_ex" (dsp convolve [cases] [pmf]);
concat "lag of the maximum correlation " [lag] " days, sum of the delay probabilities " (math round (stat sum [pmf]) 4) ", expected admissions today " (math round (math calc (concat "0.08*" (json select "rx_ex" "/[23]"))) 1) " (observed 26), tomorrow " (math round (math calc (concat "0.08*" (json select "rx_ex" "/[24]"))) 1) ", in 3 days " (math round (math calc (concat "0.08*" (json select "rx_ex" "/[26]"))) 1) " at least"

Result

mentdb

lag of the maximum correlation 4 days, sum of the delay probabilities 0.9967, expected admissions today 25.7 (observed 26), tomorrow 21.6, in 3 days 14.1 at least

Go further

Fit the hospitalized fraction by least squares between the observed and the convolved series with stat fit_model
Add the cases forecast by the growth rate to complete the admissions of the coming days
Draw observed and expected admissions with plot chart

executed on MentDB

Case-control study: odds ratio and attributable fraction

Smoking among 100 lung cancer cases and 100 controls: odds ratio with the Woolf 95 % confidence interval, chi2 and Fisher exact tests, attributable fraction in the exposed and population attributable fraction (Miettinen).

Pipeline

01 · ORad/bcmath calc→02 · CIexp(ln OR +- z.SE)stat dist_inv→03 · testschi2 and Fisher exactstat chi2_independence→04 · impactAFe and PAFmath calc

Script

admin

#1. Case-control study: smokers among 100 lung cancer cases and 100 hospital controls;
-> "[a]" 45; -> "[b]" 55; -> "[c]" 25; -> "[d]" 75;
-> "[tab]" (concat "[[" [a] ", " [b] "], [" [c] ", " [d] "]]");
#2. Odds ratio with Woolf 95 % confidence interval;
-> "[or]" (math calc "[a]*[d]/([b]*[c])");
-> "[z]" (stat dist_inv "normal" "[0, 1]" 0.975);
-> "[se]" (math calc "sqrt(1/[a]+1/[b]+1/[c]+1/[d])");
#3. Tests of association: Pearson chi2 and Fisher exact;
-> "[chi]" (json load select "rx_tu1_cc" (stat chi2_independence [tab]) "/p_value");
-> "[fis]" (json load select "rx_tu1_cf" (stat fisher_exact [tab]) "/p_value");
#4. Attributable fraction in the exposed and population attributable fraction, Miettinen formula;
-> "[afe]" (math calc "([or]-1)/[or]");
-> "[paf]" (math calc "[a]/([a]+[b])*[afe]");
concat "OR = " (math round [or] 3) ", 95 % CI " (math round (math calc "exp(ln([or])-[z]*[se])") 2) " to " (math round (math calc "exp(ln([or])+[z]*[se])") 2) "; chi2 p = " (math round_sig [chi] 3) ", Fisher p = " (math round_sig [fis] 3) "; AF exposed = " (math round (* 100 [afe]) 1) " %, PAF = " (math round (* 100 [paf]) 1) " %"

Result

mentdb

OR = 2.455, 95 % CI 1.35 to 4.47; chi2 p = 0.00303, Fisher p = 0.00468; AF exposed = 59.3 %, PAF = 26.7 %

Go further

Adjust for age with a stratified table and the Mantel-Haenszel odds ratio
Compute the sample size needed to detect OR = 2 with a power of 80 %

executed on MentDB

Diagnostic test: Wilson interval, likelihood ratios, predictive values

A rapid antigen test is validated against PCR. The script gives the sensitivity with its Wilson score interval (z from the inverse error function), the specificity, both likelihood ratios and the predictive values in a 2 % prevalence screening.

Pipeline

01 · zsqrt(2).erfinv(0.95)math erf_inv→02 · intervalWilson scoremath sqrt→03 · LRLR+ and LR-math round→04 · PPVBayes at 2 % prevalencemath calc

Script

admin

#1. Validation of a rapid antigen test against PCR on 1000 people;
-> "[tp]" 92; -> "[fn]" 8; -> "[fp]" 45; -> "[tn]" 855;
-> "[se]" (/ [tp] (+ [tp] [fn])); -> "[sp]" (/ [tn] (+ [tn] [fp]));
#2. z of a 95 % interval from the inverse error function: sqrt 2 . erfinv 0.95;
-> "[z]" (* (math sqrt 2) (math erf_inv 0.95));
#3. Wilson score interval of the sensitivity;
-> "[n]" (+ [tp] [fn]);
-> "[c]" (/ (+ [se] (/ (* [z] [z]) (* 2 [n]))) (+ 1 (/ (* [z] [z]) [n])));
-> "[h]" (/ (* [z] (math sqrt (+ (/ (* [se] (- 1 [se])) [n]) (/ (* [z] [z]) (* 4 (* [n] [n])))))) (+ 1 (/ (* [z] [z]) [n])));
#4. Likelihood ratios and predictive values at a 2 % prevalence in screening;
-> "[lrp]" (/ [se] (- 1 [sp])); -> "[lrn]" (/ (- 1 [se]) [sp]);
-> "[ppv]" (math calc "[se]*0.02/([se]*0.02+(1-[sp])*0.98)");
-> "[npv]" (math calc "[sp]*0.98/([sp]*0.98+(1-[se])*0.02)");
concat "Se = " [se] " [" (math round (- [c] [h]) 3) "; " (math round (+ [c] [h]) 3) "] with z = " (math round [z] 4) ", Sp = " (math round [sp] 4) ", LR+ = " (math round [lrp] 1) ", LR- = " (math round [lrn] 3) ", PPV at 2 % = " (math round (* 100 [ppv]) 1) " %, NPV = " (math round (* 100 [npv]) 2) " %"

Result

mentdb

Se = 0.92 [0.85; 0.959] with z = 1.96, Sp = 0.95, LR+ = 18.4, LR- = 0.084, PPV at 2 % = 27.3 %, NPV = 99.83 %

Go further

Draw the PPV against the prevalence with plot function.
Compare two tests on the same patients with the McNemar test (stat chi2_test).

executed on MentDB

Incubation period and quarantine length

From a lognormal incubation period fitted on contact-tracing data, the script gives the median and the mean, the probability that symptoms start after a 14 day quarantine and the quarantine length that covers 99 % of the cases.

Pipeline

01 · mediane^mumath e→02 · meanexp(mu + sigma^2/2)math power→03 · tailP(T > 14 d) with erfcmath erfc→04 · quantile99th percentile with erf_invmath erf_inv→05 · daysRounded upmath ceil

Script

admin

#1. Incubation period lognormal: mean and SD of ln days from a contact-tracing study;
-> "[mu]" 1.63; -> "[sg]" 0.50;
#2. Median e^mu and mean e^mu+sigma2/2;
-> "[med]" (math power (math e) [mu]);
-> "[mean]" (math exp (+ [mu] (/ (* [sg] [sg]) 2)));
#3. Probability that symptoms appear after a 14 day quarantine, with erfc;
-> "[p14]" (* 0.5 (math erfc (/ (- (math log 14) [mu]) (* [sg] (math sqrt 2)))));
#4. Quarantine length covering 99 % of the cases: quantile with erf_inv;
-> "[q99]" (math exp (+ [mu] (* [sg] (* (math sqrt 2) (math erf_inv 0.98)))));
concat "median = " (math round [med] 2) " d, mean = " (math round [mean] 2) " d, P incubation > 14 d = " (math round_sig (* 100 [p14]) 3) " %, 99th percentile = " (math round [q99] 1) " d, so a quarantine of " (math ceil [q99]) " days"

Result

mentdb

median = 5.1 d, mean = 5.78 d, P incubation > 14 d = 2.18 %, 99th percentile = 16.3 d, so a quarantine of 17 days

Go further

Fit mu and sigma on observed exposure-onset pairs with stat fit_model on the log values.
Check the quantile with stat dist_inv and the lognormal law.

executed on MentDB

Occupational cohort: SMR with exact Poisson test

27 lung cancer deaths are observed in a cohort where national rates predict 18.4. The script sums the Poisson terms in a loop with log-gamma to get the exact one-sided p-value, and gives the SMR with its exact interval from the chi2 quantiles.

Pipeline

01 · loopPoisson terms exp(k ln E - E - lnGamma(k+1))math lngamma→02 · p1 - P(X <= 26)math exp→03 · intervalChi2 quantiles / 2Estat dist_inv→04 · checkPoisson lawstat dist_cdf

Script

admin

#1. Occupational cohort: 27 lung cancer deaths observed, 18.4 expected from national rates;
-> "[o]" 27; -> "[e]" 18.4;
#2. Exact Poisson probability P of X >= 27, terms exp of k.ln e - e - lngamma k+1;
-> "[cum]" 0;
for (-> "[k]" 0) (< [k] [o]) (++ "[k]") {
	-> "[cum]" (+ [cum] (math exp (- (- (* [k] (math log [e])) [e]) (math lngamma (+ [k] 1)))));
};
-> "[p]" (- 1 [cum]);
#3. SMR and its exact 95 % interval from the chi2 quantiles;
-> "[lo]" (/ (/ (stat dist_inv "chi2" (concat "[" (* 2 [o]) "]") 0.025) 2) [e]);
-> "[hi]" (/ (/ (stat dist_inv "chi2" (concat "[" (+ (* 2 [o]) 2) "]") 0.975) 2) [e]);
concat "SMR = " (math round (/ [o] [e]) 2) " [" (math round [lo] 2) "; " (math round [hi] 2) "], one-sided exact p = P X>=27 = " (math round_sig [p] 3) ", check with the Poisson law: " (math round_sig (- 1 (stat dist_cdf "poisson" (concat "[" [e] "]") 26)) 3)

Result

mentdb

SMR = 1.47 [0.97; 2.13], one-sided exact p = P X>=27 = 0.0354, check with the Poisson law: 0.0354

Go further

Stratify by age and period and add the expected deaths with vector dot.
Use the Byar approximation and compare with the exact limits.

executed on MentDB

Vaccine effectiveness from a test-negative study

Influenza vaccine effectiveness is estimated as 1 - OR from a test-negative design, with the Woolf interval, and the length of hospital stay of breakthrough cases is compared with the unvaccinated cases by a Mann-Whitney test.

Pipeline

01 · OROdds ratiomath calc→02 · intervalWoolf log intervalstat dist_inv→03 · severityLength of stay, rank teststat mann_whitney→04 · readp-valuejson load select→05 · summaryMediansstat median

Script

admin

#1. Test-negative design: vaccinated among influenza-positive cases and test-negative controls;
-> "[a]" 48; -> "[b]" 312; -> "[c]" 151; -> "[d]" 489;
#2. Vaccine effectiveness VE = 1 - OR with the Woolf interval;
-> "[or]" (/ (* [a] [d]) (* [b] [c]));
-> "[se]" (math sqrt (+ (/ 1 [a]) (+ (/ 1 [b]) (+ (/ 1 [c]) (/ 1 [d])))));
-> "[z]" (stat dist_inv "normal" "[0, 1]" 0.975);
-> "[vel]" (- 1 (math exp (+ (math log [or]) (* [z] [se]))));
-> "[veh]" (- 1 (math exp (- (math log [or]) (* [z] [se]))));
#3. Severity among breakthrough cases: hospital length of stay in days, vaccinated vs unvaccinated;
-> "[vac]" "[2, 3, 2, 4, 3, 2, 5, 3]";
-> "[unv]" "[4, 6, 3, 8, 5, 7, 4, 9, 6]";
-> "[mw]" (json load select "rx_p4d_ve" (stat mann_whitney [vac] [unv]) "/p_value");
concat "OR = " (math round [or] 3) ", VE = " (math round (* 100 (- 1 [or])) 1) " % [" (math round (* 100 [vel]) 1) "; " (math round (* 100 [veh]) 1) "], length of stay median " (stat median [vac]) " vs " (stat median [unv]) " days, Mann-Whitney p = " (math round_sig [mw] 3)

Result

mentdb

OR = 0.498, VE = 50.2 % [29; 65], length of stay median 3 vs 6 days, Mann-Whitney p = 0.00549

Go further

Adjust the OR on age with a Mantel-Haenszel estimate over age strata.
Repeat by week of onset to see the waning of the protection.

executed on MentDB

Direct age standardization of mortality

Cardiovascular mortality of two regions with different age structures is standardized on the European Standard Population: one matrix product gives both standardized rates, a diagonal matrix gives the weight of each age group, then the comparative mortality figure.

Pipeline

01 · standardizeRates x weightsmatrix mul→02 · contributionsdiag(weights)matrix diag→03 · old85+ contributionmatrix get→04 · compareCMFjson load select→05 · crudeUnweighted meanmath avg

Script

admin

#1. Cardiovascular mortality per 100 000 by age group 45-54 ... 85+, two regions, rows = regions;
-> "[r]" "[[85, 210, 540, 1450, 4100], [70, 190, 560, 1600, 4600]]";
#2. Weights of the 2013 European Standard Population for these ages, rescaled to 1;
-> "[w]" "[0.2824, 0.2588, 0.2235, 0.1647, 0.0706]";
#3. Direct standardization: R.w, then the contribution of each age group with diag w;
-> "[asr]" (matrix mul [r] [w]);
-> "[contrib]" (matrix mul [r] (matrix diag [w]));
-> "[old]" (matrix get [contrib] 1 5);
#4. Comparative mortality figure of region 2 against region 1;
-> "[a1]" (json load select "rx_p4d_asr" [asr] "/[0]");
-> "[a2]" (json select "rx_p4d_asr" "/[1]");
concat "age-standardized rates: " (math round [a1] 1) " and " (math round [a2] 1) " per 100 000, CMF = " (math round (/ [a2] [a1]) 3) ", unweighted mean of the region 1 age rates: " (math avg (concat "[" (string replace (string replace (matrix row [r] 1) "[" "") "]" "") "]")) ", the 85+ group makes " (math round (* 100 (/ [old] [a1])) 1) " % of the region 1 rate"

Result

mentdb

age-standardized rates: 727.3 and 782.4 per 100 000, CMF = 1.076, unweighted mean of the region 1 age rates: 1277.0, the 85+ group makes 39.8 % of the region 1 rate

Go further

Add the 95 % interval of each ASR with the variance sum of w^2.r/n (vector dot).
Indirect standardization: expected deaths = population x reference rates.