kkstatfun – Kostadin Kostadinov

One package, from bedside 2×2 tables to targeted causal inference

kkstatfun is the R engine behind my research and teaching. It wraps the messy, error-prone parts of medical statistics — measures of association, standardization, diagnostic accuracy, survival, regression, causal and quasi-experimental designs, health-technology assessment, and infectious-disease modelling — into 146 tidy, pipe-friendly functions that all speak the same grammar: tidy data in, tidy data out.

Why it exists

Most applied medical analyses reach for a different package at every step — epitools for a 2×2 table, survival for Kaplan–Meier, pROC for an ROC curve, gtsummary for Table 1, something else for standardization — each with its own conventions, output shapes, and quirks. kkstatfun is my answer to that fragmentation: a single, opinionated layer on top of the best of those engines so that a whole analysis reads as one coherent story.

Consistent

Every function follows tidy data in, tidy data out. Results come back as tibbles you can pipe, filter, join, and plot — never a print-only object you have to scrape.

Genuinely stratified

The analysis functions honour dplyr::group_by(): pipe a grouped data frame in and the analysis runs once per group, with the grouping columns prefixed to the result — real stratified estimates, not a pooled number with a loop on top.

Publication-ready

Confidence intervals, effect sizes, and correct tests by default — plus themed plots — so output drops straight into a manuscript or a lecture.

New in 1.3.0 — quasi-experimental designs. Difference-in-differences (including staggered adoption), regression discontinuity and kink, interrupted time series, synthetic control, and instrumental variables. These identify a causal effect from how the intervention was assigned rather than from adjusting for measured confounders, so they survive confounders nobody measured. See the Quasi-experimental designs tab in the function reference below.

Installation

# Install devtools if you don't have it
if (!require("devtools")) install.packages("devtools")

# Install kkstatfun from GitHub
devtools::install_github("kostadinoff/kkstatfun")

library(kkstatfun)

A 60-second tour

library(kkstatfun)
library(dplyr)

# 1. Measure of association from a 2x2 table
smoking_cancer <- data.frame(
  smoking     = c(rep(1, 80), rep(0, 120)),
  lung_cancer = c(rep(1, 30), rep(0, 50), rep(1, 10), rep(0, 110))
)
smoking_cancer %>% kk_twobytwo(smoking, lung_cancer)
#> Odds Ratio 6.6 (3.0-14.5), Relative Risk 4.5 (2.3-8.7), plus RD, AF, PAF...

# 2. The same analysis, stratified — group_by() just works
cohort %>%
  group_by(site) %>%
  kk_twobytwo(exposure, outcome)
#> One row-set per site: effect modification becomes visible,
#> where a pooled estimate would hide it.

# 3. Publication-ready Table 1, stratified by treatment arm
kk_table1(cohort, by = "arm", variables = c("age", "bmi", "sex", "smoker"))

# 4. A diagnostic work-up across every site at once
cohort %>%
  group_by(site) %>%
  kk_diagnostic(gold_standard, rapid_test)

Everything below is organized by the question you’re trying to answer.

Function reference

The core epidemiology toolkit: turn exposures and outcomes into interpretable, confidence-bounded effect measures.

