--- title: "funbootband: Simultaneous Prediction and Confidence Bands for Functional Data" author: "Daniel Koska" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{funbootband: Simultaneous Prediction and Confidence Bands for Functional Data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4, message = FALSE, warning = FALSE ) # Keep the vignette fast on CRAN by using tiny B is_cran <- !identical(tolower(Sys.getenv("NOT_CRAN")), "true") B_demo <- if (is_cran) 25L else 500L # bootstrap reps k.coef_demo <- 4L set.seed(1) ``` ## Overview `funbootband` computes **simultaneous** prediction and confidence bands for dense functional data observed on a common grid. It supports both i.i.d. and clustered (hierarchical) designs and uses a fast 'Rcpp' backend. Curves are **preprocessed via finite Fourier series** (`k.coef` harmonics) and **bootstrapped** to generate empirical distributions, from which band limits are obtained as quantiles. - **Prediction bands**: target simultaneous coverage of one *future individual curve*. In the clustered analysis this is specifically one curve from a new subject drawn from the subject population. - **Confidence bands**: cover the *mean function* with probability \(1-\alpha\). The main user function is: ```r band(data, type = c("prediction","confidence"), alpha = 0.05, iid = TRUE, id = NULL, B = 1000L, k.coef = 50L) ``` This document gives a quick tour of `funbootband`, including simulated i.i.d. and hierarchical examples. The i.i.d. calibration builds on Lenhoff et al. (1999). Koska et al. (2023) motivated the hierarchical setting, but the clustered implementation in version 0.3.0 is deliberately revised: it resamples whole subjects and defines a specific new-subject/new-curve prediction target. ### Statistical targets For independent curves $Y_i(t)$, the prediction band targets one independent future curve. A confidence band instead targets the population mean function. For clustered data, let $Y_{ij}(t)$ denote repeat $j=1,\ldots,m_i$ from subject $i=1,\ldots,K$. The default clustered prediction target is \[ Y_{\mathrm{new}}(t) = \mu(t) + A_{\mathrm{new}}(t) + E_{\mathrm{new}}(t), \] that is, one curve from an independent new subject. Subjects are sampled with equal probability and curves are sampled equally within subject. The empirical weight of curve $(i,j)$ is therefore \[ q_{ij}=\frac{1}{K m_i}, \] and the fitted centre is the equally subject-weighted mean \[ \widehat\mu(t)=\frac{1}{K}\sum_{i=1}^K \frac{1}{m_i}\sum_{j=1}^{m_i}Y_{ij}(t). \] In bootstrap replicate $b$, $K$ subjects are sampled with replacement. If subject $i$ is selected $c_{bi}$ times, all of its curves are retained with weight $c_{bi}/(K m_i)$. Prediction calibration keeps one supremum statistic per possible future curve; it does not maximize jointly over the complete observed sample. ### Quick start (i.i.d.) When curves are independent and identically distributed, the function argument `iid` in `band()` should be set to `TRUE`. For this example, consider smooth periodic curves on a common grid. ```{r iid-sim} library(funbootband) set.seed(1) T <- 101L n <- 30L x <- seq(0, 1, length.out = T) mu_true <- 0.7 * sin(2 * pi * x) - 0.2 * cos(4 * pi * x) generate_iid_curve <- function() { mu_true + rnorm(1, sd = 0.35) + rnorm(1, sd = 0.30) * sin(2 * pi * x) + rnorm(1, sd = 0.20) * cos(2 * pi * x) + rnorm(1, sd = 0.15) * sin(4 * pi * x) } Y <- replicate(n, generate_iid_curve()) ``` Simultaneous prediction and confidence bands are then computed by setting the `type` argument to either `prediction` or `confidence`: ```{r iid-sim-plot} # Fit prediction and confidence bands fit_pred <- band(Y, type = "prediction", alpha = 0.10, iid = TRUE, B = B_demo, k.coef = k.coef_demo) fit_conf <- band(Y, type = "confidence", alpha = 0.10, iid = TRUE, B = B_demo, k.coef = k.coef_demo) ``` When plotting the bands alongside the original curves, we see that the shaded region is calibrated to contain *entire curves* with probability \(1-\alpha\) (simultaneous coverage). ```{r iid-plot, fig.cap="Calculated prediction (blue) and confidence (gray) bands."} ylim <- range(c(Y, fit_pred$lower, fit_pred$upper), finite = TRUE) plot(x, fit_pred$mean, type = "n", ylim = ylim, xlab = "Normalized time", ylab = "Value", main = "Simultaneous bands (i.i.d.)") matlines(x, Y, col = grDevices::adjustcolor("gray40", 0.25), lty = 1) polygon(c(x, rev(x)), c(fit_pred$lower, rev(fit_pred$upper)), col = grDevices::adjustcolor("steelblue", alpha.f = 0.25), border = NA) polygon(c(x, rev(x)), c(fit_conf$lower, rev(fit_conf$upper)), col = grDevices::adjustcolor("darkorange", alpha.f = 0.30), border = NA) lines(x, fit_pred$mean, lwd = 2) lines(x, mu_true, col = "red", lwd = 2, lty = 2) ``` ## Clustered (hierarchical) curves When the i.i.d. assumption is violated, `iid` needs to be set to `FALSE`. `band()` will then automatically detect the cluster structure from the column names. Optionally, an integer/factor vector of length `ncol(data)` giving a cluster id for each curve can be used. The clustered case is illustrated using a design where each subject contributes repeated curves. The estimand first samples a subject uniformly from the subject population and then samples one curve from that subject. Consequently, subjects receive equal weight even when they contribute different numbers of curves. ```{r clustered, eval=TRUE} library(funbootband) set.seed(2) K_subject <- 12L m <- rep(c(2L, 3L, 4L), length.out = K_subject) id <- rep(seq_len(K_subject), m) subject_effect <- sapply(seq_len(K_subject), function(i) { rnorm(1, sd = 0.35) + rnorm(1, sd = 0.30) * sin(2 * pi * x) + rnorm(1, sd = 0.20) * cos(2 * pi * x) }) within_subject_effect <- function() { rnorm(1, sd = 0.18) * sin(4 * pi * x) + rnorm(1, sd = 0.12) * cos(4 * pi * x) } Y <- sapply(seq_along(id), function(j) { mu_true + subject_effect[, id[j]] + within_subject_effect() }) trial <- ave(id, id, FUN = seq_along) colnames(Y) <- paste0("subject", id, "_trial", trial) # Fit prediction and confidence bands fit_pred <- band(Y, type = "prediction", alpha = 0.10, iid = FALSE, id = id, B = B_demo, k.coef = k.coef_demo) fit_conf <- band(Y, type = "confidence", alpha = 0.10, iid = FALSE, id = id, B = B_demo, k.coef = k.coef_demo) ``` **Important:** When `iid = FALSE`, the bootstrap samples subjects with replacement and carries **all observed curves of each selected subject intact**. There is no second-stage resampling of individual curves in this revision. If a subject is selected more than once, its entire set of curves is copied more than once. The fitted mean is the equally weighted average of subject-specific mean curves, and prediction is calibrated for one curve from a new subject using equal-subject, equal-within-subject empirical weights. The prediction target and resampling unit are recorded explicitly: ```{r clustered-meta} fit_pred$meta[c("target", "weighting", "bootstrap_unit", "n_clusters")] ``` This interpretation assumes that subjects are independent draws from the population and that the observed repeats are exchangeable representatives of a curve from their subject. It is a **marginal population prediction band**: it does not condition on data already observed for a particular subject and does not promise joint coverage of several future curves. At least two subjects and some within-subject replication are required; reliable tail calibration will usually require substantially more than the formal minimum. Finally, the target inherits the preprocessing step: it is a new curve in the same finite-Fourier representation used to reconstruct the observed curves. Coverage of raw high-frequency measurement noise is not automatically implied when that variation is removed by the chosen `k.coef`. ```{r clustered-plot} ylim <- range(c(Y, fit_pred$lower, fit_pred$upper), finite = TRUE) plot(x, fit_pred$mean, type = "n", ylim = ylim, xlab = "Normalized time", ylab = "Value", main = "Simultaneous bands (clustered)") matlines(x, Y, col = grDevices::adjustcolor("gray40", 0.20), lty = 1) polygon(c(x, rev(x)), c(fit_pred$lower, rev(fit_pred$upper)), col = grDevices::adjustcolor("steelblue", alpha.f = 0.25), border = NA) polygon(c(x, rev(x)), c(fit_conf$lower, rev(fit_conf$upper)), col = grDevices::adjustcolor("darkorange", alpha.f = 0.30), border = NA) lines(x, fit_pred$mean, lwd = 2) lines(x, mu_true, col = "red", lwd = 2, lty = 2) ``` The following optional check generates one curve from each of many independent new subjects and estimates conditional simultaneous coverage of the fitted prediction band. It is skipped during CRAN checks and is illustrative rather than a replacement for an outer-loop simulation study. ```{r clustered-future-coverage, eval = !is_cran} generate_new_subject_curve <- function() { new_subject_effect <- rnorm(1, sd = 0.35) + rnorm(1, sd = 0.30) * sin(2 * pi * x) + rnorm(1, sd = 0.20) * cos(2 * pi * x) new_curve_effect <- rnorm(1, sd = 0.18) * sin(4 * pi * x) + rnorm(1, sd = 0.12) * cos(4 * pi * x) mu_true + new_subject_effect + new_curve_effect } future_curves <- replicate(500L, generate_new_subject_curve()) covered <- apply(future_curves, 2L, function(curve) { all(curve >= fit_pred$lower & curve <= fit_pred$upper) }) mean(covered) ``` ## Choosing `k.coef` (Fourier harmonics) `k.coef` controls the number of sine/cosine harmonics (plus intercept) used to represent each curve before bootstrapping. - The code automatically **clamps** `k.coef` to the maximum meaningful value on a length-`T` periodic grid. If a larger value is supplied, it is reduced accordingly. - In practice, choose the smallest `k.coef` that **adequately reconstructs** your curves. **Heuristic for selecting `k.coef`** The appropriate value of `k.coef` depends on both the **grid length** `T` and the **smoothness** of the curves: - Longer curves (larger `T`) can support higher `k.coef` values. - Smooth, slowly varying curves typically require only a small number of harmonics. - Highly oscillatory or noisy curves may need more harmonics to capture relevant detail. As a rule of thumb, increase `k.coef` only as far as necessary. Very high values mainly fit high-frequency noise, increase runtime, and may lead to numerical instability near the Nyquist limit. One simple way to inspect the effect of `k.coef` is to evaluate the **reconstruction error** (e.g., mean squared error, MSE) for different choices of `k.coef`: ```{r kcoef-mse-insample, eval = !is_cran} # MSE vs k.coef for the i.i.d. example (uses Y from above) fourier_basis <- function(T, K) { t <- 0:(T - 1L) denom <- T - 1L if (K == 0L) return(cbind(1)) cbind( 1, sapply(1:K, function(k) cos(2*pi*k*t/denom)), sapply(1:K, function(k) sin(2*pi*k*t/denom)) ) } reconstruct <- function(B, Y) { coef <- qr.coef(qr(B), Y) # solve for all curves at once B %*% coef } mse <- function(A, B) mean((A - B)^2) Ks <- c(10L, 20L, 30L, 40L, 50L, 60L, 70L, 80L, 90L, 99L) T <- nrow(Y) tab <- do.call(rbind, lapply(Ks, function(K){ B <- fourier_basis(T, K) Yhat <- reconstruct(B, Y) data.frame(k.coef = K, mse = mse(Y, Yhat), pve = 100 * (1 - sum((Y - Yhat)^2) / sum((Y - mean(Y))^2))) })) row.names(tab) <- NULL print(tab) # Plot op <- par(mar = c(4,4,2,1)) plot(tab$k.coef, tab$mse, type = "b", xlab = "k.coef", ylab = "MSE", main = "Fourier reconstruction error vs. k.coef") ``` **Recommendation** If your curves are smooth and runtime is not a concern, moderately large `k.coef` values (e.g., 20–50 for dense grids) are often a good choice. For noisy data or when computation time matters, examine how the reconstruction error (or resulting band limits) changes with k.coef and select the smallest value beyond which improvements become negligible. Recomputing the bands for a few candidate values can also help confirm that results have converged. ## Tuning `B` (bootstrap replicates) The argument B controls the number of bootstrap replications. - Increasing `B` generally yields tighter and more stable bands, as quantile estimates become more precise. - Beyond a certain threshold, however, band limits typically stabilize and no longer change meaningfully. - For exploratory analyses, smaller values of `B` can be used to reduce runtime. For final inference, rerun with a sufficiently large B to ensure convergence of the band limits. ## API reference ```r band(data, type = c("prediction","confidence"), alpha = 0.05, iid = TRUE, id = NULL, B = 1000L, k.coef = 50L) ``` - **data**: numeric matrix `T × n` (rows: time points; cols: curves). A numeric `data.frame` is accepted. - **type**: `"prediction"` (future curve coverage) or `"confidence"` (mean). - **alpha**: e.g., `0.05` gives 95% bands. - **iid** / **id**: set `iid = FALSE` and supply `id` (length `n`) for clusters, or allow inference from column-name prefixes. - **B**: bootstrap reps. - **k.coef**: Fourier harmonics (default `50`; clamped to a grid-specific max). **Return value.** A list with numeric vectors `lower`, `mean`, `upper` (length `T`) and `meta`. The metadata records the estimand, weighting convention, bootstrap unit, curve representation, and cluster sizes where applicable. ## References Lenhoff, M. W., Santner, T. J., Otis, J. C., Peterson, M. G., Williams, B. J., & Backus, S. I. (1999). Bootstrap prediction and confidence bands: a superior statistical method for analysis of gait data. *Gait & Posture*, 9(1), 10–17. Koska, D., Oriwol, D., & Maiwald, C. (2023). Comparison of statistical models for characterizing continuous differences between two biomechanical measurement systems. *Journal of Biomechanics*, 149, 111506. ## Session info ```{r session-info} sessionInfo() ```