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, health-technology assessment, and infectious-disease modelling — into 100+ 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.

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: Brier score, Hosmer–Lemeshow test, decile table of predicted vs. observed.
`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.

Example — logistic regression returning tidy odds ratios:

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

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"))

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

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.
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.0.1},
  url     = {https://github.com/kostadinoff/kkstatfun},
  doi     = {10.5281/zenodo.21599749},
}

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)