Function	What it does
`kk_twobytwo()`	Full 2×2 work-up: OR, RR, risk difference, attributable and preventable fractions (AF/PAF/PF), phi, Yule’s Q, NNH — all with 95% CIs.
`kk_epi_stats()` / `odds_ratio()` / `risk_ratio()`	Lightweight OR and RR when you just need the point estimate and CI.
`kk_stratified_2x2()`	Stratified analysis with Cochran–Mantel–Haenszel pooled OR/RR and the Breslow–Day test for effect-modification/confounding.
`kk_reri()`	Relative Excess Risk due to Interaction (RERI), attributable proportion, and Synergy Index for additive interaction.
`kk_matched_case_control()`	Conditional logistic regression (`survival::clogit`) for 1:1, 1:m, or variable-ratio matched case-control studies.
`kk_case_crossover()`	Self-controlled case-crossover analysis — hazard vs. control windows within the same patient, controlling all time-invariant confounders.
`kk_paf()`	Crude and model-adjusted Population Attributable Fraction (Bruzzi standardization, Miettinen formula) with bootstrap CIs.
`kk_nnt()`	Number Needed to Treat / Harm with ARR and RRR, with confidence intervals.
`kk_survival_nnt()`	Time-varying ARR and NNT from Kaplan–Meier curves at specified horizons, with Greenwood-based CIs.
`kk_sensitivity_analysis()`	How strong would an unmeasured confounder have to be to explain away an observed OR/RR?
`kk_incidence_rate()`	Incidence rates per person-time with exact Poisson CIs, and rate ratios vs. a reference stratum.
`kk_trend_test()`	Cochran–Armitage test for trend in proportions across ordinal exposure levels.
`kk_std_rates()`	Direct age-standardization against a reference population (crude + standardized rates with CIs).
`kk_std_rates_ci()`	Directly standardized rates with exact Gamma CIs (Tiwari et al. 2006 / Fay–Feuer SEER method) for sparse counts.
`kk_smr()`	Indirect standardization → Standardized Mortality/Morbidity Ratio with an exact Poisson CI.

Example — a classic 2×2 exposure–outcome analysis:

df <- data.frame(
  smoking     = c(rep(1, 80), rep(0, 120)),
  lung_cancer = c(rep(1, 30), rep(0, 50), rep(1, 10), rep(0, 110))
)
df %>% kk_twobytwo(smoking, lung_cancer)

Note

Interpretation. A relative risk of 4.5 means exposed individuals have 4.5× the risk of the outcome versus unexposed. Because odds ratios overstate the risk ratio when the outcome is common, RR is generally preferred in cohort studies — kk_twobytwo() gives you both so you can report the right one.

Everything you need to validate a test, a marker, a model, or a rater.

Function	What it does
`kk_diagnostic()` / `confusion_metrics_ci()`	Sensitivity, specificity, PPV/NPV, likelihood ratios — with Wilson-score CIs.
`diagnostic_summary()` / `kk_confusion_matrix()`	Quick confusion-matrix summaries for a truth–test pair.
`kk_diagnostic_lrt()`	Bedside Bayesian update: pre-test probability + Se/Sp → post-test probability after a positive or negative result.
`kk_roc()`	ROC curve for a continuous marker: AUC with DeLong CI and the Youden-optimal cutoff.
`kk_compare_roc()`	Compare two markers’ AUCs on the same subjects (DeLong’s paired test).
`kk_calibration()`	Calibration of a risk model in one row: Hosmer–Lemeshow test, observed/expected ratio, calibration slope and intercept.
`kk_calibration_table()`	The per-group table behind a calibration curve, with the Brier score attached.
`kk_reclassification()`	Net Reclassification Improvement (NRI) and Integrated Discrimination Improvement (IDI) when adding a biomarker to a risk model.
`kk_decision_curve()`	Decision-curve analysis — net clinical benefit across threshold probabilities, vs. treat-all and treat-none.
`kk_kappa()` / `kk_agreement()`	Cohen’s κ, weighted κ (linear/quadratic), and PABAK for categorical agreement.
`kk_icc()`	Intraclass correlation in all six Shrout–Fleiss forms (continuous agreement).
`kk_bland_altman()`	Bland–Altman limits-of-agreement analysis and plot.
`kk_mcnemar()`	McNemar’s test for matched-pair (paired binary) data.
`kk_reliability()`	Cronbach’s α (raw + standardized) with alpha-if-item-dropped for scale validation.

Example — validate a rapid test against a gold standard:

df_test <- data.frame(
  gold_standard = c(rep(1, 100), rep(0, 900)),
  pcr_test      = c(rep(1, 95),  rep(0, 5), rep(1, 20), rep(0, 880))
)
df_test %>% kk_diagnostic(gold_standard, pcr_test)

# Turn Se/Sp into a bedside post-test probability
kk_diagnostic_lrt(pre_test_prob = 0.20, sensitivity = 0.95, specificity = 0.60,
                  n_diseased = 100, n_healthy = 400)

One consistent interface across the models applied research actually uses.

