R/bootstrap.R

In ERRI: Economic Resilience and Recovery Index

#' Bootstrap uncertainty for ERRI
#'
#' Uses a residual bootstrap of the pre-shock counterfactual model. The
#' observed post-shock path is held fixed while counterfactual estimation
#' uncertainty is propagated to the component measures and index.
#'
#' @param object An `erri` object.
#' @param R Number of bootstrap replications.
#' @param block_length Positive integer residual-block length. A value of one
#'   gives an ordinary residual bootstrap.
#' @param level Confidence level.
#' @param seed Optional integer seed.
#' @return An object of class `erri_bootstrap` with replicate estimates and
#'   percentile confidence intervals.
#' @export
#' @examples
#' fit <- erri(erri_example_data(), "year", "income", 2020, "region")
#' boot <- erri_bootstrap(fit, R = 49, seed = 1)
#' boot
erri_bootstrap <- function(object, R = 499L, block_length = 1L,
                           level = 0.95, seed = NULL) {
    if (!inherits(object, "erri")) .erri_stop("'object' must inherit from 'erri'.")
    R <- as.integer(R)
    block_length <- as.integer(block_length)
    if (is.na(R) || R < 20L) .erri_stop("Use at least 20 bootstrap replications.")
    if (is.na(block_length) || block_length < 1L) .erri_stop("Invalid block length.")
    if (!is.numeric(level) || length(level) != 1L || level <= 0 || level >= 1) {
        .erri_stop("'level' must be between zero and one.")
    }
    if (!is.null(seed)) set.seed(seed)
    alpha <- (1 - level) / 2
    reps <- vector("list", length(object$trajectories))
    names(reps) <- names(object$trajectories)

    for (id in names(object$trajectories)) {
        tr <- object$trajectories[[id]]
        model <- object$models[[id]]
        post <- which(tr$post_shock)
        sidx <- post[1L]
        residuals <- model$residuals[is.finite(model$residuals)]
        if (length(residuals) < 3L) .erri_stop("Too few residuals for unit '", id, "'.")
        z <- matrix(NA_real_, R, 6L,
                    dimnames = list(NULL, c("shock_depth", "cumulative_loss",
                                            "recovery_time", "recovery_strength",
                                            "stability_ratio", "transformation")))
        sc <- .scale_value(tr$observed[seq_len(sidx - 1L)], object$settings$scale)
        for (b in seq_len(R)) {
            starts <- sample(seq_along(residuals),
                             ceiling(length(tr$observed) / block_length), replace = TRUE)
            draw <- unlist(lapply(starts, function(k) {
                residuals[((k - 1L + seq_len(block_length) - 1L) %% length(residuals)) + 1L]
            }), use.names = FALSE)[seq_along(tr$observed)]
            cf <- tr$counterfactual + draw
            gap <- (cf - tr$observed) / sc
            pgap <- gap[post]
            adverse <- pmax(pgap, 0)
            depth <- max(adverse, na.rm = TRUE)
            loss <- sum(adverse, na.rm = TRUE)
            rt <- .run_recovery(pgap, object$settings$epsilon,
                                object$settings$consecutive)
            trough <- which.max(adverse)
            end <- if (is.na(rt)) length(pgap) else max(rt + 1L, trough)
            strength <- if (end <= trough) 0 else
                pmax(depth - adverse[end], 0) / (end - trough)
            pre_res <- tr$observed[seq_len(sidx - 1L)] - cf[seq_len(sidx - 1L)]
            post_res <- tr$observed[post] - cf[post]
            denom <- stats::sd(pre_res, na.rm = TRUE)
            stability <- stats::sd(post_res, na.rm = TRUE) /
                ifelse(is.finite(denom) && denom > 0, denom, sc)
            ts <- if (is.na(rt)) length(pgap) else min(rt + 1L, length(pgap))
            trans <- mean(pmax(-pgap[ts:length(pgap)], 0), na.rm = TRUE)
            if (!is.finite(trans)) trans <- 0
            one <- data.frame(shock_depth = depth, cumulative_loss = loss,
                              recovery_time = if (is.na(rt)) length(pgap) else rt,
                              recovery_strength = strength,
                              stability_ratio = stability,
                              transformation = trans)
            z[b, ] <- unlist(one, use.names = FALSE)
        }
        reps[[id]] <- as.data.frame(z)
    }
    for (b in seq_len(R)) {
        joint <- do.call(rbind, lapply(reps, function(z) z[b, , drop = FALSE]))
        joint_scores <- .component_scores(joint)
        joint_index <- 100 * as.numeric(as.matrix(joint_scores) %*% object$weights)
        for (i in seq_along(reps)) reps[[i]]$ERRI[b] <- joint_index[i]
    }
    all_rep <- do.call(rbind, lapply(names(reps), function(id) {
        cbind(unit = id, replication = seq_len(R), reps[[id]],
              stringsAsFactors = FALSE)
    }))
    intervals <- do.call(rbind, lapply(names(reps), function(id) {
        z <- reps[[id]]
        est <- unname(unlist(object$results[object$results$unit == id,
                                          names(z)[names(z) != "ERRI"],
                                          drop = FALSE], use.names = FALSE))
        data.frame(unit = id, measure = names(z),
                   estimate = c(est,
                                object$results$ERRI[object$results$unit == id]),
                   lower = vapply(z, stats::quantile, numeric(1), probs = alpha,
                                  na.rm = TRUE, names = FALSE),
                   upper = vapply(z, stats::quantile, numeric(1), probs = 1 - alpha,
                                  na.rm = TRUE, names = FALSE),
                   stringsAsFactors = FALSE)
    }))
    rownames(intervals) <- NULL
    out <- list(intervals = intervals, replicates = all_rep, level = level,
                R = R, block_length = block_length, estimate = object)
    class(out) <- "erri_bootstrap"
    out
}

#' @export
print.erri_bootstrap <- function(x, digits = 2, ...) {
    cat("ERRI residual bootstrap\n")
    cat("Replications:", x$R, "| Confidence level:", x$level, "\n\n")
    print(x$intervals, digits = digits, row.names = FALSE)
    invisible(x)
}

#' Ranking probabilities from bootstrap estimates
#'
#' Calculates pairwise probabilities that one unit's ERRI exceeds another's.
#'
#' @param object An object returned by [erri_bootstrap()].
#' @return A square matrix of pairwise probabilities.
#' @export
#' @examples
#' fit <- erri(erri_example_data(), "year", "income", 2020, "region")
#' boot <- erri_bootstrap(fit, R = 49, seed = 2)
#' rank_probability(boot)
rank_probability <- function(object) {
    if (!inherits(object, "erri_bootstrap")) {
        .erri_stop("'object' must inherit from 'erri_bootstrap'.")
    }
    ids <- unique(object$replicates$unit)
    out <- matrix(NA_real_, length(ids), length(ids), dimnames = list(ids, ids))
    for (i in seq_along(ids)) for (j in seq_along(ids)) {
        if (i == j) out[i, j] <- 0.5 else {
            a <- object$replicates$ERRI[object$replicates$unit == ids[i]]
            b <- object$replicates$ERRI[object$replicates$unit == ids[j]]
            n <- min(length(a), length(b))
            out[i, j] <- mean(a[seq_len(n)] > b[seq_len(n)], na.rm = TRUE)
        }
    }
    out
}