From c6385a2dd5f9f82f7bcb90b99250e74238013762 Mon Sep 17 00:00:00 2001 From: fabian-s Date: Mon, 28 Sep 2026 13:03:29 +0200 Subject: [PATCH 1/3] pffr: bias-aware pointwise intervals (coef bias_ref, pffr_predict_ci) coef.pffr() gains `bias_ref`: given a second fit of the same model and data (typically REML for an NCV fit), each coefficient table gets a `delta` column (estimate difference through the same linear map as the SE) and pointwise z intervals use sqrt(se^2 + delta^2). New exported pffr_predict_ci() gives pointwise intervals for the linear predictor or conditional mean at fitted points or newdata, with sandwich/freq/cluster choices as in coef.pffr() and optional bias_ref; response-scale intervals transform the link-scale endpoints. The reference fit is checked for identical coefficient layout, model frame, weights, offsets, smooth bases, ffpc/pcre metadata, family and link. Tests reproduce the pffr-ci restart study's interval (NCV estimate, exact Bayesian CL2 SE, NCV-REML difference) for Gaussian and Poisson fits, delta = 0 for a fit as its own reference, seWithMean, missing dense and sparse responses, covariance passthrough, and rejections. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01JtmaCapnpxpQhDG2q2CQny --- NAMESPACE | 1 + NEWS.md | 9 + R/pffr-bias-aware.R | 322 +++++++++++++++++++++++ R/pffr-methods.R | 115 ++++++++- man/coef.pffr.Rd | 42 ++- man/coef_get_predictions.Rd | 4 +- man/compute_coef_delta.Rd | 20 ++ man/pffr_bias_aware_se.Rd | 20 ++ man/pffr_bias_ref_difference.Rd | 28 ++ man/pffr_predict_ci.Rd | 111 ++++++++ tests/testthat/test-pffr-bias-aware.R | 357 ++++++++++++++++++++++++++ 11 files changed, 1019 insertions(+), 10 deletions(-) create mode 100644 R/pffr-bias-aware.R create mode 100644 man/compute_coef_delta.Rd create mode 100644 man/pffr_bias_aware_se.Rd create mode 100644 man/pffr_bias_ref_difference.Rd create mode 100644 man/pffr_predict_ci.Rd create mode 100644 tests/testthat/test-pffr-bias-aware.R diff --git a/NAMESPACE b/NAMESPACE index 672dfdd3..06955338 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -91,6 +91,7 @@ export(pffr_coefboot) export(pffr_dependence_check) export(pffr_gls) export(pffr_jackknife_se) +export(pffr_predict_ci) export(pffr_qq) export(pffr_simulate) export(pffr_upgrade_fit) diff --git a/NEWS.md b/NEWS.md index 56e68018..da2d32b9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,14 @@ # refund 0.1-38 +* Bias-aware pointwise intervals for NCV fits: `coef.pffr()` gains + `bias_ref`, and the new `pffr_predict_ci()` gives pointwise intervals for the + linear predictor or conditional mean (optionally bias-aware). Given a second + fit of the same model (typically REML for an NCV fit), the interval is + `est -/+ z * sqrt(se^2 + delta^2)` with `delta` the difference between the + two estimates, built on the link scale with transformed endpoints for the + response mean. Recommended: NCV with curve blocks, exact CL2 in Bayesian form + (`sandwich = "cl2", cl2_adjustment = "exact"`), REML reference. `delta` + does not capture smoothing bias that both fits share. * `pffr(method = "NCV")` now leaves out whole curves (or `cluster` groups) by default. `ncv_blocks = "point"` provides pointwise comparisons; explicit `nei` takes precedence. Dual mgcv neighbourhood names, a cached behavioural diff --git a/R/pffr-bias-aware.R b/R/pffr-bias-aware.R new file mode 100644 index 00000000..be4e257e --- /dev/null +++ b/R/pffr-bias-aware.R @@ -0,0 +1,322 @@ +# Bias-aware pointwise intervals for pffr fits. +# +# The interval for a linear estimand L theta of a fit (typically selected by +# curve-blocked NCV) is +# L theta -/+ z * sqrt(se^2 + delta^2), delta = L (theta - theta_ref), +# where theta_ref comes from a reference fit of the same model (typically REML) +# and se from the chosen covariance (recommended: exact CL2, Bayesian form). + +# Helpers ------------------------------------------------------------------------ + +#' Check that a bias reference fit matches a pffr fit +#' +#' The bias-aware interval compares two fits of the *same* model and data that +#' differ only in how the smoothing parameters were chosen, so that their +#' coefficient vectors live in the same basis. This checks the coefficient +#' layout, the model frame (response, covariates and functional-covariate +#' matrices, in the same row order), prior weights and offsets, the smooth bases +#' (labels, classes, basis dimensions, knots and penalty matrices, which carry the +#' identifiability constraints), `ffpc`/`pcre` metadata, and family and link. +#' Family parameters estimated during fitting (e.g. the negative binomial +#' \eqn{\theta}{theta}) may differ between the fits. +#' +#' @param object The fit whose intervals are computed. +#' @param bias_ref The reference fit. +#' @returns `object$coefficients - bias_ref$coefficients`. +#' @keywords internal +pffr_bias_ref_difference <- function(object, bias_ref) { + if (!inherits(bias_ref, "pffr")) { + stop("`bias_ref` must be a fitted pffr model.", call. = FALSE) + } + mismatch <- function(what) { + stop( + "`bias_ref` is not a fit of the same model and data as `object` (", + what, + " differ). Refit the reference with the same formula, bases and data, ", + "changing only `method`.", + call. = FALSE + ) + } + if (!identical(names(object$coefficients), names(bias_ref$coefficients))) { + mismatch("coefficient names") + } + if ( + !identical( + pffr_family_base(object$family), + pffr_family_base(bias_ref$family) + ) + ) { + mismatch("families or links") + } + same <- function(a, b) isTRUE(all.equal(a, b, check.attributes = FALSE)) + meta <- c( + "nobs", + "yind", + "is_sparse", + "missing_indices", + "ffpc", + "pcre_terms" + ) + if (!same(object$pffr[meta], bias_ref$pffr[meta])) { + mismatch("response grids, missing-value patterns or ffpc/pcre bases") + } + if (!same(object$model, bias_ref$model)) { + mismatch("model frames (responses or covariates)") + } + if ( + !same(object$prior.weights, bias_ref$prior.weights) || + !same(object$offset, bias_ref$offset) + ) { + mismatch("weights or offsets") + } + if (!pffr_same_smooths(object$smooth, bias_ref$smooth)) { + mismatch("smooth terms") + } + object$coefficients - bias_ref$coefficients +} + +# Family name without fitted parameters (mgcv writes e.g. "Negative +# Binomial(5.2)" into a fitted nb family), plus the link. +pffr_family_base <- function(family) { + c(tolower(sub("\\(.*$", "", family$family)), family$link) +} + +# Smooth terms match if labels, classes, coefficient ranges, basis dimensions, +# knots and penalty matrices agree. +pffr_same_smooths <- function(smooths, ref_smooths) { + if (length(smooths) != length(ref_smooths)) return(FALSE) + same <- function(a, b) { + identical(a$label, b$label) && + identical(class(a), class(b)) && + identical(c(a$first.para, a$last.para), c(b$first.para, b$last.para)) && + identical(a$bs.dim, b$bs.dim) && + isTRUE(all.equal(a$S, b$S)) && + isTRUE(all.equal(pffr_smooth_knots(a), pffr_smooth_knots(b))) + } + all(mapply(same, smooths, ref_smooths)) +} + +pffr_smooth_knots <- function(sm) { + if (!is.null(sm$margin)) return(lapply(sm$margin, pffr_smooth_knots)) + sm$knots +} + +#' Combine a variance-only standard error with a bias estimate +#' +#' @param se Standard errors (variance part). +#' @param delta Bias estimates of the same length. +#' @returns `sqrt(se^2 + delta^2)`. +#' @keywords internal +pffr_bias_aware_se <- function(se, delta) { + sqrt(se^2 + delta^2) +} + +# Pointwise intervals for the linear predictor / conditional mean ----------------- + +#' Pointwise confidence intervals for pffr predictions, optionally bias-aware +#' +#' Pointwise Wald intervals for the linear predictor, or for the conditional +#' mean \eqn{E(Y(t) \mid X)}{E(Y(t)|X)}, of a [pffr()] fit, with a choice of +#' covariance (model-based or sandwich, as in [coef.pffr()]) and an optional +#' bias-aware widening computed from a second fit of the same model. +#' +#' Intervals are always built on the link scale, +#' \deqn{\hat\eta \pm z_{(1+\mathrm{level})/2}\sqrt{\mathrm{se}^2 + +#' \delta^2},}{eta_hat -/+ z * sqrt(se^2 + delta^2),} +#' with \eqn{\delta = 0} unless `bias_ref` is supplied. For +#' `type = "response"` the estimate and both endpoints are mapped through the +#' inverse link (so the interval is not symmetric around the estimate); `se` +#' and `delta` stay on the link scale. +#' +#' @section Bias-aware intervals: +#' Smoothing-parameter selection by curve-blocked neighbourhood +#' cross-validation (`pffr(method = "NCV")`) is robust to within-curve +#' dependence but smooths more than REML, so its smoothing bias is no longer +#' negligible and a variance-only interval undercovers. The bias-aware interval +#' adds the difference between the two estimates in quadrature: +#' \deqn{\delta = L(\hat\theta_{\mathrm{NCV}} - \hat\theta_{\mathrm{REML}}),}{ +#' delta = L (theta_NCV - theta_REML),} +#' for the rows \eqn{L} of the prediction matrix. The recommended recipe is +#' an NCV fit with curve blocks (the default `ncv_blocks = "cluster"`) as +#' `object`, the exact CL2 sandwich in its Bayesian form +#' (`sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE`) and a REML fit +#' of the same model as `bias_ref`. Fit both with `sandwich = "none"` to avoid +#' computing a fit-time sandwich that is not used. +#' +#' \eqn{\delta} estimates only the part of the smoothing bias in which the two +#' fits differ. Bias that both fits share -- a basis too small for the truth, +#' or both fits oversmoothing a rough truth -- is not covered. +#' +#' @param object A fitted [pffr()] model with a single linear predictor. +#' @param newdata Optional prediction data in the format supplied to [pffr()] +#' (as in [predict.pffr()]). `NULL` (default) evaluates at the fitted +#' observation points. Fits with a model offset are supported at the fitted +#' points only. +#' @param type `"link"` (default) for the linear predictor, `"response"` for the +#' conditional mean. +#' @param level Confidence level, defaults to `0.95`. +#' @param sandwich,freq,cluster,dof_correction,edf_type,cl2_adjustment +#' Covariance choice, as in [coef.pffr()]. `sandwich = NULL` inherits the +#' fit-time choice. +#' @param bias_ref Optional reference fit of the same model and data (typically +#' the REML fit when `object` is the NCV fit), differing only in how the +#' smoothing parameters were chosen. If supplied, intervals are bias-aware +#' (see the section above). Any second fit of the same model is accepted; +#' \eqn{\delta}{delta} is then the contrast between the two estimators, which +#' is a smoothing-bias proxy only for the NCV-versus-REML pairing. +#' @returns A data frame with one row per evaluation point, in the row order of +#' `predict(object, type = "lpmatrix")` (curve-major, index fastest; rows with +#' missing responses are omitted and sparse responses keep the order of +#' `ydata` when `newdata = NULL`), with columns `.obs` (curve), `.index` +#' (value of the response index), `fit` (on the scale of `type`), `se_link` +#' (variance part, link scale), `delta_link` (link scale; only with +#' `bias_ref`), `lower` and `upper` (on the scale of `type`). Attributes +#' `type`, `level`, `crit_value` and `bias_ref_method` (the reference fit's +#' smoothing-parameter method, or `NA`). +#' @seealso [coef.pffr()] (argument `bias_ref`) for bias-aware intervals of +#' coefficient functions, [pffr_jackknife_se()]. +#' @export +#' @author Fabian Scheipl +#' @examples +#' \donttest{ +#' set.seed(1) +#' d <- pffr_simulate(Y ~ ff(X1), n = 30, nxgrid = 15, nygrid = 25) +#' yind <- attr(d, "yindex") +#' fit_ncv <- pffr(Y ~ ff(X1), yind = yind, data = d, method = "NCV", +#' sandwich = "none") +#' fit_reml <- pffr(Y ~ ff(X1), yind = yind, data = d, method = "REML", +#' sandwich = "none") +#' ci <- pffr_predict_ci(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", +#' bias_ref = fit_reml) +#' head(ci) +#' } +pffr_predict_ci <- function( + object, + newdata = NULL, + type = c("link", "response"), + level = 0.95, + sandwich = NULL, + freq = FALSE, + cluster = NULL, + dof_correction = NULL, + edf_type = NULL, + cl2_adjustment = NULL, + bias_ref = NULL +) { + if (!inherits(object, "pffr")) { + stop("`object` must be a fitted pffr model.", call. = FALSE) + } + type <- match.arg(type) + if ( + !is.numeric(level) || + length(level) != 1 || + !is.finite(level) || + level <= 0 || + level >= 1 + ) { + stop("`level` must be a single number in (0, 1).", call. = FALSE) + } + if (isTRUE((object$family$nlp %||% 1L) > 1L)) { + stop( + "pffr_predict_ci() supports single linear-predictor families only.", + call. = FALSE + ) + } + theta_diff <- if (!is.null(bias_ref)) { + pffr_bias_ref_difference(object, bias_ref) + } + + lp <- pffr_prediction_design(object, newdata) + V <- pffr_vcov( + object, + sandwich = sandwich, + freq = freq, + cluster = cluster, + dof_correction = dof_correction, + edf_type = edf_type, + cl2_adjustment = cl2_adjustment + ) + se <- sqrt(as.vector(rowSums((lp$X %*% V) * lp$X))) + delta <- if (!is.null(theta_diff)) as.vector(lp$X %*% theta_diff) + crit_value <- stats::qnorm((1 + level) / 2) + half <- crit_value * if (is.null(delta)) se else pffr_bias_aware_se(se, delta) + lower <- lp$eta - half + upper <- lp$eta + half + fit <- lp$eta + if (type == "response") { + linkinv <- object$family$linkinv + ends <- cbind(linkinv(lower), linkinv(upper)) + # A decreasing inverse link (e.g. the inverse link) swaps the endpoints. + lower <- pmin(ends[, 1], ends[, 2]) + upper <- pmax(ends[, 1], ends[, 2]) + fit <- linkinv(fit) + } + out <- lp$points + out$fit <- fit + out$se_link <- se + if (!is.null(delta)) out$delta_link <- delta + out$lower <- lower + out$upper <- upper + attr(out, "type") <- type + attr(out, "level") <- level + attr(out, "crit_value") <- crit_value + attr(out, "bias_ref_method") <- if (is.null(bias_ref)) { + NA_character_ + } else { + bias_ref$method %||% NA_character_ + } + out +} + +# Prediction matrix, linear predictor and evaluation points (.obs, .index) at +# the fitted points or at newdata. At the fitted points the stored linear +# predictor includes any offset exactly. +pffr_prediction_design <- function(object, newdata) { + meta <- object$pffr + if (is.null(newdata)) { + X <- predict(object, type = "lpmatrix", reformat = FALSE) + eta <- as.vector(object$linear.predictors) + points <- if (isTRUE(meta$is_sparse)) { + # The model frame omits ydata rows with a missing .value. + meta$ydata[!is.na(meta$ydata$.value), c(".obs", ".index")] + } else { + pffr_grid_points(meta$nobs, meta$yind, meta$missing_indices) + } + } else { + if (!is.null(object$offset) && any(object$offset != 0)) { + stop( + "This fit uses a model offset: intervals for `newdata` are not ", + "supported. Use `newdata = NULL` for the fitted points.", + call. = FALSE + ) + } + X <- predict(object, newdata = newdata, type = "lpmatrix", reformat = FALSE) + eta <- as.vector(X %*% object$coefficients) + points <- pffr_grid_points(nrow(X) / length(meta$yind), meta$yind) + } + if (!is.null(attr(X, "lpi"))) { + stop( + "pffr_predict_ci() supports single linear-predictor families only.", + call. = FALSE + ) + } + if (nrow(points) != nrow(X) || length(eta) != nrow(X)) { + stop( + "Could not align evaluation points with the prediction matrix.", + call. = FALSE + ) + } + rownames(points) <- NULL + list(X = X, eta = eta, points = points) +} + +# Evaluation points of a dense response grid in lpmatrix order (curve-major), +# without the rows of missing responses. +pffr_grid_points <- function(nobs, yind, missing_indices = NULL) { + points <- data.frame( + .obs = rep(seq_len(nobs), each = length(yind)), + .index = rep(yind, times = nobs) + ) + if (length(missing_indices)) points <- points[-missing_indices, ] + points +} diff --git a/R/pffr-methods.R b/R/pffr-methods.R index 0c955e1d..aeb92d9a 100755 --- a/R/pffr-methods.R +++ b/R/pffr-methods.R @@ -879,7 +879,9 @@ ensure_grid_axis_attributes <- function(d, trm, is_pcre, pffr_info) { #' #' @param trm Smooth term object. #' @param data_grid Data frame from coef_make_data_grid. -#' @param object_info List with: coefficients, cmX, Vp. +#' @param object_info List with: coefficients, cmX, Vp, and optionally +#' theta_diff (coefficient difference to a bias reference fit; see +#' `bias_ref` in [coef.pffr()]). #' @param pffr_info List with: yind_name. #' @param covmat Covariance matrix for SE computation. #' @param se Logical, compute standard errors? @@ -926,8 +928,11 @@ coef_get_predictions <- function( P$value <- X %*% object_info$coefficients[trmind] P$coef <- cbind(data_grid, value = P$value) - # Compute standard errors and intervals if requested - if (se) { + # Bias estimate of a bias-aware interval: the reference-fit contrast through + # the same linear map as the standard error (so it includes the mean-level + # columns when seWithMean applies). + theta_diff <- object_info$theta_diff + if (se || !is.null(theta_diff)) { linear_map <- build_coef_linear_map( X = X, trmind = trmind, @@ -935,9 +940,16 @@ coef_get_predictions <- function( object_info = object_info, seWithMean = seWithMean ) + } + if (se) { P$se <- compute_coef_se(linear_map = linear_map, covmat = covmat) P$coef <- cbind(P$coef, se = P$se) - + } + if (!is.null(theta_diff)) { + P$delta <- compute_coef_delta(linear_map, theta_diff) + P$coef <- cbind(P$coef, delta = P$delta) + } + if (se) { if (ci == "simultaneous") { crit <- compute_ci_critical( ci = ci, @@ -965,7 +977,12 @@ coef_get_predictions <- function( crit_df_const = crit_df_const ) P$crit <- pw$crit - ci_half <- pw$crit * P$se + ci_se <- if (is.null(theta_diff)) { + P$se + } else { + pffr_bias_aware_se(P$se, P$delta) + } + ci_half <- pw$crit * ci_se P$coef <- cbind( P$coef, lower = P$value - ci_half, @@ -1057,6 +1074,17 @@ compute_coef_se <- function(linear_map, covmat) { sqrt(rowSums((linear_map$X %*% covmat[trmind, trmind]) * linear_map$X)) } +#' Evaluate a coefficient difference through a coefficient linear map +#' +#' @param linear_map List returned by build_coef_linear_map(). +#' @param theta_diff Full-length coefficient difference vector. +#' @returns Numeric vector, one value per evaluation point. +#' @keywords internal +compute_coef_delta <- function(linear_map, theta_diff) { + d <- if (linear_map$use_full) theta_diff else theta_diff[linear_map$trmind] + as.vector(linear_map$X %*% d) +} + #' Draw coefficient perturbations for simultaneous intervals #' #' @param covmat Covariance matrix. @@ -1420,6 +1448,18 @@ pffr_end_undefined_df <- function(opened) { #' @param n_sim Number of simulations for simultaneous intervals, defaults to #' \code{2000}. Ignored unless \code{ci = "simultaneous"}. #' @param sim_seed Optional integer seed for simultaneous interval simulation. +#' @param bias_ref Optional reference fit for bias-aware intervals: a +#' \code{pffr} fit of the same model and data that differs from \code{object} +#' only in how the smoothing parameters were chosen (typically the REML fit +#' when \code{object} was fitted with \code{method = "NCV"}). If supplied, +#' every returned coefficient table gets a column \code{delta}, the estimate +#' of \code{object} minus that of \code{bias_ref} (evaluated through the same +#' linear map as \code{se}, so with \code{seWithMean = TRUE} it includes the +#' difference in the mean level), and pointwise intervals use +#' \eqn{\sqrt{se^2 + delta^2}}{sqrt(se^2 + delta^2)} instead of \code{se}; +#' \code{se} itself stays the variance part. Only \code{ci = "none"} and +#' \code{ci = "pointwise"} with \code{crit = "z"} are available. See the +#' section \sQuote{Bias-aware intervals}. #' @param ... other arguments, not used. #' #' @return If \code{raw==FALSE}, a list containing \itemize{ @@ -1441,7 +1481,31 @@ pffr_end_undefined_df <- function(opened) { #' \code{crit = "tG1"}, and the per-point Satterthwaite df for #' \code{crit = "satterthwaite"}/\code{"auto"}). The returned list also includes #' \code{ci_meta} with CI settings (including \code{crit} and the resolved -#' \code{crit_used}). +#' \code{crit_used}). With \code{bias_ref}, the matrices also include a column +#' \code{delta} (placed after \code{se}) and \code{ci_meta$bias_ref_method} +#' records the smoothing-parameter method of the reference fit. +#' @section Bias-aware intervals: +#' Smoothing-parameter selection by curve-blocked neighbourhood +#' cross-validation (\code{pffr(method = "NCV")}) is robust to within-curve +#' dependence but smooths more than REML, so its smoothing bias is no longer +#' negligible and variance-only intervals undercover. The bias-aware pointwise +#' interval for an estimate \eqn{L\hat\theta_{NCV}}{L theta_NCV} is +#' \deqn{L\hat\theta_{NCV} \pm z_{(1+level)/2}\sqrt{se^2 + \delta^2}, +#' \qquad \delta = L(\hat\theta_{NCV} - \hat\theta_{REML}),}{ +#' L theta_NCV -/+ z * sqrt(se^2 + delta^2), delta = L (theta_NCV - theta_REML),} +#' on the link scale. The recommended recipe fits the model twice with +#' \code{sandwich = "none"}, by NCV with curve blocks (the default +#' \code{ncv_blocks = "cluster"}) and by REML, and computes \code{se} from the +#' exact CL2 sandwich in its Bayesian form of the NCV fit: +#' \preformatted{coef(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", +#' ci = "pointwise", bias_ref = fit_reml)} +#' \eqn{\delta}{delta} estimates only the part of the smoothing bias in which +#' the two fits differ: bias that both fits share (a basis too small for the +#' truth, or both fits oversmoothing a rough truth) is not covered. Any second +#' fit of the same model is accepted as \code{bias_ref}; \eqn{\delta}{delta} +#' is then the contrast between the two estimators. For fitted values and +#' predictions see \code{\link{pffr_predict_ci}}. +#' #' @method coef pffr #' @export #' @importFrom mgcv PredictMat get.var @@ -1469,6 +1533,7 @@ coef.pffr <- function( level = 0.95, n_sim = 2000, sim_seed = NULL, + bias_ref = NULL, ... ) { # One coef() call warns at most once about undefined moment df, however many @@ -1556,6 +1621,28 @@ coef.pffr <- function( } if (!is.null(sim_seed)) sim_seed <- as.integer(sim_seed) + theta_diff <- NULL + if (!is.null(bias_ref)) { + if (raw) { + stop("`bias_ref` cannot be combined with raw = TRUE.", call. = FALSE) + } + if (ci == "simultaneous") { + stop( + "Bias-aware intervals (`bias_ref`) are pointwise only; use ", + "ci = \"pointwise\".", + call. = FALSE + ) + } + if (crit != "z") { + stop( + "Bias-aware intervals (`bias_ref`) use crit = \"z\"; the reference ", + "df of the other choices describe only the variance part.", + call. = FALSE + ) + } + theta_diff <- pffr_bias_ref_difference(object, bias_ref) + } + dots <- list(...) eval_grid <- dots$eval_grid %||% NULL @@ -1586,7 +1673,8 @@ coef.pffr <- function( object_info <- list( coefficients = object$coefficients, cmX = object$cmX, - Vp = object$Vp + Vp = object$Vp, + theta_diff = theta_diff ) getCoefs <- function(i) { @@ -1777,9 +1865,15 @@ coef.pffr <- function( })) ret$pterms <- cbind(value = object$coefficients[-smind]) if (se) ret$pterms <- cbind(ret$pterms, se = sqrt(diag(covmat)[-smind])) + if (!is.null(theta_diff)) { + ret$pterms <- cbind(ret$pterms, delta = theta_diff[-smind]) + } if (se && ci != "none") { p_se <- ret$pterms[, "se"] + if (!is.null(theta_diff)) { + p_se <- pffr_bias_aware_se(p_se, ret$pterms[, "delta"]) + } p_df <- NULL if (ci == "pointwise") { prob <- (1 + level) / 2 @@ -1863,7 +1957,12 @@ coef.pffr <- function( ci_ref_n_clusters = ci_ref_n_clusters, ci_ref_df = ci_ref_df, crit = crit, - crit_used = if (ci == "pointwise") crit_mode else NA_character_ + crit_used = if (ci == "pointwise") crit_mode else NA_character_, + bias_ref_method = if (is.null(bias_ref)) { + NA_character_ + } else { + bias_ref$method %||% NA_character_ + } ) return(ret) } diff --git a/man/coef.pffr.Rd b/man/coef.pffr.Rd index 004bf494..c8e0abfa 100644 --- a/man/coef.pffr.Rd +++ b/man/coef.pffr.Rd @@ -25,6 +25,7 @@ level = 0.95, n_sim = 2000, sim_seed = NULL, + bias_ref = NULL, ... ) } @@ -142,6 +143,19 @@ used the shortcut leverage weight \eqn{A_g=(I-H_{gg})^{-1/2}}, so \item{sim_seed}{Optional integer seed for simultaneous interval simulation.} +\item{bias_ref}{Optional reference fit for bias-aware intervals: a +\code{pffr} fit of the same model and data that differs from \code{object} +only in how the smoothing parameters were chosen (typically the REML fit +when \code{object} was fitted with \code{method = "NCV"}). If supplied, +every returned coefficient table gets a column \code{delta}, the estimate +of \code{object} minus that of \code{bias_ref} (evaluated through the same +linear map as \code{se}, so with \code{seWithMean = TRUE} it includes the +difference in the mean level), and pointwise intervals use +\eqn{\sqrt{se^2 + delta^2}}{sqrt(se^2 + delta^2)} instead of \code{se}; +\code{se} itself stays the variance part. Only \code{ci = "none"} and +\code{ci = "pointwise"} with \code{crit = "z"} are available. See the +section \sQuote{Bias-aware intervals}.} + \item{...}{other arguments, not used.} } \value{ @@ -164,7 +178,9 @@ critical value (\code{Inf} for \code{crit = "z"}, \eqn{G-1} for \code{crit = "tG1"}, and the per-point Satterthwaite df for \code{crit = "satterthwaite"}/\code{"auto"}). The returned list also includes \code{ci_meta} with CI settings (including \code{crit} and the resolved -\code{crit_used}). +\code{crit_used}). With \code{bias_ref}, the matrices also include a column +\code{delta} (placed after \code{se}) and \code{ci_meta$bias_ref_method} +records the smoothing-parameter method of the reference fit. } \description{ Returns estimated coefficient functions/surfaces \eqn{\beta(t), \beta(s,t)} @@ -182,6 +198,30 @@ With \code{sandwich="hc"}, mgcv's observation-level HC sandwich is used. If the model was fitted with a matching sandwich option in \code{\link{pffr}}, the pre-computed covariance matrices are used directly. } +\section{Bias-aware intervals}{ + +Smoothing-parameter selection by curve-blocked neighbourhood +cross-validation (\code{pffr(method = "NCV")}) is robust to within-curve +dependence but smooths more than REML, so its smoothing bias is no longer +negligible and variance-only intervals undercover. The bias-aware pointwise +interval for an estimate \eqn{L\hat\theta_{NCV}}{L theta_NCV} is +\deqn{L\hat\theta_{NCV} \pm z_{(1+level)/2}\sqrt{se^2 + \delta^2}, + \qquad \delta = L(\hat\theta_{NCV} - \hat\theta_{REML}),}{ + L theta_NCV -/+ z * sqrt(se^2 + delta^2), delta = L (theta_NCV - theta_REML),} +on the link scale. The recommended recipe fits the model twice with +\code{sandwich = "none"}, by NCV with curve blocks (the default +\code{ncv_blocks = "cluster"}) and by REML, and computes \code{se} from the +exact CL2 sandwich in its Bayesian form of the NCV fit: +\preformatted{coef(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", + ci = "pointwise", bias_ref = fit_reml)} +\eqn{\delta}{delta} estimates only the part of the smoothing bias in which +the two fits differ: bias that both fits share (a basis too small for the +truth, or both fits oversmoothing a rough truth) is not covered. Any second +fit of the same model is accepted as \code{bias_ref}; \eqn{\delta}{delta} +is then the contrast between the two estimators. For fitted values and +predictions see \code{\link{pffr_predict_ci}}. +} + \seealso{ \code{\link[mgcv]{plot.gam}}, \code{\link[mgcv]{predict.gam}} which this routine is based on. diff --git a/man/coef_get_predictions.Rd b/man/coef_get_predictions.Rd index 10ae182b..123b94a3 100644 --- a/man/coef_get_predictions.Rd +++ b/man/coef_get_predictions.Rd @@ -27,7 +27,9 @@ coef_get_predictions( \item{data_grid}{Data frame from coef_make_data_grid.} -\item{object_info}{List with: coefficients, cmX, Vp.} +\item{object_info}{List with: coefficients, cmX, Vp, and optionally +theta_diff (coefficient difference to a bias reference fit; see +`bias_ref` in [coef.pffr()]).} \item{pffr_info}{List with: yind_name.} diff --git a/man/compute_coef_delta.Rd b/man/compute_coef_delta.Rd new file mode 100644 index 00000000..7b8e43c1 --- /dev/null +++ b/man/compute_coef_delta.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-methods.R +\name{compute_coef_delta} +\alias{compute_coef_delta} +\title{Evaluate a coefficient difference through a coefficient linear map} +\usage{ +compute_coef_delta(linear_map, theta_diff) +} +\arguments{ +\item{linear_map}{List returned by build_coef_linear_map().} + +\item{theta_diff}{Full-length coefficient difference vector.} +} +\value{ +Numeric vector, one value per evaluation point. +} +\description{ +Evaluate a coefficient difference through a coefficient linear map +} +\keyword{internal} diff --git a/man/pffr_bias_aware_se.Rd b/man/pffr_bias_aware_se.Rd new file mode 100644 index 00000000..9f25e085 --- /dev/null +++ b/man/pffr_bias_aware_se.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-bias-aware.R +\name{pffr_bias_aware_se} +\alias{pffr_bias_aware_se} +\title{Combine a variance-only standard error with a bias estimate} +\usage{ +pffr_bias_aware_se(se, delta) +} +\arguments{ +\item{se}{Standard errors (variance part).} + +\item{delta}{Bias estimates of the same length.} +} +\value{ +`sqrt(se^2 + delta^2)`. +} +\description{ +Combine a variance-only standard error with a bias estimate +} +\keyword{internal} diff --git a/man/pffr_bias_ref_difference.Rd b/man/pffr_bias_ref_difference.Rd new file mode 100644 index 00000000..0e6eb764 --- /dev/null +++ b/man/pffr_bias_ref_difference.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-bias-aware.R +\name{pffr_bias_ref_difference} +\alias{pffr_bias_ref_difference} +\title{Check that a bias reference fit matches a pffr fit} +\usage{ +pffr_bias_ref_difference(object, bias_ref) +} +\arguments{ +\item{object}{The fit whose intervals are computed.} + +\item{bias_ref}{The reference fit.} +} +\value{ +`object$coefficients - bias_ref$coefficients`. +} +\description{ +The bias-aware interval compares two fits of the *same* model and data that +differ only in how the smoothing parameters were chosen, so that their +coefficient vectors live in the same basis. This checks the coefficient +layout, the model frame (response, covariates and functional-covariate +matrices, in the same row order), prior weights and offsets, the smooth bases +(labels, classes, basis dimensions, knots and penalty matrices, which carry the +identifiability constraints), `ffpc`/`pcre` metadata, and family and link. +Family parameters estimated during fitting (e.g. the negative binomial +\eqn{\theta}{theta}) may differ between the fits. +} +\keyword{internal} diff --git a/man/pffr_predict_ci.Rd b/man/pffr_predict_ci.Rd new file mode 100644 index 00000000..a1dd3474 --- /dev/null +++ b/man/pffr_predict_ci.Rd @@ -0,0 +1,111 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-bias-aware.R +\name{pffr_predict_ci} +\alias{pffr_predict_ci} +\title{Pointwise confidence intervals for pffr predictions, optionally bias-aware} +\usage{ +pffr_predict_ci( + object, + newdata = NULL, + type = c("link", "response"), + level = 0.95, + sandwich = NULL, + freq = FALSE, + cluster = NULL, + dof_correction = NULL, + edf_type = NULL, + cl2_adjustment = NULL, + bias_ref = NULL +) +} +\arguments{ +\item{object}{A fitted [pffr()] model with a single linear predictor.} + +\item{newdata}{Optional prediction data in the format supplied to [pffr()] +(as in [predict.pffr()]). `NULL` (default) evaluates at the fitted +observation points. Fits with a model offset are supported at the fitted +points only.} + +\item{type}{`"link"` (default) for the linear predictor, `"response"` for the +conditional mean.} + +\item{level}{Confidence level, defaults to `0.95`.} + +\item{sandwich, freq, cluster, dof_correction, edf_type, cl2_adjustment}{Covariance choice, as in [coef.pffr()]. `sandwich = NULL` inherits the +fit-time choice.} + +\item{bias_ref}{Optional reference fit of the same model and data (typically +the REML fit when `object` is the NCV fit), differing only in how the +smoothing parameters were chosen. If supplied, intervals are bias-aware +(see the section above). Any second fit of the same model is accepted; +\eqn{\delta}{delta} is then the contrast between the two estimators, which +is a smoothing-bias proxy only for the NCV-versus-REML pairing.} +} +\value{ +A data frame with one row per evaluation point, in the row order of + `predict(object, type = "lpmatrix")` (curve-major, index fastest; rows with + missing responses are omitted and sparse responses keep the order of + `ydata` when `newdata = NULL`), with columns `.obs` (curve), `.index` + (value of the response index), `fit` (on the scale of `type`), `se_link` + (variance part, link scale), `delta_link` (link scale; only with + `bias_ref`), `lower` and `upper` (on the scale of `type`). Attributes + `type`, `level`, `crit_value` and `bias_ref_method` (the reference fit's + smoothing-parameter method, or `NA`). +} +\description{ +Pointwise Wald intervals for the linear predictor, or for the conditional +mean \eqn{E(Y(t) \mid X)}{E(Y(t)|X)}, of a [pffr()] fit, with a choice of +covariance (model-based or sandwich, as in [coef.pffr()]) and an optional +bias-aware widening computed from a second fit of the same model. +} +\details{ +Intervals are always built on the link scale, +\deqn{\hat\eta \pm z_{(1+\mathrm{level})/2}\sqrt{\mathrm{se}^2 + + \delta^2},}{eta_hat -/+ z * sqrt(se^2 + delta^2),} +with \eqn{\delta = 0} unless `bias_ref` is supplied. For +`type = "response"` the estimate and both endpoints are mapped through the +inverse link (so the interval is not symmetric around the estimate); `se` +and `delta` stay on the link scale. +} +\section{Bias-aware intervals}{ + +Smoothing-parameter selection by curve-blocked neighbourhood +cross-validation (`pffr(method = "NCV")`) is robust to within-curve +dependence but smooths more than REML, so its smoothing bias is no longer +negligible and a variance-only interval undercovers. The bias-aware interval +adds the difference between the two estimates in quadrature: +\deqn{\delta = L(\hat\theta_{\mathrm{NCV}} - \hat\theta_{\mathrm{REML}}),}{ + delta = L (theta_NCV - theta_REML),} +for the rows \eqn{L} of the prediction matrix. The recommended recipe is +an NCV fit with curve blocks (the default `ncv_blocks = "cluster"`) as +`object`, the exact CL2 sandwich in its Bayesian form +(`sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE`) and a REML fit +of the same model as `bias_ref`. Fit both with `sandwich = "none"` to avoid +computing a fit-time sandwich that is not used. + +\eqn{\delta} estimates only the part of the smoothing bias in which the two +fits differ. Bias that both fits share -- a basis too small for the truth, +or both fits oversmoothing a rough truth -- is not covered. +} + +\examples{ +\donttest{ +set.seed(1) +d <- pffr_simulate(Y ~ ff(X1), n = 30, nxgrid = 15, nygrid = 25) +yind <- attr(d, "yindex") +fit_ncv <- pffr(Y ~ ff(X1), yind = yind, data = d, method = "NCV", + sandwich = "none") +fit_reml <- pffr(Y ~ ff(X1), yind = yind, data = d, method = "REML", + sandwich = "none") +ci <- pffr_predict_ci(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", + bias_ref = fit_reml) +head(ci) +} +} +\seealso{ +[coef.pffr()] (argument `bias_ref`) for bias-aware intervals of + coefficient functions, [pffr_jackknife_se()]. +} +\author{ +Fabian Scheipl +} diff --git a/tests/testthat/test-pffr-bias-aware.R b/tests/testthat/test-pffr-bias-aware.R new file mode 100644 index 00000000..ccd9645e --- /dev/null +++ b/tests/testthat/test-pffr-bias-aware.R @@ -0,0 +1,357 @@ +# Bias-aware pointwise intervals: coef.pffr(bias_ref = ) and pffr_predict_ci(). + +bias_aware_data <- function(family = gaussian()) { + dat <- if (family$family == "gaussian") { + ncv_test_data() + } else { + ncv_test_glm_data(family) + } + dat$z <- rnorm(nrow(dat)) + dat +} + +bias_aware_fit <- function(dat, method, family = gaussian()) { + pffr( + Y ~ + ff( + X, + splinepars = list(bs = "ps", m = list(c(2, 1), c(2, 1)), k = c(5, 5)) + ) + + c(z), + data = dat, + yind = seq(0, 1, length.out = ncol(dat$Y)), + family = family, + method = method, + bs.yindex = list(bs = "ps", k = 6, m = c(2, 1)), + bs.int = list(bs = "ps", k = 6, m = c(2, 1)), + sandwich = "none" + ) +} + +# Exact CL2 covariance in Bayesian form, as scored by the study. +bias_aware_cl2 <- function(fit) { + V <- pffr_vcov(fit, sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE) + expect_identical(attr(V, "cl2_adjustment"), "exact") + V +} + +# The study's interval from rows L of an estimand: NCV estimate, CL2 SE of the +# NCV fit and the NCV - REML difference in quadrature, z critical value. +bias_aware_direct <- function(L, fit_ncv, fit_reml, V, level = 0.95) { + est <- as.vector(L %*% fit_ncv$coefficients) + se <- unname(sqrt(rowSums((L %*% V) * L))) + delta <- as.vector(L %*% (fit_ncv$coefficients - fit_reml$coefficients)) + half <- qnorm((1 + level) / 2) * sqrt(se^2 + delta^2) + list( + est = est, + se = se, + delta = delta, + lower = est - half, + upper = est + half + ) +} + +# Rows of the ff coefficient surface on the grid coef.pffr() used (by = 1, no +# intercept columns: the seWithMean = FALSE convention). +bias_aware_ff_rows <- function(fit, grid) { + sm <- fit$smooth[[grep("X.smat", names(fit$smooth), fixed = TRUE)]] + L <- matrix(0, nrow(grid), length(fit$coefficients)) + L[, sm$first.para:sm$last.para] <- mgcv::PredictMat(sm, grid) + L +} + +bias_aware_newdata <- function(dat, rows = 1:4) { + list(X = I(dat$X[rows, ]), z = dat$z[rows]) +} + +expect_study_intervals <- function(family) { + set.seed(4211) + dat <- bias_aware_data(family) + fit_ncv <- bias_aware_fit(dat, "NCV", family) + fit_reml <- bias_aware_fit(dat, "REML", family) + V <- bias_aware_cl2(fit_ncv) + + # Coefficient surface of the ff term. + cf <- coef( + fit_ncv, + sandwich = "cl2", + cl2_adjustment = "exact", + ci = "pointwise", + seWithMean = FALSE, + bias_ref = fit_reml + ) + ff_coef <- cf$smterms[["ff(X)"]]$coef + ref <- bias_aware_direct( + bias_aware_ff_rows(fit_ncv, ff_coef[, c("X.smat", "X.tmat", "L.X")]), + fit_ncv, + fit_reml, + V + ) + expect_gt(max(abs(ref$delta)), 1e-4) + expect_equal(ff_coef$value, ref$est, tolerance = 1e-10) + expect_equal(ff_coef$se, ref$se, tolerance = 1e-10) + expect_equal(ff_coef$delta, ref$delta, tolerance = 1e-10) + expect_equal(ff_coef$lower, ref$lower, tolerance = 1e-10) + expect_equal(ff_coef$upper, ref$upper, tolerance = 1e-10) + expect_identical(cf$ci_meta$bias_ref_method, "REML") + + # Parametric coefficient of z. + j <- which(names(fit_ncv$coefficients) == "z") + ref_z <- bias_aware_direct( + diag(length(fit_ncv$coefficients))[j, , drop = FALSE], + fit_ncv, + fit_reml, + V + ) + expect_equal( + unname(cf$pterms["z", c("value", "se", "delta", "lower", "upper")]), + unlist(ref_z, use.names = FALSE), + tolerance = 1e-10 + ) + + # Conditional mean for new curves: link scale, then transformed endpoints. + newdata <- bias_aware_newdata(dat) + L <- predict(fit_ncv, newdata, type = "lpmatrix", reformat = FALSE) + ref_mean <- bias_aware_direct(L, fit_ncv, fit_reml, V) + args <- list( + fit_ncv, + newdata = newdata, + sandwich = "cl2", + cl2_adjustment = "exact", + bias_ref = fit_reml + ) + link <- do.call(pffr_predict_ci, c(args, type = "link")) + response <- do.call(pffr_predict_ci, c(args, type = "response")) + expect_equal(link$fit, ref_mean$est, tolerance = 1e-10) + expect_equal(link$se_link, ref_mean$se, tolerance = 1e-10) + expect_equal(link$delta_link, ref_mean$delta, tolerance = 1e-10) + expect_equal(link$lower, ref_mean$lower, tolerance = 1e-10) + expect_equal(link$upper, ref_mean$upper, tolerance = 1e-10) + linkinv <- fit_ncv$family$linkinv + expect_equal(response$fit, linkinv(ref_mean$est), tolerance = 1e-10) + expect_equal(response$lower, linkinv(ref_mean$lower), tolerance = 1e-10) + expect_equal(response$upper, linkinv(ref_mean$upper), tolerance = 1e-10) + expect_equal(response$se_link, link$se_link) + expect_equal(nrow(link), 4 * ncol(dat$Y)) + expect_equal(link$.obs, rep(1:4, each = ncol(dat$Y))) + expect_equal(link$.index, rep(fit_ncv$pffr$yind, 4)) + expect_identical(attr(link, "bias_ref_method"), "REML") +} + +test_that("bias-aware intervals reproduce the study's computation (Gaussian)", { + skip_if_not_installed("mgcv", "1.9.0") + expect_study_intervals(gaussian()) +}) + +test_that("bias-aware intervals reproduce the study's computation (Poisson)", { + skip_if_not_installed("mgcv", "1.9.0") + expect_study_intervals(poisson()) +}) + +test_that("a fit as its own bias reference gives delta = 0 and plain intervals", { + skip_if_not_installed("mgcv", "1.9.0") + set.seed(4212) + dat <- bias_aware_data() + fit <- bias_aware_fit(dat, "NCV") + args <- list( + fit, + sandwich = "cl2", + cl2_adjustment = "exact", + ci = "pointwise" + ) + plain <- suppressMessages(do.call(coef, args)) + self <- suppressMessages(do.call(coef, c(args, list(bias_ref = fit)))) + for (term in names(plain$smterms)) { + expect_true(all(self$smterms[[term]]$coef$delta == 0)) + expect_equal( + self$smterms[[term]]$coef[, c("value", "se", "lower", "upper")], + plain$smterms[[term]]$coef[, c("value", "se", "lower", "upper")] + ) + } + expect_true(all(self$pterms[, "delta"] == 0)) + expect_equal( + self$pterms[, c("lower", "upper")], + plain$pterms[, c("lower", "upper")] + ) + + pred_plain <- pffr_predict_ci(fit, sandwich = "cl2", cl2_adjustment = "exact") + pred_self <- pffr_predict_ci( + fit, + sandwich = "cl2", + cl2_adjustment = "exact", + bias_ref = fit + ) + expect_true(all(pred_self$delta_link == 0)) + for (column in names(pred_plain)) { + expect_equal(pred_self[[column]], pred_plain[[column]]) + } + expect_false("delta_link" %in% names(pred_plain)) +}) + +test_that("with seWithMean = TRUE delta uses the same linear map as the SE", { + skip_if_not_installed("mgcv", "1.9.0") + set.seed(4213) + # Poisson: for a Gaussian fit the mean level (mean fitted value) is the same + # for both fits, so the seWithMean shift would be zero. + dat <- bias_aware_data(poisson()) + fit_ncv <- bias_aware_fit(dat, "NCV", poisson()) + fit_reml <- bias_aware_fit(dat, "REML", poisson()) + get_intercept <- function(seWithMean) { + cf <- suppressMessages(coef( + fit_ncv, + sandwich = "cl2", + cl2_adjustment = "exact", + ci = "pointwise", + seWithMean = seWithMean, + bias_ref = fit_reml + )) + cf$smterms[["Intercept(yindex)"]]$coef + } + with_mean <- get_intercept(TRUE) + without <- get_intercept(FALSE) + sm <- fit_ncv$smooth[["s(yindex.vec)"]] + expect_gt(attr(sm, "nCons"), 0) + # The mean-level contrast: column means of the other model-matrix columns + # (incl. the scalar intercept) times the coefficient difference. + others <- setdiff(seq_along(fit_ncv$coefficients), sm$first.para:sm$last.para) + dtheta <- fit_ncv$coefficients - fit_reml$coefficients + shift <- sum(fit_ncv$cmX[others] * dtheta[others]) / (sm$meanL1 %||% 1) + expect_gt(abs(shift), 1e-6) + expect_equal(with_mean$delta - without$delta, rep(shift, nrow(with_mean))) + expect_equal( + with_mean$upper - with_mean$value, + qnorm(0.975) * sqrt(with_mean$se^2 + with_mean$delta^2) + ) +}) + +test_that("pffr_predict_ci passes covariance options through", { + skip_if_not_installed("mgcv", "1.9.0") + set.seed(4214) + dat <- bias_aware_data() + fit <- bias_aware_fit(dat, "NCV") + L <- predict(fit, type = "lpmatrix", reformat = FALSE) + for (freq in c(FALSE, TRUE)) { + V <- pffr_vcov(fit, sandwich = "cl2", cl2_adjustment = "exact", freq = freq) + ci <- pffr_predict_ci( + fit, + sandwich = "cl2", + cl2_adjustment = "exact", + freq = freq, + level = 0.9 + ) + se <- unname(sqrt(rowSums((L %*% V) * L))) + expect_equal(ci$se_link, se, tolerance = 1e-10) + expect_equal(ci$upper - ci$fit, qnorm(0.95) * se, tolerance = 1e-10) + expect_equal(ci$fit, as.vector(fit$linear.predictors)) + } + # Model-based covariance of an NCV fit is Vp. + ci_model <- pffr_predict_ci(fit, sandwich = "none") + expect_equal(ci_model$se_link, unname(sqrt(rowSums((L %*% fit$Vp) * L)))) +}) + +test_that("pffr_predict_ci aligns fitted points with missing responses", { + skip_if_not_installed("mgcv", "1.9.0") + set.seed(4215) + dat <- bias_aware_data() + dat$Y[cbind(c(2, 5, 5), c(3, 1, 16))] <- NA + fit_ncv <- bias_aware_fit(dat, "NCV") + fit_reml <- bias_aware_fit(dat, "REML") + ci <- pffr_predict_ci(fit_ncv, bias_ref = fit_reml, sandwich = "cl2") + n_grid <- ncol(dat$Y) + expect_equal(nrow(ci), length(dat$Y) - 3) + observed <- which(!is.na(t(dat$Y))) + expect_equal(ci$.obs, (observed - 1) %/% n_grid + 1) + expect_equal(ci$.index, fit_ncv$pffr$yind[(observed - 1) %% n_grid + 1]) + expect_equal(ci$fit, as.vector(fit_ncv$linear.predictors)) + expect_equal( + ci$delta_link, + as.vector(fit_ncv$linear.predictors - fit_reml$linear.predictors), + tolerance = 1e-8 + ) +}) + +test_that("pffr_predict_ci aligns sparse responses, skipping missing values", { + skip_if_not_installed("mgcv", "1.9.0") + set.seed(4217) + n <- 16 + dat <- data.frame(z = rnorm(n)) + yd <- data.frame( + .obs = rep(seq_len(n), times = rep(c(6, 8, 10, 12), n / 4)), + .index = runif(n * 9) + ) + yd$.value <- dat$z[yd$.obs] * sin(2 * pi * yd$.index) + rnorm(nrow(yd)) + yd <- yd[sample(nrow(yd)), ] + yd$.value[c(3, 40)] <- NA + fit_sparse <- function(method) { + pffr( + Y ~ z, + data = dat, + ydata = yd, + method = method, + bs.yindex = list(bs = "ps", k = 5, m = c(2, 1)), + bs.int = list(bs = "ps", k = 5, m = c(2, 1)), + sandwich = "none" + ) + } + # pffr's NCV does not accept missing sparse responses, so an ML fit stands in + # for the second fit here: the alignment does not depend on the method. + fit_ncv <- fit_sparse("REML") + fit_reml <- fit_sparse("ML") + # Model-based SEs: pffr_vcov() cannot yet build cluster ids for sparse fits + # with missing .value rows. + ci <- pffr_predict_ci(fit_ncv, bias_ref = fit_reml, sandwich = "none") + kept <- yd[!is.na(yd$.value), ] + expect_equal(nrow(ci), nrow(kept)) + expect_equal(ci$.obs, kept$.obs) + expect_equal(ci$.index, kept$.index) + expect_equal(ci$fit, as.vector(fit_ncv$linear.predictors)) + expect_gt(max(abs(ci$delta_link)), 1e-6) + expect_equal( + ci$delta_link, + as.vector(fit_ncv$linear.predictors - fit_reml$linear.predictors), + tolerance = 1e-8 + ) + expect_identical(attr(ci, "bias_ref_method"), "ML") +}) + +test_that("bias references of a different model or data are rejected", { + skip_if_not_installed("mgcv", "1.9.0") + set.seed(4216) + dat <- bias_aware_data() + fit <- bias_aware_fit(dat, "NCV") + other <- dat + other$X[1, ] <- other$X[1, ] + 1 + fit_other <- bias_aware_fit(other, "REML") + expect_error(coef(fit, bias_ref = fit_other), "same model and data") + expect_error( + pffr_predict_ci(fit, bias_ref = fit_other), + "same model and data" + ) + fit_k <- pffr( + Y ~ + ff( + X, + splinepars = list(bs = "ps", m = list(c(2, 1), c(2, 1)), k = c(5, 5)) + ) + + c(z), + data = dat, + yind = seq(0, 1, length.out = ncol(dat$Y)), + method = "REML", + bs.yindex = list(bs = "ps", k = 6, m = c(2, 1)), + bs.int = list(bs = "ps", k = 6, m = c(2, 1)), + knots = list(yindex.vec = seq(-1.2, 2.2, length.out = 10)), + sandwich = "none" + ) + expect_error(coef(fit, bias_ref = fit_k), "same model and data") + expect_error(coef(fit, bias_ref = fit$coefficients), "fitted pffr model") + expect_error( + coef(fit, ci = "simultaneous", bias_ref = fit), + "pointwise only" + ) + expect_error( + coef(fit, ci = "pointwise", crit = "tG1", bias_ref = fit), + "crit = \"z\"" + ) + expect_error(coef(fit, crit = "tG1", bias_ref = fit), "crit = \"z\"") + expect_error(coef(fit, raw = TRUE, bias_ref = fit), "raw = TRUE") +}) From 7a09e872b333d2b79faad36b17b7d73fb74f1056 Mon Sep 17 00:00:00 2001 From: fabian-s Date: Mon, 28 Sep 2026 17:48:33 +0200 Subject: [PATCH 2/3] pffr bias-aware intervals: exact-CL2 default, full intercept docs, test fixes - With `bias_ref`, coef.pffr() and pffr_predict_ci() default to the covariance the interval was evaluated with: sandwich = NULL -> "cl2", cl2_adjustment = NULL -> "exact" (Bayesian form via freq = FALSE). Explicit choices are respected; without bias_ref nothing changes. - Document the full functional intercept alpha(t): Intercept(yindex) is centred and the level sits in "(Intercept)"; pffr_predict_ci() at covariate values where all other terms vanish gives alpha(t) with its bias-aware interval (exact intercept rows; tested and in the example). - Document seWithMean for the bias-aware path: keep the default TRUE. It acts only on constrained terms (the intercept; ff() and varying coefficients are unconstrained by-variable smooths). In the pffr-ci study cells TRUE matched the exact full-intercept rows, FALSE undercovered alpha by 7-15 pp. - Add pffr-bias-aware.R to Collate (it was missing, so an installed build lacked pffr_predict_ci()). - test-pffr-ncv: the NCV covariance test no longer asserts mgcv's buggy Vc < Vp; it checks Vc != Vp and that refund returns Vp, under both unpatched and patched mgcv 1.9-5. - test-pffr-ar: the mgcv 1.9-5 skip helper called the nonexistent utils::package_version(); use base package_version(). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01JtmaCapnpxpQhDG2q2CQny --- DESCRIPTION | 1 + NEWS.md | 6 +- R/pffr-bias-aware.R | 54 ++++++++++++++--- R/pffr-methods.R | 54 +++++++++++++++-- man/coef.pffr.Rd | 49 ++++++++++++++-- man/pffr_bias_aware_cov_defaults.Rd | 22 +++++++ man/pffr_predict_ci.Rd | 27 +++++++-- tests/testthat/test-pffr-ar.R | 2 +- tests/testthat/test-pffr-bias-aware.R | 84 +++++++++++++++++++++++++++ tests/testthat/test-pffr-ncv.R | 7 ++- 10 files changed, 276 insertions(+), 30 deletions(-) create mode 100644 man/pffr_bias_aware_cov_defaults.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 7146172d..8d02a698 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -107,6 +107,7 @@ Collate: 'pffr-ff.R' 'pffr-ffpc.R' 'pffr-methods.R' + 'pffr-bias-aware.R' 'pffr-pcre.R' 'pffr-robust.R' 'pffr-sff.R' diff --git a/NEWS.md b/NEWS.md index da2d32b9..a21e57e8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,8 +7,10 @@ `est -/+ z * sqrt(se^2 + delta^2)` with `delta` the difference between the two estimates, built on the link scale with transformed endpoints for the response mean. Recommended: NCV with curve blocks, exact CL2 in Bayesian form - (`sandwich = "cl2", cl2_adjustment = "exact"`), REML reference. `delta` - does not capture smoothing bias that both fits share. + (`sandwich = "cl2", cl2_adjustment = "exact"`, the default covariance when + `bias_ref` is given), REML reference, default `seWithMean = TRUE`. `delta` + does not capture smoothing bias that both fits share. The docs show how to + get the full functional intercept (level included) and its interval. * `pffr(method = "NCV")` now leaves out whole curves (or `cluster` groups) by default. `ncv_blocks = "point"` provides pointwise comparisons; explicit `nei` takes precedence. Dual mgcv neighbourhood names, a cached behavioural diff --git a/R/pffr-bias-aware.R b/R/pffr-bias-aware.R index be4e257e..085c29dd 100644 --- a/R/pffr-bias-aware.R +++ b/R/pffr-bias-aware.R @@ -101,6 +101,25 @@ pffr_smooth_knots <- function(sm) { sm$knots } +#' Default covariance of bias-aware intervals +#' +#' Bias-aware intervals were evaluated with the exact CL2 sandwich in its +#' Bayesian form, so that is what they use unless the caller chose otherwise: +#' `sandwich = NULL` becomes `"cl2"`, and `cl2_adjustment = NULL` becomes +#' `"exact"` whenever the resolved sandwich is CL2. `freq` keeps its default +#' `FALSE` (Bayesian form) in the callers. +#' +#' @param sandwich,cl2_adjustment As supplied by the caller. +#' @returns A list with the resolved `sandwich` and `cl2_adjustment`. +#' @keywords internal +pffr_bias_aware_cov_defaults <- function(sandwich, cl2_adjustment) { + sandwich <- sandwich %||% "cl2" + if (identical(normalize_sandwich_type(sandwich), "cl2")) { + cl2_adjustment <- cl2_adjustment %||% "exact" + } + list(sandwich = sandwich, cl2_adjustment = cl2_adjustment) +} + #' Combine a variance-only standard error with a bias estimate #' #' @param se Standard errors (variance part). @@ -139,9 +158,15 @@ pffr_bias_aware_se <- function(se, delta) { #' for the rows \eqn{L} of the prediction matrix. The recommended recipe is #' an NCV fit with curve blocks (the default `ncv_blocks = "cluster"`) as #' `object`, the exact CL2 sandwich in its Bayesian form -#' (`sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE`) and a REML fit -#' of the same model as `bias_ref`. Fit both with `sandwich = "none"` to avoid -#' computing a fit-time sandwich that is not used. +#' (`sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE`, the default +#' covariance whenever `bias_ref` is supplied) and a REML fit of the same model +#' as `bias_ref`. Fit both with `sandwich = "none"` to avoid computing a +#' fit-time sandwich that is not used. +#' +#' The full functional intercept \eqn{\alpha(t)}{alpha(t)} (level included) +#' is the linear predictor at covariate values at which all other terms vanish, +#' e.g. `X = 0` for `ff(X)` and `z = 0` for a linear effect of `z`; see the +#' examples and [coef.pffr()]. #' #' \eqn{\delta} estimates only the part of the smoothing bias in which the two #' fits differ. Bias that both fits share -- a basis too small for the truth, @@ -157,7 +182,10 @@ pffr_bias_aware_se <- function(se, delta) { #' @param level Confidence level, defaults to `0.95`. #' @param sandwich,freq,cluster,dof_correction,edf_type,cl2_adjustment #' Covariance choice, as in [coef.pffr()]. `sandwich = NULL` inherits the -#' fit-time choice. +#' fit-time choice; with `bias_ref`, `sandwich = NULL` and +#' `cl2_adjustment = NULL` instead default to the exact CL2 sandwich (with the +#' default `freq = FALSE`, its Bayesian form), the covariance the bias-aware +#' interval was evaluated with. #' @param bias_ref Optional reference fit of the same model and data (typically #' the REML fit when `object` is the NCV fit), differing only in how the #' smoothing parameters were chosen. If supplied, intervals are bias-aware @@ -186,9 +214,15 @@ pffr_bias_aware_se <- function(se, delta) { #' sandwich = "none") #' fit_reml <- pffr(Y ~ ff(X1), yind = yind, data = d, method = "REML", #' sandwich = "none") -#' ci <- pffr_predict_ci(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", -#' bias_ref = fit_reml) +#' # exact CL2 (Bayesian form) is the default covariance with bias_ref +#' ci <- pffr_predict_ci(fit_ncv, bias_ref = fit_reml) #' head(ci) +#' +#' # Full functional intercept alpha(t) with its bias-aware interval: the +#' # linear predictor of a curve with X1 = 0, where the ff() term vanishes. +#' zero <- data.frame(X1 = I(matrix(0, 1, ncol(d$X1)))) +#' alpha <- pffr_predict_ci(fit_ncv, newdata = zero, bias_ref = fit_reml) +#' head(alpha[, c(".index", "fit", "lower", "upper")]) #' } pffr_predict_ci <- function( object, @@ -222,8 +256,12 @@ pffr_predict_ci <- function( call. = FALSE ) } - theta_diff <- if (!is.null(bias_ref)) { - pffr_bias_ref_difference(object, bias_ref) + theta_diff <- NULL + if (!is.null(bias_ref)) { + theta_diff <- pffr_bias_ref_difference(object, bias_ref) + cov_defaults <- pffr_bias_aware_cov_defaults(sandwich, cl2_adjustment) + sandwich <- cov_defaults$sandwich + cl2_adjustment <- cov_defaults$cl2_adjustment } lp <- pffr_prediction_design(object, newdata) diff --git a/R/pffr-methods.R b/R/pffr-methods.R index aeb92d9a..95657748 100755 --- a/R/pffr-methods.R +++ b/R/pffr-methods.R @@ -1399,7 +1399,12 @@ pffr_end_undefined_df <- function(opened) { #' applies the documented feasibility rule, while \code{"exact"} and #' \code{"shortcut"} force either CL2 variant. Supplying a value forces a #' covariance recomputation. -#' @param seWithMean logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in standard errors returned for smooth components? +#' @param seWithMean logical, defaults to TRUE. Include uncertainty about the +#' intercept/overall mean in standard errors returned for smooth components? +#' Acts only on terms with an identifiability constraint (e.g. +#' \code{Intercept(yindex)}); \code{ff()} terms and varying coefficients of +#' scalar covariates are unconstrained by-variable smooths and are unaffected. +#' With \code{bias_ref}, \code{delta} uses the same linear map as the SE. #' @param n1 see below #' @param n2 see below #' @param n3 \code{n1, n2, n3} give the number of gridpoints for 1-/2-/3-dimensional smooth terms @@ -1458,8 +1463,13 @@ pffr_end_undefined_df <- function(opened) { #' difference in the mean level), and pointwise intervals use #' \eqn{\sqrt{se^2 + delta^2}}{sqrt(se^2 + delta^2)} instead of \code{se}; #' \code{se} itself stays the variance part. Only \code{ci = "none"} and -#' \code{ci = "pointwise"} with \code{crit = "z"} are available. See the -#' section \sQuote{Bias-aware intervals}. +#' \code{ci = "pointwise"} with \code{crit = "z"} are available. With +#' \code{bias_ref}, \code{sandwich = NULL} and \code{cl2_adjustment = NULL} +#' default to the exact CL2 sandwich (\code{"cl2"}, \code{"exact"}) instead +#' of the fit-time choice; together with the default \code{freq = FALSE} +#' (Bayesian form) this is the covariance the interval was evaluated with. +#' Explicit values are respected. See the section +#' \sQuote{Bias-aware intervals}. #' @param ... other arguments, not used. #' #' @return If \code{raw==FALSE}, a list containing \itemize{ @@ -1496,9 +1506,38 @@ pffr_end_undefined_df <- function(opened) { #' on the link scale. The recommended recipe fits the model twice with #' \code{sandwich = "none"}, by NCV with curve blocks (the default #' \code{ncv_blocks = "cluster"}) and by REML, and computes \code{se} from the -#' exact CL2 sandwich in its Bayesian form of the NCV fit: -#' \preformatted{coef(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", +#' exact CL2 sandwich in its Bayesian form of the NCV fit, which is the default +#' covariance whenever \code{bias_ref} is supplied: +#' \preformatted{coef(fit_ncv, ci = "pointwise", bias_ref = fit_reml) +#' # same as +#' coef(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE, #' ci = "pointwise", bias_ref = fit_reml)} +#' +#' Keep the default \code{seWithMean = TRUE} for bias-aware intervals. It +#' changes only constrained terms, of which the functional intercept is the +#' main one: in the simulation study behind this recipe (full intercept +#' scored, dependent errors, Gaussian, Poisson and binary responses) the +#' \code{seWithMean = TRUE} intervals were practically identical to those from +#' the exact full-intercept rows below (root-mean-square undercoverage 0 to 3 +#' percentage points), while \code{seWithMean = FALSE} omits the level's +#' uncertainty and undercovered by 7 to 15 points, with 11 to 22\% larger +#' interval scores. +#' +#' \strong{Functional intercept.} The \code{Intercept(yindex)} term is centred +#' (sum-to-zero over the observed response grid); the level of +#' \eqn{\alpha(t)}{alpha(t)} is the scalar \code{"(Intercept)"} in +#' \code{pterms}. The full intercept is their sum, +#' \preformatted{alpha_hat <- cf$smterms[["Intercept(yindex)"]]$coef[, "value"] + +#' cf$pterms["(Intercept)", "value"]} +#' but the \code{se}, \code{delta} and interval of the centred term do not +#' describe it: they omit the level (\code{seWithMean = FALSE}) or add the +#' column means of all other terms (\code{seWithMean = TRUE}). For an interval +#' for the full \eqn{\alpha(t)}{alpha(t)} use \code{\link{pffr_predict_ci}} at +#' covariate values at which every other term vanishes (e.g. \code{X = 0} for +#' \code{ff(X)} and \code{z = 0} for a linear effect of a scalar \code{z}); its +#' prediction rows are then exactly the intercept rows +#' \eqn{(1, B(t))}{(1, B(t))}. See the examples of +#' \code{\link{pffr_predict_ci}}. #' \eqn{\delta}{delta} estimates only the part of the smoothing bias in which #' the two fits differ: bias that both fits share (a basis too small for the #' truth, or both fits oversmoothing a rough truth) is not covered. Any second @@ -1542,6 +1581,11 @@ coef.pffr <- function( on.exit(pffr_end_undefined_df(df_warn_window), add = TRUE) sandwich_missing <- missing(sandwich) + if (!is.null(bias_ref)) { + cov_defaults <- pffr_bias_aware_cov_defaults(sandwich, cl2_adjustment) + sandwich <- cov_defaults$sandwich + cl2_adjustment <- cov_defaults$cl2_adjustment + } # Backward compat: TRUE -> "cluster", FALSE -> "none" if (is.logical(sandwich)) sandwich <- if (sandwich) "cluster" else "none" if (is.null(sandwich)) sandwich <- pffr_canonicalize_cov(object)$fit_type diff --git a/man/coef.pffr.Rd b/man/coef.pffr.Rd index c8e0abfa..96d3baf3 100644 --- a/man/coef.pffr.Rd +++ b/man/coef.pffr.Rd @@ -83,7 +83,12 @@ applies the documented feasibility rule, while \code{"exact"} and \code{"shortcut"} force either CL2 variant. Supplying a value forces a covariance recomputation.} -\item{seWithMean}{logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in standard errors returned for smooth components?} +\item{seWithMean}{logical, defaults to TRUE. Include uncertainty about the +intercept/overall mean in standard errors returned for smooth components? +Acts only on terms with an identifiability constraint (e.g. +\code{Intercept(yindex)}); \code{ff()} terms and varying coefficients of +scalar covariates are unconstrained by-variable smooths and are unaffected. +With \code{bias_ref}, \code{delta} uses the same linear map as the SE.} \item{n1}{see below} @@ -153,8 +158,13 @@ linear map as \code{se}, so with \code{seWithMean = TRUE} it includes the difference in the mean level), and pointwise intervals use \eqn{\sqrt{se^2 + delta^2}}{sqrt(se^2 + delta^2)} instead of \code{se}; \code{se} itself stays the variance part. Only \code{ci = "none"} and -\code{ci = "pointwise"} with \code{crit = "z"} are available. See the -section \sQuote{Bias-aware intervals}.} +\code{ci = "pointwise"} with \code{crit = "z"} are available. With +\code{bias_ref}, \code{sandwich = NULL} and \code{cl2_adjustment = NULL} +default to the exact CL2 sandwich (\code{"cl2"}, \code{"exact"}) instead +of the fit-time choice; together with the default \code{freq = FALSE} +(Bayesian form) this is the covariance the interval was evaluated with. +Explicit values are respected. See the section +\sQuote{Bias-aware intervals}.} \item{...}{other arguments, not used.} } @@ -211,9 +221,38 @@ interval for an estimate \eqn{L\hat\theta_{NCV}}{L theta_NCV} is on the link scale. The recommended recipe fits the model twice with \code{sandwich = "none"}, by NCV with curve blocks (the default \code{ncv_blocks = "cluster"}) and by REML, and computes \code{se} from the -exact CL2 sandwich in its Bayesian form of the NCV fit: -\preformatted{coef(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", +exact CL2 sandwich in its Bayesian form of the NCV fit, which is the default +covariance whenever \code{bias_ref} is supplied: +\preformatted{coef(fit_ncv, ci = "pointwise", bias_ref = fit_reml) +# same as +coef(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE, ci = "pointwise", bias_ref = fit_reml)} + +Keep the default \code{seWithMean = TRUE} for bias-aware intervals. It +changes only constrained terms, of which the functional intercept is the +main one: in the simulation study behind this recipe (full intercept +scored, dependent errors, Gaussian, Poisson and binary responses) the +\code{seWithMean = TRUE} intervals were practically identical to those from +the exact full-intercept rows below (root-mean-square undercoverage 0 to 3 +percentage points), while \code{seWithMean = FALSE} omits the level's +uncertainty and undercovered by 7 to 15 points, with 11 to 22\% larger +interval scores. + +\strong{Functional intercept.} The \code{Intercept(yindex)} term is centred +(sum-to-zero over the observed response grid); the level of +\eqn{\alpha(t)}{alpha(t)} is the scalar \code{"(Intercept)"} in +\code{pterms}. The full intercept is their sum, +\preformatted{alpha_hat <- cf$smterms[["Intercept(yindex)"]]$coef[, "value"] + + cf$pterms["(Intercept)", "value"]} +but the \code{se}, \code{delta} and interval of the centred term do not +describe it: they omit the level (\code{seWithMean = FALSE}) or add the +column means of all other terms (\code{seWithMean = TRUE}). For an interval +for the full \eqn{\alpha(t)}{alpha(t)} use \code{\link{pffr_predict_ci}} at +covariate values at which every other term vanishes (e.g. \code{X = 0} for +\code{ff(X)} and \code{z = 0} for a linear effect of a scalar \code{z}); its +prediction rows are then exactly the intercept rows +\eqn{(1, B(t))}{(1, B(t))}. See the examples of +\code{\link{pffr_predict_ci}}. \eqn{\delta}{delta} estimates only the part of the smoothing bias in which the two fits differ: bias that both fits share (a basis too small for the truth, or both fits oversmoothing a rough truth) is not covered. Any second diff --git a/man/pffr_bias_aware_cov_defaults.Rd b/man/pffr_bias_aware_cov_defaults.Rd new file mode 100644 index 00000000..9e4a62b1 --- /dev/null +++ b/man/pffr_bias_aware_cov_defaults.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-bias-aware.R +\name{pffr_bias_aware_cov_defaults} +\alias{pffr_bias_aware_cov_defaults} +\title{Default covariance of bias-aware intervals} +\usage{ +pffr_bias_aware_cov_defaults(sandwich, cl2_adjustment) +} +\arguments{ +\item{sandwich, cl2_adjustment}{As supplied by the caller.} +} +\value{ +A list with the resolved `sandwich` and `cl2_adjustment`. +} +\description{ +Bias-aware intervals were evaluated with the exact CL2 sandwich in its +Bayesian form, so that is what they use unless the caller chose otherwise: +`sandwich = NULL` becomes `"cl2"`, and `cl2_adjustment = NULL` becomes +`"exact"` whenever the resolved sandwich is CL2. `freq` keeps its default +`FALSE` (Bayesian form) in the callers. +} +\keyword{internal} diff --git a/man/pffr_predict_ci.Rd b/man/pffr_predict_ci.Rd index a1dd3474..ceda877d 100644 --- a/man/pffr_predict_ci.Rd +++ b/man/pffr_predict_ci.Rd @@ -32,7 +32,10 @@ conditional mean.} \item{level}{Confidence level, defaults to `0.95`.} \item{sandwich, freq, cluster, dof_correction, edf_type, cl2_adjustment}{Covariance choice, as in [coef.pffr()]. `sandwich = NULL` inherits the -fit-time choice.} +fit-time choice; with `bias_ref`, `sandwich = NULL` and +`cl2_adjustment = NULL` instead default to the exact CL2 sandwich (with the +default `freq = FALSE`, its Bayesian form), the covariance the bias-aware +interval was evaluated with.} \item{bias_ref}{Optional reference fit of the same model and data (typically the REML fit when `object` is the NCV fit), differing only in how the @@ -79,9 +82,15 @@ adds the difference between the two estimates in quadrature: for the rows \eqn{L} of the prediction matrix. The recommended recipe is an NCV fit with curve blocks (the default `ncv_blocks = "cluster"`) as `object`, the exact CL2 sandwich in its Bayesian form -(`sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE`) and a REML fit -of the same model as `bias_ref`. Fit both with `sandwich = "none"` to avoid -computing a fit-time sandwich that is not used. +(`sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE`, the default +covariance whenever `bias_ref` is supplied) and a REML fit of the same model +as `bias_ref`. Fit both with `sandwich = "none"` to avoid computing a +fit-time sandwich that is not used. + +The full functional intercept \eqn{\alpha(t)}{alpha(t)} (level included) +is the linear predictor at covariate values at which all other terms vanish, +e.g. `X = 0` for `ff(X)` and `z = 0` for a linear effect of `z`; see the +examples and [coef.pffr()]. \eqn{\delta} estimates only the part of the smoothing bias in which the two fits differ. Bias that both fits share -- a basis too small for the truth, @@ -97,9 +106,15 @@ fit_ncv <- pffr(Y ~ ff(X1), yind = yind, data = d, method = "NCV", sandwich = "none") fit_reml <- pffr(Y ~ ff(X1), yind = yind, data = d, method = "REML", sandwich = "none") -ci <- pffr_predict_ci(fit_ncv, sandwich = "cl2", cl2_adjustment = "exact", - bias_ref = fit_reml) +# exact CL2 (Bayesian form) is the default covariance with bias_ref +ci <- pffr_predict_ci(fit_ncv, bias_ref = fit_reml) head(ci) + +# Full functional intercept alpha(t) with its bias-aware interval: the +# linear predictor of a curve with X1 = 0, where the ff() term vanishes. +zero <- data.frame(X1 = I(matrix(0, 1, ncol(d$X1)))) +alpha <- pffr_predict_ci(fit_ncv, newdata = zero, bias_ref = fit_reml) +head(alpha[, c(".index", "fit", "lower", "upper")]) } } \seealso{ diff --git a/tests/testthat/test-pffr-ar.R b/tests/testthat/test-pffr-ar.R index 97be927e..854a5404 100644 --- a/tests/testthat/test-pffr-ar.R +++ b/tests/testthat/test-pffr-ar.R @@ -3,7 +3,7 @@ ############################################################################### skip_if_mgcv_1_9_5_binomial_ar <- function() { - if (utils::packageVersion("mgcv") == utils::package_version("1.9-5")) { + if (utils::packageVersion("mgcv") == package_version("1.9-5")) { skip("Skipping known mgcv 1.9-5 binomial+rho segfault case.") } } diff --git a/tests/testthat/test-pffr-bias-aware.R b/tests/testthat/test-pffr-bias-aware.R index ccd9645e..bbe8422e 100644 --- a/tests/testthat/test-pffr-bias-aware.R +++ b/tests/testthat/test-pffr-bias-aware.R @@ -249,6 +249,90 @@ test_that("pffr_predict_ci passes covariance options through", { expect_equal(ci_model$se_link, unname(sqrt(rowSums((L %*% fit$Vp) * L)))) }) +test_that("bias_ref defaults to exact Bayesian CL2 unless overridden", { + skip_if_not_installed("mgcv", "1.9.0") + set.seed(4218) + dat <- bias_aware_data() + fit_ncv <- bias_aware_fit(dat, "NCV") + fit_reml <- bias_aware_fit(dat, "REML") + explicit <- list(sandwich = "cl2", cl2_adjustment = "exact", freq = FALSE) + ff_se <- function(...) { + cf <- suppressMessages(coef( + fit_ncv, + ci = "pointwise", + seWithMean = FALSE, + bias_ref = fit_reml, + ... + )) + cf$smterms[["ff(X)"]]$coef$se + } + expect_equal(ff_se(), do.call(ff_se, explicit)) + expect_equal( + pffr_predict_ci(fit_ncv, bias_ref = fit_reml), + do.call(pffr_predict_ci, c(list(fit_ncv, bias_ref = fit_reml), explicit)) + ) + # Explicit choices are respected: the fit-time model-based covariance, the + # CL2 shortcut and the frequentist form each change the SE. + for (override in list( + list(sandwich = "none"), + list(sandwich = "cl2", cl2_adjustment = "shortcut"), + list(freq = TRUE) + )) { + expect_false(isTRUE(all.equal(do.call(ff_se, override), ff_se()))) + } + L <- predict(fit_ncv, type = "lpmatrix", reformat = FALSE) + model_based <- pffr_predict_ci( + fit_ncv, + bias_ref = fit_reml, + sandwich = "none" + ) + expect_equal( + model_based$se_link, + unname(sqrt(rowSums((L %*% fit_ncv$Vp) * L))) + ) + # Without bias_ref, NULL still inherits the fit-time (here model-based) choice. + expect_equal(pffr_predict_ci(fit_ncv)$se_link, model_based$se_link) +}) + +test_that("pffr_predict_ci at zero covariates gives the full intercept alpha(t)", { + skip_if_not_installed("mgcv", "1.9.0") + set.seed(4219) + dat <- bias_aware_data() + fit_ncv <- bias_aware_fit(dat, "NCV") + fit_reml <- bias_aware_fit(dat, "REML") + zero <- list(X = I(matrix(0, 1, ncol(dat$X))), z = 0) + alpha <- pffr_predict_ci(fit_ncv, newdata = zero, bias_ref = fit_reml) + # The full intercept's rows: the scalar intercept plus the centred + # Intercept(yindex) basis at the response grid (the study's alpha rows). + sm <- fit_ncv$smooth[["s(yindex.vec)"]] + yind <- fit_ncv$pffr$yind + L <- matrix(0, length(yind), length(fit_ncv$coefficients)) + L[, sm$first.para:sm$last.para] <- mgcv::PredictMat( + sm, + data.frame(yindex.vec = yind) + ) + L[, names(fit_ncv$coefficients) == "(Intercept)"] <- 1 + ref <- bias_aware_direct(L, fit_ncv, fit_reml, bias_aware_cl2(fit_ncv)) + expect_equal(alpha$.index, yind) + expect_equal(alpha$fit, ref$est, tolerance = 1e-10) + expect_equal(alpha$se_link, ref$se, tolerance = 1e-10) + expect_equal(alpha$delta_link, ref$delta, tolerance = 1e-10) + expect_equal(alpha$lower, ref$lower, tolerance = 1e-10) + expect_equal(alpha$upper, ref$upper, tolerance = 1e-10) + # Its estimate is the centred term plus the scalar intercept from coef(). + cf <- suppressMessages(coef( + fit_ncv, + bias_ref = fit_reml, + eval_grid = list(`Intercept(yindex)` = data.frame(yindex.vec = yind)) + )) + expect_equal( + alpha$fit, + cf$smterms[["Intercept(yindex)"]]$coef$value + + cf$pterms["(Intercept)", "value"], + tolerance = 1e-10 + ) +}) + test_that("pffr_predict_ci aligns fitted points with missing responses", { skip_if_not_installed("mgcv", "1.9.0") set.seed(4215) diff --git a/tests/testthat/test-pffr-ncv.R b/tests/testthat/test-pffr-ncv.R index 8fb8d471..0b931f40 100644 --- a/tests/testthat/test-pffr-ncv.R +++ b/tests/testthat/test-pffr-ncv.R @@ -248,9 +248,10 @@ test_that("model-based covariance of NCV fits is Vp, not mgcv's Vc", { set.seed(102) dat <- ncv_test_data() ncv <- ncv_test_fit(dat) - # mgcv's Vc for NCV fits is far smaller than Vp, so it cannot be a covariance of the - # estimate that adds smoothing-parameter uncertainty to Vp - expect_lt(median(sqrt(diag(ncv$Vc)) / sqrt(diag(ncv$Vp))), 0.9) + # Unpatched mgcv 1.9-5 returns only the smoothing-parameter correction as Vc for + # NCV fits (far smaller than Vp); a fixed mgcv returns Vp plus that correction. + # Either way Vc differs from Vp, and refund must return Vp. + expect_false(isTRUE(all.equal(ncv$Vc, ncv$Vp))) expect_identical(pffr_vcov(ncv, sandwich = "none"), ncv$Vp) expect_identical(pffr_vcov(ncv, sandwich = "none", freq = TRUE), ncv$Ve) reml <- pffr( From 3b422c6f394ca2bc7a01a0746db6556c30c7f5fa Mon Sep 17 00:00:00 2001 From: fabian-s Date: Tue, 29 Sep 2026 09:40:37 +0200 Subject: [PATCH 3/3] pffr bias_ref: reject multi-linear-predictor families in coef() too The check now lives in pffr_bias_ref_difference(), shared by coef.pffr() and pffr_predict_ci(). Addresses Copilot review on #132. Co-Authored-By: Claude Opus 5.5 --- R/pffr-bias-aware.R | 10 +++++++++- man/pffr_bias_ref_difference.Rd | 3 ++- tests/testthat/test-pffr-bias-aware.R | 7 +++++++ 3 files changed, 18 insertions(+), 2 deletions(-) diff --git a/R/pffr-bias-aware.R b/R/pffr-bias-aware.R index 085c29dd..a95a2c9d 100644 --- a/R/pffr-bias-aware.R +++ b/R/pffr-bias-aware.R @@ -18,7 +18,8 @@ #' (labels, classes, basis dimensions, knots and penalty matrices, which carry the #' identifiability constraints), `ffpc`/`pcre` metadata, and family and link. #' Family parameters estimated during fitting (e.g. the negative binomial -#' \eqn{\theta}{theta}) may differ between the fits. +#' \eqn{\theta}{theta}) may differ between the fits. Fits with several linear +#' predictors (e.g. `gaulss`) are rejected. #' #' @param object The fit whose intervals are computed. #' @param bias_ref The reference fit. @@ -28,6 +29,13 @@ pffr_bias_ref_difference <- function(object, bias_ref) { if (!inherits(bias_ref, "pffr")) { stop("`bias_ref` must be a fitted pffr model.", call. = FALSE) } + if (isTRUE((object$family$nlp %||% 1L) > 1L)) { + stop( + "Bias-aware intervals (`bias_ref`) support single linear-predictor ", + "families only.", + call. = FALSE + ) + } mismatch <- function(what) { stop( "`bias_ref` is not a fit of the same model and data as `object` (", diff --git a/man/pffr_bias_ref_difference.Rd b/man/pffr_bias_ref_difference.Rd index 0e6eb764..25dca252 100644 --- a/man/pffr_bias_ref_difference.Rd +++ b/man/pffr_bias_ref_difference.Rd @@ -23,6 +23,7 @@ matrices, in the same row order), prior weights and offsets, the smooth bases (labels, classes, basis dimensions, knots and penalty matrices, which carry the identifiability constraints), `ffpc`/`pcre` metadata, and family and link. Family parameters estimated during fitting (e.g. the negative binomial -\eqn{\theta}{theta}) may differ between the fits. +\eqn{\theta}{theta}) may differ between the fits. Fits with several linear +predictors (e.g. `gaulss`) are rejected. } \keyword{internal} diff --git a/tests/testthat/test-pffr-bias-aware.R b/tests/testthat/test-pffr-bias-aware.R index bbe8422e..10131a1c 100644 --- a/tests/testthat/test-pffr-bias-aware.R +++ b/tests/testthat/test-pffr-bias-aware.R @@ -439,3 +439,10 @@ test_that("bias references of a different model or data are rejected", { expect_error(coef(fit, crit = "tG1", bias_ref = fit), "crit = \"z\"") expect_error(coef(fit, raw = TRUE, bias_ref = fit), "raw = TRUE") }) + +test_that("bias references are rejected for multi-linear-predictor families", { + skip_if_not_installed("mgcv", "1.9.0") + fit <- get_gaulss_model() + expect_error(coef(fit, bias_ref = fit), "single linear-predictor") + expect_error(pffr_predict_ci(fit, bias_ref = fit), "single linear-predictor") +})