Function	What it does
`kk_reg()`	Unified wrapper: auto-detects binary outcomes (logistic, with ORs + ROC) vs. continuous (linear, with diagnostic plots).
`kk_firth()`	Firth-penalized logistic regression — finite, bias-reduced ORs under separation / rare events where ordinary logistic fails.
`kk_rr_reg()`	Adjusted risk-ratio regression (modified Poisson with robust SEs, Zou 2004) — RRs instead of ORs for common binary outcomes.
`kk_rate_reg()` / `kk_poisson()`	Poisson / negative-binomial rate regression with a person-time offset; overdispersion detected and handled automatically.
`kk_trend_nonlinearity()`	Test whether an exposure–outcome dose–response departs from linearity.
`kk_dose_response()`	Estimate dose–response effects relative to a reference exposure value.
`kk_compare_independent_correlations()` / `kk_compare_dependent_correlations()`	Compare two correlation coefficients — across independent samples or within the same sample.
`kk_apc()`	Age–Period–Cohort analysis for trends in rates over time.
`kk_model()` / `kk_model_data()`	Retrieve the fitted `lm`/`glm`/`coxph` object (and its data) that a modelling function kept, so nothing has to be refitted by hand to go further.
`kk_emmeans()`	Estimated marginal means and contrasts via `emmeans`, back-transformed to ORs, RRs or probabilities — model-adjusted means for the results table instead of raw group means.

Example — logistic regression returning tidy odds ratios:

fit <- kk_reg(cohort, outcome = "event", predictors = c("age", "smoker", "arm"))

# The fitted model is kept, so carry on with it rather than refitting
kk_emmeans(fit, ~ arm, contrast = "pairwise", type = "response")
emmeans::joint_tests(kk_model(fit))

Sample size, power, and trial-specific analyses alongside proportion comparisons.

Function	What it does
`kk_sample_size_epi()`	Sample-size calculations for common epidemiological designs (given baseline risk, target RR/OR, power, α).
`power_proportions()`	Power analysis for comparing two proportions.
`compare_proportions()` / `compare_proportions_by()`	Compare proportions across groups, optionally stratified.
`pcit()`	Proportion confidence intervals and tests.
`compare_proportions_kk_glm()`	GLM-based proportion comparison for grouped count data.
`plot_proportion_comparisons()`	Ready-made visualization of proportion-comparison results.
`kk_cluster_trial()`	Analysis of cluster-randomized trials, accounting for intra-cluster correlation.
`kk_crossover_trial()`	Analysis of two-period crossover trial data (period and sequence effects).

Example — sample size for a cohort study:

kk_sample_size_epi(design = "cohort", p0 = 0.10, rr_or = 2.0, power = 0.8, alpha = 0.05)

From Kaplan–Meier to competing risks and RMST.

Function	What it does
`kk_survival_plot()`	Publication-quality Kaplan–Meier curves with CIs and a risk table.
`kk_coxph()`	Cox proportional-hazards model with tidy hazard ratios and PH diagnostics.
`kk_logrank()`	Log-rank and related survival tests between groups.
`kk_rmst()`	Restricted Mean Survival Time (area under the KM curve to a horizon τ) with group differences — robust when proportional hazards fails.
`kk_cuminc()`	Competing-risks cumulative incidence functions.
`kk_finegray()`	Fine–Gray subdistribution hazards model for competing risks.
`kk_survival_nnt()`	Time-varying NNT/ARR from survival curves at chosen follow-up times.
`kk_partsa()`	Partitioned survival analysis (oncology cost-effectiveness modelling).

Example:

library(survival)
kk_survival_plot(lung, time, status, sex)
kk_rmst(lung, time, status, sex, tau = 365)

The descriptive backbone of every paper.

Function	What it does
`kk_table1()` / `table1_summary()`	The standard “Table 1” of baseline characteristics by group, with correct parametric/non-parametric tests chosen automatically.
`kk_compare_groups_table()`	Detailed comparison: per-group summaries, differences, CIs, effect sizes, p-values — fully tidyselect and `group_by()` aware.
`kk_summary()`	Extensive numeric summary — central tendency, dispersion, robust (Huber-M) estimators, skewness/kurtosis, normality tests.
`kk_chisq_test()`	Chi-square test for r×c contingency tables with effect sizes.
`kk_median_test()`	Median test across two or more groups.
`kk_jonckheere_test()`	Jonckheere–Terpstra trend test for ordered groups.
`kk_vdw_test()`	van der Waerden normal-scores test — a robust non-parametric alternative to one-way ANOVA.

