In medical diagnostics, a central task is to evaluate how accurately a continuous biomarker distinguishes between diseased and non-diseased individuals. The receiver operating characteristic (ROC) curve and the area under the curve (AUC) are standard tools for this purpose. The ROC curve summarizes the trade-off between sensitivity and specificity over all possible decision thresholds, while the AUC provides a single-number summary of overall discriminative ability (the probability that a randomly chosen diseased subject has a higher biomarker value than a randomly chosen non-diseased subject).
Let \(X\) and \(Y\) denote biomarker values from non-diseased and diseased subjects with cumulative distribution functions (CDFs) \(F\) and \(G\), and survival functions \(\bar F = 1 - F\) and \(\bar G = 1 - G\). The ROC curve can be written as \[ \mathrm{ROC}(p) = \bar G\{\bar F^{-1}(p)\}, \quad p \in [0,1], \] and the AUC as \[ \mathrm{AUC} = \int_{-\infty}^{\infty} F(x)\,dG(x) = P(Y > X). \]
Empirical ROC curves, constructed directly from the empirical CDFs of \(X\) and \(Y\), are fully nonparametric but stepwise and potentially jagged, especially in small or moderate samples. Parametric ROC models are smooth but require strong distributional assumptions that may be unrealistic in practice.
The smoothROC package implements a
kernel-based, distribution-free approach that produces smooth
ROC curves, a kernel-based AUC estimator with confidence intervals, and
a smooth Youden index summary with an associated optimal cutoff. The
core function is smoothROC(), which this vignette
introduces and illustrates.
smoothROC() functionThe main user-facing function is:
smoothROC(
data,
biomarker,
status,
diseased,
kernel = c("gaussian", "biweight", "epanechnikov"),
bw_method = c("pdf", "AL", "PB", "BHP", "AR"),
alpha = 0.05,
logtrans = FALSE,
grid_n = 1000
)It estimates a kernel-smoothed ROC curve for a continuous biomarker using kernel CDF estimators in the non-diseased and diseased groups. It returns:
data
A data frame containing the biomarker and status variables.
biomarker
Character string; name of the numeric column containing biomarker
values.
status
Character string; name of the column containing binary disease
status.
diseased
The value in status indicating the diseased class
(e.g. "carrier").
kernel
Character string; kernel function for smoothing. One of
"gaussian", "biweight", or
"epanechnikov".
bw_method
Character string; bandwidth selection method. One of:
"pdf" – density-based rule-of-thumb (Silverman),"BHP" – CDF-based normal-reference
(Bowman–Hall–Prvan),"AR" – adjusted CDF reference bandwidth,"AL" – Altman–Leger CDF plug-in,"PB" – Polansky–Baker multistage plug-in.alpha
Numeric; significance level for \((1 -
\alpha)\) confidence intervals (default:
0.05).
logtrans
Logical; if TRUE, applies a log-transformation to biomarker
values prior to ROC estimation (useful for right-skewed biomarkers).
Default: FALSE.
grid_n
Integer; number of grid points for evaluating the ROC curve (default:
1000).
smoothROC() returns an object of class
"smoothROC" with components including:
curve
Data frame with columns FPR, TPR,
threshold, and Youden index \(J\), one row per unique false-positive rate
(deduplicated, keeping the maximum TPR at each
FPR), with the corner points \((0,0)\) and \((1,1)\) always included.
AUC, AUC_ci, AUC_ci_lo,
AUC_ci_hi
Kernel-smoothed AUC estimate and its confidence interval.
J, J_ci, J_ci_lo,
J_ci_hi
Kernel-smoothed Youden index estimate and its confidence
interval.
t0
Estimated optimal cutoff associated with the Youden index.
sensitivity, specificity
Sensitivity and specificity at the Youden cutoff.
kernel, bandwidth_method
The chosen kernel and bandwidth selection method.
hX, hY
Selected CDF bandwidths for the non-diseased and diseased
groups.
plot
A ggplot2 ROC plot object including the Youden point and a
textual annotation.
Print, summary, and plot methods are available:
and are invoked automatically via print(),
summary(), and plot().
To obtain a smooth ROC curve while remaining nonparametric,
smoothROC() applies kernel-based CDF estimators of the form
\[
\hat F(x)
= \frac{1}{m} \sum_{i=1}^m K\!\left(\frac{x - X_i}{h_m}\right),
\quad
\hat G(x)
= \frac{1}{n} \sum_{j=1}^n K\!\left(\frac{x - Y_j}{h_n}\right),
\] where \(K(u) = \int_{-\infty}^u
k(v)\,dv\) is the integrated kernel, and \(h_m\), \(h_n\) are bandwidths for the non-diseased
and diseased groups, respectively. The smooth ROC curve is then obtained
by plugging \(\hat F\) and \(\hat G\) into the ROC functional.
Three univariate kernels are implemented:
Gaussian kernel
\(k(u) = (2\pi)^{-1/2}\exp(-u^2/2)\)
with CDF \(K(u)\) equal to the standard
normal distribution function. This kernel has infinite support and is a
default choice in many smoothing problems.
Epanechnikov kernel
\(k(u) =
\tfrac{3}{4}(1-u^2)\mathbf{1}_{\{|u|\le 1\}}\), with compact
support on \([-1,1]\) and optimal
second-order efficiency under many mean squared error criteria.
Biweight kernel
\(k(u) =
\tfrac{15}{16}(1-u^2)^2\mathbf{1}_{\{|u|\le 1\}}\), a
higher-order, compactly supported kernel that produces more rounded
estimates near the boundaries.
Bandwidth selection is critical for balancing bias and variance in
the smoothed CDFs and the resulting ROC curve. smoothROC()
focuses on bandwidths that are optimal for CDF estimation
(rather than densities), which aligns more directly with ROC
functionals.
The bw_method argument implements several
strategies:
"pdf"
A density-based rule-of-thumb (Silverman) using a kernel density
bandwidth. Convenient and widely used, but it does not satisfy the usual
asymptotic conditions for CDF estimation.
"BHP"
A CDF-based normal-reference bandwidth that approximately minimizes the
integrated mean squared error of \(\hat
F\). It uses a robust scale estimate based on \(\min(\mathrm{SD},
\mathrm{IQR}/1.34)\).
"AR"
An adjusted CDF reference bandwidth for the Gaussian kernel, obtained by
shrinking the normal-reference constant to reduce oversmoothing for
non-Gaussian data while preserving the \(m^{-1/3}\) CDF rate.
"AL"
A fully data-driven CDF-based bandwidth in the spirit of Altman and
Leger, where the unknown roughness functional is estimated using an
auxiliary kernel estimator.
"PB"
A multistage plug-in CDF-based bandwidth (a two-stage version of
Polansky and Baker) that uses an initial normal-reference pilot followed
by a data-driven refinement.
In simulation studies (not shown here), the "PB" method
often provides stable performance across a range of underlying
distributions, especially when sample sizes are small to moderate.
Given the kernel CDFs, smoothROC() computes the AUC by
trapezoidal integration of the smoothed ROC curve, \[
\hat \delta
= \int_{-\infty}^{\infty} \hat F(x)\, d\hat G(x)
\approx \sum_i \tfrac{1}{2}\big(\mathrm{TPR}_i +
\mathrm{TPR}_{i-1}\big)\big(\mathrm{FPR}_i - \mathrm{FPR}_{i-1}\big),
\] which is asymptotically equivalent to the empirical AUC based
on the Mann–Whitney statistic. This link ensures that classical
large-sample results for the empirical AUC remain valid in the smoothed
setting.
To quantify uncertainty, smoothROC() implements a
kernel-smoothed analogue of DeLong’s variance estimator. Instead of
using empirical placement values, the variance expression replaces them
by their kernel-smoothed counterparts. This typically yields a more
stable variance estimate in small samples while retaining the
large-sample properties of the original DeLong method. The confidence
interval is then centered at the trapezoidal AUC estimate above, using
this kernel-smoothed DeLong-type standard error.
The resulting AUC estimate and confidence interval are available via:
The Youden index \[ J = \max_t \{\mathrm{sensitivity}(t) + \mathrm{specificity}(t) - 1\} = \max_t \{F(t) - G(t)\} \] provides a summary of the optimal trade-off between sensitivity and specificity, with \(J \in [0,1]\). The corresponding optimal cutoff is \[ t_0 = \operatorname*{arg\,max}_t \{F(t) - G(t)\}. \]
Using the kernel CDFs, smoothROC() computes a smoothed
Youden index \(\hat J\) and the
maximizing cutoff \(\hat t_0\) on a
grid of thresholds. When multiple cutoffs achieve the same maximum,
secondary rules (favoring higher sensitivity or specificity) or
median-based summaries can be used. A Delta-method approximation
provides the variance of \(\hat J\),
from which a Wald-type confidence interval is constructed.
These quantities are returned as:
The confidence intervals returned by smoothROC() for the
AUC and the Youden index are kernel-smoothed DeLong-type variance and
Delta-method, respectively. These asymptotic approximations can be less
reliable in small samples. As an alternative,
smoothROCboot() computes bootstrap percentile confidence
intervals for the AUC, the Youden index, and its optimal cutoff, by
independently resampling the non-diseased and diseased groups with
replacement and recomputing the full smoothed analysis on each
replicate.
boot <- smoothROCboot(
data = dystrophy,
biomarker = "CK",
status = "Class",
diseased = "carrier",
kernel = "biweight",
bw_method = "PB",
alpha = 0.05,
logtrans = TRUE,
B = 1000,
grid_n = 1000,
seed = 1691
)
print(boot)smoothROCboot() shares the same
data/biomarker/status/diseased/
kernel/bw_method/alpha/logtrans/grid_n
arguments as smoothROC(), plus:
B
Number of bootstrap replicates (default: 1000).
seed
Optional integer seed, for reproducible bootstrap resampling.
It returns an object of class "smoothROCboot" with:
boot$J # Youden index estimate
boot$J_ci # Bootstrap percentile CI for J
boot$t0 # Optimal cutoff estimate
boot$t0_ci # Bootstrap percentile CI for the cutoff
boot$AUC # AUC estimate (trapezoidal)
boot$AUC_ci # Bootstrap percentile CI for AUC
boot$sensitivity # Sensitivity at the Youden point (observed data)
boot$specificity # Specificity at the Youden point (observed data)along with the raw bootstrap replicate vectors
boot$boot_J, boot$boot_t0, and
boot$boot_AUC, for users who want to inspect the bootstrap
distribution directly (for example, plotting a histogram of
boot$boot_AUC).
Because each replicate refits bandwidths and re-evaluates the
smoothed ROC curve from scratch, smoothROCboot() is more
computationally intensive than smoothROC();
B = 1000 with grid_n = 1000 is typically a
reasonable default, but both can be reduced for exploratory work.
The package includes an example dataset, dystrophy, with
biomarker measurements for Duchenne muscular dystrophy (DMD) carriers
and non-carriers. We treat the serum marker CK as the
primary biomarker and Class as the disease status.
data(dystrophy)
str(dystrophy)
#> spc_tbl_ [209 × 10] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
#> $ OBS : num [1:209] 1 1 1 1 1 1 1 1 1 1 ...
#> $ HospID: num [1:209] 1007 786 778 1306 895 ...
#> $ AGE : num [1:209] 22 32 36 22 23 30 27 30 25 26 ...
#> $ M : num [1:209] 6 8 7 11 1 5 8 11 10 2 ...
#> $ Y : num [1:209] 79 78 78 79 78 79 78 78 79 79 ...
#> $ CK : num [1:209] 52 20 28 30 40 24 15 22 42 130 ...
#> $ H : num [1:209] 83.5 77 86.5 104 83 78.8 87 91 65.5 80.3 ...
#> $ PK : num [1:209] 10.9 11 13.2 22.6 15.2 9.6 13.5 17.5 13.3 17.1 ...
#> $ LD : num [1:209] 176 200 171 230 205 151 232 198 216 211 ...
#> $ Class : Factor w/ 2 levels "normal","carrier": 1 1 1 1 1 1 1 1 1 1 ...
#> - attr(*, "spec")=List of 3
#> ..$ cols :List of 10
#> .. ..$ OBS : list()
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
#> .. ..$ HospID: list()
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
#> .. ..$ AGE : list()
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
#> .. ..$ M : list()
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
#> .. ..$ Y : list()
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
#> .. ..$ CK : list()
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
#> .. ..$ H : list()
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
#> .. ..$ PK : list()
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
#> .. ..$ LD : list()
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_double" "collector"
#> .. ..$ Class :List of 3
#> .. .. ..$ levels : chr [1:2] "normal" "carrier"
#> .. .. ..$ ordered : logi FALSE
#> .. .. ..$ include_na: logi FALSE
#> .. .. ..- attr(*, "class")= chr [1:2] "collector_factor" "collector"
#> ..$ default: list()
#> .. ..- attr(*, "class")= chr [1:2] "collector_guess" "collector"
#> ..$ delim : chr ","
#> ..- attr(*, "class")= chr "col_spec"
#> - attr(*, "problems")=<externalptr>A basic smooth ROC analysis is:
roc <- smoothROC(
data = dystrophy,
biomarker = "CK",
status = "Class",
diseased = "carrier",
kernel = "biweight",
bw_method = "PB",
alpha = 0.05,
logtrans = TRUE,
grid_n = 1000
)The print and summary methods provide a concise summary:
roc
#> Kernel-smoothed ROC
#> Kernel: biweight
#> Bandwidth: PB
#> AUC = 0.8456 (95.0% CI: 0.7915, 0.8997)
#> Youden Index, J = 0.5724 (95.0% CI: 0.4543, 0.6905)
#> Youden point: (FPR = 0.1305, TPR = 0.7029)
#> Threshold (cutoff) = 4.0707
#> At Youden point Sensitivity = 0.7029, Specificity = 0.8695
summary(roc)
#> Kernel-smoothed ROC
#> Kernel: biweight
#> Bandwidth: PB
#> AUC = 0.8456 (95.0% CI: 0.7915, 0.8997)
#> Youden Index, J = 0.5724 (95.0% CI: 0.4543, 0.6905)
#> Youden point: (FPR = 0.1305, TPR = 0.7029)
#> Threshold (cutoff) = 4.0707
#> At Youden point Sensitivity = 0.7029, Specificity = 0.8695They report the kernel and bandwidth method, AUC with confidence interval, the Youden index and its confidence interval, and the Youden point (FPR, TPR, cutoff, sensitivity, specificity).
We can visualize the ROC curve. By default, plot() shows
the full ROC curve, including:
ROC curve with Youden point and annotation
The label and youden arguments toggle the
annotation and the Youden point (with its guide segments)
independently:
ROC curve with Youden point, no annotation
ROC curve with annotation, no Youden point
ROC curve only
The underlying ROC data and key summaries can be accessed directly:
head(roc$curve) # FPR, TPR, threshold, J (one row per unique FPR)
#> FPR TPR threshold J
#> 1 0.000000e+00 0.3107176 5.187889 0.3107176
#> 2 3.008076e-08 0.3134597 5.177544 0.3134597
#> 3 9.766935e-07 0.3162222 5.167200 0.3162212
#> 4 4.589869e-06 0.3190053 5.156855 0.3190007
#> 5 1.248845e-05 0.3218091 5.146511 0.3217966
#> 6 2.611099e-05 0.3246338 5.136166 0.3246077
roc$AUC # AUC estimate
#> [1] 0.8455996
roc$AUC_ci # AUC confidence interval
#> [1] 0.7914592 0.8997401
roc$J # Youden index estimate
#> [1] 0.5723763
roc$J_ci # Youden index CI
#> [1] 0.4542678 0.6904848
roc$t0 # Youden cutoff
#> [1] 4.07068
roc$sensitivity # Sensitivity at Youden point
#> [1] 0.7028534
roc$specificity # Specificity at Youden point
#> [1] 0.8695229
roc$hX # Bandwidth for non-diseased CDF
#> [1] 0.3146762
roc$hY # Bandwidth for diseased CDF
#> [1] 0.9802306This section summarizes the more technical aspects of
smoothROC() that may be useful for advanced users.
Setting logtrans = TRUE applies a natural log
transformation to the biomarker prior to ROC estimation. This is often
appropriate for biomarkers with strong right skew or multiplicative
variability (e.g. enzyme concentrations, cytokines). When
logtrans = TRUE, biomarker values must be strictly
positive.
The argument grid_n controls the resolution of the ROC
curve. Larger values produce a smoother ROC curve and more precise
localization of the Youden cutoff, at the cost of increased computation.
Reasonable values include:
grid_n = 500 – fast and adequate for exploratory
work;grid_n = 1000 – default, smoother curve and better
stability;grid_n = 3000 – more refined.For most applications, a good starting point is:
kernel = "biweight",bw_method = "PB".The "PB" method tends to perform well across a range of
scenarios and sample sizes. When computation time is a concern, the
"AR" method offers a simple, robust alternative that is
easy to compute.
In large samples, the differences between bandwidth methods may be minor. However, in small or moderate samples, bandwidth selection can substantially affect ROC shape, AUC estimates, and the stability of the Youden index.
Altman, N., & Leger, C. (1995). Bandwidth selection for kernel distribution function estimation. Journal of Statistical Planning and Inference, 46(2), 195–214.
Andrews, D. F., & Herzberg, A. M. (2012). Data: A Collection of Problems from Many Fields for the Student and Research Worker. Springer.
Bowman, A., Hall, P., & Prvan, T. (1998). Bandwidth selection for the smoothing of distribution functions. Biometrika, 85(4), 799–808.
DeLong, E. R., DeLong, D. M., & Clarke-Pearson, D. L. (1988). Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach. Biometrics, 44(3), 837–845.
Khan, R. A., & Ghebremichael, M. (2025). Smooth ROC Curve Estimation. Journal Name (preprint).
Lloyd, C. J. (1998). Using smoothed receiver operating characteristic curves to summarize and compare diagnostic systems. Journal of the American Statistical Association, 93(444), 1356–1364.
Polansky, A. M., & Baker, E. R. (2000). Multistage plug-in bandwidth selection for kernel distribution function estimates. Journal of Statistical Computation and Simulation, 65(1–4), 63–80.
Silverman, B. W. (1986). Density Estimation for Statistics and Data Analysis. Chapman & Hall, London.
Youden, W. J. (1950). Index for rating diagnostic tests. Cancer, 3(1), 32–35.
Zhou, X.-H., & Harezlak, J. (2002). Comparison of bandwidth selection methods for kernel smoothing of ROC curves. Statistics in Medicine, 21(14), 2045–2055.
Zou, K. H., Hall, W. J., & Shapiro, D. E. (1997). Smooth nonparametric receiver operating characteristic (ROC) curves for continuous diagnostic tests. Statistics in Medicine, 16(19), 2143–2156.