---
title: "Preterm Bicarbonate Reference Interval"
author: "Millen et al."
date: "30 June 2026"
output:
  pdf_document:
    toc: true
  html_document:
    toc: true
    toc_float: true
    df_print: kable
subtitle: "Data-Driven Physiologic Bicarbonate Range for Preterm Infants"
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
```

Place `bicarbonate_data_request.xlsx` in this directory and knit. Each row is one
bicarbonate sample, with a flag to designate which analyses the sample was a
part of.

| column | meaning |
|---|---|
| `Patient` | randomly assigned number |
| `GA_group` | `< 32 weeks` or `32 - 35 weeks` |
| `HOL` | hours of life at which the bicarbonate was drawn |
| `Bicarbonate` | measured serum bicarbonate, mEq/L |
| `Gas_pH`, `Gas_pCO2` | paired blood gas within ±1 hour, blank if none |
| `Reference_sample` | `Yes` for the one sample per patient used for the reference interval |
| `Gas_confirmed_sample` | `Yes` for the 51 samples in the sensitivity analysis |

```{r setup-data}
library(Hmisc)
library(dplyr)
library(boot)

# Harrell-Davis quantile (Harrell & Davis, Biometrika 1982;69(3):635-40).
hd_quantile <- function(x, p) unname(hdquantile(x, probs = p))

dat <- as.data.frame(readxl::read_excel("bicarbonate_data_request.xlsx",
                                        sheet = "bicarbonate_data_request"))
dat <- dat[!is.na(dat$Patient), ]
ref <- filter(dat, Reference_sample == "Yes")
sub <- filter(dat, Gas_confirmed_sample == "Yes")

c(samples = nrow(dat), patients = n_distinct(dat$Patient),
  reference = nrow(ref), gas_confirmed = nrow(sub))
```

# 1. Estimators

Limits are `mean ± 1.96 SD`; each limit's 90% CI uses the standard error of a
reference limit, `SD * sqrt(1/n + 1.96^2 / (2n))`. 

```{r estimators}
param_row <- function(x, label) {
  x <- as.numeric(na.omit(x)); n <- length(x)
  mu <- mean(x); s <- sd(x)
  se_mean <- s / sqrt(n)
  se_lim  <- s * sqrt(1 / n + 1.96^2 / (2 * n))
  lo <- mu - 1.96 * s; hi <- mu + 1.96 * s
  data.frame(
    Cohort = label, Method = "Parametric", N = n, Mean = mu, SD = s,
    MeanLo  = mu - 1.645 * se_mean, MeanHi  = mu + 1.645 * se_mean,
    Lower   = lo, LowerLo = lo - 1.645 * se_lim, LowerHi = lo + 1.645 * se_lim,
    Upper   = hi, UpperLo = hi - 1.645 * se_lim, UpperHi = hi + 1.645 * se_lim,
    stringsAsFactors = FALSE)
}

# Harrell-Davis limits with 90% BCa bootstrap CIs (10,000 resamples)
hd_row <- function(x, label, R = 10000, seed = 42) {
  x <- as.numeric(na.omit(x))
  stat <- function(d, i) { s <- d[i]
    c(mean(s), hd_quantile(s, 0.025), hd_quantile(s, 0.975)) }
  set.seed(seed)
  b <- boot(data = x, statistic = stat, R = R)
  ci <- function(k) boot.ci(b, conf = 0.90, type = "bca", index = k)$bca[4:5]
  cm <- ci(1); cl <- ci(2); cu <- ci(3)
  data.frame(
    Cohort = label, Method = "Harrell-Davis", N = length(x),
    Mean = b$t0[1], SD = sd(x), MeanLo = cm[1], MeanHi = cm[2],
    Lower = b$t0[2], LowerLo = cl[1], LowerHi = cl[2],
    Upper = b$t0[3], UpperLo = cu[1], UpperHi = cu[2],
    stringsAsFactors = FALSE)
}

fmt <- function(tbl, dp = 1) {
  f <- function(v) sprintf(paste0("%.", dp, "f"), v)
  transmute(tbl, Cohort, Method, n = N,
    `Mean (90% CI)`        = sprintf("%s (%s – %s)", f(Mean),  f(MeanLo),  f(MeanHi)),
    `Lower Limit (90% CI)` = sprintf("%s (%s – %s)", f(Lower), f(LowerLo), f(LowerHi)),
    `Upper Limit (90% CI)` = sprintf("%s (%s – %s)", f(Upper), f(UpperLo), f(UpperHi)))
}
```

# 2. Supplemental Table 3

```{r suppl3}
by_ga <- bind_rows(
  param_row(ref$Bicarbonate,                                  "Reference Cohort (< 35 wk)"),
  param_row(ref$Bicarbonate[ref$GA_group == "32 - 35 weeks"], "Reference Cohort (32-35 wk)"),
  param_row(ref$Bicarbonate[ref$GA_group == "< 32 weeks"],    "Reference Cohort (< 32 wk)"))
primary <- by_ga[1, ]

sensitivity <- bind_rows(
  param_row(sub$Bicarbonate,                                  "Gas-Confirmed Subset (< 35 wk)"),
  hd_row(   sub$Bicarbonate,                                  "Gas-Confirmed Subset (< 35 wk)"),
  param_row(sub$Bicarbonate[sub$GA_group == "32 - 35 weeks"], "Gas-Confirmed Subset (32-35 wk)"),
  param_row(sub$Bicarbonate[sub$GA_group == "< 32 weeks"],    "Gas-Confirmed Subset (< 32 wk)"))

suppl_table_3 <- bind_rows(by_ga, sensitivity) %>% fmt(dp = 1) %>%
  mutate(`Gestational Age` = sub("wk$", "weeks", sub(".*\\((.*)\\)$", "\\1", Cohort)),
         Cohort            = trimws(sub("\\s*\\([^)]*\\)$", "", Cohort))) %>%
  select(Cohort, `Gestational Age`, Method, n, everything())
suppl_table_3
```

# 3. Comparisons

Gestational-age groups within the reference cohort, then the subset against the
reference cohort excluding those infants.

```{r comparisons}
a <- ref$Bicarbonate[ref$GA_group == "< 32 weeks"]
b <- ref$Bicarbonate[ref$GA_group == "32 - 35 weeks"]
tt <- t.test(a, b)
ga_cmp <- data.frame(
  N_lt32 = length(a), Mean_lt32 = round(mean(a), 1),
  N_32_35 = length(b), Mean_32_35 = round(mean(b), 1),
  Mean_Diff = round(mean(a) - mean(b), 2),
  t = round(unname(tt$statistic), 2),
  p = base::format.pval(tt$p.value, digits = 3))

ref_excl <- ref$Bicarbonate[!ref$Patient %in% unique(sub$Patient)]
td <- t.test(sub$Bicarbonate, ref_excl)
disjoint <- data.frame(
  n_subset = nrow(sub), Mean_subset = round(mean(sub$Bicarbonate), 2),
  n_reference = length(ref_excl), Mean_reference = round(mean(ref_excl), 2),
  Mean_Diff = round(mean(sub$Bicarbonate) - mean(ref_excl), 2),
  t = round(unname(td$statistic), 2), p = signif(td$p.value, 3))

list(gestational_age = ga_cmp, disjoint = disjoint)
```