Example:

kk_table1(cohort, by = "arm", variables = c("age", "bmi", "sex", "smoker"))

cohort %>%
  group_by(sex) %>%
  kk_compare_groups_table(arm, c(sbp, dbp, ldl))

Tests for serial independence and randomness in ordered sequences.

Function	What it does
`kk_runs_test()` / `kk_random_seq()`	Runs test for randomness in a binary or numeric sequence.
`kk_frequency_test()`	Frequency (monobit) test — do the categories occur in expected proportions?
`kk_mssd_test()`	Mean successive squared difference test for serial independence.

Modern tools for estimating treatment effects from observational data.

Function	What it does
`kk_iptw()`	Inverse-probability-of-treatment weighting (ATE/ATT) with stabilized weights and robust CIs.
`kk_tmle()`	Targeted maximum likelihood estimation of the ATE — double-robust, with efficient-influence-curve inference.
`kk_smd()` / `kk_balance_table()`	Standardized mean differences and a covariate-balance table to check weighting/matching.
`kk_causal_mediation()`	Causal mediation analysis: natural direct and indirect effects.
`kk_simex()`	SIMEX correction for classical measurement error in a continuous covariate.
`kk_butler_ks()`	Butler’s symmetry-based KS randomization test of the sharp null in matched/stratified experiments.
`kk_gsf()`	Generalized sensitivity functions — how model output responds to perturbations in parameters over time.

Example — double-robust average treatment effect:

kk_tmle(cohort,
        outcome    = "death",
        treatment  = "new_drug",
        covariates = c("age", "sex", "comorbidity", "severity"))

Important

Every function on this tab assumes you measured the confounders that matter — an assumption the data can never verify. When the treatment was switched on by a rule instead of by a clinician, the next tab gets you an effect estimate that does not need it.

New in 1.3.0. When an intervention was assigned by a rule — a policy start date, an eligibility threshold, a staggered rollout, a lottery — that rule can identify the causal effect on its own. No claim about having measured every confounder is required. What each design does require is its own identifying assumption, so every function here returns the diagnostic for it.

Function	What it does
`kk_did()`	Difference-in-differences: the change in the treated group net of the change in the control group. Canonical 2×2, or two-way fixed effects with cluster-robust SEs when `unit`/`time` are supplied. Returns the raw 2×2 cell means beside the regression estimate.
`kk_event_study()`	One coefficient per period relative to adoption — the parallel-trends check a 2×2 design cannot provide. Flat leads are the evidence; the joint Wald test over them comes attached.
`kk_did_staggered()`	Callaway–Sant’Anna group-time ATT for staggered adoption, where two-way fixed effects uses already-treated units as controls and can flip sign. Bootstrap SEs, simultaneous bands, and dynamic/cohort/calendar aggregations.
`kk_rdd()`	Regression discontinuity, sharp or fuzzy: the jump in the outcome at a threshold, by local polynomial regression with the MSE-optimal Imbens–Kalyanaraman bandwidth and a bandwidth-sensitivity table.
`kk_rkd()`	Regression kink design — the effect read off a change in slope at the threshold rather than a jump, as when a benefit formula bends at an income limit.
`kk_rd_density()`	McCrary manipulation test: did units sort themselves across the cutoff? A null result is a precondition for reading an RD estimate causally.
`kk_its()`	Interrupted time series by segmented regression — immediate level change separated from the change in trend, with Newey–West standard errors, seasonality, count models, and a counterfactual series with its cumulative impact.
`kk_synth()`	Synthetic control: builds the comparator for a single treated unit as a non-negative, sum-to-one weighted average of donors fitted to the pre-intervention period, with placebo permutation inference.
`kk_iv()`	Two-stage least squares — the complier effect (LATE), with the first-stage F, Sargan and Wu–Hausman tests, and an Anderson–Rubin interval that stays valid when the instrument is weak.
`kk_ancova()`	Pretest–posttest baseline adjustment; beats both the follow-up-means comparison and the change score, and reports them alongside for contrast.

