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:
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.
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:
# 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).
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)
Calculated prediction (blue) and confidence
(gray) bands.
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.
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:
fit_pred$meta[c("target", "weighting", "bootstrap_unit", "n_clusters")]
#> $target
#> [1] "new_subject_new_curve"
#>
#> $weighting
#> [1] "equal_subject_then_equal_curve_within_subject"
#>
#> $bootstrap_unit
#> [1] "intact_subject"
#>
#> $n_clusters
#> [1] 12
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.
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.
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:
# 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.
API reference
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. doi:10.1016/S0966-6362(98)00043-5
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. doi:10.1016/j.jbiomech.2023.111506