---
title: "Kernel-Smoothed ROC Curves with smoothROC"
author: "Ruhul Ali Khan & Musie Ghebremichael"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: |
  %\VignetteIndexEntry{Kernel-Smoothed ROC Curves with smoothROC}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.align = "center"
)
library(smoothROC)
```

# Introduction

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.

# The `smoothROC()` function

The main user-facing function is:

```r
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:

- a smooth ROC curve on a fine grid of thresholds,
- a kernel-based estimator of the AUC with a confidence interval, and
- a kernel-smoothed Youden index with its optimal cutoff and confidence
  interval.

## Arguments

- `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`).

## Returned value

`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:

```r
print.smoothROC()
summary.smoothROC()
plot.smoothROC()
```

and are invoked automatically via `print()`, `summary()`, and `plot()`.

# Method overview

## Kernel CDF estimators and the smooth ROC curve

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 strategies

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.

## AUC estimation and kernel DeLong-type variance

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:

```r
roc$AUC
roc$AUC_ci
```

## Youden index and optimal cutoff

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:

```r
roc$J
roc$J_ci
roc$t0
roc$sensitivity
roc$specificity
```

## Bootstrap percentile confidence intervals

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.

```r
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:

```r
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.

# Example: Duchenne muscular dystrophy dataset

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.

```{r}
data(dystrophy)
str(dystrophy)
```

A basic smooth ROC analysis is:

```{r}
roc <- smoothROC(
  data      = dystrophy,
  biomarker = "CK",
  status    = "Class",
  diseased  = "carrier",
  kernel    = "biweight",
  bw_method = "PB",
  alpha     = 0.05,
  logtrans  = TRUE,
  grid_n    = 1000
)
```

## Displaying the result

The print and summary methods provide a concise summary:

```{r}
roc
summary(roc)
```

They 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:

- the 45-degree reference line (no-discrimination),
- the smooth ROC curve,
- the Youden point marked in red, and
- a label showing AUC, its confidence interval, the Youden index, and the
  Youden cutoff.

```{r roc-plot-full, fig.cap="ROC curve with Youden point and annotation"}
plot(roc)
```

The `label` and `youden` arguments toggle the annotation and the Youden point (with its guide segments) independently:

```{r roc-plot-nolabel, fig.cap="ROC curve with Youden point, no annotation"}
plot(roc, label = FALSE)
```

```{r roc-plot-noyouden, fig.cap="ROC curve with annotation, no Youden point"}
plot(roc, youden = FALSE)
```

```{r roc-plot-clean, fig.cap="ROC curve only"}
plot(roc, label = FALSE, youden = FALSE)
```

The underlying ROC data and key summaries can be accessed directly:

```{r}
head(roc$curve)     # FPR, TPR, threshold, J (one row per unique FPR)
roc$AUC             # AUC estimate
roc$AUC_ci          # AUC confidence interval
roc$J               # Youden index estimate
roc$J_ci            # Youden index CI
roc$t0              # Youden cutoff
roc$sensitivity     # Sensitivity at Youden point
roc$specificity     # Specificity at Youden point
roc$hX              # Bandwidth for non-diseased CDF
roc$hY              # Bandwidth for diseased CDF
```

# Advanced options

This section summarizes the more technical aspects of `smoothROC()` that may
be useful for advanced users.

## Log transformation

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.

## Controlling the ROC grid

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.

## Choosing a bandwidth method

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.

# References

- 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.