Choosing the design. Treated and control groups, before and after → kk_did(). Units adopting at different times → kk_did_staggered(). One series, no control group → kk_its(). One treated region, several comparators → kk_synth(). Assignment by a threshold → kk_rdd(). Non-compliance, or an encouragement or lottery → kk_iv().

Example — a policy evaluated three ways:

# 1. Difference-in-differences, clustered on the panel unit
kk_did(panel, cases, treated, post, unit = region, time = year)

# 2. Is the parallel-trends assumption credible? Look at the leads.
es <- kk_event_study(panel, cases, region, year, adoption_year)
attr(es, "pretrend")     # joint Wald test over all pre-periods

# 3. Staggered rollout: one clean 2x2 per cohort-period, then aggregated
res <- kk_did_staggered(panel, cases, region, year, adoption_year)
res$dynamic              # effect by time since adoption
res$overall              # single summary

Note

Interpretation. A difference-in-differences estimate is only as good as the parallel-trends assumption behind it, and that assumption is not testable from a 2×2 table — which is why kk_event_study() exists and why its pre-period coefficients belong in the paper, not in a drawer. The same discipline applies across the tab: report the bandwidth-sensitivity table with an RD estimate, the manipulation test with the design, the first-stage F with an IV, and the placebo distribution with a synthetic control.

A full health-technology-assessment workflow, from ICERs to value of information.

Function	What it does
`kk_icer()`	Incremental cost-effectiveness analysis over strategies: ranks by effect, flags (extended) dominance, traces the efficiency frontier.
`kk_nmb()`	Net monetary and net health benefit at one or more willingness-to-pay thresholds, flagging the optimal strategy.
`kk_ceac()`	Cost-effectiveness acceptability curve from PSA draws.
`kk_markov()`	Markov cohort model for cost-utility analysis (time-varying transitions, discounting, half-cycle correction).
`kk_partsa()`	Partitioned survival model for oncology CEA.
`kk_evpi()`	Expected Value of Perfect Information from PSA.
`kk_discount()`	Discounting of future costs/effects.

Example — rank strategies and find the frontier:

strategies %>% kk_icer(cost, qaly, strategy)
strategies %>% kk_nmb(cost, qaly, wtp = 25000)

Compartmental models, R₀, serological incidence, and outbreak detection.

Function	What it does
`kk_seir()`	Simulate SIR / SEIR dynamics with optional vital dynamics and vaccination; returns compartments, incidence, and implied R₀.
`kk_r0()`	Estimate the basic reproduction number from transmission parameters, early growth rate, or final attack rate.
`kk_final_size()`	Final epidemic size from R₀.
`kk_sero_incidence()`	Force of infection from two or more seroprevalence surveys (birth-cohort catalytic model), or from a single cross-sectional survey.
`kk_nb_scan()`	Negative-binomial space–time scan statistic for detecting localized disease clusters under overdispersion.

Example — simulate an outbreak and read off R₀:

out <- kk_seir(beta = 0.4, gamma = 0.2, sigma = 0.25,
               S0 = 1e6, I0 = 10, times = 0:180)
attr(out, "R0")

The theming layer, quick plots, and the Bulgarian-registry helpers.

Function	What it does
`kkplot()` / `set_plot_font()` / `set_plot_colors()`	The house ggplot2 theme and font/colour setup used across my figures.
`scale_color_kk()` / `scale_fill_kk()` (+ `_c` continuous)	Discrete and continuous colour/fill scales in the kkstatfun palette.
`kk_pal()` / `kk_gen_palettes()` / `kk_show_palettes()`	The raw palette, palette generation, and an interactive palette gallery.
`kk_risk_plot()`	Forest plot of OR/RR estimates with CIs straight from `kk_twobytwo()` output.
`kk_fullcorplot()`	Full correlation-matrix plot for a data frame.
`univariate_plot()`	Quick distribution plot for a single variable (auto categorical/continuous).
`kk_time_series()`	End-to-end time-series analysis: decomposition, ADF/KPSS stationarity, autocorrelation, auto-ARIMA.
`kk_time_metrics()`	Summary metrics for a time series.
`kkonehot()` / `one_hot_encode()`	One-hot encoding of categorical columns.
`mutate_round()` / `format_tibble()`	Column-wise rounding and consistent tibble formatting.
`kk_setup()`	Session setup: parallel cores and sane scientific-notation defaults.
`extract_egn_info()` / `extract_age_from_egn()`	Parse Bulgarian EGN identifiers → date of birth, sex, age, birth region.

