Dynamic Diastolic Loading Index (DDLI): Reproducible Analysis
Canonical p = 1 form, DDLI = 4 * E^2 / (PEWAT * SVI)
Davis JW et al.
June 22, 2026
This document reproduces every quantitative analysis, table, and
figure in the DDLI manuscript from the single source-of-truth CSV. It is
fully self-contained: it reads only the raw analytic columns
(ev, PEWAT_ms, svi,
ee, grade, lap) and recomputes
DDLI from scratch under the canonical p = 1 form. The
stored DDLI/logDDLI columns in the CSV encode
the superseded p = 0.5 sqrt(PEWAT) form and are deliberately
ignored.
Conventions enforced throughout: seed = 42; bootstrap
B = 2000 with percentile CIs, each iteration refit and
evaluated on the same resample; AUROC and c-index to 3 decimals, ORs to
2, E/e’ and SVI to 1, medians and raw DDLI as integers, log values and
Spearman rho to 2.
Descriptive, non-quantitative figures (Figure 1, the PEWAT measurement schematic; Figure 2, the pressure-volume loop) are not reproduced here.
Setup and helper functions
library(dplyr)
library(tidyr)
library(ggplot2)
library(pROC) # ROC curve geometry and plotting
library(MASS) # polr (proportional-odds ordinal regression)
library(brant) # proportional-odds (Brant) test
library(logistf) # Firth-penalized logistic regression
library(DiagrammeR) # cohort-flow diagram
library(knitr)
library(kableExtra)
SEED <- 42
B <- 2000
# --- formatting helpers (project significant-figure rules) ---
f3 <- function(x) formatC(x, format = "f", digits = 3)
f2 <- function(x) formatC(x, format = "f", digits = 2)
f1 <- function(x) formatC(x, format = "f", digits = 1)
ci3 <- function(e, lo, hi) sprintf("%s (%s-%s)", f3(e), f3(lo), f3(hi))
ci2 <- function(e, lo, hi) sprintf("%s (%s-%s)", f2(e), f2(lo), f2(hi))
# --- AUROC via the Mann-Whitney U statistic (ties handled by mid-ranks) ---
fast_auc <- function(score, y) {
y <- as.integer(y); r <- rank(score)
n1 <- sum(y == 1); n0 <- sum(y == 0)
(sum(r[y == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}
# percentile bootstrap CI for a single-predictor AUROC (resample, re-score)
auc_boot <- function(score, y, B = 2000, seed = SEED) {
pt <- fast_auc(score, y)
set.seed(seed); n <- length(y); out <- numeric(B)
for (b in seq_len(B)) {
idx <- sample.int(n, n, replace = TRUE)
out[b] <- if (length(unique(y[idx])) < 2) NA_real_ else fast_auc(score[idx], y[idx])
}
out <- out[is.finite(out)]
c(est = pt, lo = unname(quantile(out, .025)), hi = unname(quantile(out, .975)))
}
# combined-model AUROC: refit logistic on each resample, evaluate on that resample
comb_auc_boot <- function(data, outcome, preds = c("ee", "logDDLI"),
B = 2000, seed = SEED) {
f <- reformulate(preds, response = outcome)
m <- glm(f, data = data, family = binomial)
pt <- fast_auc(predict(m, type = "response"), data[[outcome]])
set.seed(seed); n <- nrow(data); out <- numeric(B)
for (b in seq_len(B)) {
idx <- sample.int(n, n, replace = TRUE); dd <- data[idx, ]
if (length(unique(dd[[outcome]])) < 2) { out[b] <- NA_real_; next }
mb <- tryCatch(suppressWarnings(glm(f, data = dd, family = binomial)),
error = function(e) NULL)
out[b] <- if (is.null(mb)) NA_real_ else
fast_auc(predict(mb, newdata = dd, type = "response"), dd[[outcome]])
}
out <- out[is.finite(out)]
c(est = pt, lo = unname(quantile(out, .025)), hi = unname(quantile(out, .975)))
}
# Youden-optimal threshold (>= rule) on the raw predictor scale
youden_thr <- function(score, y) {
cuts <- sort(unique(score)); bj <- -Inf; best <- cuts[1]
np <- sum(y == 1); nn <- sum(y == 0)
for (c in cuts) {
pred <- as.integer(score >= c)
j <- sum(pred == 1 & y == 1) / np + sum(pred == 0 & y == 0) / nn - 1
if (j > bj) { bj <- j; best <- c }
}
best
}
# Wilson score interval for a binomial proportion
wilson <- function(x, n, conf = 0.95) {
if (n == 0) return(c(NA, NA))
z <- qnorm(1 - (1 - conf) / 2); p <- x / n; d <- 1 + z^2 / n
ctr <- (p + z^2 / (2 * n)) / d
hw <- z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2)) / d
c(max(0, ctr - hw), min(1, ctr + hw))
}
# rank-concordance (c-index) of a linear predictor against an ordinal outcome
c_index_ord <- function(lp, y) {
y <- as.integer(y); n <- length(y); num <- 0; den <- 0
for (i in 1:(n - 1)) for (j in (i + 1):n) {
if (y[i] == y[j]) next
den <- den + 1
if (lp[i] == lp[j]) num <- num + 0.5
else if ((y[i] < y[j]) == (lp[i] < lp[j])) num <- num + 1
}
num / den
}Data and cohort derivation
DDLI is computed under the canonical p = 1 form. E/e’ is kept on the raw scale. DDLI is log-transformed for all regression modeling. Cohort membership is taken from the pre-built flags in the CSV; counts are asserted against the locked design (82 valid DDLI -> 80 LAP / 66 Grade).
# Set this to wherever the CSV lives. By default the document looks for it
# beside this .Rmd and in the current working directory.
csv_name <- "Diastology_FINAL_ANALYTIC_corrected_v4_sqrtAT.csv"
rmd_dir <- tryCatch(dirname(knitr::current_input(dir = TRUE)), error = function(e) NA)
candidates <- c(csv_name, if (!is.na(rmd_dir)) file.path(rmd_dir, csv_name))
csv_path <- candidates[file.exists(candidates)][1]
if (is.na(csv_path))
stop("CSV not found. Set 'csv_name' to the full path of ", csv_name,
" (looked in: ", paste(candidates, collapse = ", "), ")")
raw <- read.csv(csv_path, check.names = FALSE, stringsAsFactors = FALSE)
is_true <- function(x) x %in% c("True", "TRUE", TRUE, "1", 1)
d <- raw %>%
mutate(
DDLI = 4 * ev^2 / (PEWAT_ms * svi), # canonical p = 1
logDDLI = log(DDLI),
in_DDLI = is_true(valid_DDLI),
in_LAP = is_true(LAP_cohort),
in_Grade = is_true(Grade_cohort)
)
lap <- d %>% filter(in_LAP)
grd <- d %>% filter(in_Grade) %>%
mutate(g2 = as.integer(grade >= 2),
g3 = as.integer(grade == 3),
gc = factor(ifelse(grade <= 1, "0-1", ifelse(grade == 2, "2", "3")),
levels = c("0-1", "2", "3")))
n_valid <- sum(d$in_DDLI); n_lap <- nrow(lap); n_grd <- nrow(grd)
stopifnot(n_valid == 82, n_lap == 80, n_grd == 66)
c(valid_DDLI = n_valid, LAP = n_lap, LAP_elevated = sum(lap$lap == 1),
Grade = n_grd, grade_ge2 = sum(grd$g2), grade3 = sum(grd$g3))## valid_DDLI LAP LAP_elevated Grade grade_ge2 grade3
## 82 80 33 66 22 11table(grd$grade)##
## 0 1 2 3
## 24 20 11 11Supplementary Figure S1. Cohort derivation
Counts only, per the project’s standing rule (no exclusion mechanisms in summary output). Screening-level counts (261 tercile-sampled echocardiograms, 20 enrichment, 199 not abstracted) are documented study-flow constants; the 82 / 80 / 66 analytic counts are derived above and substituted into the diagram.
flow <- sprintf('
digraph cohort {
graph [layout = dot, rankdir = TB, fontname = "Times New Roman"]
node [shape = box, style = filled, fontname = "Times New Roman",
fontsize = 11, color = "#55667F", penwidth = 1.2]
edge [color = "#55667F"]
tercile [label = "E/e\u2032 Tercile-Based Echos\\nN = 261", fillcolor = "#DCE3EE"]
enrich [label = "Enrichment Sample\\nN = 20", fillcolor = "#E3EEDC"]
valid [label = "Valid DDLI\\nn = %d", fillcolor = "#FFFFFF"]
lapc [label = "Valid LAP\\nn = %d", fillcolor = "#1C3354", fontcolor = white]
grdc [label = "Valid ASE DD Grade\\nn = %d", fillcolor = "#8B1A1A", fontcolor = white]
excl [label = "Not abstracted /\\nineligible\\nn = %d", fillcolor = "#F2F2F2"]
tercile -> valid
enrich -> valid
valid -> lapc
valid -> grdc
valid -> excl [style = dashed]
{ rank = same; lapc; grdc }
}', n_valid, n_lap, n_grd, as.integer(281 - n_valid))
grViz(flow)Table 1. Baseline clinical and echocardiographic characteristics (LAP cohort)
ckd_stage <- suppressWarnings(as.numeric(substr(trimws(lap[["CKD Stage"]]), 1, 1)))
bnp <- suppressWarnings(as.numeric(as.character(lap[["Last Pro-BNP"]])))
med_iqr_i <- function(x) { x <- x[!is.na(x)]; sprintf("%d (%d-%d)",
as.integer(round(median(x))), as.integer(round(quantile(x, .25))),
as.integer(round(quantile(x, .75)))) }
med_iqr_1 <- function(x) { x <- x[!is.na(x)]; sprintf("%s (%s-%s)",
f1(median(x)), f1(quantile(x, .25)), f1(quantile(x, .75))) }
np <- function(k) sprintf("%d (%d)", as.integer(k), as.integer(round(100 * k / n_lap)))
t1 <- tribble(
~Characteristic, ~`Value (n = 80)`,
"Age, y, median (IQR)", med_iqr_i(lap$Age),
"Male, n (%)", np(sum(lap$Sex == "Male")),
"Hypertension, n (%)", np(sum(lap[["History HTN?"]] == "Yes", na.rm = TRUE)),
"Ischemic heart disease, n (%)", np(sum(is_true(lap[["History of Ischemic Heart Disease? CAD, PCI, CABG History"]]))),
"CKD stage >= 3, n (%)", np(sum(ckd_stage >= 3, na.rm = TRUE)),
"Atrial fibrillation history, n (%)", np(sum(is_true(lap[["Atrial Fibrillation History?"]]))),
"LVEF, %, median (IQR)", med_iqr_i(lap[["TTE Ejection Fraction"]]),
"E/e' (average), median (IQR)", med_iqr_1(lap$ee),
"DDLI, median (IQR)", med_iqr_i(lap$DDLI),
"NT-proBNP, pg/mL, median (IQR), n = 49", med_iqr_i(bnp),
"Elevated LAP, n (%)", sprintf("%d (%s)", sum(lap$lap == 1), f1(100 * mean(lap$lap == 1)))
)
kable(t1, align = c("l", "r")) %>% kable_styling(full_width = FALSE)| Characteristic | Value (n = 80) |
|---|---|
| Age, y, median (IQR) | 74 (62-82) |
| Male, n (%) | 37 (46) |
| Hypertension, n (%) | 65 (81) |
| Ischemic heart disease, n (%) | 45 (56) |
| CKD stage >= 3, n (%) | 31 (39) |
| Atrial fibrillation history, n (%) | 32 (40) |
| LVEF, %, median (IQR) | 60 (54-65) |
| E/e’ (average), median (IQR) | 13.0 (10.0-16.4) |
| DDLI, median (IQR) | 18 (10-33) |
| NT-proBNP, pg/mL, median (IQR), n = 49 | 1327 (369-5275) |
| Elevated LAP, n (%) | 33 (41.2) |
Distribution, normality, and correlation
sw <- shapiro.test(lap$logDDLI)
mw <- wilcox.test(logDDLI ~ lap, data = lap)
rho_ee_lap <- cor(lap$logDDLI, lap$ee, method = "spearman")
rho_ee_grd <- cor(grd$logDDLI, grd$ee, method = "spearman")
rho_lap <- cor(lap$logDDLI, lap$lap, method = "spearman")
rho_grade <- cor(grd$logDDLI, as.integer(grd$grade), method = "spearman")
data.frame(
Statistic = c("Shapiro-Wilk W (log-DDLI, LAP)", "Shapiro-Wilk P",
"Mann-Whitney P (log-DDLI by LAP)",
"Spearman rho: DDLI~E/e' (LAP)", "Spearman rho: DDLI~E/e' (Grade)",
"Spearman rho: DDLI~LAP", "Spearman rho: DDLI~ordinal grade"),
Value = c(f2(sw$statistic), formatC(sw$p.value, format = "f", digits = 2),
formatC(mw$p.value, format = "e", digits = 1),
f2(rho_ee_lap), f2(rho_ee_grd), f2(rho_lap), f2(rho_grade))
)grd %>% group_by(grade) %>%
summarise(n = n(), median_DDLI = round(median(DDLI)), .groups = "drop")theme_ddli <- theme_minimal(base_size = 12) +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", size = 12))
pal <- c(navy = "#1C3354", maroon = "#8B1A1A", gold = "#C68B00", slate = "#55667F")Supplementary Figure S2. Distribution of DDLI by LAP
meds <- lap %>% group_by(lap) %>% summarise(m = median(logDDLI), .groups = "drop")
ggplot(lap, aes(factor(lap, labels = c("Non-elevated", "Elevated")), logDDLI)) +
geom_boxplot(width = .55, outlier.shape = NA, fill = "grey92") +
geom_jitter(width = .12, alpha = .6, color = pal["navy"]) +
geom_segment(data = meds, aes(x = lap + 0.6, xend = lap + 1.4,
y = m, yend = m), linetype = "dashed", color = pal["maroon"]) +
labs(x = "Adjudicated LAP", y = "log-DDLI",
title = "DDLI by LAP status",
subtitle = sprintf("Mann-Whitney P = %s", formatC(mw$p.value, format = "e", digits = 1))) +
theme_ddli
Supplementary Figure S3. Correlation of DDLI with E/e’ (LAP cohort)
ggplot(lap, aes(ee, logDDLI)) +
geom_point(alpha = .65, color = pal["navy"]) +
geom_smooth(method = "lm", se = FALSE, color = pal["maroon"], linewidth = .7) +
labs(x = "E/e' (average)", y = "log-DDLI",
title = "DDLI vs E/e', LAP cohort",
subtitle = sprintf("Spearman rho = %s, P < 0.001", f2(rho_ee_lap))) +
theme_ddli
Supplementary Figure S4. DDLI by ASE diastolic dysfunction grade
ggplot(grd, aes(factor(grade), DDLI)) +
geom_boxplot(width = .6, outlier.shape = NA, fill = "grey92") +
geom_jitter(width = .12, alpha = .6, color = pal["navy"]) +
scale_y_log10() +
labs(x = "ASE diastolic dysfunction grade", y = "DDLI (log scale)",
title = "DDLI across ASE grade strata",
subtitle = sprintf("Spearman rho vs ordinal grade = %s, P < 0.001", f2(rho_grade))) +
theme_ddli
Table 2. Binary discrimination (AUROC, bootstrap 95% CI)
lap_ee <- auc_boot(lap$ee, lap$lap)
lap_dd <- auc_boot(lap$logDDLI, lap$lap)
lap_cb <- comb_auc_boot(lap, "lap")
grd_ee <- auc_boot(grd$ee, grd$g2)
grd_dd <- auc_boot(grd$logDDLI, grd$g2)
grd_cb <- comb_auc_boot(grd, "g2")
t2 <- tribble(
~Model, ~`Elevated LAP (n = 80)`, ~`ASE grade >= 2 (n = 66)`,
"E/e' alone", ci3(lap_ee[1], lap_ee[2], lap_ee[3]), ci3(grd_ee[1], grd_ee[2], grd_ee[3]),
"DDLI alone", ci3(lap_dd[1], lap_dd[2], lap_dd[3]), ci3(grd_dd[1], grd_dd[2], grd_dd[3]),
"E/e' + DDLI (combined)", ci3(lap_cb[1], lap_cb[2], lap_cb[3]), ci3(grd_cb[1], grd_cb[2], grd_cb[3])
)
kable(t2, align = c("l", "c", "c")) %>% kable_styling(full_width = FALSE)| Model | Elevated LAP (n = 80) | ASE grade >= 2 (n = 66) |
|---|---|---|
| E/e’ alone | 0.949 (0.895-0.988) | 0.942 (0.868-0.991) |
| DDLI alone | 0.954 (0.898-0.993) | 0.941 (0.854-0.994) |
| E/e’ + DDLI (combined) | 0.968 (0.922-0.999) | 0.958 (0.889-1.000) |
Figure 3. ROC curves for elevated LAP
roc_df <- function(score, y, label, panel) {
r <- roc(y, score, quiet = TRUE, levels = c(0, 1), direction = "<")
data.frame(fpr = 1 - r$specificities, tpr = r$sensitivities,
model = label, panel = panel)
}
comb_score <- function(data, outcome) {
m <- glm(reformulate(c("ee", "logDDLI"), outcome), data, family = binomial)
predict(m, type = "response")
}
roc_figure <- function(df, title, sub) {
ggplot(df, aes(fpr, tpr, color = model)) +
geom_abline(linetype = "dotted", color = "grey60") +
geom_path(linewidth = .8) +
facet_wrap(~ panel) +
scale_color_manual(values = unname(pal[c("slate", "navy", "maroon")])) +
coord_equal() +
labs(x = "1 - specificity", y = "Sensitivity", color = NULL,
title = title, subtitle = sub) +
theme_ddli + theme(legend.position = "bottom")
}df3 <- rbind(
roc_df(lap$ee, lap$lap, sprintf("E/e' (%s)", f3(lap_ee[1])), "Individual predictors"),
roc_df(lap$logDDLI, lap$lap, sprintf("DDLI (%s)", f3(lap_dd[1])), "Individual predictors"),
roc_df(lap$ee, lap$lap, sprintf("E/e' (%s)", f3(lap_ee[1])), "Combined vs E/e'"),
roc_df(comb_score(lap,"lap"), lap$lap, sprintf("Combined (%s)", f3(lap_cb[1])), "Combined vs E/e'")
)
roc_figure(df3, "ROC: elevated LAP", "n = 80 (33 cases, 47 controls); 2,000-iteration bootstrap CIs")
Figure 4. ROC curves for ASE grade >= 2
df4 <- rbind(
roc_df(grd$ee, grd$g2, sprintf("E/e' (%s)", f3(grd_ee[1])), "Individual predictors"),
roc_df(grd$logDDLI, grd$g2, sprintf("DDLI (%s)", f3(grd_dd[1])), "Individual predictors"),
roc_df(grd$ee, grd$g2, sprintf("E/e' (%s)", f3(grd_ee[1])), "Combined vs E/e'"),
roc_df(comb_score(grd,"g2"), grd$g2, sprintf("Combined (%s)", f3(grd_cb[1])), "Combined vs E/e'")
)
roc_figure(df4, "ROC: ASE grade >= 2", "n = 66 (22 cases, 44 controls); 2,000-iteration bootstrap CIs")
Supplementary Table S1. Threshold-specific operating characteristics
Youden-optimal thresholds are derived in-sample; apparent estimates carry Wilson 95% CIs. Optimism-corrected estimates subtract the mean bootstrap optimism (threshold re-derived on each resample, evaluated on the original sample) over 2000 resamples, and carry a Wilson 95% CI computed on the optimism-corrected proportion.
op_chars <- function(score, y, thr) {
pred <- as.integer(score >= thr)
TP <- sum(pred == 1 & y == 1); FP <- sum(pred == 1 & y == 0)
TN <- sum(pred == 0 & y == 0); FN <- sum(pred == 0 & y == 1)
c(Sensitivity = TP / (TP + FN), Specificity = TN / (TN + FP),
PPV = TP / (TP + FP), NPV = TN / (TN + FN))
}
op_counts <- function(score, y, thr) {
pred <- as.integer(score >= thr)
list(TP = sum(pred == 1 & y == 1), FP = sum(pred == 1 & y == 0),
TN = sum(pred == 0 & y == 0), FN = sum(pred == 0 & y == 1))
}
# mean optimism for each metric via Efron-Gong bootstrap
optimism <- function(score, y, B = 2000, seed = SEED) {
set.seed(seed); n <- length(y)
acc <- matrix(0, B, 4); colnames(acc) <- c("Sensitivity","Specificity","PPV","NPV")
for (b in seq_len(B)) {
idx <- sample.int(n, n, replace = TRUE)
if (length(unique(y[idx])) < 2) { acc[b, ] <- NA; next }
thr_b <- youden_thr(score[idx], y[idx])
acc[b, ] <- op_chars(score[idx], y[idx], thr_b) - op_chars(score, y, thr_b)
}
colMeans(acc, na.rm = TRUE)
}
build_S1 <- function(score, y, thr) {
app <- op_chars(score, y, thr); cnt <- op_counts(score, y, thr)
opt <- optimism(score, y, B)
denom <- list(Sensitivity = cnt$TP + cnt$FN, Specificity = cnt$TN + cnt$FP,
PPV = cnt$TP + cnt$FP, NPV = cnt$TN + cnt$FN)
num <- list(Sensitivity = cnt$TP, Specificity = cnt$TN, PPV = cnt$TP, NPV = cnt$TN)
out <- lapply(names(app), function(m) {
wa <- wilson(num[[m]], denom[[m]]) # apparent Wilson CI
corr <- max(0, min(1, app[m] - opt[m])) # optimism-corrected point
wc <- wilson(corr * denom[[m]], denom[[m]]) # Wilson CI on corrected proportion
data.frame(Statistic = m,
Apparent = sprintf("%s (%s-%s)", f2(app[m]), f2(wa[1]), f2(wa[2])),
`Optimism-corrected` = sprintf("%s (%s-%s)", f2(corr), f2(wc[1]), f2(wc[2])),
check.names = FALSE)
})
do.call(rbind, out)
}
thr_dd_lap <- youden_thr(lap$logDDLI, lap$lap)
thr_dd_grd <- youden_thr(grd$logDDLI, grd$g2)
thr_ee_lap <- youden_thr(lap$ee, lap$lap)
thr_ee_grd <- youden_thr(grd$ee, grd$g2)
s1 <- rbind(
cbind(Predictor = sprintf("DDLI >= %d (log %s) | LAP", as.integer(round(exp(thr_dd_lap))), f2(thr_dd_lap)),
build_S1(lap$logDDLI, lap$lap, thr_dd_lap)),
cbind(Predictor = sprintf("DDLI >= %d (log %s) | Grade>=2", as.integer(round(exp(thr_dd_grd))), f2(thr_dd_grd)),
build_S1(grd$logDDLI, grd$g2, thr_dd_grd)),
cbind(Predictor = sprintf("E/e' >= %s | LAP", f1(thr_ee_lap)),
build_S1(lap$ee, lap$lap, thr_ee_lap)),
cbind(Predictor = sprintf("E/e' >= %s | Grade>=2", f1(thr_ee_grd)),
build_S1(grd$ee, grd$g2, thr_ee_grd))
)
kable(s1, row.names = FALSE) %>% kable_styling(full_width = FALSE) %>%
collapse_rows(columns = 1, valign = "top")| Predictor | Statistic | Apparent | Optimism-corrected |
|---|---|---|---|
| DDLI >= 21 (log 3.05) | LAP | Sensitivity | 0.88 (0.73-0.95) | 0.85 (0.69-0.93) |
| Specificity | 0.96 (0.86-0.99) | 0.94 (0.84-0.98) | |
| PPV | 0.94 (0.79-0.98) | 0.92 (0.77-0.97) | |
| NPV | 0.92 (0.81-0.97) | 0.90 (0.78-0.96) | |
| DDLI >= 21 (log 3.05) | Grade>=2 | Sensitivity | 0.82 (0.61-0.93) | 0.77 (0.57-0.90) |
| Specificity | 0.98 (0.88-1.00) | 0.96 (0.86-0.99) | |
| PPV | 0.95 (0.75-0.99) | 0.92 (0.72-0.98) | |
| NPV | 0.91 (0.80-0.97) | 0.89 (0.78-0.95) | |
| E/e’ >= 14.0 | LAP | Sensitivity | 0.88 (0.73-0.95) | 0.86 (0.70-0.94) |
| Specificity | 0.94 (0.83-0.98) | 0.93 (0.82-0.97) | |
| PPV | 0.91 (0.76-0.97) | 0.89 (0.74-0.96) | |
| NPV | 0.92 (0.80-0.97) | 0.90 (0.78-0.96) | |
| E/e’ >= 14.0 | Grade>=2 | Sensitivity | 0.86 (0.67-0.95) | 0.83 (0.63-0.94) |
| Specificity | 0.95 (0.85-0.99) | 0.95 (0.84-0.98) | |
| PPV | 0.90 (0.71-0.97) | 0.89 (0.70-0.97) | |
| NPV | 0.93 (0.82-0.98) | 0.92 (0.80-0.97) |
Supplementary Table S2. DDLI misclassifications at the Youden threshold
Cases where the DDLI threshold call disagrees with the adjudicated LAP and/or ASE grade >= 2 outcome, at the full-sample LAP threshold. Generated programmatically under p = 1, so the case set reflects the current formula.
thr <- thr_dd_lap
lap2 <- lap %>%
mutate(call = as.integer(logDDLI >= thr),
miss_lap = call != lap,
g2_self = ifelse(in_Grade, as.integer(grade >= 2), NA_integer_),
miss_g2 = !is.na(g2_self) & call != g2_self,
AF = ifelse(is_true(`Atrial Fibrillation History?`), "Hx AF", "-")) %>%
filter(miss_lap | (miss_g2 %in% TRUE)) %>%
transmute(
LVEF = round(`TTE Ejection Fraction`), E = round(ev), AT = round(PEWAT_ms),
SVI = round(svi, 1), `E/e'` = round(ee, 1), DDLI = round(DDLI),
Grade = ifelse(in_Grade, as.character(grade), "-"),
LAP = ifelse(lap == 1, "Elev.", "Norm."), `AF status` = AF,
`Missed outcome` = case_when(
miss_lap & (miss_g2 %in% TRUE) ~ "LAP + Grade >= 2",
miss_lap ~ "LAP",
TRUE ~ "Grade >= 2")
) %>% arrange(desc(LAP == "Elev."), DDLI)
kable(lap2, row.names = FALSE) %>% kable_styling(full_width = FALSE)| LVEF | E | AT | SVI | E/e’ | DDLI | Grade | LAP | AF status | Missed outcome |
|---|---|---|---|---|---|---|---|---|---|
| 68 | 90 | 65 | 65.2 | 10.0 | 8 | 2 | Elev. |
|
LAP + Grade >= 2 |
| 36 | 109 | 97 | 28.0 | 23.0 | 17 | 3 | Elev. | Hx AF | LAP + Grade >= 2 |
| 64 | 95 | 65 | 33.0 | 14.0 | 17 | 3 | Elev. |
|
LAP + Grade >= 2 |
| 65 | 110 | 56 | 44.0 | 12.2 | 20 | 2 | Elev. | Hx AF | LAP + Grade >= 2 |
| 64 | 95 | 43 | 31.4 | 13.6 | 27 | 0 | Norm. | Hx AF | LAP + Grade >= 2 |
| 68 | 94 | 63 | 20.0 | 14.4 | 28 |
|
Norm. | Hx AF | LAP |
cat(sprintf("LAP Youden threshold: raw DDLI >= %d (log-DDLI >= %s)\n",
as.integer(round(exp(thr))), f2(thr)))## LAP Youden threshold: raw DDLI >= 21 (log-DDLI >= 3.05)Table 3. Ordinal proportional-odds logistic regression (grades 0-1, 2, 3)
Predictors are standardized to z-scores; ORs express the change in odds of a higher grade per SD. The c-index is rank concordance of the linear predictor against observed grade; the Brant test assesses proportional odds.
grd$z_ee <- as.numeric(scale(grd$ee))
grd$z_dd <- as.numeric(scale(grd$logDDLI))
m_ee <- polr(gc ~ z_ee, data = grd, Hess = TRUE, method = "logistic")
m_dd <- polr(gc ~ z_dd, data = grd, Hess = TRUE, method = "logistic")
m_cb <- polr(gc ~ z_dd + z_ee, data = grd, Hess = TRUE, method = "logistic")
or_ci <- function(m, term) {
b <- coef(m)[term]; se <- sqrt(diag(vcov(m)))[term]
sprintf("%s (%s-%s)", f2(exp(b)), f2(exp(b - 1.96 * se)), f2(exp(b + 1.96 * se)))
}
p_term <- function(m, term) {
b <- coef(m)[term]; se <- sqrt(diag(vcov(m)))[term]
pv <- 2 * pnorm(-abs(b / se))
if (pv < 0.001) "< 0.001" else f2(pv)
}
c_ee <- c_index_ord(grd$z_ee * coef(m_ee)["z_ee"], grd$gc)
c_dd <- c_index_ord(grd$z_dd * coef(m_dd)["z_dd"], grd$gc)
c_cb <- c_index_ord(as.numeric(model.matrix(~ z_dd + z_ee, grd)[, -1] %*% coef(m_cb)), grd$gc)
brant_p <- function(m) {
bt <- tryCatch(suppressWarnings(brant::brant(m)), error = function(e) NULL)
if (is.null(bt)) return(NA_character_)
row <- if ("Omnibus" %in% rownames(bt)) "Omnibus" else 1
p <- bt[row, "probability"]
if (p < 0.001) "< 0.001" else f2(p)
}
t3 <- tribble(
~Model, ~`OR per SD (95% CI)`, ~`P value`, ~`c-index`, ~`Brant P`,
"E/e' alone", or_ci(m_ee, "z_ee"), p_term(m_ee, "z_ee"), f3(c_ee), brant_p(m_ee),
"DDLI alone", or_ci(m_dd, "z_dd"), p_term(m_dd, "z_dd"), f3(c_dd), brant_p(m_dd),
"Combined: DDLI", or_ci(m_cb, "z_dd"), p_term(m_cb, "z_dd"), f3(c_cb), brant_p(m_cb),
"Combined: E/e'", or_ci(m_cb, "z_ee"), p_term(m_cb, "z_ee"), "", ""
)## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 11.96 1 0
## z_ee 11.96 1 0
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 5.37 1 0.02
## z_dd 5.37 1 0.02
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 7.59 2 0.02
## z_dd 1.15 1 0.28
## z_ee 2.9 1 0.09
## --------------------------------------------
##
## H0: Parallel Regression Assumption holdskable(t3, align = c("l", "c", "c", "c", "c")) %>% kable_styling(full_width = FALSE)| Model | OR per SD (95% CI) | P value | c-index | Brant P |
|---|---|---|---|---|
| E/e’ alone | 7.58 (2.92-19.67) | < 0.001 | 0.912 | < 0.001 |
| DDLI alone | 10.45 (4.03-27.08) | < 0.001 | 0.904 | 0.02 |
| Combined: DDLI | 6.82 (2.25-20.74) | < 0.001 | 0.908 | 0.02 |
| Combined: E/e’ | 1.75 (0.70-4.38) | 0.23 |
# incremental value and model comparison (likelihood-ratio and AIC)
lr <- function(full, reduced) {
stat <- 2 * (logLik(full) - logLik(reduced))
df <- attr(logLik(full), "df") - attr(logLik(reduced), "df")
c(chisq = as.numeric(stat), df = df, p = pchisq(as.numeric(stat), df, lower.tail = FALSE))
}
data.frame(
Comparison = c("Add E/e' to DDLI", "Add DDLI to E/e'"),
ChiSq = c(f2(lr(m_cb, m_dd)["chisq"]), f2(lr(m_cb, m_ee)["chisq"])),
P = c(f2(lr(m_cb, m_dd)["p"]),
ifelse(lr(m_cb, m_ee)["p"] < 0.001, "< 0.001", f2(lr(m_cb, m_ee)["p"]))),
AIC = c(sprintf("DDLI alone %s", f1(AIC(m_dd))),
sprintf("E/e' alone %s", f1(AIC(m_ee))))
)Supplementary Table S3. Firth-penalized logistic regression (grade 3 vs 0-2)
firth_row <- function(form, label, predcol) {
m <- logistf(form, data = grd)
i <- which(names(coef(m)) == predcol)
lp <- as.numeric(model.matrix(form, grd) %*% coef(m))
auc <- fast_auc(lp, grd$g3)
data.frame(
Model = label,
`OR per SD (95% CI)` = sprintf("%s (%s-%s)", f2(exp(coef(m)[i])),
f2(exp(m$ci.lower[i])), f2(exp(m$ci.upper[i]))),
`P value` = ifelse(m$prob[i] < 0.001, "< 0.001", f3(m$prob[i])),
`AUROC (95% CI)` = f3(auc), check.names = FALSE
)
}
s3 <- rbind(
firth_row(g3 ~ z_ee, "E/e' alone", "z_ee"),
firth_row(g3 ~ z_dd, "DDLI alone", "z_dd"),
firth_row(g3 ~ z_dd + z_ee, "Combined: DDLI", "z_dd"),
firth_row(g3 ~ z_dd + z_ee, "Combined: E/e'", "z_ee")
)
kable(s3, row.names = FALSE) %>% kable_styling(full_width = FALSE)| Model | OR per SD (95% CI) | P value | AUROC (95% CI) |
|---|---|---|---|
| E/e’ alone | 3.45 (1.63-9.45) | < 0.001 | 0.922 |
| DDLI alone | 4.89 (2.20-14.08) | < 0.001 | 0.893 |
| Combined: DDLI | 3.67 (1.33-12.23) | 0.012 | 0.901 |
| Combined: E/e’ | 1.37 (0.66-3.72) | 0.389 | 0.901 |
Sensitivity analysis: deceleration time substituted for PEWAT
Canonical p = 1 form with deceleration time in place of PEWAT, DDLI_DT = 4 * E^2 / (DT * SVI). The “Deceleration time (ms)” column is scaled x100 to true ms and DT = 0 records are set to missing.
d <- d %>% mutate(
DT_ms = ifelse(`Deceleration time (ms)` == 0, NA, `Deceleration time (ms)` * 100),
DDLI_dt = 4 * ev^2 / (DT_ms * svi),
logDD_dt = log(DDLI_dt)
)
lap_dt <- d %>% filter(in_LAP, is.finite(logDD_dt))
grd_dt <- d %>% filter(in_Grade, is.finite(logDD_dt)) %>%
mutate(g2 = as.integer(grade >= 2),
gc = factor(ifelse(grade <= 1, "0-1", ifelse(grade == 2, "2", "3")),
levels = c("0-1", "2", "3")),
z = as.numeric(scale(logDD_dt)))
m_dt <- polr(gc ~ z, data = grd_dt, Hess = TRUE, method = "logistic")
data.frame(
Metric = c("LAP AUROC", "Grade >= 2 AUROC", "Ordinal c-index", "LAP n", "Grade n"),
Value = c(f3(fast_auc(lap_dt$logDD_dt, lap_dt$lap)),
f3(fast_auc(grd_dt$logDD_dt, grd_dt$g2)),
f3(c_index_ord(grd_dt$z * coef(m_dt)["z"], grd_dt$gc)),
nrow(lap_dt), nrow(grd_dt))
)Reproducibility
cat("Seed:", SEED, " Bootstrap B:", B, "\n")## Seed: 42 Bootstrap B: 2000cat("DDLI formula (p=1): 4 * E^2 / (PEWAT_ms * SVI), log-transformed for modeling\n")## DDLI formula (p=1): 4 * E^2 / (PEWAT_ms * SVI), log-transformed for modelingsessionInfo()## R version 4.5.3 (2026-03-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] kableExtra_1.4.0 knitr_1.51 DiagrammeR_1.0.12 logistf_1.26.1
## [5] brant_0.3-0 MASS_7.3-65 pROC_1.19.0.1 ggplot2_4.0.3
## [9] tidyr_1.3.2 dplyr_1.2.1
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 shape_1.4.6.1 xfun_0.57
## [4] bslib_0.10.0 visNetwork_2.1.4 htmlwidgets_1.6.4
## [7] formula.tools_1.7.1 lattice_0.22-9 vctrs_0.7.3
## [10] tools_4.5.3 Rdpack_2.6.6 generics_0.1.4
## [13] tibble_3.3.1 pan_1.9 pkgconfig_2.0.3
## [16] jomo_2.7-6 Matrix_1.7-4 RColorBrewer_1.1-3
## [19] S7_0.2.2 lifecycle_1.0.5 stringr_1.6.0
## [22] compiler_4.5.3 farver_2.1.2 textshaping_1.0.5
## [25] codetools_0.2-20 htmltools_0.5.9 sass_0.4.10
## [28] yaml_2.3.12 glmnet_5.0 mice_3.19.0
## [31] pillar_1.11.1 nloptr_2.2.1 jquerylib_0.1.4
## [34] cachem_1.1.0 reformulas_0.4.4 iterators_1.0.14
## [37] rpart_4.1.24 boot_1.3-32 foreach_1.5.2
## [40] mitml_0.4-5 nlme_3.1-168 tidyselect_1.2.1
## [43] digest_0.6.39 stringi_1.8.7 purrr_1.2.2
## [46] labeling_0.4.3 splines_4.5.3 operator.tools_1.6.3.1
## [49] fastmap_1.2.0 grid_4.5.3 cli_3.6.6
## [52] magrittr_2.0.5 survival_3.8-6 broom_1.0.12
## [55] withr_3.0.2 scales_1.4.0 backports_1.5.1
## [58] rmarkdown_2.31 otel_0.2.0 nnet_7.3-20
## [61] lme4_2.0-1 evaluate_1.0.5 rbibutils_2.4.1
## [64] viridisLite_0.4.3 mgcv_1.9-4 rlang_1.2.0
## [67] Rcpp_1.1.1-1 glue_1.8.1 xml2_1.5.2
## [70] svglite_2.2.2 rstudioapi_0.18.0 minqa_1.2.8
## [73] jsonlite_2.0.0 R6_2.6.1 systemfonts_1.3.2