Estimating the average annual healthcare expenditure attributable to diagnosed diabetes among US adults, using the Medical Expenditure Panel Survey (MEPS) 2023 Full Year Consolidated file.
The analysis is written in R with cohort construction in SQL (DuckDB), respects the MEPS complex survey design throughout, and is reproducible from a clean checkout with one command.
Among 14,756 US adults (2,083 with diagnosed diabetes), representing 259.4 million civilian non-institutionalized adults:
Diagnosed diabetes is associated with $3,466 in additional annual healthcare expenditure per person (95% CI $2,350 – $4,591; Rao–Wu subbootstrap, 200 replicates).
Adjusted mean annual expenditure: $11,542 with diabetes vs $8,075 without.
This is an adjusted association, not a causal effect, and not the saving available from preventing a case of diabetes. See Limitations.
Healthcare cost data breaks ordinary regression in two specific ways, and both are present here:
| Feature | In this cohort |
|---|---|
| A point mass at zero | 12.8% of adults spent nothing at all |
| Severe right skew | Mean ($9,763) is 3.8× the median ($2,550); the top 1% exceeds $106,000 and the maximum is $574,675 |
A two-part model addresses both: part one models whether someone spends (logistic), part two models how much given they spend (Gamma GLM with a log link), and the parts multiply:
E[Y | X] = Pr(Y > 0 | X) × E[Y | Y > 0, X]
The estimand is a marginal effect computed by recycled predictions — predict every person twice, once under each exposure value, and average the weighted difference. The exponentiated coefficient is not the answer: in a nonlinear two-part model it is a conditional effect on the log scale, and it is not what a payer is asking for.
All four rows below answer the same question on the same data.
| Estimator | Incremental cost | vs. preferred | Handles zeros | Handles skew | Adjusted |
|---|---|---|---|---|---|
| Unadjusted weighted difference | $10,940 | 3.16× | ✗ | ✗ | ✗ |
| OLS on raw dollars (adjusted) | $3,967 | 1.14× | ✗ | ✗ | ✓ |
| Two-part: logit + log-OLS w/ Duan smearing | $6,716 | 1.94× | ✓ | ✓ | ✓ |
| Two-part: logit + Gamma GLM (preferred) | $3,466 | 1.00× | ✓ | ✓ | ✓ |
Two things worth noting:
- Confounding is most of the raw gap. The unadjusted difference is 3.2× larger than the adjusted estimate. People with diabetes are older, sicker, and more likely to be insured — each independently raising cost.
- Log-OLS nearly doubles the estimate. Retransforming from the log scale with a single Duan smearing factor (here 3.02) assumes homoscedastic log-scale errors. That assumption is rejected in these data (p ≈ 6×10⁻³⁵), and the resulting estimate is inflated to $6,716. This is the concrete argument for modelling the mean on the dollar scale with a GLM.
The Gamma family is not assumed. A modified Park test regresses squared
residuals on the fitted mean to recover λ in Var(Y|X) ∝ μ^λ:
λ = 1.52 → between Poisson (1) and Gamma (2), closer to Gamma
λ = 1.52 sits between the two, so Gamma is the nearest standard family rather than a decisive winner. A Poisson-family (quasi-Poisson) GLM is a reasonable sensitivity analysis and is flagged as planned future work.
Source. MEPS Household Component, Full Year Consolidated File 2023 (AHRQ HC-251).
Access. Public use file. No data use agreement, no IRB, no identifiers.
R/01_download.R fetches it directly; nothing is committed to the repository.
data/raw/MANIFEST.csv records the source URL, file size, MD5, and retrieval
date so the exact input is identifiable later.
Validation. Two independent checks that the survey weights are applied correctly, both enforced as tests:
- Weighted population totals 259.4M adults, against roughly 262M US civilian non-institutionalized adults.
- Weighted diagnosed-diabetes prevalence is 11.65%, against a national rate of roughly 11–12%.
Neither is a modelling result — both are checks that would have caught a wrong weight variable, a dropped stratum, or a miscoded exposure flag.
├── protocol.md Pre-specified study protocol
├── run_all.R Full pipeline, one command
├── sql/
│ └── 01_analytic_table.sql Cohort construction, recoding, comorbidity index
├── R/
│ ├── 00_config.R Year, paths, variable manifest
│ ├── 01_download.R Fetch PUF + provenance manifest
│ ├── 02_build_analytic_file.R PUF → DuckDB → analytic parquet
│ ├── 03_survey_design.R svydesign, covariate set, data prep
│ ├── 04_two_part_model.R Model fits, Park test, comparators
│ ├── 05_incremental_cost.R Recycled predictions + bootstrap CI
│ └── estimator.R Reusable, study-agnostic estimator
└── tests/testthat/ 43 tests
R/estimator.R is deliberately independent of this study: the exposure
variable, its two levels, and both model formulas are arguments, so the same
estimator serves any binary-exposure cost analysis.
Requires R ≥ 4.4.
install.packages(c("survey", "haven", "duckdb", "DBI", "dplyr",
"tidyr", "testthat", "ggplot2"))Rscript run_all.R # full pipeline, ~6 minutes
N_REPLICATES=20 Rscript run_all.R # quick check, ~1 minuteResults land in output/tables/.
Rscript -e 'testthat::test_dir("tests/testthat")'The suite targets failure modes that produce plausible but wrong numbers rather than errors:
- MEPS negative missing codes (
-1,-7,-8,-9,-15) leaking into arithmetic and being averaged as if they were dollars - SQL three-valued logic silently NULLing the cohort flag, so
WHERE in_cohortdrops rows nobody intended to drop - The exposure appearing inside the comorbidity index, partialling out the very effect being estimated
- Weighted population and prevalence drifting away from national benchmarks
- Weight-scale invariance of the estimator, and equality of the normalised-
weight GLM coefficients with
svyglm— the property that justifies the numerical fix described below
Fitting part one with raw MEPS person weights diverges: coefficients reach order 10¹⁴ and fitted probabilities pin to 0 and 1. Binomial IRLS reads prior weights as trial counts, and weights up to ~130,000 across 14,756 rows imply an effective sample size near 2.6×10⁸.
Rescaling the weights to mean 1 fixes it. Scaling every weight by a constant
multiplies the log-likelihood by that constant and leaves its maximiser
unchanged, so this is exact rather than approximate — and it reproduces the
svyglm coefficients to every printed digit. Both claims are asserted in
tests/testthat/test-estimator.R rather than left as a comment.
The full pre-specified list is in protocol.md. The four that
most constrain interpretation:
-
Cross-sectional and non-randomised. Exposure and outcome are measured over the same period. Reverse causation cannot be excluded — costly comorbid illness may precipitate a diagnosis as readily as the reverse.
-
Residual confounding, biasing away from the null. Obesity, smoking, diet, and physical activity cause both diabetes and cost and are not adjusted for. MEPS collects smoking only in a self-administered subsample requiring a different weight, so including it would change the estimand's population.
-
Self-reported, diagnosed diabetes only, biasing toward the null. Undiagnosed cases — roughly a fifth of all cases nationally — are misclassified into the comparison group.
-
MEPS excludes institutionalized people and is known to fall below the National Health Expenditure Accounts benchmark in aggregate. Nursing-home residents are older, sicker, and costlier, so the absolute dollar figure is an underestimate.
Not yet done, and honestly labelled as such: this estimate has not been formally benchmarked against the published cost-of-diabetes literature. Direct comparison requires matching population, cost definition, and price year, and published figures using broader cost definitions are not directly comparable to a MEPS direct-medical-cost estimate.
- Quasi-Poisson sensitivity analysis, given λ = 1.52
- Propensity score weighting (IPTW) as an alternative to regression adjustment, with standardized mean difference balance diagnostics
- Decomposition of the incremental cost by service category (inpatient, office-based, emergency, prescription drugs)
- Repeat across MEPS years for a trend series
Analysis code released under the MIT License. MEPS data is public domain, produced by the US Agency for Healthcare Research and Quality.