Example — parse Bulgarian personal identifiers:

extract_egn_info(c("9201014321", "8812128765"))
#> Date of birth, sex, age, and birth region for each EGN

New in 1.2.0. Roughly one man in twelve has a red–green colour deficiency, journals and university accessibility policies increasingly ask for the check, and a figure still has to survive being photocopied. These audit a palette before it reaches a manuscript — no extra dependencies, base grDevices only.

Function	What it does
`kk_pal_check()`	The one to reach for: for every colour and every vision type, the perceptual distance (OKLab ΔE) to the nearest other colour plus the WCAG contrast against the background. An empty result means the palette passes.
`kk_show_cvd()`	Swatch grid, one row per vision type — categories that collapse into each other become obvious at a glance.
`kk_cvd()`	Simulate deuteranopia, protanopia, tritanopia or achromatopsia (Machado et al. 2009, applied in linear-light RGB), with interpolated severity for the common mild forms.
`kk_pal_safe()`	Build a categorical palette that stays distinguishable under normal, deuteranopic, protanopic and tritanopic vision at once — deterministic, so figures stay reproducible, and anchorable on institutional colours.
`kk_contrast()`	WCAG 2.1 contrast ratios with the AA/AAA flags.
`kk_color_convert()`	OKLab / OKLCh coordinates and relative luminance.

kk_gen_palettes() and set_plot_colors() also take space = "oklch", which builds palettes in perceptually uniform space: hue rotations hold their apparent lightness instead of lurching between a muddy green and a fluorescent one, and out-of-gamut colours are mapped by reducing chroma rather than clipping RGB, which shifts hue.

Example — audit a palette before the figures ship:

pal <- kk_pal(5)
kk_pal_check(pal) %>% filter(!distinct | !graphic)   # empty = passes

# Or build one that is safe by construction
kk_pal_safe(5, seed_colors = "#003049")

Technical notes

Philosophy — tidy data in, tidy data out; full integration with dplyr verbs and genuinely stratified group_by() analyses.
Built on — the best-in-class engines under a single interface: survival, gtsummary, broom, rstatix, pROC, epiR/epitools, easystats, marginaleffects, sandwich, and more.
Methods — implemented according to standard frameworks: Rothman, Greenland & Lash’s Modern Epidemiology (3rd ed.), Szklo & Nieto’s Epidemiology: Beyond the Basics, Jewell’s Statistics for Epidemiology, Zar’s Biostatistical Analysis, and Sheskin’s Handbook of Parametric and Nonparametric Statistical Procedures. The quasi-experimental designs follow Angrist & Pischke, Callaway & Sant’Anna (2021), Imbens & Kalyanaraman (2012), McCrary (2008), Abadie et al. (2010) and Bernal et al. (2017).
Tested — every exported function is covered by the test suite, and the causal estimators are tested against a known simulated effect rather than merely for running. R CMD check passes with 0 errors and 0 warnings.
Reproducible — versioned and archived on Zenodo with a citable concept DOI that always resolves to the latest release.

Citation

If kkstatfun supports your work, please cite it:

@Software{kkstatfun,
  title   = {kkstatfun: R Statistical Analysis Toolkit for Medical
             Statistics and Epidemiology},
  author  = {Kostadinov, Kostadin},
  year    = {2026},
  version = {1.3.0},
  url     = {https://github.com/kostadinoff/kkstatfun},
  doi     = {10.5281/zenodo.18936019},
}

The DOI above is the Zenodo concept DOI — it always resolves to the most recent release. To cite the exact version you used, take the version-specific DOI from the Zenodo record, or run citation("kkstatfun") in R.

View the source on GitHub Archived release (Zenodo)