Nothing
R/network.interpolation.R
In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data
Defines functions
.network_interp_empty_plot
.network_interp_validation_seasonal_subset
.network_interp_validation_points_subset
.network_interp_validation_subset
.network_interp_diag_subset
.plot_network_seasonal_error
.plot_network_compare
.plot_network_coverage
.plot_network_validation
.plot_network_summary
.plot_network_availability
.plot_network_uncertainty
.plot_network_interpolation
plot.network_interpolation
network.interpolation
Documented in
network.interpolation
plot.network_interpolation
#' @title Interpolate missing dendrometer values using a site network
#'
#' @description
#' Fills \code{NA} gaps in one or more focal dendrometer series by borrowing
#' information from a reference network recorded at the same site and on the
#' same time grid.
#'
#' The focal dataset \code{df} is treated as the master time axis. The
#' reference network is aligned to that axis, so missing timestamps already
#' inserted into \code{df} are preserved during interpolation.
#'
#' Two interpolation modes are available:
#' \itemize{
#' \item \code{niMethod = "proportional"} uses the average relative change in
#' the reference network between consecutive timestamps and propagates that
#' change to the focal series.
#' \item \code{niMethod = "linear"} fits a cross-sensor regression between
#' reference values at time \eqn{t-1} and \eqn{t}, and predicts the focal
#' value at time \eqn{t} from the focal value at time \eqn{t-1}.
#' }
#'
#' Optional bootstrap-based uncertainty limits can be attached, synthetic-gap
#' validation can be used to assess recovery performance, and an optional
#' post-processing step can detect and correct abrupt post-gap jumps in the
#' interpolated series.
#'
#' @details
#' The function assumes that missing timestamps have already been inserted into
#' \code{df} and corresponding measurements set to \code{NA}, for example with
#' \code{\link{dm.na.interpolation}}.
#'
#' For a focal series \eqn{a}, let \eqn{a_{t-1}} be the last available focal
#' value before the current timestamp. Let \eqn{b_{t-1}} and \eqn{b_t} be the
#' vectors of reference values at the previous and current timestamps.
#'
#' In proportional mode, the relative reference change is
#' \deqn{
#' \delta = \mathrm{mean}\left(\frac{b_t - b_{t-1}}{\max(b_{t-1}, \epsilon)}\right)
#' }
#' and the focal estimate is
#' \deqn{
#' \hat{a}_t = a_{t-1}(1 + \delta).
#' }
#' Bootstrap resampling of reference sensors is used to estimate a 95\%
#' interval.
#'
#' In linear mode, a simple regression of \eqn{b_t} on \eqn{b_{t-1}} is fit
#' across reference sensors at each step, and the focal value is predicted from
#' \eqn{a_{t-1}} with a 95\% prediction interval.
#'
#' If fewer than two valid reference sensors are available at a timestamp pair,
#' or if the focal previous value is missing, the focal value is left unchanged.
#'
#' If \code{correct_gap_jumps = TRUE}, the function inspects the first observed
#' point immediately after each imputed run. If a jump larger than
#' \code{jump_threshold} is detected, the estimated offset is removed from that
#' point onward. This is useful when a dendrometer resumes from a shifted
#' baseline after a gap.
#'
#' If \code{assess = TRUE}, the function inserts artificial gaps into observed
#' sections of each focal series, interpolates them with the chosen method, and
#' compares predictions with the true values. This produces both per-gap
#' summaries and seasonal diagnostics of interpolation error.
#'
#' The returned object is a data frame with class
#' \code{"network_interpolation"}.
#'
#' @param df Data frame containing focal dendrometer data. The first column must
#' be datetime (\code{POSIXct}, \code{Date}, or character convertible by
#' \code{lubridate::ymd_hms}). Remaining columns are focal series.
#' @param referenceDF Data frame containing reference dendrometer data. The
#' first column must be datetime, and remaining columns are reference series.
#' @param niMethod Character. Interpolation method. One of
#' \code{"proportional"} or \code{"linear"}.
#' @param n_boot Integer. Number of bootstrap resamples used in proportional
#' mode.
#' @param return_flags Logical. If \code{TRUE}, add a logical column
#' \code{<series>_interp} indicating which focal values were imputed.
#' @param return_pi Logical. If \code{TRUE}, add 95\% prediction-limit columns
#' \code{<series>_pi_lo} and \code{<series>_pi_hi}.
#' @param pi_for_all Logical. If \code{TRUE}, prediction limits are also
#' attached for non-imputed rows where the model could be formed. If
#' \code{FALSE}, limits are returned only for imputed rows.
#' @param return_fit Logical. If \code{TRUE}, add \code{<series>_fit} columns
#' containing the model-implied fitted value.
#' @param assess Logical. If \code{TRUE}, run synthetic-gap validation.
#' @param assess_lengths_steps Integer vector. Artificial gap lengths,
#' expressed in number of rows/time steps, to test during validation.
#' @param assess_samples_per_length Integer. Number of sampled artificial gaps
#' per gap length and per focal series.
#' @param assess_buffer_steps Integer. Number of observed steps required before
#' and after each artificial gap when
#' \code{assess_use_only_observed = TRUE}.
#' @param assess_seed Integer or \code{NULL}. Random seed for reproducible
#' validation sampling.
#' @param assess_use_only_observed Logical. If \code{TRUE}, synthetic gaps are
#' sampled only from windows fully observed in the focal series, including
#' the requested buffer.
#' @param progress Logical. If \code{TRUE}, display a single-line progress
#' indicator.
#' @param correct_gap_jumps Logical. If \code{TRUE}, inspect the first observed
#' point after each imputed run and optionally remove a detected jump.
#' @param jump_threshold Numeric. Minimum absolute jump size required before a
#' post-gap correction is applied. Required if
#' \code{correct_gap_jumps = TRUE}.
#' @param jump_adjustment_method Character. Method used to estimate the
#' post-gap offset. One of \code{"window_median"} or \code{"point_diff"}.
#' @param jump_adjust_window Integer. Window size used by
#' \code{jump_adjustment_method = "window_median"}.
#' @param jump_set_point_na Logical. If \code{TRUE}, the jump point itself is
#' set to \code{NA} after correction.
#'
#' @return
#' A data frame with class \code{"network_interpolation"}.
#'
#' The data frame contains the original \code{TIME} column, focal series
#' columns, and optionally:
#' \describe{
#' \item{\code{<series>_interp}}{Logical flag for imputed rows.}
#' \item{\code{<series>_pi_lo}, \code{<series>_pi_hi}}{95\% prediction
#' limits.}
#' \item{\code{<series>_fit}}{Model-implied fitted values.}
#' \item{\code{<series>_jump_correction}}{Numeric jump offset removed at the
#' first observed point after an imputed run. \code{NA} means no correction
#' was applied.}
#' }
#'
#' In addition, the following attributes are attached:
#' \describe{
#' \item{\code{network_original}}{Original focal data on the aligned time
#' grid.}
#' \item{\code{network_reference}}{Reference data aligned to the focal time
#' grid.}
#' \item{\code{network_diagnostics}}{Long-format diagnostics for plotting.}
#' \item{\code{network_summary}}{Per-series interpolation summary table.}
#' \item{\code{network_jump_corrections}}{Table of inspected and corrected
#' post-gap jumps.}
#' \item{\code{network_validation_raw}}{Per-gap validation summary table, or
#' \code{NULL}.}
#' \item{\code{network_validation_points}}{Point-level validation results, or
#' \code{NULL}.}
#' \item{\code{network_validation_summary}}{Summarized validation metrics by
#' series and gap length, or \code{NULL}.}
#' \item{\code{network_validation_seasonal}}{Seasonal validation diagnostics
#' by day of year, or \code{NULL}.}
#' \item{\code{network_method}}{Interpolation method used.}
#' \item{\code{network_n_boot}}{Number of bootstrap resamples used.}
#' \item{\code{network_assessed}}{Logical indicating whether validation was
#' run.}
#' \item{\code{network_jump_settings}}{List of jump-correction settings.}
#' }
#'
#' @section Interpretation:
#' Smaller validation metrics such as \code{MAE}, \code{RMSE}, \code{MAPE}, and
#' \code{MdAPE} indicate better gap recovery. \code{PI_coverage_95} indicates
#' how well the nominal 95\% uncertainty interval captures the true values
#' under synthetic-gap validation. Seasonal diagnostics help identify periods
#' of the year when network interpolation is more or less reliable.
#'
#' If many post-gap corrections are applied or \code{max_abs_gap_jump} is
#' large, it may indicate re-zeroing, sensor reset, or baseline
#' discontinuities after gaps.
#'
#' @section Notes:
#' \itemize{
#' \item The focal dataset \code{df} defines the output time grid.
#' \item Reference series are aligned to that grid using the datetime column.
#' \item At least two valid reference sensors are required at consecutive
#' timestamps for interpolation.
#' \item Jump correction is a post-processing step and does not alter the
#' core interpolation logic.
#' }
#'
#' @examples
#' \donttest{
#' #library(dendRoAnalyst)
#' #data("gf_nepa17")
#'
#' ## Create an artificial focal gap
#' #df1 <- gf_nepa17
#' #df1[40:50, "T2"] <- NA
#'
#' ## Build a small reference network
#' #ref <- cbind(gf_nepa17, gf_nepa17[, 2:3], gf_nepa17[, 2:3])
#' #colnames(ref) <- c("Time", "T1", "T2", "T3", "T4", "T5", "T6")
#'
#' ## Interpolate with uncertainty limits
#' #out <- network.interpolation(
#' # df1, ref,
#' # niMethod = "proportional",
#' # n_boot = 100,
#' # return_flags = TRUE,
#' # return_pi = TRUE
#' #)
#'
#' #head(out, 10)
#'
#' # Interpolate with jump correction and validation
#' #out2 <- network.interpolation(
#' # df1, ref,
#' # niMethod = "proportional",
#' # n_boot = 100,
#' # return_flags = TRUE,
#' # return_pi = TRUE,
#' # assess = TRUE,
#' # assess_lengths_steps = c(1, 2, 4),
#' # correct_gap_jumps = TRUE,
#' # jump_threshold = 0.05,
#' # jump_adjustment_method = "window_median",
#' # jump_adjust_window = 5
#' #)
#' # Plotting
#' #plot(out2)
#' #plot(out2, type = "compare")
#' #plot(out2, type = "seasonal_error")
#'
#' # Extracting the information using attr()
#' #attr(out, "network_validation_summary")
#' #attr(out, "network_validation_points")
#' #attr(out, "network_validation_seasonal")
#' }
#'
#' @seealso
#' \code{\link{plot.network_interpolation}},
#' \code{\link{dm.na.interpolation}}
#'
#' @export
network.interpolation <- function(df, referenceDF,
niMethod = c("proportional", "linear"),
n_boot = 1000,
return_flags = TRUE,
return_pi = TRUE,
pi_for_all = FALSE,
return_fit = FALSE,
assess = FALSE,
assess_lengths_steps = c(1, 2, 3, 6, 12, 24),
assess_samples_per_length = 20,
assess_buffer_steps = 1,
assess_seed = NULL,
assess_use_only_observed = TRUE,
progress = interactive(),
correct_gap_jumps = FALSE,
jump_threshold = NULL,
jump_adjustment_method = c("window_median", "point_diff"),
jump_adjust_window = 5,
jump_set_point_na = FALSE) {
niMethod <- match.arg(niMethod)
jump_adjustment_method <- match.arg(jump_adjustment_method)
if (isTRUE(correct_gap_jumps) &&
(is.null(jump_threshold) || !is.finite(jump_threshold) || jump_threshold <= 0)) {
stop("When correct_gap_jumps = TRUE, provide a positive numeric 'jump_threshold'.")
}
Bias <- End_value_abs_diff <- MAE <- MAPE <- Max_abs_diff_within_gap <- MdAPE <- Mean_PI_width <- Mean_ref_n <- PI_coverage_95 <- RMSE <- Success_rate <- TIME <- abs_pct_error <- applied_offset <- case_id <- covered <- doy <- error <- filled <- gap_steps <- jump_correction <- jump_removed <- jump_time <- max_abs_gap_jump <- mean_abs_error <- mean_abs_pct_error <- mean_error <- mean_pred <- mean_truth <- n_gap_jumps_removed <- pred <- raw_jump <- ref_n <- rmse <- series <- start_time <- end_time <- truth <- NULL
to_posixct <- function(x) {
if (inherits(x, "POSIXct")) return(x)
if (inherits(x, "Date")) return(as.POSIXct(x))
lubridate::ymd_hms(x)
}
mean_or_na <- function(x) {
if (!length(x) || all(is.na(x))) return(NA_real_)
mean(x, na.rm = TRUE)
}
median_or_na <- function(x) {
if (!length(x) || all(is.na(x))) return(NA_real_)
stats::median(x, na.rm = TRUE)
}
boot_mean_delta_fast <- function(delta, n_boot) {
delta <- delta[is.finite(delta)]
n <- length(delta)
if (n < 2L) {
return(list(
mean = if (n) mean(delta) else NA_real_,
lo = NA_real_,
hi = NA_real_
))
}
idx <- matrix(sample.int(n, size = n * n_boot, replace = TRUE), nrow = n)
means <- colMeans(matrix(delta[idx], nrow = n))
list(
mean = mean(means, na.rm = TRUE),
lo = stats::quantile(means, 0.025, na.rm = TRUE, names = FALSE),
hi = stats::quantile(means, 0.975, na.rm = TRUE, names = FALSE)
)
}
fast_linear_pred <- function(x, y, x0, level = 0.95) {
ok <- is.finite(x) & is.finite(y)
x <- x[ok]
y <- y[ok]
n <- length(x)
if (n < 2L || !is.finite(x0)) {
return(c(fit = NA_real_, lwr = NA_real_, upr = NA_real_))
}
xbar <- mean(x)
ybar <- mean(y)
sxx <- sum((x - xbar)^2)
if (!is.finite(sxx) || sxx <= 0) {
return(c(fit = NA_real_, lwr = NA_real_, upr = NA_real_))
}
beta <- sum((x - xbar) * (y - ybar)) / sxx
alpha <- ybar - beta * xbar
fit <- alpha + beta * x0
if (n <= 2L) {
return(c(fit = fit, lwr = NA_real_, upr = NA_real_))
}
resid <- y - (alpha + beta * x)
mse <- sum(resid^2) / (n - 2L)
if (!is.finite(mse)) {
return(c(fit = fit, lwr = NA_real_, upr = NA_real_))
}
se_pred <- sqrt(mse * (1 + 1 / n + ((x0 - xbar)^2 / sxx)))
tval <- stats::qt((1 + level) / 2, df = n - 2L)
c(
fit = fit,
lwr = fit - tval * se_pred,
upr = fit + tval * se_pred
)
}
estimate_jump_offset_gap <- function(x, j,
adjustment_method = c("window_median", "point_diff"),
adjust_window = 5) {
adjustment_method <- match.arg(adjustment_method)
if (length(x) == 0 || j <= 1 || j > length(x)) {
return(NA_real_)
}
if (adjustment_method == "point_diff") {
if (is.na(x[j]) || is.na(x[j - 1])) {
return(NA_real_)
}
return(x[j] - x[j - 1])
}
n <- length(x)
pre_idx <- seq.int(max(1, j - adjust_window), j - 1)
post_idx <- seq.int(j, min(n, j + adjust_window - 1))
pre_vals <- x[pre_idx]
post_vals <- x[post_idx]
pre_vals <- pre_vals[!is.na(pre_vals)]
post_vals <- post_vals[!is.na(post_vals)]
if (length(pre_vals) < 2 || length(post_vals) < 2) {
if (is.na(x[j]) || is.na(x[j - 1])) {
return(NA_real_)
}
return(x[j] - x[j - 1])
}
stats::median(post_vals) - stats::median(pre_vals)
}
apply_jump_correction_gap <- function(x, j, offset, set_jump_point_na = FALSE) {
if (length(x) == 0 || is.na(offset) || j > length(x)) {
return(x)
}
idx <- j:length(x)
idx <- idx[!is.na(x[idx])]
if (length(idx) > 0) {
x[idx] <- x[idx] - offset
}
if (isTRUE(set_jump_point_na) && j <= length(x)) {
x[j] <- NA_real_
}
x
}
correct_gap_end_jumps <- function(value, fit, pi_lo, pi_hi,
interp_flag, original_values, time_vec,
jump_threshold,
adjustment_method = c("window_median", "point_diff"),
adjust_window = 5,
set_jump_point_na = FALSE) {
adjustment_method <- match.arg(adjustment_method)
n <- length(value)
jump_correction <- rep(NA_real_, n)
if (!any(interp_flag, na.rm = TRUE)) {
return(list(
value = value,
fit = fit,
pi_lo = pi_lo,
pi_hi = pi_hi,
jump_correction = jump_correction,
table = tibble::tibble()
))
}
r <- rle(interp_flag)
ends <- cumsum(r$lengths)
gap_end_idx <- ends[r$values]
candidate_idx <- gap_end_idx + 1L
candidate_idx <- candidate_idx[candidate_idx <= n]
if (!length(candidate_idx)) {
return(list(
value = value,
fit = fit,
pi_lo = pi_lo,
pi_hi = pi_hi,
jump_correction = jump_correction,
table = tibble::tibble()
))
}
rows <- vector("list", length(candidate_idx))
rr <- 0L
for (j in candidate_idx) {
if (!is.finite(original_values[j])) next
raw_jump <- if (j > 1 && is.finite(value[j]) && is.finite(value[j - 1])) {
value[j] - value[j - 1]
} else {
NA_real_
}
offset <- estimate_jump_offset_gap(
x = value,
j = j,
adjustment_method = adjustment_method,
adjust_window = adjust_window
)
removed <- is.finite(offset) && abs(offset) >= jump_threshold
if (removed) {
value <- apply_jump_correction_gap(value, j, offset, set_jump_point_na = set_jump_point_na)
fit <- apply_jump_correction_gap(fit, j, offset, set_jump_point_na = FALSE)
pi_lo <- apply_jump_correction_gap(pi_lo, j, offset, set_jump_point_na = FALSE)
pi_hi <- apply_jump_correction_gap(pi_hi, j, offset, set_jump_point_na = FALSE)
jump_correction[j] <- offset
}
rr <- rr + 1L
rows[[rr]] <- tibble::tibble(
jump_index = j,
jump_time = time_vec[j],
raw_jump = raw_jump,
estimated_offset = offset,
applied_offset = if (removed) offset else NA_real_,
jump_removed = removed
)
}
jump_table <- if (rr > 0L) dplyr::bind_rows(rows[seq_len(rr)]) else tibble::tibble()
list(
value = value,
fit = fit,
pi_lo = pi_lo,
pi_hi = pi_hi,
jump_correction = jump_correction,
table = jump_table
)
}
interpolate_one_series <- function(a_raw, ref_mat, ref_ok, niMethod, n_boot, tiny = 1e-6) {
n <- length(a_raw)
a <- as.numeric(a_raw)
interp_flag <- rep(FALSE, n)
ref_n <- rep(NA_integer_, n)
fit <- rep(NA_real_, n)
pi_lo <- rep(NA_real_, n)
pi_hi <- rep(NA_real_, n)
for (i in 2:n) {
a_prev <- a[i - 1L]
if (!is.finite(a_prev)) next
ok <- ref_ok[i - 1L, ] & ref_ok[i, ]
n_ok <- sum(ok)
ref_n[i] <- n_ok
if (n_ok < 2L) next
b_prev <- ref_mat[i - 1L, ok]
b_cur <- ref_mat[i, ok]
if (niMethod == "proportional") {
rel <- (b_cur - b_prev) / pmax(b_prev, tiny)
bm <- boot_mean_delta_fast(rel, n_boot)
if (!is.finite(bm$mean)) next
fit_i <- a_prev * (1 + bm$mean)
lo_i <- if (is.na(bm$lo)) NA_real_ else a_prev * (1 + bm$lo)
hi_i <- if (is.na(bm$hi)) NA_real_ else a_prev * (1 + bm$hi)
} else {
pred <- fast_linear_pred(b_prev, b_cur, a_prev, level = 0.95)
fit_i <- unname(pred["fit"])
lo_i <- unname(pred["lwr"])
hi_i <- unname(pred["upr"])
if (!is.finite(fit_i)) next
}
fit[i] <- fit_i
pi_lo[i] <- lo_i
pi_hi[i] <- hi_i
if (is.na(a[i]) && is.finite(fit_i)) {
a[i] <- fit_i
interp_flag[i] <- TRUE
}
}
list(
value = a,
interp_flag = interp_flag,
ref_n = ref_n,
fit = fit,
pi_lo = pi_lo,
pi_hi = pi_hi
)
}
df <- tibble::as_tibble(df)
referenceDF <- tibble::as_tibble(referenceDF)
names(df)[1] <- "TIME"
names(referenceDF)[1] <- "TIME"
df$TIME <- to_posixct(df$TIME)
referenceDF$TIME <- to_posixct(referenceDF$TIME)
if (any(is.na(df$TIME))) {
stop("Could not convert the first column of 'df' to POSIXct.")
}
if (any(is.na(referenceDF$TIME))) {
stop("Could not convert the first column of 'referenceDF' to POSIXct.")
}
df <- df[order(df$TIME), , drop = FALSE]
referenceDF <- referenceDF[order(referenceDF$TIME), , drop = FALSE]
time_min <- min(df$TIME, na.rm = TRUE)
time_max <- max(df$TIME, na.rm = TRUE)
referenceDF <- referenceDF[
referenceDF$TIME >= time_min & referenceDF$TIME <= time_max,
,
drop = FALSE
]
if (nrow(referenceDF) == 0) {
stop("'referenceDF' has no overlap with the time span of 'df'.")
}
referenceDF <- dplyr::left_join(df["TIME"], referenceDF, by = "TIME")
if (ncol(referenceDF) < 2L) {
stop("'referenceDF' must contain at least one reference series.")
}
ref_mat <- as.matrix(referenceDF[, -1, drop = FALSE])
suppressWarnings(storage.mode(ref_mat) <- "double")
ref_ok <- is.finite(ref_mat)
if (!any(ref_ok)) {
stop("No usable reference values remain after aligning 'referenceDF' to the time grid of 'df'.")
}
original_df <- df
out <- df
add_cols <- list()
diag_list <- vector("list", max(1L, ncol(df) - 1L))
sum_list <- vector("list", max(1L, ncol(df) - 1L))
jump_tables <- vector("list", max(1L, ncol(df) - 1L))
n_series <- max(1L, ncol(df) - 1L)
progress_total <- n_series + if (assess) n_series * length(assess_lengths_steps) else 0L
progress_value <- 0L
.progress_last_width <- 0L
render_progress_line <- function(note = "") {
if (!isTRUE(progress) || progress_total <= 0L) return(invisible(NULL))
width_bar <- 30L
frac <- progress_value / progress_total
n_done <- floor(width_bar * frac)
bar <- paste0(
"[",
paste(rep("=", n_done), collapse = ""),
paste(rep("-", width_bar - n_done), collapse = ""),
"]"
)
txt <- sprintf(
"\r%s %3d%% (%d/%d) %s",
bar,
round(100 * frac),
progress_value,
progress_total,
note
)
pad <- max(0L, .progress_last_width - nchar(txt, type = "width"))
cat(txt, paste(rep(" ", pad), collapse = ""))
utils::flush.console()
.progress_last_width <<- nchar(txt, type = "width")
invisible(NULL)
}
update_progress <- function(step = 1L, note = "") {
progress_value <<- min(progress_total, progress_value + step)
render_progress_line(note)
invisible(NULL)
}
progress_note <- function(note = "") {
render_progress_line(note)
invisible(NULL)
}
finish_progress <- function() {
if (!isTRUE(progress)) return(invisible(NULL))
cat("\n")
utils::flush.console()
invisible(NULL)
}
if (isTRUE(progress) && progress_total > 0L) {
on.exit(finish_progress(), add = TRUE)
render_progress_line("Starting...")
}
for (j in seq.int(2, ncol(df))) {
colname <- names(df)[j]
progress_note(paste0("Interpolating series ", colname, "..."))
raw_a <- df[[j]]
res <- interpolate_one_series(
a_raw = raw_a,
ref_mat = ref_mat,
ref_ok = ref_ok,
niMethod = niMethod,
n_boot = n_boot
)
jump_correction <- rep(NA_real_, length(raw_a))
jump_tbl <- tibble::tibble()
if (isTRUE(correct_gap_jumps)) {
jc <- correct_gap_end_jumps(
value = res$value,
fit = res$fit,
pi_lo = res$pi_lo,
pi_hi = res$pi_hi,
interp_flag = res$interp_flag,
original_values = raw_a,
time_vec = df$TIME,
jump_threshold = jump_threshold,
adjustment_method = jump_adjustment_method,
adjust_window = jump_adjust_window,
set_jump_point_na = jump_set_point_na
)
res$value <- jc$value
res$fit <- jc$fit
res$pi_lo <- jc$pi_lo
res$pi_hi <- jc$pi_hi
jump_correction <- jc$jump_correction
if (nrow(jc$table) > 0) {
jump_tbl <- jc$table
jump_tbl$series <- colname
}
}
jump_tables[[j - 1L]] <- jump_tbl
out[[j]] <- res$value
if (return_flags) {
add_cols[[paste0(colname, "_interp")]] <- res$interp_flag
}
if (return_pi) {
out_pi_lo <- ifelse(res$interp_flag | pi_for_all, res$pi_lo, NA_real_)
out_pi_hi <- ifelse(res$interp_flag | pi_for_all, res$pi_hi, NA_real_)
add_cols[[paste0(colname, "_pi_lo")]] <- out_pi_lo
add_cols[[paste0(colname, "_pi_hi")]] <- out_pi_hi
}
if (return_fit) {
out_fit <- ifelse(res$interp_flag | pi_for_all, res$fit, NA_real_)
add_cols[[paste0(colname, "_fit")]] <- out_fit
}
if (isTRUE(correct_gap_jumps)) {
add_cols[[paste0(colname, "_jump_correction")]] <- jump_correction
}
diag_list[[j - 1L]] <- tibble::tibble(
TIME = df$TIME,
series = colname,
original = raw_a,
value = res$value,
observed = is.finite(raw_a),
interp = res$interp_flag,
remaining_na = is.na(res$value),
ref_n = res$ref_n,
fit = res$fit,
pi_lo = res$pi_lo,
pi_hi = res$pi_hi,
pi_width = res$pi_hi - res$pi_lo,
jump_correction = jump_correction
)
sum_list[[j - 1L]] <- tibble::tibble(
series = colname,
n_rows = length(raw_a),
n_missing_input = sum(is.na(raw_a)),
n_imputed = sum(res$interp_flag),
n_remaining_na = sum(is.na(res$value)),
mean_ref_n = mean_or_na(res$ref_n),
median_ref_n = median_or_na(res$ref_n),
mean_pi_width = mean_or_na((res$pi_hi - res$pi_lo)[res$interp_flag]),
median_pi_width = median_or_na((res$pi_hi - res$pi_lo)[res$interp_flag]),
n_gap_jumps_removed = sum(!is.na(jump_correction)),
max_abs_gap_jump = if (all(is.na(jump_correction))) NA_real_ else max(abs(jump_correction), na.rm = TRUE)
)
update_progress(note = paste0("Finished series ", colname))
}
if (length(add_cols)) {
out <- dplyr::bind_cols(out, tibble::as_tibble(add_cols))
}
diagnostics <- dplyr::bind_rows(diag_list)
summary_tbl <- dplyr::bind_rows(sum_list)
jump_corrections <- dplyr::bind_rows(jump_tables)
validation_raw <- NULL
validation_points <- NULL
validation_summary <- NULL
validation_seasonal <- NULL
if (assess) {
if (!is.null(assess_seed)) {
set.seed(assess_seed)
}
max_cases <- n_series * length(assess_lengths_steps) * assess_samples_per_length
max_points <- max_cases * max(assess_lengths_steps)
val_rows <- vector("list", max(1L, max_cases))
point_rows <- vector("list", max(1L, max_points))
val_id <- 0L
point_id <- 0L
push_val_row <- function(row) {
val_id <<- val_id + 1L
if (val_id > length(val_rows)) {
length(val_rows) <<- max(length(val_rows) * 2L, val_id)
}
val_rows[[val_id]] <<- row
}
push_point_row <- function(row) {
point_id <<- point_id + 1L
if (point_id > length(point_rows)) {
length(point_rows) <<- max(length(point_rows) * 2L, point_id)
}
point_rows[[point_id]] <<- row
}
for (j in seq.int(2, ncol(df))) {
colname <- names(df)[j]
raw_a <- as.numeric(df[[j]])
n <- length(raw_a)
if (n < 2L) {
update_progress(length(assess_lengths_steps), note = paste0("Skipping validation for ", colname))
next
}
finite_raw <- is.finite(raw_a)
cs_finite <- c(0L, cumsum(finite_raw))
for (gap_steps in as.integer(assess_lengths_steps)) {
progress_note(paste0("Validating ", colname, " | gap=", gap_steps, " step(s)"))
if (!is.finite(gap_steps) || gap_steps < 1L) {
update_progress(note = paste0("Skipped ", colname, " | gap=", gap_steps))
next
}
min_start <- 2L + assess_buffer_steps
max_start <- n - gap_steps - assess_buffer_steps + 1L
if (max_start >= min_start) {
starts <- seq.int(min_start, max_start)
gap_counts <- cs_finite[starts + gap_steps] - cs_finite[starts]
prev_ok <- finite_raw[starts - 1L]
valid <- (gap_counts == gap_steps) & prev_ok
if (assess_use_only_observed) {
win_left <- starts - assess_buffer_steps
win_right <- starts + gap_steps - 1L + assess_buffer_steps
win_counts <- cs_finite[win_right + 1L] - cs_finite[win_left]
valid <- valid & (win_counts == (gap_steps + 2L * assess_buffer_steps))
}
valid_starts <- starts[valid]
if (length(valid_starts)) {
starts_to_use <- if (length(valid_starts) > assess_samples_per_length) {
sample(valid_starts, assess_samples_per_length, replace = FALSE)
} else {
valid_starts
}
for (s in starts_to_use) {
idx_gap <- s:(s + gap_steps - 1L)
case_id <- paste(colname, gap_steps, s, sep = "_")
a_mask <- raw_a
a_mask[idx_gap] <- NA_real_
res_case <- interpolate_one_series(
a_raw = a_mask,
ref_mat = ref_mat,
ref_ok = ref_ok,
niMethod = niMethod,
n_boot = n_boot
)
if (isTRUE(correct_gap_jumps)) {
jc_case <- correct_gap_end_jumps(
value = res_case$value,
fit = res_case$fit,
pi_lo = res_case$pi_lo,
pi_hi = res_case$pi_hi,
interp_flag = res_case$interp_flag,
original_values = a_mask,
time_vec = df$TIME,
jump_threshold = jump_threshold,
adjustment_method = jump_adjustment_method,
adjust_window = jump_adjust_window,
set_jump_point_na = jump_set_point_na
)
res_case$value <- jc_case$value
res_case$fit <- jc_case$fit
res_case$pi_lo <- jc_case$pi_lo
res_case$pi_hi <- jc_case$pi_hi
}
truth <- raw_a[idx_gap]
pred <- res_case$value[idx_gap]
pi_lo_gap <- res_case$pi_lo[idx_gap]
pi_hi_gap <- res_case$pi_hi[idx_gap]
ref_n_gap <- res_case$ref_n[idx_gap]
filled <- res_case$interp_flag[idx_gap] & is.finite(pred)
n_points <- length(idx_gap)
n_filled <- sum(filled)
success_rate <- n_filled / n_points
for (k in seq_along(idx_gap)) {
err_k <- if (is.finite(pred[k]) && is.finite(truth[k])) pred[k] - truth[k] else NA_real_
ape_k <- if (is.finite(err_k) && is.finite(truth[k])) {
100 * abs(err_k) / pmax(abs(truth[k]), 1e-6)
} else {
NA_real_
}
cov_k <- if (is.finite(truth[k]) && is.finite(pi_lo_gap[k]) && is.finite(pi_hi_gap[k])) {
truth[k] >= pi_lo_gap[k] && truth[k] <= pi_hi_gap[k]
} else {
NA
}
push_point_row(tibble::tibble(
series = colname,
gap_steps = gap_steps,
case_id = case_id,
TIME = df$TIME[idx_gap[k]],
truth = truth[k],
pred = pred[k],
filled = filled[k],
error = err_k,
abs_pct_error = ape_k,
pi_lo = pi_lo_gap[k],
pi_hi = pi_hi_gap[k],
covered = cov_k,
ref_n = ref_n_gap[k]
))
}
if (n_filled > 0L) {
err <- pred[filled] - truth[filled]
denom <- pmax(abs(truth[filled]), 1e-6)
pe <- 100 * err / denom
pi_ok <- filled & is.finite(pi_lo_gap) & is.finite(pi_hi_gap)
PI_coverage_95 <- if (any(pi_ok)) {
mean(truth[pi_ok] >= pi_lo_gap[pi_ok] & truth[pi_ok] <= pi_hi_gap[pi_ok], na.rm = TRUE)
} else {
NA_real_
}
Mean_PI_width <- if (any(pi_ok)) {
mean(pi_hi_gap[pi_ok] - pi_lo_gap[pi_ok], na.rm = TRUE)
} else {
NA_real_
}
End_value_abs_diff <- if (is.finite(pred[length(pred)])) {
abs(pred[length(pred)] - truth[length(truth)])
} else {
NA_real_
}
Max_abs_diff_within_gap <- max(abs(err), na.rm = TRUE)
row <- tibble::tibble(
series = colname,
gap_steps = gap_steps,
case_id = case_id,
start_time = df$TIME[min(idx_gap)],
end_time = df$TIME[max(idx_gap)],
n_points = n_points,
n_filled = n_filled,
Success_rate = success_rate,
MAE = mean(abs(err), na.rm = TRUE),
MAPE = mean(abs(pe), na.rm = TRUE),
MdAPE = stats::median(abs(pe), na.rm = TRUE),
RMSE = sqrt(mean(err^2, na.rm = TRUE)),
Bias = mean(err, na.rm = TRUE),
PI_coverage_95 = PI_coverage_95,
Mean_PI_width = Mean_PI_width,
Mean_ref_n = mean(ref_n_gap, na.rm = TRUE),
End_value_abs_diff = End_value_abs_diff,
Max_abs_diff_within_gap = Max_abs_diff_within_gap
)
} else {
row <- tibble::tibble(
series = colname,
gap_steps = gap_steps,
case_id = case_id,
start_time = df$TIME[min(idx_gap)],
end_time = df$TIME[max(idx_gap)],
n_points = n_points,
n_filled = n_filled,
Success_rate = success_rate,
MAE = NA_real_,
MAPE = NA_real_,
MdAPE = NA_real_,
RMSE = NA_real_,
Bias = NA_real_,
PI_coverage_95 = NA_real_,
Mean_PI_width = NA_real_,
Mean_ref_n = mean(ref_n_gap, na.rm = TRUE),
End_value_abs_diff = NA_real_,
Max_abs_diff_within_gap = NA_real_
)
}
push_val_row(row)
}
}
}
update_progress(note = paste0("Validated ", colname, " | gap=", gap_steps))
}
}
if (val_id > 0L) {
validation_raw <- dplyr::bind_rows(val_rows[seq_len(val_id)])
validation_points <- dplyr::bind_rows(point_rows[seq_len(point_id)])
validation_summary <- validation_raw %>%
dplyr::group_by(series, gap_steps) %>%
dplyr::summarise(
n_cases = dplyr::n(),
n_points = sum(n_points, na.rm = TRUE),
n_filled = sum(n_filled, na.rm = TRUE),
Success_rate = mean_or_na(Success_rate),
MAE = mean_or_na(MAE),
MAPE = mean_or_na(MAPE),
MdAPE = mean_or_na(MdAPE),
RMSE = mean_or_na(RMSE),
Bias = mean_or_na(Bias),
PI_coverage_95 = mean_or_na(PI_coverage_95),
Mean_PI_width = mean_or_na(Mean_PI_width),
Mean_ref_n = mean_or_na(Mean_ref_n),
End_value_abs_diff = mean_or_na(End_value_abs_diff),
Max_abs_diff_within_gap = mean_or_na(Max_abs_diff_within_gap),
.groups = "drop"
) %>%
dplyr::arrange(series, gap_steps)
validation_seasonal <- validation_points %>%
dplyr::filter(filled, is.finite(truth), is.finite(pred)) %>%
dplyr::mutate(doy = as.integer(format(TIME, "%j"))) %>%
dplyr::group_by(series, doy) %>%
dplyr::summarise(
n_points = dplyr::n(),
mean_error = mean(error, na.rm = TRUE),
mean_abs_error = mean(abs(error), na.rm = TRUE),
mean_abs_pct_error = mean(abs_pct_error, na.rm = TRUE),
rmse = sqrt(mean(error^2, na.rm = TRUE)),
mean_truth = mean(truth, na.rm = TRUE),
mean_pred = mean(pred, na.rm = TRUE),
coverage = mean(covered, na.rm = TRUE),
.groups = "drop"
) %>%
dplyr::arrange(series, doy)
} else {
validation_raw <- tibble::tibble()
validation_points <- tibble::tibble()
validation_summary <- tibble::tibble()
validation_seasonal <- tibble::tibble()
warning("Assessment produced no results (no valid synthetic windows found).")
}
}
attr(out, "network_original") <- original_df
attr(out, "network_reference") <- referenceDF
attr(out, "network_diagnostics") <- diagnostics
attr(out, "network_summary") <- summary_tbl
attr(out, "network_jump_corrections") <- jump_corrections
attr(out, "network_validation_raw") <- validation_raw
attr(out, "network_validation_points") <- validation_points
attr(out, "network_validation_summary") <- validation_summary
attr(out, "network_validation_seasonal") <- validation_seasonal
attr(out, "network_method") <- niMethod
attr(out, "network_n_boot") <- n_boot
attr(out, "network_return_flags") <- return_flags
attr(out, "network_return_pi") <- return_pi
attr(out, "network_pi_for_all") <- pi_for_all
attr(out, "network_return_fit") <- return_fit
attr(out, "network_assessed") <- assess
attr(out, "network_jump_settings") <- list(
correct_gap_jumps = correct_gap_jumps,
jump_threshold = jump_threshold,
jump_adjustment_method = jump_adjustment_method,
jump_adjust_window = jump_adjust_window,
jump_set_point_na = jump_set_point_na
)
class(out) <- c("network_interpolation", class(out))
out
}
#' @title Plot method for network interpolation output
#'
#' @description
#' S3 plot method for objects returned by
#' \code{\link{network.interpolation}}.
#'
#' Depending on \code{type}, the method can display:
#' \itemize{
#' \item interpolated series through time,
#' \item uncertainty width,
#' \item reference-network availability,
#' \item summary diagnostics,
#' \item synthetic-gap validation metrics,
#' \item interval coverage,
#' \item actual versus interpolated validation points,
#' \item seasonal validation diagnostics.
#' }
#'
#' Jump-corrected post-gap points are marked in the interpolation plot when
#' jump correction was enabled in \code{\link{network.interpolation}}.
#'
#' @param x Object returned by \code{\link{network.interpolation}}.
#' @param type Character. Plot type. One of \code{"interpolation"},
#' \code{"uncertainty"}, \code{"availability"}, \code{"summary"},
#' \code{"validation"}, \code{"coverage"}, \code{"compare"}, or
#' \code{"seasonal_error"}.
#' @param series Optional character vector of focal series names to plot.
#' @param start Optional start datetime for time-based plots.
#' @param end Optional end datetime for time-based plots.
#' @param show_pi Logical. If \code{TRUE}, show 95\% prediction intervals in
#' the interpolation plot where available.
#' @param show_fit Logical. If \code{TRUE}, show model-implied fitted values in
#' the interpolation plot where available.
#' @param free_y Logical. If \code{TRUE}, facets use free y scales where
#' applicable.
#' @param ncol Integer. Number of facet columns.
#' @param summary_metric Character metric used when
#' \code{type = "summary"}. One of \code{"n_imputed"},
#' \code{"n_remaining_na"}, \code{"n_missing_input"},
#' \code{"mean_pi_width"}, \code{"median_pi_width"}, \code{"mean_ref_n"},
#' \code{"median_ref_n"}, \code{"n_gap_jumps_removed"}, or
#' \code{"max_abs_gap_jump"}.
#' @param validation_metric Character metric used when
#' \code{type = "validation"}. One of \code{"MAE"}, \code{"MAPE"},
#' \code{"MdAPE"}, \code{"RMSE"}, \code{"Bias"},
#' \code{"Success_rate"}, \code{"Mean_PI_width"}, \code{"Mean_ref_n"},
#' \code{"End_value_abs_diff"}, or \code{"Max_abs_diff_within_gap"}.
#' @param seasonal_metric Character metric used when
#' \code{type = "seasonal_error"}. One of \code{"mean_error"},
#' \code{"mean_abs_error"}, \code{"mean_abs_pct_error"}, \code{"rmse"},
#' or \code{"coverage"}.
#' @param compare_colour Character. Colouring used when
#' \code{type = "compare"}. One of \code{"series"} or
#' \code{"gap_steps"}.
#' @param facet Logical. If \code{TRUE}, use faceting where supported.
#' @param empty Character. Behavior when no suitable data are available:
#' \code{"plot"} returns an informative empty plot and \code{"error"} throws
#' an error.
#' @param ... Further arguments passed through the generic.
#'
#' @return A \code{ggplot2} object.
#'
#' @details
#' Plot types:
#' \describe{
#' \item{\code{"interpolation"}}{Shows original and interpolated focal series
#' through time. Imputed points are highlighted. If jump correction was
#' applied, corrected post-gap points are marked in green.}
#' \item{\code{"uncertainty"}}{Shows the width of the prediction interval
#' through time.}
#' \item{\code{"availability"}}{Shows the number of valid reference sensors
#' available at each step.}
#' \item{\code{"summary"}}{Shows per-series summary diagnostics produced by
#' the interpolation run.}
#' \item{\code{"validation"}}{Shows synthetic-gap recovery metrics across
#' tested gap lengths.}
#' \item{\code{"coverage"}}{Shows how often true values fall inside the
#' nominal 95\% prediction interval across gap lengths.}
#' \item{\code{"compare"}}{Scatter plot of actual versus interpolated values
#' from synthetic-gap validation.}
#' \item{\code{"seasonal_error"}}{Shows how validation error varies through
#' the season, summarized by day of year.}
#' }
#'
#' Validation-based plots require that
#' \code{network.interpolation(..., assess = TRUE)} was used to create
#' \code{x}.
#'
#' @examples
#' \donttest{
#' #library(dendRoAnalyst)
#' #data("gf_nepa17")
#'
#' #df1 <- gf_nepa17
#' #df1[40:50, "T2"] <- NA
#'
#' #ref <- cbind(gf_nepa17, gf_nepa17[, 2:3], gf_nepa17[, 2:3])
#' #colnames(ref) <- c("Time", "T1", "T2", "T3", "T4", "T5", "T6")
#'
#' #out <- network.interpolation(
#' # df1, ref,
#' # niMethod = "proportional",
#' # n_boot = 100,
#' # assess = TRUE,
#' # assess_lengths_steps = c(1, 2, 4),
#' # correct_gap_jumps = TRUE,
#' # jump_threshold = 0.05
#' #)
#'
#' #plot(out)
#' #plot(out, type = "uncertainty")
#' #plot(out, type = "availability")
#' #plot(out, type = "summary", summary_metric = "n_gap_jumps_removed")
#' #plot(out, type = "validation", validation_metric = "MAPE")
#' #plot(out, type = "coverage")
#' #plot(out, type = "compare")
#' #plot(out, type = "seasonal_error", seasonal_metric = "mean_abs_error")
#' }
#'
#' @seealso
#' \code{\link{network.interpolation}}
#'
#' @method plot network_interpolation
#' @export
plot.network_interpolation <- function(x,
type = c("interpolation", "uncertainty", "availability", "summary", "validation", "coverage", "compare", "seasonal_error"),
series = NULL,
start = NULL,
end = NULL,
show_pi = TRUE,
show_fit = FALSE,
free_y = TRUE,
ncol = 1,
summary_metric = c("n_imputed", "n_remaining_na", "n_missing_input", "mean_pi_width", "median_pi_width", "mean_ref_n", "median_ref_n", "n_gap_jumps_removed", "max_abs_gap_jump"),
validation_metric = c("MAE", "MAPE", "MdAPE", "RMSE", "Bias", "Success_rate", "Mean_PI_width", "Mean_ref_n", "End_value_abs_diff", "Max_abs_diff_within_gap"),
seasonal_metric = c("mean_error", "mean_abs_error", "mean_abs_pct_error", "rmse", "coverage"),
compare_colour = c("series", "gap_steps"),
facet = TRUE,
empty = c("plot", "error"),
...) {
type <- match.arg(type)
summary_metric <- match.arg(summary_metric)
validation_metric <- match.arg(validation_metric)
seasonal_metric <- match.arg(seasonal_metric)
compare_colour <- match.arg(compare_colour)
empty <- match.arg(empty)
if (!inherits(x, "network_interpolation")) {
stop("x must be an object returned by network.interpolation().")
}
if (type == "interpolation") {
return(.plot_network_interpolation(
x = x,
series = series,
start = start,
end = end,
show_pi = show_pi,
show_fit = show_fit,
free_y = free_y,
ncol = ncol,
empty = empty
))
}
if (type == "uncertainty") {
return(.plot_network_uncertainty(
x = x,
series = series,
start = start,
end = end,
free_y = free_y,
ncol = ncol,
empty = empty
))
}
if (type == "availability") {
return(.plot_network_availability(
x = x,
series = series,
start = start,
end = end,
free_y = free_y,
ncol = ncol,
empty = empty
))
}
if (type == "summary") {
return(.plot_network_summary(
x = x,
metric = summary_metric,
empty = empty
))
}
if (type == "validation") {
return(.plot_network_validation(
x = x,
metric = validation_metric,
series = series,
facet = facet,
empty = empty
))
}
if (type == "coverage") {
return(.plot_network_coverage(
x = x,
series = series,
facet = facet,
empty = empty
))
}
if (type == "seasonal_error") {
return(.plot_network_seasonal_error(
x = x,
metric = seasonal_metric,
series = series,
facet = facet,
empty = empty
))
}
.plot_network_compare(
x = x,
series = series,
facet = facet,
colour_by = compare_colour,
empty = empty
)
}
#' @keywords internal
#' @noRd
.plot_network_interpolation <- function(x, series = NULL, start = NULL, end = NULL,
show_pi = TRUE, show_fit = FALSE,
free_y = TRUE, ncol = 1,
empty = c("plot", "error")) {
empty <- match.arg(empty)
TIME <- fit <- original <- pi_hi <- pi_lo <- value <- jump_correction <- NULL
d <- .network_interp_diag_subset(x, series = series, start = start, end = end)
if (nrow(d) == 0) {
if (empty == "error") stop("No rows available for plotting.")
return(.network_interp_empty_plot("No rows available for plotting."))
}
method <- attr(x, "network_method", exact = TRUE)
p <- ggplot2::ggplot(d, ggplot2::aes(x = TIME)) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Time",
y = "Dendrometer value",
title = paste("Network interpolation:", method)
)
if (show_pi && any(is.finite(d$pi_lo) & is.finite(d$pi_hi))) {
p <- p +
ggplot2::geom_ribbon(
ggplot2::aes(ymin = pi_lo, ymax = pi_hi),
fill = "grey70",
alpha = 0.30,
na.rm = TRUE
)
}
p <- p +
ggplot2::geom_line(
ggplot2::aes(y = original),
linewidth = 0.35,
colour = "grey45",
na.rm = TRUE
) +
ggplot2::geom_line(
ggplot2::aes(y = value),
linewidth = 0.45,
colour = "#1F78B4",
na.rm = TRUE
)
if (show_fit && any(is.finite(d$fit))) {
p <- p +
ggplot2::geom_line(
ggplot2::aes(y = fit),
linewidth = 0.35,
linetype = 2,
colour = "black",
na.rm = TRUE
)
}
if (any(d$interp, na.rm = TRUE)) {
p <- p +
ggplot2::geom_point(
data = d[d$interp, , drop = FALSE],
ggplot2::aes(x = TIME, y = value),
colour = "#D7301F",
size = 1.2,
inherit.aes = FALSE,
na.rm = TRUE
)
}
jump_df <- d[!is.na(d$jump_correction), , drop = FALSE]
if (nrow(jump_df) > 0) {
p <- p +
ggplot2::geom_point(
data = jump_df,
ggplot2::aes(x = TIME, y = value),
colour = "darkgreen",
size = 2,
inherit.aes = FALSE,
na.rm = TRUE
)
}
p +
ggplot2::facet_wrap(
~series,
scales = if (free_y) "free_y" else "fixed",
ncol = ncol
)
}
#' @keywords internal
#' @noRd
.plot_network_uncertainty <- function(x, series = NULL, start = NULL, end = NULL,
free_y = TRUE, ncol = 1,
empty = c("plot", "error")) {
empty <- match.arg(empty)
TIME <- pi_width <- NULL
d <- .network_interp_diag_subset(x, series = series, start = start, end = end)
d <- d[is.finite(d$pi_width), , drop = FALSE]
if (nrow(d) == 0) {
if (empty == "error") stop("No prediction intervals available for plotting.")
return(.network_interp_empty_plot("No prediction intervals available for plotting."))
}
ggplot2::ggplot(d, ggplot2::aes(x = TIME, y = pi_width)) +
ggplot2::geom_line(linewidth = 0.45, na.rm = TRUE) +
ggplot2::geom_point(
data = d[d$interp, , drop = FALSE],
size = 1.1,
na.rm = TRUE
) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Time",
y = "95% interval width",
title = "Interpolation uncertainty over time"
) +
ggplot2::facet_wrap(
~series,
scales = if (free_y) "free_y" else "fixed",
ncol = ncol
)
}
#' @keywords internal
#' @noRd
.plot_network_availability <- function(x, series = NULL, start = NULL, end = NULL,
free_y = TRUE, ncol = 1,
empty = c("plot", "error")) {
TIME <- ref_n <- NULL
empty <- match.arg(empty)
d <- .network_interp_diag_subset(x, series = series, start = start, end = end)
if (nrow(d) == 0 || !any(is.finite(d$ref_n))) {
if (empty == "error") stop("No reference availability data available for plotting.")
return(.network_interp_empty_plot("No reference availability data available for plotting."))
}
ggplot2::ggplot(d, ggplot2::aes(x = TIME, y = ref_n)) +
ggplot2::geom_line(linewidth = 0.45, na.rm = TRUE) +
ggplot2::geom_point(
data = d[d$interp, , drop = FALSE],
colour = "#D7301F",
size = 1.2,
na.rm = TRUE
) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Time",
y = "Available reference sensors",
title = "Reference network availability"
) +
ggplot2::facet_wrap(
~series,
scales = if (free_y) "free_y" else "fixed",
ncol = ncol
)
}
#' @keywords internal
#' @noRd
.plot_network_summary <- function(x,
metric = c("n_imputed", "n_remaining_na", "n_missing_input", "mean_pi_width", "median_pi_width", "mean_ref_n", "median_ref_n", "n_gap_jumps_removed", "max_abs_gap_jump"),
empty = c("plot", "error")) {
metric <- match.arg(metric)
empty <- match.arg(empty)
series <- value <- NULL
s <- attr(x, "network_summary", exact = TRUE)
if (is.null(s) || nrow(s) == 0) {
if (empty == "error") stop("No summary information available for plotting.")
return(.network_interp_empty_plot("No summary information available for plotting."))
}
s <- tibble::as_tibble(s)
s$value <- s[[metric]]
ggplot2::ggplot(s, ggplot2::aes(x = series, y = value)) +
ggplot2::geom_col() +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Series",
y = metric,
title = paste("Network interpolation summary:", metric)
)
}
#' @keywords internal
#' @noRd
.plot_network_validation <- function(x,
metric = c("MAE", "MAPE", "MdAPE", "RMSE", "Bias", "Success_rate", "Mean_PI_width", "Mean_ref_n", "End_value_abs_diff", "Max_abs_diff_within_gap"),
series = NULL,
facet = TRUE,
empty = c("plot", "error")) {
metric <- match.arg(metric)
empty <- match.arg(empty)
gap_steps <- value <- NULL
s <- .network_interp_validation_subset(x, series = series)
if (nrow(s) == 0) {
if (empty == "error") stop("No validation results available for plotting.")
return(.network_interp_empty_plot("No validation results available for plotting."))
}
s$value <- s[[metric]]
p <- ggplot2::ggplot(
s,
ggplot2::aes(x = gap_steps, y = value, group = series, colour = series)
) +
ggplot2::geom_line(linewidth = 0.5, na.rm = TRUE) +
ggplot2::geom_point(size = 1.8, na.rm = TRUE) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Artificial gap length (steps)",
y = metric,
colour = "Series",
title = paste("Validation:", metric)
)
if (facet) {
p + ggplot2::facet_wrap(~series, scales = "free_y")
} else {
p
}
}
#' @keywords internal
#' @noRd
.plot_network_coverage <- function(x,
series = NULL,
facet = TRUE,
empty = c("plot", "error")) {
gap_steps <- PI_coverage_95 <- NULL
empty <- match.arg(empty)
s <- .network_interp_validation_subset(x, series = series)
if (nrow(s) == 0 || !any(is.finite(s$PI_coverage_95))) {
if (empty == "error") stop("No validation coverage results available for plotting.")
return(.network_interp_empty_plot("No validation coverage results available for plotting."))
}
p <- ggplot2::ggplot(
s,
ggplot2::aes(x = gap_steps, y = PI_coverage_95, group = series, colour = series)
) +
ggplot2::geom_hline(yintercept = 0.95, linetype = 2) +
ggplot2::geom_line(linewidth = 0.5, na.rm = TRUE) +
ggplot2::geom_point(size = 1.8, na.rm = TRUE) +
ggplot2::coord_cartesian(ylim = c(0, 1)) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Artificial gap length (steps)",
y = "95% PI coverage",
colour = "Series",
title = "Validation of interval coverage"
)
if (facet) {
p + ggplot2::facet_wrap(~series, scales = "fixed")
} else {
p
}
}
#' @keywords internal
#' @noRd
.plot_network_compare <- function(x,
series = NULL,
facet = TRUE,
colour_by = c("series", "gap_steps"),
empty = c("plot", "error")) {
truth <- pred <- gap_steps <- gap_steps_f <- NULL
empty <- match.arg(empty)
colour_by <- match.arg(colour_by)
pdat <- .network_interp_validation_points_subset(x, series = series)
pdat <- pdat[pdat$filled & is.finite(pdat$truth) & is.finite(pdat$pred), , drop = FALSE]
if (nrow(pdat) == 0) {
if (empty == "error") stop("No point-level validation comparison data available for plotting.")
return(.network_interp_empty_plot("No point-level validation comparison data available for plotting."))
}
pdat$gap_steps_f <- factor(pdat$gap_steps)
aes_colour <- if (colour_by == "series") {
ggplot2::aes(colour = series)
} else {
ggplot2::aes(colour = gap_steps_f)
}
p <- ggplot2::ggplot(pdat, ggplot2::aes(x = truth, y = pred)) +
ggplot2::geom_abline(slope = 1, intercept = 0, linetype = 2) +
ggplot2::geom_point(aes_colour, alpha = 0.75, size = 1.8, na.rm = TRUE) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Actual value",
y = "Interpolated value",
colour = if (colour_by == "series") "Series" else "Gap steps",
title = "Actual vs interpolated values from validation"
)
if (facet) {
p + ggplot2::facet_wrap(~series, scales = "free")
} else {
p
}
}
#' @keywords internal
#' @noRd
.plot_network_seasonal_error <- function(x,
metric = c("mean_error", "mean_abs_error", "mean_abs_pct_error", "rmse", "coverage"),
series = NULL,
facet = TRUE,
empty = c("plot", "error")) {
metric <- match.arg(metric)
empty <- match.arg(empty)
doy <- value <- NULL
s <- .network_interp_validation_seasonal_subset(x, series = series)
if (nrow(s) == 0) {
if (empty == "error") stop("No seasonal validation results available for plotting.")
return(.network_interp_empty_plot("No seasonal validation results available for plotting."))
}
s$value <- s[[metric]]
p <- ggplot2::ggplot(
s,
ggplot2::aes(x = doy, y = value, colour = series, group = series)
) +
ggplot2::geom_line(linewidth = 0.45, na.rm = TRUE) +
ggplot2::geom_point(size = 1.2, na.rm = TRUE) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Day of year",
y = metric,
colour = "Series",
title = paste("Seasonal validation diagnostic:", metric)
)
if (facet) {
p + ggplot2::facet_wrap(~series, scales = "free_y")
} else {
p
}
}
#' @keywords internal
#' @noRd
.network_interp_diag_subset <- function(x, series = NULL, start = NULL, end = NULL) {
d <- attr(x, "network_diagnostics", exact = TRUE)
if (is.null(d)) {
stop("No diagnostics found on this object.")
}
d <- tibble::as_tibble(d)
if (!is.null(series)) {
d <- d[d$series %in% series, , drop = FALSE]
}
if (!is.null(start)) {
d <- d[d$TIME >= as.POSIXct(start), , drop = FALSE]
}
if (!is.null(end)) {
d <- d[d$TIME <= as.POSIXct(end), , drop = FALSE]
}
d
}
#' @keywords internal
#' @noRd
.network_interp_validation_subset <- function(x, series = NULL) {
s <- attr(x, "network_validation_summary", exact = TRUE)
if (is.null(s)) {
return(tibble::tibble())
}
s <- tibble::as_tibble(s)
if (!is.null(series)) {
s <- s[s$series %in% series, , drop = FALSE]
}
s
}
#' @keywords internal
#' @noRd
.network_interp_validation_points_subset <- function(x, series = NULL) {
s <- attr(x, "network_validation_points", exact = TRUE)
if (is.null(s)) {
return(tibble::tibble())
}
s <- tibble::as_tibble(s)
if (!is.null(series)) {
s <- s[s$series %in% series, , drop = FALSE]
}
s
}
#' @keywords internal
#' @noRd
.network_interp_validation_seasonal_subset <- function(x, series = NULL) {
s <- attr(x, "network_validation_seasonal", exact = TRUE)
if (is.null(s)) {
return(tibble::tibble())
}
s <- tibble::as_tibble(s)
if (!is.null(series)) {
s <- s[s$series %in% series, , drop = FALSE]
}
s
}
#' @keywords internal
#' @noRd
.network_interp_empty_plot <- function(label) {
ggplot2::ggplot() +
ggplot2::annotate("text", x = 0, y = 0, label = label, size = 5) +
ggplot2::xlim(-1, 1) +
ggplot2::ylim(-1, 1) +
ggplot2::labs(x = NULL, y = NULL) +
ggplot2::theme_void()
}
Try the dendRoAnalyst package in your browser
Any scripts or data that you put into this service are public.
dendRoAnalyst documentation built on Oct. 2, 2026, 9:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .network_interp_empty_plot .network_interp_validation_seasonal_subset .network_interp_validation_points_subset .network_interp_validation_subset .network_interp_diag_subset .plot_network_seasonal_error .plot_network_compare .plot_network_coverage .plot_network_validation .plot_network_summary .plot_network_availability .plot_network_uncertainty .plot_network_interpolation plot.network_interpolation network.interpolation
Documented in network.interpolation plot.network_interpolation
#' @title Interpolate missing dendrometer values using a site network
#'
#' @description
#' Fills \code{NA} gaps in one or more focal dendrometer series by borrowing
#' information from a reference network recorded at the same site and on the
#' same time grid.
#'
#' The focal dataset \code{df} is treated as the master time axis. The
#' reference network is aligned to that axis, so missing timestamps already
#' inserted into \code{df} are preserved during interpolation.
#'
#' Two interpolation modes are available:
#' \itemize{
#' \item \code{niMethod = "proportional"} uses the average relative change in
#' the reference network between consecutive timestamps and propagates that
#' change to the focal series.
#' \item \code{niMethod = "linear"} fits a cross-sensor regression between
#' reference values at time \eqn{t-1} and \eqn{t}, and predicts the focal
#' value at time \eqn{t} from the focal value at time \eqn{t-1}.
#' }
#'
#' Optional bootstrap-based uncertainty limits can be attached, synthetic-gap
#' validation can be used to assess recovery performance, and an optional
#' post-processing step can detect and correct abrupt post-gap jumps in the
#' interpolated series.
#'
#' @details
#' The function assumes that missing timestamps have already been inserted into
#' \code{df} and corresponding measurements set to \code{NA}, for example with
#' \code{\link{dm.na.interpolation}}.
#'
#' For a focal series \eqn{a}, let \eqn{a_{t-1}} be the last available focal
#' value before the current timestamp. Let \eqn{b_{t-1}} and \eqn{b_t} be the
#' vectors of reference values at the previous and current timestamps.
#'
#' In proportional mode, the relative reference change is
#' \deqn{
#' \delta = \mathrm{mean}\left(\frac{b_t - b_{t-1}}{\max(b_{t-1}, \epsilon)}\right)
#' }
#' and the focal estimate is
#' \deqn{
#' \hat{a}_t = a_{t-1}(1 + \delta).
#' }
#' Bootstrap resampling of reference sensors is used to estimate a 95\%
#' interval.
#'
#' In linear mode, a simple regression of \eqn{b_t} on \eqn{b_{t-1}} is fit
#' across reference sensors at each step, and the focal value is predicted from
#' \eqn{a_{t-1}} with a 95\% prediction interval.
#'
#' If fewer than two valid reference sensors are available at a timestamp pair,
#' or if the focal previous value is missing, the focal value is left unchanged.
#'
#' If \code{correct_gap_jumps = TRUE}, the function inspects the first observed
#' point immediately after each imputed run. If a jump larger than
#' \code{jump_threshold} is detected, the estimated offset is removed from that
#' point onward. This is useful when a dendrometer resumes from a shifted
#' baseline after a gap.
#'
#' If \code{assess = TRUE}, the function inserts artificial gaps into observed
#' sections of each focal series, interpolates them with the chosen method, and
#' compares predictions with the true values. This produces both per-gap
#' summaries and seasonal diagnostics of interpolation error.
#'
#' The returned object is a data frame with class
#' \code{"network_interpolation"}.
#'
#' @param df Data frame containing focal dendrometer data. The first column must
#' be datetime (\code{POSIXct}, \code{Date}, or character convertible by
#' \code{lubridate::ymd_hms}). Remaining columns are focal series.
#' @param referenceDF Data frame containing reference dendrometer data. The
#' first column must be datetime, and remaining columns are reference series.
#' @param niMethod Character. Interpolation method. One of
#' \code{"proportional"} or \code{"linear"}.
#' @param n_boot Integer. Number of bootstrap resamples used in proportional
#' mode.
#' @param return_flags Logical. If \code{TRUE}, add a logical column
#' \code{<series>_interp} indicating which focal values were imputed.
#' @param return_pi Logical. If \code{TRUE}, add 95\% prediction-limit columns
#' \code{<series>_pi_lo} and \code{<series>_pi_hi}.
#' @param pi_for_all Logical. If \code{TRUE}, prediction limits are also
#' attached for non-imputed rows where the model could be formed. If
#' \code{FALSE}, limits are returned only for imputed rows.
#' @param return_fit Logical. If \code{TRUE}, add \code{<series>_fit} columns
#' containing the model-implied fitted value.
#' @param assess Logical. If \code{TRUE}, run synthetic-gap validation.
#' @param assess_lengths_steps Integer vector. Artificial gap lengths,
#' expressed in number of rows/time steps, to test during validation.
#' @param assess_samples_per_length Integer. Number of sampled artificial gaps
#' per gap length and per focal series.
#' @param assess_buffer_steps Integer. Number of observed steps required before
#' and after each artificial gap when
#' \code{assess_use_only_observed = TRUE}.
#' @param assess_seed Integer or \code{NULL}. Random seed for reproducible
#' validation sampling.
#' @param assess_use_only_observed Logical. If \code{TRUE}, synthetic gaps are
#' sampled only from windows fully observed in the focal series, including
#' the requested buffer.
#' @param progress Logical. If \code{TRUE}, display a single-line progress
#' indicator.
#' @param correct_gap_jumps Logical. If \code{TRUE}, inspect the first observed
#' point after each imputed run and optionally remove a detected jump.
#' @param jump_threshold Numeric. Minimum absolute jump size required before a
#' post-gap correction is applied. Required if
#' \code{correct_gap_jumps = TRUE}.
#' @param jump_adjustment_method Character. Method used to estimate the
#' post-gap offset. One of \code{"window_median"} or \code{"point_diff"}.
#' @param jump_adjust_window Integer. Window size used by
#' \code{jump_adjustment_method = "window_median"}.
#' @param jump_set_point_na Logical. If \code{TRUE}, the jump point itself is
#' set to \code{NA} after correction.
#'
#' @return
#' A data frame with class \code{"network_interpolation"}.
#'
#' The data frame contains the original \code{TIME} column, focal series
#' columns, and optionally:
#' \describe{
#' \item{\code{<series>_interp}}{Logical flag for imputed rows.}
#' \item{\code{<series>_pi_lo}, \code{<series>_pi_hi}}{95\% prediction
#' limits.}
#' \item{\code{<series>_fit}}{Model-implied fitted values.}
#' \item{\code{<series>_jump_correction}}{Numeric jump offset removed at the
#' first observed point after an imputed run. \code{NA} means no correction
#' was applied.}
#' }
#'
#' In addition, the following attributes are attached:
#' \describe{
#' \item{\code{network_original}}{Original focal data on the aligned time
#' grid.}
#' \item{\code{network_reference}}{Reference data aligned to the focal time
#' grid.}
#' \item{\code{network_diagnostics}}{Long-format diagnostics for plotting.}
#' \item{\code{network_summary}}{Per-series interpolation summary table.}
#' \item{\code{network_jump_corrections}}{Table of inspected and corrected
#' post-gap jumps.}
#' \item{\code{network_validation_raw}}{Per-gap validation summary table, or
#' \code{NULL}.}
#' \item{\code{network_validation_points}}{Point-level validation results, or
#' \code{NULL}.}
#' \item{\code{network_validation_summary}}{Summarized validation metrics by
#' series and gap length, or \code{NULL}.}
#' \item{\code{network_validation_seasonal}}{Seasonal validation diagnostics
#' by day of year, or \code{NULL}.}
#' \item{\code{network_method}}{Interpolation method used.}
#' \item{\code{network_n_boot}}{Number of bootstrap resamples used.}
#' \item{\code{network_assessed}}{Logical indicating whether validation was
#' run.}
#' \item{\code{network_jump_settings}}{List of jump-correction settings.}
#' }
#'
#' @section Interpretation:
#' Smaller validation metrics such as \code{MAE}, \code{RMSE}, \code{MAPE}, and
#' \code{MdAPE} indicate better gap recovery. \code{PI_coverage_95} indicates
#' how well the nominal 95\% uncertainty interval captures the true values
#' under synthetic-gap validation. Seasonal diagnostics help identify periods
#' of the year when network interpolation is more or less reliable.
#'
#' If many post-gap corrections are applied or \code{max_abs_gap_jump} is
#' large, it may indicate re-zeroing, sensor reset, or baseline
#' discontinuities after gaps.
#'
#' @section Notes:
#' \itemize{
#' \item The focal dataset \code{df} defines the output time grid.
#' \item Reference series are aligned to that grid using the datetime column.
#' \item At least two valid reference sensors are required at consecutive
#' timestamps for interpolation.
#' \item Jump correction is a post-processing step and does not alter the
#' core interpolation logic.
#' }
#'
#' @examples
#' \donttest{
#' #library(dendRoAnalyst)
#' #data("gf_nepa17")
#'
#' ## Create an artificial focal gap
#' #df1 <- gf_nepa17
#' #df1[40:50, "T2"] <- NA
#'
#' ## Build a small reference network
#' #ref <- cbind(gf_nepa17, gf_nepa17[, 2:3], gf_nepa17[, 2:3])
#' #colnames(ref) <- c("Time", "T1", "T2", "T3", "T4", "T5", "T6")
#'
#' ## Interpolate with uncertainty limits
#' #out <- network.interpolation(
#' # df1, ref,
#' # niMethod = "proportional",
#' # n_boot = 100,
#' # return_flags = TRUE,
#' # return_pi = TRUE
#' #)
#'
#' #head(out, 10)
#'
#' # Interpolate with jump correction and validation
#' #out2 <- network.interpolation(
#' # df1, ref,
#' # niMethod = "proportional",
#' # n_boot = 100,
#' # return_flags = TRUE,
#' # return_pi = TRUE,
#' # assess = TRUE,
#' # assess_lengths_steps = c(1, 2, 4),
#' # correct_gap_jumps = TRUE,
#' # jump_threshold = 0.05,
#' # jump_adjustment_method = "window_median",
#' # jump_adjust_window = 5
#' #)
#' # Plotting
#' #plot(out2)
#' #plot(out2, type = "compare")
#' #plot(out2, type = "seasonal_error")
#'
#' # Extracting the information using attr()
#' #attr(out, "network_validation_summary")
#' #attr(out, "network_validation_points")
#' #attr(out, "network_validation_seasonal")
#' }
#'
#' @seealso
#' \code{\link{plot.network_interpolation}},
#' \code{\link{dm.na.interpolation}}
#'
#' @export
network.interpolation <- function(df, referenceDF,
niMethod = c("proportional", "linear"),
n_boot = 1000,
return_flags = TRUE,
return_pi = TRUE,
pi_for_all = FALSE,
return_fit = FALSE,
assess = FALSE,
assess_lengths_steps = c(1, 2, 3, 6, 12, 24),
assess_samples_per_length = 20,
assess_buffer_steps = 1,
assess_seed = NULL,
assess_use_only_observed = TRUE,
progress = interactive(),
correct_gap_jumps = FALSE,
jump_threshold = NULL,
jump_adjustment_method = c("window_median", "point_diff"),
jump_adjust_window = 5,
jump_set_point_na = FALSE) {
niMethod <- match.arg(niMethod)
jump_adjustment_method <- match.arg(jump_adjustment_method)
if (isTRUE(correct_gap_jumps) &&
(is.null(jump_threshold) || !is.finite(jump_threshold) || jump_threshold <= 0)) {
stop("When correct_gap_jumps = TRUE, provide a positive numeric 'jump_threshold'.")
}
Bias <- End_value_abs_diff <- MAE <- MAPE <- Max_abs_diff_within_gap <- MdAPE <- Mean_PI_width <- Mean_ref_n <- PI_coverage_95 <- RMSE <- Success_rate <- TIME <- abs_pct_error <- applied_offset <- case_id <- covered <- doy <- error <- filled <- gap_steps <- jump_correction <- jump_removed <- jump_time <- max_abs_gap_jump <- mean_abs_error <- mean_abs_pct_error <- mean_error <- mean_pred <- mean_truth <- n_gap_jumps_removed <- pred <- raw_jump <- ref_n <- rmse <- series <- start_time <- end_time <- truth <- NULL
to_posixct <- function(x) {
if (inherits(x, "POSIXct")) return(x)
if (inherits(x, "Date")) return(as.POSIXct(x))
lubridate::ymd_hms(x)
}
mean_or_na <- function(x) {
if (!length(x) || all(is.na(x))) return(NA_real_)
mean(x, na.rm = TRUE)
}
median_or_na <- function(x) {
if (!length(x) || all(is.na(x))) return(NA_real_)
stats::median(x, na.rm = TRUE)
}
boot_mean_delta_fast <- function(delta, n_boot) {
delta <- delta[is.finite(delta)]
n <- length(delta)
if (n < 2L) {
return(list(
mean = if (n) mean(delta) else NA_real_,
lo = NA_real_,
hi = NA_real_
))
}
idx <- matrix(sample.int(n, size = n * n_boot, replace = TRUE), nrow = n)
means <- colMeans(matrix(delta[idx], nrow = n))
list(
mean = mean(means, na.rm = TRUE),
lo = stats::quantile(means, 0.025, na.rm = TRUE, names = FALSE),
hi = stats::quantile(means, 0.975, na.rm = TRUE, names = FALSE)
)
}
fast_linear_pred <- function(x, y, x0, level = 0.95) {
ok <- is.finite(x) & is.finite(y)
x <- x[ok]
y <- y[ok]
n <- length(x)
if (n < 2L || !is.finite(x0)) {
return(c(fit = NA_real_, lwr = NA_real_, upr = NA_real_))
}
xbar <- mean(x)
ybar <- mean(y)
sxx <- sum((x - xbar)^2)
if (!is.finite(sxx) || sxx <= 0) {
return(c(fit = NA_real_, lwr = NA_real_, upr = NA_real_))
}
beta <- sum((x - xbar) * (y - ybar)) / sxx
alpha <- ybar - beta * xbar
fit <- alpha + beta * x0
if (n <= 2L) {
return(c(fit = fit, lwr = NA_real_, upr = NA_real_))
}
resid <- y - (alpha + beta * x)
mse <- sum(resid^2) / (n - 2L)
if (!is.finite(mse)) {
return(c(fit = fit, lwr = NA_real_, upr = NA_real_))
}
se_pred <- sqrt(mse * (1 + 1 / n + ((x0 - xbar)^2 / sxx)))
tval <- stats::qt((1 + level) / 2, df = n - 2L)
c(
fit = fit,
lwr = fit - tval * se_pred,
upr = fit + tval * se_pred
)
}
estimate_jump_offset_gap <- function(x, j,
adjustment_method = c("window_median", "point_diff"),
adjust_window = 5) {
adjustment_method <- match.arg(adjustment_method)
if (length(x) == 0 || j <= 1 || j > length(x)) {
return(NA_real_)
}
if (adjustment_method == "point_diff") {
if (is.na(x[j]) || is.na(x[j - 1])) {
return(NA_real_)
}
return(x[j] - x[j - 1])
}
n <- length(x)
pre_idx <- seq.int(max(1, j - adjust_window), j - 1)
post_idx <- seq.int(j, min(n, j + adjust_window - 1))
pre_vals <- x[pre_idx]
post_vals <- x[post_idx]
pre_vals <- pre_vals[!is.na(pre_vals)]
post_vals <- post_vals[!is.na(post_vals)]
if (length(pre_vals) < 2 || length(post_vals) < 2) {
if (is.na(x[j]) || is.na(x[j - 1])) {
return(NA_real_)
}
return(x[j] - x[j - 1])
}
stats::median(post_vals) - stats::median(pre_vals)
}
apply_jump_correction_gap <- function(x, j, offset, set_jump_point_na = FALSE) {
if (length(x) == 0 || is.na(offset) || j > length(x)) {
return(x)
}
idx <- j:length(x)
idx <- idx[!is.na(x[idx])]
if (length(idx) > 0) {
x[idx] <- x[idx] - offset
}
if (isTRUE(set_jump_point_na) && j <= length(x)) {
x[j] <- NA_real_
}
x
}
correct_gap_end_jumps <- function(value, fit, pi_lo, pi_hi,
interp_flag, original_values, time_vec,
jump_threshold,
adjustment_method = c("window_median", "point_diff"),
adjust_window = 5,
set_jump_point_na = FALSE) {
adjustment_method <- match.arg(adjustment_method)
n <- length(value)
jump_correction <- rep(NA_real_, n)
if (!any(interp_flag, na.rm = TRUE)) {
return(list(
value = value,
fit = fit,
pi_lo = pi_lo,
pi_hi = pi_hi,
jump_correction = jump_correction,
table = tibble::tibble()
))
}
r <- rle(interp_flag)
ends <- cumsum(r$lengths)
gap_end_idx <- ends[r$values]
candidate_idx <- gap_end_idx + 1L
candidate_idx <- candidate_idx[candidate_idx <= n]
if (!length(candidate_idx)) {
return(list(
value = value,
fit = fit,
pi_lo = pi_lo,
pi_hi = pi_hi,
jump_correction = jump_correction,
table = tibble::tibble()
))
}
rows <- vector("list", length(candidate_idx))
rr <- 0L
for (j in candidate_idx) {
if (!is.finite(original_values[j])) next
raw_jump <- if (j > 1 && is.finite(value[j]) && is.finite(value[j - 1])) {
value[j] - value[j - 1]
} else {
NA_real_
}
offset <- estimate_jump_offset_gap(
x = value,
j = j,
adjustment_method = adjustment_method,
adjust_window = adjust_window
)
removed <- is.finite(offset) && abs(offset) >= jump_threshold
if (removed) {
value <- apply_jump_correction_gap(value, j, offset, set_jump_point_na = set_jump_point_na)
fit <- apply_jump_correction_gap(fit, j, offset, set_jump_point_na = FALSE)
pi_lo <- apply_jump_correction_gap(pi_lo, j, offset, set_jump_point_na = FALSE)
pi_hi <- apply_jump_correction_gap(pi_hi, j, offset, set_jump_point_na = FALSE)
jump_correction[j] <- offset
}
rr <- rr + 1L
rows[[rr]] <- tibble::tibble(
jump_index = j,
jump_time = time_vec[j],
raw_jump = raw_jump,
estimated_offset = offset,
applied_offset = if (removed) offset else NA_real_,
jump_removed = removed
)
}
jump_table <- if (rr > 0L) dplyr::bind_rows(rows[seq_len(rr)]) else tibble::tibble()
list(
value = value,
fit = fit,
pi_lo = pi_lo,
pi_hi = pi_hi,
jump_correction = jump_correction,
table = jump_table
)
}
interpolate_one_series <- function(a_raw, ref_mat, ref_ok, niMethod, n_boot, tiny = 1e-6) {
n <- length(a_raw)
a <- as.numeric(a_raw)
interp_flag <- rep(FALSE, n)
ref_n <- rep(NA_integer_, n)
fit <- rep(NA_real_, n)
pi_lo <- rep(NA_real_, n)
pi_hi <- rep(NA_real_, n)
for (i in 2:n) {
a_prev <- a[i - 1L]
if (!is.finite(a_prev)) next
ok <- ref_ok[i - 1L, ] & ref_ok[i, ]
n_ok <- sum(ok)
ref_n[i] <- n_ok
if (n_ok < 2L) next
b_prev <- ref_mat[i - 1L, ok]
b_cur <- ref_mat[i, ok]
if (niMethod == "proportional") {
rel <- (b_cur - b_prev) / pmax(b_prev, tiny)
bm <- boot_mean_delta_fast(rel, n_boot)
if (!is.finite(bm$mean)) next
fit_i <- a_prev * (1 + bm$mean)
lo_i <- if (is.na(bm$lo)) NA_real_ else a_prev * (1 + bm$lo)
hi_i <- if (is.na(bm$hi)) NA_real_ else a_prev * (1 + bm$hi)
} else {
pred <- fast_linear_pred(b_prev, b_cur, a_prev, level = 0.95)
fit_i <- unname(pred["fit"])
lo_i <- unname(pred["lwr"])
hi_i <- unname(pred["upr"])
if (!is.finite(fit_i)) next
}
fit[i] <- fit_i
pi_lo[i] <- lo_i
pi_hi[i] <- hi_i
if (is.na(a[i]) && is.finite(fit_i)) {
a[i] <- fit_i
interp_flag[i] <- TRUE
}
}
list(
value = a,
interp_flag = interp_flag,
ref_n = ref_n,
fit = fit,
pi_lo = pi_lo,
pi_hi = pi_hi
)
}
df <- tibble::as_tibble(df)
referenceDF <- tibble::as_tibble(referenceDF)
names(df)[1] <- "TIME"
names(referenceDF)[1] <- "TIME"
df$TIME <- to_posixct(df$TIME)
referenceDF$TIME <- to_posixct(referenceDF$TIME)
if (any(is.na(df$TIME))) {
stop("Could not convert the first column of 'df' to POSIXct.")
}
if (any(is.na(referenceDF$TIME))) {
stop("Could not convert the first column of 'referenceDF' to POSIXct.")
}
df <- df[order(df$TIME), , drop = FALSE]
referenceDF <- referenceDF[order(referenceDF$TIME), , drop = FALSE]
time_min <- min(df$TIME, na.rm = TRUE)
time_max <- max(df$TIME, na.rm = TRUE)
referenceDF <- referenceDF[
referenceDF$TIME >= time_min & referenceDF$TIME <= time_max,
,
drop = FALSE
]
if (nrow(referenceDF) == 0) {
stop("'referenceDF' has no overlap with the time span of 'df'.")
}
referenceDF <- dplyr::left_join(df["TIME"], referenceDF, by = "TIME")
if (ncol(referenceDF) < 2L) {
stop("'referenceDF' must contain at least one reference series.")
}
ref_mat <- as.matrix(referenceDF[, -1, drop = FALSE])
suppressWarnings(storage.mode(ref_mat) <- "double")
ref_ok <- is.finite(ref_mat)
if (!any(ref_ok)) {
stop("No usable reference values remain after aligning 'referenceDF' to the time grid of 'df'.")
}
original_df <- df
out <- df
add_cols <- list()
diag_list <- vector("list", max(1L, ncol(df) - 1L))
sum_list <- vector("list", max(1L, ncol(df) - 1L))
jump_tables <- vector("list", max(1L, ncol(df) - 1L))
n_series <- max(1L, ncol(df) - 1L)
progress_total <- n_series + if (assess) n_series * length(assess_lengths_steps) else 0L
progress_value <- 0L
.progress_last_width <- 0L
render_progress_line <- function(note = "") {
if (!isTRUE(progress) || progress_total <= 0L) return(invisible(NULL))
width_bar <- 30L
frac <- progress_value / progress_total
n_done <- floor(width_bar * frac)
bar <- paste0(
"[",
paste(rep("=", n_done), collapse = ""),
paste(rep("-", width_bar - n_done), collapse = ""),
"]"
)
txt <- sprintf(
"\r%s %3d%% (%d/%d) %s",
bar,
round(100 * frac),
progress_value,
progress_total,
note
)
pad <- max(0L, .progress_last_width - nchar(txt, type = "width"))
cat(txt, paste(rep(" ", pad), collapse = ""))
utils::flush.console()
.progress_last_width <<- nchar(txt, type = "width")
invisible(NULL)
}
update_progress <- function(step = 1L, note = "") {
progress_value <<- min(progress_total, progress_value + step)
render_progress_line(note)
invisible(NULL)
}
progress_note <- function(note = "") {
render_progress_line(note)
invisible(NULL)
}
finish_progress <- function() {
if (!isTRUE(progress)) return(invisible(NULL))
cat("\n")
utils::flush.console()
invisible(NULL)
}
if (isTRUE(progress) && progress_total > 0L) {
on.exit(finish_progress(), add = TRUE)
render_progress_line("Starting...")
}
for (j in seq.int(2, ncol(df))) {
colname <- names(df)[j]
progress_note(paste0("Interpolating series ", colname, "..."))
raw_a <- df[[j]]
res <- interpolate_one_series(
a_raw = raw_a,
ref_mat = ref_mat,
ref_ok = ref_ok,
niMethod = niMethod,
n_boot = n_boot
)
jump_correction <- rep(NA_real_, length(raw_a))
jump_tbl <- tibble::tibble()
if (isTRUE(correct_gap_jumps)) {
jc <- correct_gap_end_jumps(
value = res$value,
fit = res$fit,
pi_lo = res$pi_lo,
pi_hi = res$pi_hi,
interp_flag = res$interp_flag,
original_values = raw_a,
time_vec = df$TIME,
jump_threshold = jump_threshold,
adjustment_method = jump_adjustment_method,
adjust_window = jump_adjust_window,
set_jump_point_na = jump_set_point_na
)
res$value <- jc$value
res$fit <- jc$fit
res$pi_lo <- jc$pi_lo
res$pi_hi <- jc$pi_hi
jump_correction <- jc$jump_correction
if (nrow(jc$table) > 0) {
jump_tbl <- jc$table
jump_tbl$series <- colname
}
}
jump_tables[[j - 1L]] <- jump_tbl
out[[j]] <- res$value
if (return_flags) {
add_cols[[paste0(colname, "_interp")]] <- res$interp_flag
}
if (return_pi) {
out_pi_lo <- ifelse(res$interp_flag | pi_for_all, res$pi_lo, NA_real_)
out_pi_hi <- ifelse(res$interp_flag | pi_for_all, res$pi_hi, NA_real_)
add_cols[[paste0(colname, "_pi_lo")]] <- out_pi_lo
add_cols[[paste0(colname, "_pi_hi")]] <- out_pi_hi
}
if (return_fit) {
out_fit <- ifelse(res$interp_flag | pi_for_all, res$fit, NA_real_)
add_cols[[paste0(colname, "_fit")]] <- out_fit
}
if (isTRUE(correct_gap_jumps)) {
add_cols[[paste0(colname, "_jump_correction")]] <- jump_correction
}
diag_list[[j - 1L]] <- tibble::tibble(
TIME = df$TIME,
series = colname,
original = raw_a,
value = res$value,
observed = is.finite(raw_a),
interp = res$interp_flag,
remaining_na = is.na(res$value),
ref_n = res$ref_n,
fit = res$fit,
pi_lo = res$pi_lo,
pi_hi = res$pi_hi,
pi_width = res$pi_hi - res$pi_lo,
jump_correction = jump_correction
)
sum_list[[j - 1L]] <- tibble::tibble(
series = colname,
n_rows = length(raw_a),
n_missing_input = sum(is.na(raw_a)),
n_imputed = sum(res$interp_flag),
n_remaining_na = sum(is.na(res$value)),
mean_ref_n = mean_or_na(res$ref_n),
median_ref_n = median_or_na(res$ref_n),
mean_pi_width = mean_or_na((res$pi_hi - res$pi_lo)[res$interp_flag]),
median_pi_width = median_or_na((res$pi_hi - res$pi_lo)[res$interp_flag]),
n_gap_jumps_removed = sum(!is.na(jump_correction)),
max_abs_gap_jump = if (all(is.na(jump_correction))) NA_real_ else max(abs(jump_correction), na.rm = TRUE)
)
update_progress(note = paste0("Finished series ", colname))
}
if (length(add_cols)) {
out <- dplyr::bind_cols(out, tibble::as_tibble(add_cols))
}
diagnostics <- dplyr::bind_rows(diag_list)
summary_tbl <- dplyr::bind_rows(sum_list)
jump_corrections <- dplyr::bind_rows(jump_tables)
validation_raw <- NULL
validation_points <- NULL
validation_summary <- NULL
validation_seasonal <- NULL
if (assess) {
if (!is.null(assess_seed)) {
set.seed(assess_seed)
}
max_cases <- n_series * length(assess_lengths_steps) * assess_samples_per_length
max_points <- max_cases * max(assess_lengths_steps)
val_rows <- vector("list", max(1L, max_cases))
point_rows <- vector("list", max(1L, max_points))
val_id <- 0L
point_id <- 0L
push_val_row <- function(row) {
val_id <<- val_id + 1L
if (val_id > length(val_rows)) {
length(val_rows) <<- max(length(val_rows) * 2L, val_id)
}
val_rows[[val_id]] <<- row
}
push_point_row <- function(row) {
point_id <<- point_id + 1L
if (point_id > length(point_rows)) {
length(point_rows) <<- max(length(point_rows) * 2L, point_id)
}
point_rows[[point_id]] <<- row
}
for (j in seq.int(2, ncol(df))) {
colname <- names(df)[j]
raw_a <- as.numeric(df[[j]])
n <- length(raw_a)
if (n < 2L) {
update_progress(length(assess_lengths_steps), note = paste0("Skipping validation for ", colname))
next
}
finite_raw <- is.finite(raw_a)
cs_finite <- c(0L, cumsum(finite_raw))
for (gap_steps in as.integer(assess_lengths_steps)) {
progress_note(paste0("Validating ", colname, " | gap=", gap_steps, " step(s)"))
if (!is.finite(gap_steps) || gap_steps < 1L) {
update_progress(note = paste0("Skipped ", colname, " | gap=", gap_steps))
next
}
min_start <- 2L + assess_buffer_steps
max_start <- n - gap_steps - assess_buffer_steps + 1L
if (max_start >= min_start) {
starts <- seq.int(min_start, max_start)
gap_counts <- cs_finite[starts + gap_steps] - cs_finite[starts]
prev_ok <- finite_raw[starts - 1L]
valid <- (gap_counts == gap_steps) & prev_ok
if (assess_use_only_observed) {
win_left <- starts - assess_buffer_steps
win_right <- starts + gap_steps - 1L + assess_buffer_steps
win_counts <- cs_finite[win_right + 1L] - cs_finite[win_left]
valid <- valid & (win_counts == (gap_steps + 2L * assess_buffer_steps))
}
valid_starts <- starts[valid]
if (length(valid_starts)) {
starts_to_use <- if (length(valid_starts) > assess_samples_per_length) {
sample(valid_starts, assess_samples_per_length, replace = FALSE)
} else {
valid_starts
}
for (s in starts_to_use) {
idx_gap <- s:(s + gap_steps - 1L)
case_id <- paste(colname, gap_steps, s, sep = "_")
a_mask <- raw_a
a_mask[idx_gap] <- NA_real_
res_case <- interpolate_one_series(
a_raw = a_mask,
ref_mat = ref_mat,
ref_ok = ref_ok,
niMethod = niMethod,
n_boot = n_boot
)
if (isTRUE(correct_gap_jumps)) {
jc_case <- correct_gap_end_jumps(
value = res_case$value,
fit = res_case$fit,
pi_lo = res_case$pi_lo,
pi_hi = res_case$pi_hi,
interp_flag = res_case$interp_flag,
original_values = a_mask,
time_vec = df$TIME,
jump_threshold = jump_threshold,
adjustment_method = jump_adjustment_method,
adjust_window = jump_adjust_window,
set_jump_point_na = jump_set_point_na
)
res_case$value <- jc_case$value
res_case$fit <- jc_case$fit
res_case$pi_lo <- jc_case$pi_lo
res_case$pi_hi <- jc_case$pi_hi
}
truth <- raw_a[idx_gap]
pred <- res_case$value[idx_gap]
pi_lo_gap <- res_case$pi_lo[idx_gap]
pi_hi_gap <- res_case$pi_hi[idx_gap]
ref_n_gap <- res_case$ref_n[idx_gap]
filled <- res_case$interp_flag[idx_gap] & is.finite(pred)
n_points <- length(idx_gap)
n_filled <- sum(filled)
success_rate <- n_filled / n_points
for (k in seq_along(idx_gap)) {
err_k <- if (is.finite(pred[k]) && is.finite(truth[k])) pred[k] - truth[k] else NA_real_
ape_k <- if (is.finite(err_k) && is.finite(truth[k])) {
100 * abs(err_k) / pmax(abs(truth[k]), 1e-6)
} else {
NA_real_
}
cov_k <- if (is.finite(truth[k]) && is.finite(pi_lo_gap[k]) && is.finite(pi_hi_gap[k])) {
truth[k] >= pi_lo_gap[k] && truth[k] <= pi_hi_gap[k]
} else {
NA
}
push_point_row(tibble::tibble(
series = colname,
gap_steps = gap_steps,
case_id = case_id,
TIME = df$TIME[idx_gap[k]],
truth = truth[k],
pred = pred[k],
filled = filled[k],
error = err_k,
abs_pct_error = ape_k,
pi_lo = pi_lo_gap[k],
pi_hi = pi_hi_gap[k],
covered = cov_k,
ref_n = ref_n_gap[k]
))
}
if (n_filled > 0L) {
err <- pred[filled] - truth[filled]
denom <- pmax(abs(truth[filled]), 1e-6)
pe <- 100 * err / denom
pi_ok <- filled & is.finite(pi_lo_gap) & is.finite(pi_hi_gap)
PI_coverage_95 <- if (any(pi_ok)) {
mean(truth[pi_ok] >= pi_lo_gap[pi_ok] & truth[pi_ok] <= pi_hi_gap[pi_ok], na.rm = TRUE)
} else {
NA_real_
}
Mean_PI_width <- if (any(pi_ok)) {
mean(pi_hi_gap[pi_ok] - pi_lo_gap[pi_ok], na.rm = TRUE)
} else {
NA_real_
}
End_value_abs_diff <- if (is.finite(pred[length(pred)])) {
abs(pred[length(pred)] - truth[length(truth)])
} else {
NA_real_
}
Max_abs_diff_within_gap <- max(abs(err), na.rm = TRUE)
row <- tibble::tibble(
series = colname,
gap_steps = gap_steps,
case_id = case_id,
start_time = df$TIME[min(idx_gap)],
end_time = df$TIME[max(idx_gap)],
n_points = n_points,
n_filled = n_filled,
Success_rate = success_rate,
MAE = mean(abs(err), na.rm = TRUE),
MAPE = mean(abs(pe), na.rm = TRUE),
MdAPE = stats::median(abs(pe), na.rm = TRUE),
RMSE = sqrt(mean(err^2, na.rm = TRUE)),
Bias = mean(err, na.rm = TRUE),
PI_coverage_95 = PI_coverage_95,
Mean_PI_width = Mean_PI_width,
Mean_ref_n = mean(ref_n_gap, na.rm = TRUE),
End_value_abs_diff = End_value_abs_diff,
Max_abs_diff_within_gap = Max_abs_diff_within_gap
)
} else {
row <- tibble::tibble(
series = colname,
gap_steps = gap_steps,
case_id = case_id,
start_time = df$TIME[min(idx_gap)],
end_time = df$TIME[max(idx_gap)],
n_points = n_points,
n_filled = n_filled,
Success_rate = success_rate,
MAE = NA_real_,
MAPE = NA_real_,
MdAPE = NA_real_,
RMSE = NA_real_,
Bias = NA_real_,
PI_coverage_95 = NA_real_,
Mean_PI_width = NA_real_,
Mean_ref_n = mean(ref_n_gap, na.rm = TRUE),
End_value_abs_diff = NA_real_,
Max_abs_diff_within_gap = NA_real_
)
}
push_val_row(row)
}
}
}
update_progress(note = paste0("Validated ", colname, " | gap=", gap_steps))
}
}
if (val_id > 0L) {
validation_raw <- dplyr::bind_rows(val_rows[seq_len(val_id)])
validation_points <- dplyr::bind_rows(point_rows[seq_len(point_id)])
validation_summary <- validation_raw %>%
dplyr::group_by(series, gap_steps) %>%
dplyr::summarise(
n_cases = dplyr::n(),
n_points = sum(n_points, na.rm = TRUE),
n_filled = sum(n_filled, na.rm = TRUE),
Success_rate = mean_or_na(Success_rate),
MAE = mean_or_na(MAE),
MAPE = mean_or_na(MAPE),
MdAPE = mean_or_na(MdAPE),
RMSE = mean_or_na(RMSE),
Bias = mean_or_na(Bias),
PI_coverage_95 = mean_or_na(PI_coverage_95),
Mean_PI_width = mean_or_na(Mean_PI_width),
Mean_ref_n = mean_or_na(Mean_ref_n),
End_value_abs_diff = mean_or_na(End_value_abs_diff),
Max_abs_diff_within_gap = mean_or_na(Max_abs_diff_within_gap),
.groups = "drop"
) %>%
dplyr::arrange(series, gap_steps)
validation_seasonal <- validation_points %>%
dplyr::filter(filled, is.finite(truth), is.finite(pred)) %>%
dplyr::mutate(doy = as.integer(format(TIME, "%j"))) %>%
dplyr::group_by(series, doy) %>%
dplyr::summarise(
n_points = dplyr::n(),
mean_error = mean(error, na.rm = TRUE),
mean_abs_error = mean(abs(error), na.rm = TRUE),
mean_abs_pct_error = mean(abs_pct_error, na.rm = TRUE),
rmse = sqrt(mean(error^2, na.rm = TRUE)),
mean_truth = mean(truth, na.rm = TRUE),
mean_pred = mean(pred, na.rm = TRUE),
coverage = mean(covered, na.rm = TRUE),
.groups = "drop"
) %>%
dplyr::arrange(series, doy)
} else {
validation_raw <- tibble::tibble()
validation_points <- tibble::tibble()
validation_summary <- tibble::tibble()
validation_seasonal <- tibble::tibble()
warning("Assessment produced no results (no valid synthetic windows found).")
}
}
attr(out, "network_original") <- original_df
attr(out, "network_reference") <- referenceDF
attr(out, "network_diagnostics") <- diagnostics
attr(out, "network_summary") <- summary_tbl
attr(out, "network_jump_corrections") <- jump_corrections
attr(out, "network_validation_raw") <- validation_raw
attr(out, "network_validation_points") <- validation_points
attr(out, "network_validation_summary") <- validation_summary
attr(out, "network_validation_seasonal") <- validation_seasonal
attr(out, "network_method") <- niMethod
attr(out, "network_n_boot") <- n_boot
attr(out, "network_return_flags") <- return_flags
attr(out, "network_return_pi") <- return_pi
attr(out, "network_pi_for_all") <- pi_for_all
attr(out, "network_return_fit") <- return_fit
attr(out, "network_assessed") <- assess
attr(out, "network_jump_settings") <- list(
correct_gap_jumps = correct_gap_jumps,
jump_threshold = jump_threshold,
jump_adjustment_method = jump_adjustment_method,
jump_adjust_window = jump_adjust_window,
jump_set_point_na = jump_set_point_na
)
class(out) <- c("network_interpolation", class(out))
out
}
#' @title Plot method for network interpolation output
#'
#' @description
#' S3 plot method for objects returned by
#' \code{\link{network.interpolation}}.
#'
#' Depending on \code{type}, the method can display:
#' \itemize{
#' \item interpolated series through time,
#' \item uncertainty width,
#' \item reference-network availability,
#' \item summary diagnostics,
#' \item synthetic-gap validation metrics,
#' \item interval coverage,
#' \item actual versus interpolated validation points,
#' \item seasonal validation diagnostics.
#' }
#'
#' Jump-corrected post-gap points are marked in the interpolation plot when
#' jump correction was enabled in \code{\link{network.interpolation}}.
#'
#' @param x Object returned by \code{\link{network.interpolation}}.
#' @param type Character. Plot type. One of \code{"interpolation"},
#' \code{"uncertainty"}, \code{"availability"}, \code{"summary"},
#' \code{"validation"}, \code{"coverage"}, \code{"compare"}, or
#' \code{"seasonal_error"}.
#' @param series Optional character vector of focal series names to plot.
#' @param start Optional start datetime for time-based plots.
#' @param end Optional end datetime for time-based plots.
#' @param show_pi Logical. If \code{TRUE}, show 95\% prediction intervals in
#' the interpolation plot where available.
#' @param show_fit Logical. If \code{TRUE}, show model-implied fitted values in
#' the interpolation plot where available.
#' @param free_y Logical. If \code{TRUE}, facets use free y scales where
#' applicable.
#' @param ncol Integer. Number of facet columns.
#' @param summary_metric Character metric used when
#' \code{type = "summary"}. One of \code{"n_imputed"},
#' \code{"n_remaining_na"}, \code{"n_missing_input"},
#' \code{"mean_pi_width"}, \code{"median_pi_width"}, \code{"mean_ref_n"},
#' \code{"median_ref_n"}, \code{"n_gap_jumps_removed"}, or
#' \code{"max_abs_gap_jump"}.
#' @param validation_metric Character metric used when
#' \code{type = "validation"}. One of \code{"MAE"}, \code{"MAPE"},
#' \code{"MdAPE"}, \code{"RMSE"}, \code{"Bias"},
#' \code{"Success_rate"}, \code{"Mean_PI_width"}, \code{"Mean_ref_n"},
#' \code{"End_value_abs_diff"}, or \code{"Max_abs_diff_within_gap"}.
#' @param seasonal_metric Character metric used when
#' \code{type = "seasonal_error"}. One of \code{"mean_error"},
#' \code{"mean_abs_error"}, \code{"mean_abs_pct_error"}, \code{"rmse"},
#' or \code{"coverage"}.
#' @param compare_colour Character. Colouring used when
#' \code{type = "compare"}. One of \code{"series"} or
#' \code{"gap_steps"}.
#' @param facet Logical. If \code{TRUE}, use faceting where supported.
#' @param empty Character. Behavior when no suitable data are available:
#' \code{"plot"} returns an informative empty plot and \code{"error"} throws
#' an error.
#' @param ... Further arguments passed through the generic.
#'
#' @return A \code{ggplot2} object.
#'
#' @details
#' Plot types:
#' \describe{
#' \item{\code{"interpolation"}}{Shows original and interpolated focal series
#' through time. Imputed points are highlighted. If jump correction was
#' applied, corrected post-gap points are marked in green.}
#' \item{\code{"uncertainty"}}{Shows the width of the prediction interval
#' through time.}
#' \item{\code{"availability"}}{Shows the number of valid reference sensors
#' available at each step.}
#' \item{\code{"summary"}}{Shows per-series summary diagnostics produced by
#' the interpolation run.}
#' \item{\code{"validation"}}{Shows synthetic-gap recovery metrics across
#' tested gap lengths.}
#' \item{\code{"coverage"}}{Shows how often true values fall inside the
#' nominal 95\% prediction interval across gap lengths.}
#' \item{\code{"compare"}}{Scatter plot of actual versus interpolated values
#' from synthetic-gap validation.}
#' \item{\code{"seasonal_error"}}{Shows how validation error varies through
#' the season, summarized by day of year.}
#' }
#'
#' Validation-based plots require that
#' \code{network.interpolation(..., assess = TRUE)} was used to create
#' \code{x}.
#'
#' @examples
#' \donttest{
#' #library(dendRoAnalyst)
#' #data("gf_nepa17")
#'
#' #df1 <- gf_nepa17
#' #df1[40:50, "T2"] <- NA
#'
#' #ref <- cbind(gf_nepa17, gf_nepa17[, 2:3], gf_nepa17[, 2:3])
#' #colnames(ref) <- c("Time", "T1", "T2", "T3", "T4", "T5", "T6")
#'
#' #out <- network.interpolation(
#' # df1, ref,
#' # niMethod = "proportional",
#' # n_boot = 100,
#' # assess = TRUE,
#' # assess_lengths_steps = c(1, 2, 4),
#' # correct_gap_jumps = TRUE,
#' # jump_threshold = 0.05
#' #)
#'
#' #plot(out)
#' #plot(out, type = "uncertainty")
#' #plot(out, type = "availability")
#' #plot(out, type = "summary", summary_metric = "n_gap_jumps_removed")
#' #plot(out, type = "validation", validation_metric = "MAPE")
#' #plot(out, type = "coverage")
#' #plot(out, type = "compare")
#' #plot(out, type = "seasonal_error", seasonal_metric = "mean_abs_error")
#' }
#'
#' @seealso
#' \code{\link{network.interpolation}}
#'
#' @method plot network_interpolation
#' @export
plot.network_interpolation <- function(x,
type = c("interpolation", "uncertainty", "availability", "summary", "validation", "coverage", "compare", "seasonal_error"),
series = NULL,
start = NULL,
end = NULL,
show_pi = TRUE,
show_fit = FALSE,
free_y = TRUE,
ncol = 1,
summary_metric = c("n_imputed", "n_remaining_na", "n_missing_input", "mean_pi_width", "median_pi_width", "mean_ref_n", "median_ref_n", "n_gap_jumps_removed", "max_abs_gap_jump"),
validation_metric = c("MAE", "MAPE", "MdAPE", "RMSE", "Bias", "Success_rate", "Mean_PI_width", "Mean_ref_n", "End_value_abs_diff", "Max_abs_diff_within_gap"),
seasonal_metric = c("mean_error", "mean_abs_error", "mean_abs_pct_error", "rmse", "coverage"),
compare_colour = c("series", "gap_steps"),
facet = TRUE,
empty = c("plot", "error"),
...) {
type <- match.arg(type)
summary_metric <- match.arg(summary_metric)
validation_metric <- match.arg(validation_metric)
seasonal_metric <- match.arg(seasonal_metric)
compare_colour <- match.arg(compare_colour)
empty <- match.arg(empty)
if (!inherits(x, "network_interpolation")) {
stop("x must be an object returned by network.interpolation().")
}
if (type == "interpolation") {
return(.plot_network_interpolation(
x = x,
series = series,
start = start,
end = end,
show_pi = show_pi,
show_fit = show_fit,
free_y = free_y,
ncol = ncol,
empty = empty
))
}
if (type == "uncertainty") {
return(.plot_network_uncertainty(
x = x,
series = series,
start = start,
end = end,
free_y = free_y,
ncol = ncol,
empty = empty
))
}
if (type == "availability") {
return(.plot_network_availability(
x = x,
series = series,
start = start,
end = end,
free_y = free_y,
ncol = ncol,
empty = empty
))
}
if (type == "summary") {
return(.plot_network_summary(
x = x,
metric = summary_metric,
empty = empty
))
}
if (type == "validation") {
return(.plot_network_validation(
x = x,
metric = validation_metric,
series = series,
facet = facet,
empty = empty
))
}
if (type == "coverage") {
return(.plot_network_coverage(
x = x,
series = series,
facet = facet,
empty = empty
))
}
if (type == "seasonal_error") {
return(.plot_network_seasonal_error(
x = x,
metric = seasonal_metric,
series = series,
facet = facet,
empty = empty
))
}
.plot_network_compare(
x = x,
series = series,
facet = facet,
colour_by = compare_colour,
empty = empty
)
}
#' @keywords internal
#' @noRd
.plot_network_interpolation <- function(x, series = NULL, start = NULL, end = NULL,
show_pi = TRUE, show_fit = FALSE,
free_y = TRUE, ncol = 1,
empty = c("plot", "error")) {
empty <- match.arg(empty)
TIME <- fit <- original <- pi_hi <- pi_lo <- value <- jump_correction <- NULL
d <- .network_interp_diag_subset(x, series = series, start = start, end = end)
if (nrow(d) == 0) {
if (empty == "error") stop("No rows available for plotting.")
return(.network_interp_empty_plot("No rows available for plotting."))
}
method <- attr(x, "network_method", exact = TRUE)
p <- ggplot2::ggplot(d, ggplot2::aes(x = TIME)) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Time",
y = "Dendrometer value",
title = paste("Network interpolation:", method)
)
if (show_pi && any(is.finite(d$pi_lo) & is.finite(d$pi_hi))) {
p <- p +
ggplot2::geom_ribbon(
ggplot2::aes(ymin = pi_lo, ymax = pi_hi),
fill = "grey70",
alpha = 0.30,
na.rm = TRUE
)
}
p <- p +
ggplot2::geom_line(
ggplot2::aes(y = original),
linewidth = 0.35,
colour = "grey45",
na.rm = TRUE
) +
ggplot2::geom_line(
ggplot2::aes(y = value),
linewidth = 0.45,
colour = "#1F78B4",
na.rm = TRUE
)
if (show_fit && any(is.finite(d$fit))) {
p <- p +
ggplot2::geom_line(
ggplot2::aes(y = fit),
linewidth = 0.35,
linetype = 2,
colour = "black",
na.rm = TRUE
)
}
if (any(d$interp, na.rm = TRUE)) {
p <- p +
ggplot2::geom_point(
data = d[d$interp, , drop = FALSE],
ggplot2::aes(x = TIME, y = value),
colour = "#D7301F",
size = 1.2,
inherit.aes = FALSE,
na.rm = TRUE
)
}
jump_df <- d[!is.na(d$jump_correction), , drop = FALSE]
if (nrow(jump_df) > 0) {
p <- p +
ggplot2::geom_point(
data = jump_df,
ggplot2::aes(x = TIME, y = value),
colour = "darkgreen",
size = 2,
inherit.aes = FALSE,
na.rm = TRUE
)
}
p +
ggplot2::facet_wrap(
~series,
scales = if (free_y) "free_y" else "fixed",
ncol = ncol
)
}
#' @keywords internal
#' @noRd
.plot_network_uncertainty <- function(x, series = NULL, start = NULL, end = NULL,
free_y = TRUE, ncol = 1,
empty = c("plot", "error")) {
empty <- match.arg(empty)
TIME <- pi_width <- NULL
d <- .network_interp_diag_subset(x, series = series, start = start, end = end)
d <- d[is.finite(d$pi_width), , drop = FALSE]
if (nrow(d) == 0) {
if (empty == "error") stop("No prediction intervals available for plotting.")
return(.network_interp_empty_plot("No prediction intervals available for plotting."))
}
ggplot2::ggplot(d, ggplot2::aes(x = TIME, y = pi_width)) +
ggplot2::geom_line(linewidth = 0.45, na.rm = TRUE) +
ggplot2::geom_point(
data = d[d$interp, , drop = FALSE],
size = 1.1,
na.rm = TRUE
) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Time",
y = "95% interval width",
title = "Interpolation uncertainty over time"
) +
ggplot2::facet_wrap(
~series,
scales = if (free_y) "free_y" else "fixed",
ncol = ncol
)
}
#' @keywords internal
#' @noRd
.plot_network_availability <- function(x, series = NULL, start = NULL, end = NULL,
free_y = TRUE, ncol = 1,
empty = c("plot", "error")) {
TIME <- ref_n <- NULL
empty <- match.arg(empty)
d <- .network_interp_diag_subset(x, series = series, start = start, end = end)
if (nrow(d) == 0 || !any(is.finite(d$ref_n))) {
if (empty == "error") stop("No reference availability data available for plotting.")
return(.network_interp_empty_plot("No reference availability data available for plotting."))
}
ggplot2::ggplot(d, ggplot2::aes(x = TIME, y = ref_n)) +
ggplot2::geom_line(linewidth = 0.45, na.rm = TRUE) +
ggplot2::geom_point(
data = d[d$interp, , drop = FALSE],
colour = "#D7301F",
size = 1.2,
na.rm = TRUE
) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Time",
y = "Available reference sensors",
title = "Reference network availability"
) +
ggplot2::facet_wrap(
~series,
scales = if (free_y) "free_y" else "fixed",
ncol = ncol
)
}
#' @keywords internal
#' @noRd
.plot_network_summary <- function(x,
metric = c("n_imputed", "n_remaining_na", "n_missing_input", "mean_pi_width", "median_pi_width", "mean_ref_n", "median_ref_n", "n_gap_jumps_removed", "max_abs_gap_jump"),
empty = c("plot", "error")) {
metric <- match.arg(metric)
empty <- match.arg(empty)
series <- value <- NULL
s <- attr(x, "network_summary", exact = TRUE)
if (is.null(s) || nrow(s) == 0) {
if (empty == "error") stop("No summary information available for plotting.")
return(.network_interp_empty_plot("No summary information available for plotting."))
}
s <- tibble::as_tibble(s)
s$value <- s[[metric]]
ggplot2::ggplot(s, ggplot2::aes(x = series, y = value)) +
ggplot2::geom_col() +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Series",
y = metric,
title = paste("Network interpolation summary:", metric)
)
}
#' @keywords internal
#' @noRd
.plot_network_validation <- function(x,
metric = c("MAE", "MAPE", "MdAPE", "RMSE", "Bias", "Success_rate", "Mean_PI_width", "Mean_ref_n", "End_value_abs_diff", "Max_abs_diff_within_gap"),
series = NULL,
facet = TRUE,
empty = c("plot", "error")) {
metric <- match.arg(metric)
empty <- match.arg(empty)
gap_steps <- value <- NULL
s <- .network_interp_validation_subset(x, series = series)
if (nrow(s) == 0) {
if (empty == "error") stop("No validation results available for plotting.")
return(.network_interp_empty_plot("No validation results available for plotting."))
}
s$value <- s[[metric]]
p <- ggplot2::ggplot(
s,
ggplot2::aes(x = gap_steps, y = value, group = series, colour = series)
) +
ggplot2::geom_line(linewidth = 0.5, na.rm = TRUE) +
ggplot2::geom_point(size = 1.8, na.rm = TRUE) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Artificial gap length (steps)",
y = metric,
colour = "Series",
title = paste("Validation:", metric)
)
if (facet) {
p + ggplot2::facet_wrap(~series, scales = "free_y")
} else {
p
}
}
#' @keywords internal
#' @noRd
.plot_network_coverage <- function(x,
series = NULL,
facet = TRUE,
empty = c("plot", "error")) {
gap_steps <- PI_coverage_95 <- NULL
empty <- match.arg(empty)
s <- .network_interp_validation_subset(x, series = series)
if (nrow(s) == 0 || !any(is.finite(s$PI_coverage_95))) {
if (empty == "error") stop("No validation coverage results available for plotting.")
return(.network_interp_empty_plot("No validation coverage results available for plotting."))
}
p <- ggplot2::ggplot(
s,
ggplot2::aes(x = gap_steps, y = PI_coverage_95, group = series, colour = series)
) +
ggplot2::geom_hline(yintercept = 0.95, linetype = 2) +
ggplot2::geom_line(linewidth = 0.5, na.rm = TRUE) +
ggplot2::geom_point(size = 1.8, na.rm = TRUE) +
ggplot2::coord_cartesian(ylim = c(0, 1)) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Artificial gap length (steps)",
y = "95% PI coverage",
colour = "Series",
title = "Validation of interval coverage"
)
if (facet) {
p + ggplot2::facet_wrap(~series, scales = "fixed")
} else {
p
}
}
#' @keywords internal
#' @noRd
.plot_network_compare <- function(x,
series = NULL,
facet = TRUE,
colour_by = c("series", "gap_steps"),
empty = c("plot", "error")) {
truth <- pred <- gap_steps <- gap_steps_f <- NULL
empty <- match.arg(empty)
colour_by <- match.arg(colour_by)
pdat <- .network_interp_validation_points_subset(x, series = series)
pdat <- pdat[pdat$filled & is.finite(pdat$truth) & is.finite(pdat$pred), , drop = FALSE]
if (nrow(pdat) == 0) {
if (empty == "error") stop("No point-level validation comparison data available for plotting.")
return(.network_interp_empty_plot("No point-level validation comparison data available for plotting."))
}
pdat$gap_steps_f <- factor(pdat$gap_steps)
aes_colour <- if (colour_by == "series") {
ggplot2::aes(colour = series)
} else {
ggplot2::aes(colour = gap_steps_f)
}
p <- ggplot2::ggplot(pdat, ggplot2::aes(x = truth, y = pred)) +
ggplot2::geom_abline(slope = 1, intercept = 0, linetype = 2) +
ggplot2::geom_point(aes_colour, alpha = 0.75, size = 1.8, na.rm = TRUE) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Actual value",
y = "Interpolated value",
colour = if (colour_by == "series") "Series" else "Gap steps",
title = "Actual vs interpolated values from validation"
)
if (facet) {
p + ggplot2::facet_wrap(~series, scales = "free")
} else {
p
}
}
#' @keywords internal
#' @noRd
.plot_network_seasonal_error <- function(x,
metric = c("mean_error", "mean_abs_error", "mean_abs_pct_error", "rmse", "coverage"),
series = NULL,
facet = TRUE,
empty = c("plot", "error")) {
metric <- match.arg(metric)
empty <- match.arg(empty)
doy <- value <- NULL
s <- .network_interp_validation_seasonal_subset(x, series = series)
if (nrow(s) == 0) {
if (empty == "error") stop("No seasonal validation results available for plotting.")
return(.network_interp_empty_plot("No seasonal validation results available for plotting."))
}
s$value <- s[[metric]]
p <- ggplot2::ggplot(
s,
ggplot2::aes(x = doy, y = value, colour = series, group = series)
) +
ggplot2::geom_line(linewidth = 0.45, na.rm = TRUE) +
ggplot2::geom_point(size = 1.2, na.rm = TRUE) +
ggplot2::theme_bw() +
ggplot2::labs(
x = "Day of year",
y = metric,
colour = "Series",
title = paste("Seasonal validation diagnostic:", metric)
)
if (facet) {
p + ggplot2::facet_wrap(~series, scales = "free_y")
} else {
p
}
}
#' @keywords internal
#' @noRd
.network_interp_diag_subset <- function(x, series = NULL, start = NULL, end = NULL) {
d <- attr(x, "network_diagnostics", exact = TRUE)
if (is.null(d)) {
stop("No diagnostics found on this object.")
}
d <- tibble::as_tibble(d)
if (!is.null(series)) {
d <- d[d$series %in% series, , drop = FALSE]
}
if (!is.null(start)) {
d <- d[d$TIME >= as.POSIXct(start), , drop = FALSE]
}
if (!is.null(end)) {
d <- d[d$TIME <= as.POSIXct(end), , drop = FALSE]
}
d
}
#' @keywords internal
#' @noRd
.network_interp_validation_subset <- function(x, series = NULL) {
s <- attr(x, "network_validation_summary", exact = TRUE)
if (is.null(s)) {
return(tibble::tibble())
}
s <- tibble::as_tibble(s)
if (!is.null(series)) {
s <- s[s$series %in% series, , drop = FALSE]
}
s
}
#' @keywords internal
#' @noRd
.network_interp_validation_points_subset <- function(x, series = NULL) {
s <- attr(x, "network_validation_points", exact = TRUE)
if (is.null(s)) {
return(tibble::tibble())
}
s <- tibble::as_tibble(s)
if (!is.null(series)) {
s <- s[s$series %in% series, , drop = FALSE]
}
s
}
#' @keywords internal
#' @noRd
.network_interp_validation_seasonal_subset <- function(x, series = NULL) {
s <- attr(x, "network_validation_seasonal", exact = TRUE)
if (is.null(s)) {
return(tibble::tibble())
}
s <- tibble::as_tibble(s)
if (!is.null(series)) {
s <- s[s$series %in% series, , drop = FALSE]
}
s
}
#' @keywords internal
#' @noRd
.network_interp_empty_plot <- function(label) {
ggplot2::ggplot() +
ggplot2::annotate("text", x = 0, y = 0, label = label, size = 5) +
ggplot2::xlim(-1, 1) +
ggplot2::ylim(-1, 1) +
ggplot2::labs(x = NULL, y = NULL) +
ggplot2::theme_void()
}
Try the dendRoAnalyst package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.