Nothing
if (getRversion() >= "2.15.1") utils::globalVariables("Group")
# TABTS_PUBLICATION_READY_2026_09_02
.r4vn_tabts_revision <- "2026-09-02-v6"
#' Comprehensive time-series analysis with publication-ready output
#'
#' `tabts()` provides a single R4VN-style interface for descriptive time-series
#' analysis, decomposition, stationarity assessment, ACF/PACF, ARIMA/SARIMA/
#' ARIMAX, ETS, model comparison, validation, forecasting, interrupted time
#' series (ITS), controlled ITS, count ITS, residual diagnostics, flexible
#' graphics, and evidence-linked interpretation.
#'
#' The function is intentionally simple for routine use:
#'
#' `tabts(cases, time = month)`
#'
#' while advanced behavior can be changed through `model_options`,
#' `diagnostic_options`, `forecast_options`, `its_options`, `table_options`,
#' `interpret_options`, and `plot_options`. This design keeps the public API
#' stable while allowing future extensions.
#'
#' @param outcome Outcome variable. A bare column name or a single character
#' column name. The outcome must be numeric.
#' @param time Time variable. A bare column name or a single character column
#' name. `Date`, `POSIXct`, integer/numeric index, and ordered time values are
#' supported. Date-like variables are recommended.
#' @param data Data frame. If `NULL`, `tabts()` attempts to use the active R4VN
#' data frame.
#' @param model Analysis engine. One of `"auto"`, `"arima"`, `"ets"`,
#' `"compare"`, `"its"`, or `"regression"`. A vector such as
#' `c("arima","ets")` is treated as a candidate-model comparison.
#' @param period Seasonal period. Use `"auto"` to infer common periods from the
#' time variable, or supply a positive integer such as 12 for monthly annual
#' seasonality, 4 for quarterly data, or 52 for weekly data.
#' @param order ARIMA order `c(p,d,q)`. If `NULL`, `forecast::auto.arima()` is
#' used when `forecast` is installed; otherwise `tabts()` performs a compact
#' AICc-based search with `stats::arima()` so the default command still works.
#' @param seasonal Seasonal ARIMA order `c(P,D,Q)`. The seasonal period comes
#' from `period`. Ignored for non-ARIMA engines.
#' @param xreg Optional external regressors. Accepts `vars(x1, x2)`, `c("x1",
#' "x2")`, or a single bare column name.
#' @param group Optional grouping variable. For ordinary time-series models,
#' separate models are fitted for each group. For ITS, a group variable
#' produces a controlled ITS when `its_options$controlled = TRUE` (the
#' default when a group is supplied).
#' @param intervention Intervention time(s) for ITS. May be a Date/POSIXct,
#' a value comparable with `time`, a vector of intervention times, or the
#' name/bare name of a 0/1 intervention indicator column.
#' @param population Optional population/exposure variable for count ITS.
#' When supplied with a count family, `log(population)` is used as an offset.
#' @param rate Display rate multiplier used for descriptive rate calculations,
#' e.g. 100000. It does not change the count-model offset.
#' @param family ITS family: `"auto"`, `"gaussian"`, `"poisson"`,
#' `"quasipoisson"`, or `"negativebinomial"`.
#' @param correlation ITS residual correlation handling: `"auto"`, `"none"`,
#' `"ar1"`, or `"nw"`. AR(1) via `nlme::gls()` is available for Gaussian ITS.
#' Newey-West robust covariance requires `sandwich`.
#' @param decompose Decomposition: `"auto"`, `"stl"`, `"classical"`, or
#' `"none"`. `"auto"` uses STL when at least two seasonal cycles are present.
#' @param stationarity Logical; run stationarity assessment when possible.
#' ADF and KPSS require `tseries`; suggested differencing additionally uses
#' `forecast` when available.
#' @param acf,pacf Logical; calculate ACF/PACF tables and plots.
#' @param diagnostic Logical; calculate residual diagnostics.
#' @param forecast Number of future periods to forecast. Zero disables
#' forecasting. Forecast intervals are prediction intervals for ARIMA/ETS.
#' @param level Forecast interval levels, e.g. `c(.80,.95)`.
#' @param future_xreg Optional data frame/matrix of future xreg values for
#' ARIMAX or time-series regression forecasts. It must contain at least
#' `forecast` rows and the same xreg variables used for fitting.
#' @param test Optional holdout size for validation. An integer means the
#' number of final observations; a value between 0 and 1 means that
#' proportion of observations.
#' @param criterion Model-selection criterion: `"aicc"`, `"rmse"`, `"mae"`,
#' or `"mape"`. Out-of-sample criteria require `test`.
#' @param missing_time Handling of missing time points: `"warn"`, `"error"`,
#' `"NA"`, `"zero"`, or `"interpolate"`. Missing time points are never
#' silently converted to zero.
#' @param duplicate_time Handling of duplicate time values within a series:
#' `"error"`, `"mean"`, `"sum"`, or `"first"`.
#' @param plot Logical; create plots. Publication-ready base-R plots are always
#' available; if `ggplot2` is installed, `tabts()` uses ggplot2 automatically.
#' @param plots Character vector selecting plots. `"auto"` creates the
#' relevant set; `"all"` is an explicit synonym and `"none"` suppresses
#' plot creation. Other values include `"series"`, `"decomposition"`,
#' `"acf"`, `"pacf"`, `"residual"`, `"forecast"`, `"its"`, and
#' `"counterfactual"`.
#' @param theme Plot theme: `"r4vn"`, `"minimal"`, `"classic"`, or `"bw"`.
#' @param title,subtitle,xlab,ylab Common plot labels.
#' @param legend Logical; show legends where relevant.
#' @param legend_position Legend position such as `"bottom"`, `"top"`,
#' `"left"`, `"right"`, or `"none"`.
#' @param observed_color,fitted_color,forecast_color,counterfactual_color Optional common layer colors. Leave `NULL` to use R4VN defaults.
#' @param intervention_color Optional intervention-layer color. Leave `NULL` to use the R4VN default.
#' @param linewidth Default line width for main series layers.
#' @param point Logical; display observed points on the main series plot.
#' @param point_size Default observed point size.
#' @param pi Logical; display prediction-interval ribbons on forecast plots.
#' @param pi_alpha Prediction-interval ribbon transparency.
#' @param width,height,dpi Default figure width/height in inches and raster
#' resolution used when `plot(result, file = ...)` saves a figure.
#' @param interpret Logical or one of `"brief"`/`"full"`. Generate
#' deterministic rule-based interpretation linked to exact numerical
#' evidence from the output. The default is `FALSE`, keeping the routine
#' report concise; request `TRUE`, `"brief"`, or `"full"` when needed.
#' @param language Output language. Currently only `"en"` is supported; all tables, plots, diagnostics, warnings, and interpretations are produced in English.
#' @param detail Interpretation detail: `"full"` or `"brief"`.
#' @param show Logical; show a publication-style HTML report. In RStudio the
#' report opens in the Viewer and contains all tables and every generated
#' plot; the primary plot is also sent to the Plots pane. No HTML package is
#' required for this Viewer report.
#' @param console Logical; also print tables to the console.
#' @param digits Number of digits for estimates.
#' @param p_digits Number of digits for p-values.
#' @param model_options Named list of advanced model controls. Important
#' entries include `auto_arima`, `ets`, `diagnostic_gate`, `candidate`,
#' `include_drift`, and `include_mean`. For fixed ARIMA/SARIMA models,
#' `tabts()` automatically disables a drift term when `d + D != 1` and
#' disables a mean term when differencing is present.
#' @param diagnostic_options Named list controlling diagnostics. Important
#' entries include `ljung_lag`, `normality`, `acf_lag`, and `alpha`.
#' @param forecast_options Named list controlling forecasts. Entries may
#' include `bootstrap`, `biasadj`, and `history`.
#' @param its_options Named list controlling ITS. Entries include
#' `controlled`, `reference_group`, `time_scale`, `effect_at`,
#' `counterfactual`, `include_season`, `seasonal_harmonics`, `nw_lag`,
#' and `post_min`. Seasonal adjustment uses sine/cosine Fourier pairs;
#' `seasonal_harmonics` controls how many pairs are included. For controlled
#' ITS, the first factor level is the reference unless `reference_group` is
#' supplied.
#' @param table_options Named list controlling output tables. Entries include
#' `ci_level`, `show_model_selection`, `show_stationarity`,
#' `show_diagnostics`, and `show_interpretation`.
#' @param interpret_options Named list controlling interpretation, including
#' `alpha`, `include_assumptions`, `include_model`, `include_forecast`,
#' `include_limitations`, and `evidence`.
#' @param plot_options Deep list controlling plots. This is the main extension
#' point for colors, line types/widths, point shapes/sizes, interval ribbons,
#' axes, date breaks/labels, limits, legends, fonts, grids, reference lines,
#' annotations, intervention/counterfactual layers, facets, and component
#' plots. See Details and examples.
#' @param ai `FALSE`, `TRUE`, or a named R4VN AI endpoint. AI is an optional
#' additional layer and is used only when `aiask()` is available. For a
#' deterministic evidence-linked interpretation, set `interpret = TRUE`.
#' @param ... Reserved for future compatible extensions.
#'
#' @details
#' ## Core workflow
#'
#' `tabts()` follows the R4VN workflow:
#'
#' 1. validate and regularize the time index;
#' 2. describe the series;
#' 3. assess seasonality/stationarity;
#' 4. fit candidate models;
#' 5. validate/select the model;
#' 6. diagnose residuals;
#' 7. forecast when requested;
#' 8. produce publication-ready tables and plots;
#' 9. generate evidence-linked interpretation at the end.
#'
#' ## Optional packages and dependency-light defaults
#'
#' Routine ARIMA/SARIMA, forecasting from fixed/base-selected ARIMA models,
#' regression, decomposition, ACF/PACF, Gaussian ITS, Poisson/quasi-Poisson
#' ITS, Viewer tables, and Viewer figures can run without extra analysis or
#' reporting packages. Optional packages add specialized methods: `forecast`
#' for Hyndman-Khandakar auto-ARIMA and ETS; `tseries` for ADF/KPSS; `nlme`
#' for Gaussian AR(1) ITS; `sandwich` for Newey-West covariance; `MASS` for
#' negative-binomial ITS; and `ggplot2` for editable ggplot objects.
#'
#' ## Flexible plot options
#'
#' Advanced graphics are changed with nested `plot_options`. For example:
#'
#' ```
#' plot_options = list(
#' observed = list(color = "black", linewidth = .8,
#' linetype = "solid", point = TRUE,
#' point_shape = 16, point_size = 2),
#' fitted = list(color = "steelblue", linewidth = 1,
#' linetype = "dashed"),
#' forecast = list(color = "firebrick", linewidth = 1.1),
#' pi = list(show = TRUE, alpha = .15, border = FALSE),
#' axis = list(
#' xlim = NULL, ylim = NULL,
#' date_breaks = "6 months", date_labels = "%b %Y",
#' y_breaks = NULL, y_log = FALSE
#' ),
#' legend = list(position = "bottom", title = NULL),
#' grid = list(major = TRUE, minor = FALSE),
#' intervention = list(line = TRUE, label = TRUE,
#' linetype = "dashed", linewidth = .8),
#' counterfactual = list(show = TRUE, linetype = "dotted"),
#' reference = list(xline = NULL, yline = NULL),
#' facet = list(show = TRUE, ncol = NULL, scales = "fixed"),
#' font = list(family = NULL, base_size = 11),
#' annotation = NULL
#' )
#' ```
#'
#' When `ggplot2` is installed, plots are returned as ordinary `ggplot`
#' objects and may be edited with standard ggplot2 syntax. Without ggplot2,
#' `tabts()` returns lightweight `r4vn_tabts_plot` objects drawn with base R;
#' this keeps the default analysis and Viewer graphics dependency-light.
#'
#' ## Interpretation
#'
#' Every rule-based interpretation row contains a `finding`, the exact
#' `evidence` used to create it, a `status`, and a `source`. Thus statements
#' about stationarity, model selection, residual adequacy, ITS effects, or
#' forecast uncertainty can always be traced back to a specific result.
#'
#' @return An object of class `r4vn_tabts`. Important components include
#' `data`, `summary`, `stationarity`, `decomposition`, `acf`, `pacf`,
#' `model`, `models`, `model_info`, `model_selection`, `coefficients`,
#' `diagnostics`, `validation`, `forecast`, `counterfactual`, `its`,
#' `interpretation`, `tables`, `plots`, and `metadata`. The `its` component
#' stores the family, correlation
#' structure, intervention timing, and pre/post counts used by ITS. Grouped
#' analyses return class `r4vn_tabts_grouped`.
#'
#' @examples
#' \donttest{
#' data(dengue_ts)
#'
#' # 1. Simplest command. Automatic ARIMA works with base R; when forecast is
#' # installed, ARIMA/ETS comparison becomes available automatically.
#' m1 <- tabts(cases, time = month, data = dengue_ts)
#' m1$tables
#' m1$plots$series
#'
#' # 2. Forecast six future months with 80% and 95% prediction intervals.
#' m2 <- tabts(cases, time = month, data = dengue_ts, forecast = 6)
#' m2$forecast
#' plot(m2, "forecast")
#'
#' # 3. Fixed ARIMA using base R only.
#' m3 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1), forecast = 6)
#' m3$coefficients
#' m3$diagnostics
#'
#' # 4. Seasonal ARIMA/SARIMA.
#' m4 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 1, 1),
#' seasonal = c(0, 1, 1), period = 12, forecast = 12)
#'
#' # 5. ARIMAX with external regressors.
#' future_weather <- tail(dengue_ts[c("rainfall", "temperature")], 6)
#' m5 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1),
#' xreg = vars(rainfall, temperature), forecast = 6,
#' future_xreg = future_weather)
#'
#' # 6. Time-series regression with trend and seasonal terms; base R only.
#' m6 <- tabts(cases, time = month, data = dengue_ts,
#' model = "regression", forecast = 6)
#'
#' # 7. Compare ARIMA and ETS when forecast is installed.
#' if (requireNamespace("forecast", quietly = TRUE)) {
#' m7 <- tabts(cases, time = month, data = dengue_ts,
#' model = c("arima", "ets"), test = 12,
#' criterion = "rmse", forecast = 12)
#' m7$model_selection
#' m7$validation
#' }
#'
#' # 8. Decomposition, ACF and PACF are available in the result.
#' m8 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1))
#' m8$decomposition
#' m8$acf
#' m8$pacf
#' plot(m8, "decomposition")
#' plot(m8, "acf")
#' plot(m8, "pacf")
#'
#' # 9. Select only the plots needed in a report.
#' m9 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1), forecast = 6,
#' plots = c("series", "forecast", "residual"))
#'
#' # 10. Missing time points are never silently converted to zero.
#' dmiss <- dengue_ts[-20, ]
#' m10 <- tabts(cases, time = month, data = dmiss,
#' missing_time = "interpolate", model = "arima",
#' order = c(1, 0, 1), forecast = 6)
#'
#' # 11. Duplicate time points can be handled explicitly.
#' ddup <- rbind(dengue_ts, dengue_ts[1, ])
#' m11 <- tabts(cases, time = month, data = ddup,
#' duplicate_time = "mean", model = "arima",
#' order = c(1, 0, 1))
#'
#' # 12. Gaussian interrupted time series using base R.
#' m12 <- tabts(cases, time = month, data = dengue_ts,
#' model = "its", intervention = as.Date("2023-01-01"),
#' family = "gaussian", correlation = "none")
#' m12$coefficients
#' m12$effect_at
#' plot(m12, "counterfactual")
#'
#' # 13. Poisson ITS with a population offset; also base R.
#' m13 <- tabts(cases, time = month, data = dengue_ts,
#' model = "its", intervention = as.Date("2023-01-01"),
#' family = "poisson", population = population,
#' rate = 100000, correlation = "none")
#'
#' # 14. Request effects at clinically meaningful post-intervention times.
#' m14 <- tabts(cases, time = month, data = dengue_ts,
#' model = "its", intervention = as.Date("2023-01-01"),
#' family = "poisson", population = population,
#' correlation = "none",
#' its_options = list(effect_at = c(1, 3, 6, 12, 24)))
#' m14$effect_at
#'
#' # 15. Intervention may be supplied as a 0/1 indicator column.
#' dind <- dengue_ts
#' dind$policy <- as.integer(dind$month >= as.Date("2023-01-01"))
#' m15 <- tabts(cases, time = month, data = dind, model = "its",
#' intervention = policy, family = "poisson",
#' population = population, correlation = "none")
#'
#' # 16. Controlled ITS.
#' data(dengue_its_control)
#' m16 <- tabts(cases, time = month, group = group,
#' data = dengue_its_control, model = "its",
#' intervention = as.Date("2023-01-01"), family = "poisson",
#' population = population, correlation = "none",
#' its_options = list(reference_group = "Control"))
#' m16$coefficients
#' m16$effect_at
#'
#' # 17. Fit separate ordinary time-series models by group.
#' m17 <- tabts(cases, time = month, group = group,
#' data = dengue_its_control, model = "arima",
#' order = c(1, 0, 1), forecast = 3)
#' names(m17$results)
#'
#' # 18. Interpretation is OFF by default.
#' m18 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1))
#' m18$interpretation
#'
#' # 19. Turn on evidence-linked interpretation when desired.
#' m19 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1),
#' forecast = 6, interpret = TRUE)
#' m19$interpretation
#'
#' # 20. Brief interpretation.
#' m20 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1),
#' interpret = "brief")
#'
#' # 21. Publication-ready plot customization.
#' m21 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1), forecast = 12,
#' title = "Monthly dengue cases",
#' subtitle = "Observed, fitted and forecast values",
#' xlab = "Month", ylab = "Cases",
#' plot_options = list(
#' observed = list(point = TRUE, point_size = 1.6),
#' forecast = list(linewidth = 1),
#' pi = list(show = TRUE, alpha = .12),
#' axis = list(date_breaks = "1 year", date_labels = "%Y"),
#' legend = list(position = "bottom")
#' ))
#'
#' # 22. The Viewer report contains every generated table and plot when show=TRUE.
#' m22 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1),
#' forecast = 6, show = TRUE)
#'
#' # 23. Save a publication figure without another export package.
#' plot(m22, "forecast",
#' file = file.path(tempdir(), "tabts_forecast_300dpi.png"),
#' width = 8, height = 5, dpi = 300)
#'
#' # 24. Use the active R4VN data frame.
#' usedf(dengue_ts, quiet = TRUE)
#' m24 <- tabts(cases, time = month, model = "arima",
#' order = c(1, 0, 1), forecast = 3)
#' usedf(clear = TRUE, quiet = TRUE)
#'
#' # 25. Optional enhancements only when needed:
#' # tseries -> ADF/KPSS; nlme -> Gaussian AR(1) ITS;
#' # sandwich -> Newey-West ITS; MASS -> negative-binomial ITS;
#' # forecast -> auto.arima/ETS; ggplot2 -> editable ggplot objects.
#'
#' # 26. Explicit ETS when forecast is available.
#' if (requireNamespace("forecast", quietly = TRUE)) {
#' m26 <- tabts(cases, time = month, data = dengue_ts,
#' model = "ets", forecast = 6)
#' m26$model_info
#' m26$forecast
#' }
#'
#' # 27. ADF/KPSS stationarity tests are added when tseries is installed.
#' m27 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1),
#' stationarity = TRUE)
#' m27$stationarity
#'
#' # 28. Hold out the final 20% of observations for validation.
#' m28 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1), test = .20)
#' m28$validation
#'
#' # 29. Keep tables but suppress plot creation completely.
#' m29 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1),
#' plot = FALSE, show = TRUE)
#' m29$tables
#'
#' # 30. Request all relevant figures explicitly.
#' m30 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1),
#' forecast = 6, plots = "all")
#' names(m30$plots)
#'
#' # 31. Gaussian ITS: correlation = "auto" uses AR(1) only when nlme is
#' # available and the correlated model improves AIC sufficiently.
#' m31 <- tabts(cases, time = month, data = dengue_ts,
#' model = "its", intervention = as.Date("2023-01-01"),
#' family = "gaussian", correlation = "auto")
#' m31$its
#'
#' # 32. Newey-West covariance is an optional enhancement via sandwich.
#' m32 <- tabts(cases, time = month, data = dengue_ts,
#' model = "its", intervention = as.Date("2023-01-01"),
#' family = "poisson", population = population,
#' correlation = "nw")
#' m32$coefficients
#'
#' # 33. Automatic count-family choice: Poisson, quasi-Poisson, or
#' # negative-binomial when MASS is available and overdispersion is marked.
#' m33 <- tabts(cases, time = month, data = dengue_ts,
#' model = "its", intervention = as.Date("2023-01-01"),
#' family = "auto", population = population,
#' correlation = "none")
#' m33$its$family
#'
#' # 34. Explicit quasi-Poisson ITS requires no additional package.
#' m34 <- tabts(cases, time = month, data = dengue_ts,
#' model = "its", intervention = as.Date("2023-01-01"),
#' family = "quasipoisson", population = population,
#' correlation = "none")
#' m34$coefficients
#'
#' # 35. Multiple intervention dates in one segmented model.
#' m35 <- tabts(cases, time = month, data = dengue_ts,
#' model = "its",
#' intervention = as.Date(c("2022-01-01", "2023-01-01")),
#' family = "poisson", population = population,
#' correlation = "none")
#' m35$coefficients
#' m35$its
#'
#' # 36. Add dependency-light Fourier seasonal terms to ITS when required.
#' m36 <- tabts(cases, time = month, data = dengue_ts,
#' model = "its", intervention = as.Date("2023-01-01"),
#' family = "poisson", population = population,
#' correlation = "none", period = 12,
#' its_options = list(include_season = TRUE, seasonal_harmonics = 2))
#' m36$coefficients
#'
#' # 37. Customize which tables are shown in the Viewer without deleting the
#' # underlying result components.
#' m37 <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1),
#' table_options = list(show_stationarity = FALSE,
#' show_validation = FALSE))
#' names(m37$tables)
#' m37$stationarity
#'
#' # 38. If ggplot2 is installed, edit a returned plot as an ordinary ggplot.
#' if (requireNamespace("ggplot2", quietly = TRUE)) {
#' p38 <- m22$plots$forecast + ggplot2::labs(caption = "R4VN tabts")
#' print(p38)
#' }
#'
#' # 39. Save TIFF or vector PDF directly through plot().
#' plot(m22, "forecast",
#' file = file.path(tempdir(), "tabts_forecast.tiff"),
#' width = 8, height = 5, dpi = 300)
#' plot(m22, "forecast",
#' file = file.path(tempdir(), "tabts_forecast.pdf"),
#' width = 8, height = 5)
#'
#' # 40. Inspect reusable components for a custom manuscript/report workflow.
#' names(m22)
#' names(m22$tables)
#' names(m22$plots)
#' m22$metadata
#' summary(m22)
#' }
#'
#' @export
tabts <- function(
outcome,
time,
data = NULL,
model = "auto",
period = "auto",
order = NULL,
seasonal = NULL,
xreg = NULL,
group = NULL,
intervention = NULL,
population = NULL,
rate = 100000,
family = c("auto", "gaussian", "poisson", "quasipoisson", "negativebinomial"),
correlation = c("auto", "none", "ar1", "nw"),
decompose = c("auto", "stl", "classical", "none"),
stationarity = TRUE,
acf = TRUE,
pacf = TRUE,
diagnostic = TRUE,
forecast = 0,
level = c(.80, .95),
future_xreg = NULL,
test = NULL,
criterion = c("aicc", "rmse", "mae", "mape"),
missing_time = c("warn", "error", "NA", "zero", "interpolate"),
duplicate_time = c("error", "mean", "sum", "first"),
plot = TRUE,
plots = "auto",
theme = c("r4vn", "minimal", "classic", "bw"),
title = NULL,
subtitle = NULL,
xlab = NULL,
ylab = NULL,
legend = TRUE,
legend_position = "bottom",
observed_color = NULL,
fitted_color = NULL,
forecast_color = NULL,
counterfactual_color = NULL,
intervention_color = NULL,
linewidth = .8,
point = TRUE,
point_size = 1.8,
pi = TRUE,
pi_alpha = .15,
width = 8,
height = 5,
dpi = 300,
interpret = FALSE,
language = "en",
detail = c("full", "brief"),
show = TRUE,
console = FALSE,
digits = 2,
p_digits = 3,
model_options = list(),
diagnostic_options = list(),
forecast_options = list(),
its_options = list(),
table_options = list(),
interpret_options = list(),
plot_options = list(),
ai = FALSE,
...) {
env <- parent.frame()
outcome_expr <- substitute(outcome)
time_expr <- substitute(time)
xreg_expr <- substitute(xreg)
group_expr <- substitute(group)
pop_expr <- substitute(population)
int_expr <- substitute(intervention)
data <- .r4vn_ts_resolve_data(data, env)
outcome_name <- .r4vn_ts_one_name(outcome_expr, data, env, "outcome")
time_name <- .r4vn_ts_one_name(time_expr, data, env, "time")
xreg_names <- .r4vn_ts_many_names(xreg_expr, data, env)
group_name <- .r4vn_ts_optional_name(group_expr, data, env)
pop_name <- .r4vn_ts_optional_name(pop_expr, data, env)
if (!is.numeric(data[[outcome_name]])) {
stop("`outcome` must be numeric.", call. = FALSE)
}
family <- match.arg(family)
correlation <- match.arg(correlation)
decompose <- match.arg(decompose)
criterion <- match.arg(criterion)
missing_time <- match.arg(missing_time)
duplicate_time <- match.arg(duplicate_time)
theme <- match.arg(theme)
language <- tolower(as.character(language)[1L])
if (!identical(language, "en")) {
stop("`language` currently supports only \"en\".", call. = FALSE)
}
detail <- match.arg(detail)
if (is.character(interpret) && length(interpret) == 1L) {
if (interpret %in% c("brief", "full")) {
detail <- interpret
interpret <- TRUE
}
}
interpret <- isTRUE(interpret)
if (!is.numeric(forecast) || length(forecast) != 1L || is.na(forecast) || forecast < 0) {
stop("`forecast` must be a non-negative number of future periods.", call. = FALSE)
}
forecast <- as.integer(forecast)
level <- sort(unique(as.numeric(level)))
if (length(level) == 0L || any(!is.finite(level)) || any(level <= 0 | level >= 1)) {
stop("`level` must contain probabilities between 0 and 1.", call. = FALSE)
}
defaults <- .r4vn_ts_defaults()
model_options <- .r4vn_ts_merge(defaults$model, model_options)
diagnostic_options <- .r4vn_ts_merge(defaults$diagnostic, diagnostic_options)
forecast_options <- .r4vn_ts_merge(defaults$forecast, forecast_options)
its_options <- .r4vn_ts_merge(defaults$its, its_options)
table_options <- .r4vn_ts_merge(defaults$table, table_options)
interpret_options <- .r4vn_ts_merge(defaults$interpret, interpret_options)
plot_options <- .r4vn_ts_merge(defaults$plot, plot_options)
# `xline` denotes a vertical line on the x axis and `yline` a horizontal
# line on the y axis. Keep the former names as compatibility aliases only.
reference <- plot_options$reference
if (!is.null(reference$xline) && !is.null(reference$vline)) {
stop("Use only `plot_options$reference$xline`; `vline` is its deprecated alias.", call. = FALSE)
}
if (!is.null(reference$yline) && !is.null(reference$hline)) {
stop("Use only `plot_options$reference$yline`; `hline` is its deprecated alias.", call. = FALSE)
}
if (is.null(reference$xline)) reference$xline <- reference$vline
if (is.null(reference$yline)) reference$yline <- reference$hline
if (is.null(reference$xline_label)) reference$xline_label <- reference$vline_label
if (is.null(reference$yline_label)) reference$yline_label <- reference$hline_label
reference[c("vline", "hline", "vline_label", "hline_label")] <- NULL
plot_options$reference <- reference
# Common plot arguments override defaults, but advanced plot_options remain
# the preferred extension point.
plot_options$theme <- theme
plot_options$title <- title %||% plot_options$title
plot_options$subtitle <- subtitle %||% plot_options$subtitle
plot_options$xlab <- xlab %||% plot_options$xlab
plot_options$ylab <- ylab %||% plot_options$ylab
plot_options$legend$show <- isTRUE(legend)
plot_options$legend$position <- legend_position
plot_options$observed$linewidth <- linewidth
plot_options$observed$point <- isTRUE(point)
plot_options$observed$point_size <- point_size
plot_options$fitted$linewidth <- linewidth
plot_options$forecast$linewidth <- linewidth
plot_options$pi$show <- isTRUE(pi)
plot_options$pi$alpha <- pi_alpha
if (!is.null(observed_color)) plot_options$observed$color <- observed_color
if (!is.null(fitted_color)) plot_options$fitted$color <- fitted_color
if (!is.null(forecast_color)) plot_options$forecast$color <- forecast_color
if (!is.null(counterfactual_color)) plot_options$counterfactual$color <- counterfactual_color
if (!is.null(intervention_color)) plot_options$intervention$color <- intervention_color
model <- tolower(as.character(model))
model[model %in% c("sarima", "arimax")] <- "arima"
if (length(model) > 1L) {
model_options$candidate <- unique(model)
model <- "compare"
}
model <- match.arg(model, c("auto", "arima", "ets", "compare", "its", "regression"))
# Resolve intervention only after data have been resolved. A bare column
# name is interpreted as an intervention indicator.
intervention_value <- .r4vn_ts_resolve_intervention(
int_expr, data, env, time_name
)
# For ordinary grouped time series, fit each group separately.
controlled <- isTRUE(its_options$controlled)
if (!is.null(group_name) && model != "its") {
out <- .r4vn_ts_grouped(
data = data,
group_name = group_name,
outcome_name = outcome_name,
time_name = time_name,
xreg_names = xreg_names,
pop_name = pop_name,
model = model,
period = period,
order = order,
seasonal = seasonal,
intervention_value = intervention_value,
rate = rate,
family = family,
correlation = correlation,
decompose = decompose,
stationarity = stationarity,
do_acf = acf,
do_pacf = pacf,
diagnostic = diagnostic,
h = forecast,
level = level,
future_xreg = future_xreg,
test = test,
criterion = criterion,
missing_time = missing_time,
duplicate_time = duplicate_time,
plot = plot,
plots = plots,
interpret = interpret,
language = language,
detail = detail,
digits = digits,
p_digits = p_digits,
model_options = model_options,
diagnostic_options = diagnostic_options,
forecast_options = forecast_options,
its_options = its_options,
table_options = table_options,
interpret_options = interpret_options,
plot_options = plot_options,
export_size = list(width = width, height = height, dpi = dpi)
)
if (console) print(out)
if (show) .r4vn_ts_show_grouped(out)
return(invisible(out))
}
out <- .r4vn_ts_core(
data = data,
outcome_name = outcome_name,
time_name = time_name,
xreg_names = xreg_names,
group_name = if (model == "its" && controlled) group_name else NULL,
pop_name = pop_name,
model = model,
period = period,
order = order,
seasonal = seasonal,
intervention_value = intervention_value,
rate = rate,
family = family,
correlation = correlation,
decompose = decompose,
stationarity = stationarity,
do_acf = acf,
do_pacf = pacf,
diagnostic = diagnostic,
h = forecast,
level = level,
future_xreg = future_xreg,
test = test,
criterion = criterion,
missing_time = missing_time,
duplicate_time = duplicate_time,
plot = plot,
plots = plots,
interpret = interpret,
language = language,
detail = detail,
digits = digits,
p_digits = p_digits,
model_options = model_options,
diagnostic_options = diagnostic_options,
forecast_options = forecast_options,
its_options = its_options,
table_options = table_options,
interpret_options = interpret_options,
plot_options = plot_options,
export_size = list(width = width, height = height, dpi = dpi)
)
out$call <- match.call()
out$ai <- .r4vn_ts_ai(out, ai)
if (console) print(out)
if (show) .r4vn_ts_show(out)
invisible(out)
}
`%||%` <- function(x, y) if (is.null(x)) y else x
.r4vn_ts_defaults <- function() {
list(
model = list(
candidate = c("arima", "ets"),
diagnostic_gate = TRUE,
include_drift = TRUE,
include_mean = TRUE,
approximation = FALSE,
stepwise = TRUE,
auto_arima = list(),
arima = list(),
ets = list(),
regression = list()
),
diagnostic = list(
alpha = .05,
ljung_lag = NULL,
normality = FALSE,
acf_lag = NULL,
residual_standardize = FALSE
),
forecast = list(
bootstrap = FALSE,
biasadj = FALSE,
history = NULL
),
its = list(
controlled = TRUE,
reference_group = NULL,
counterfactual = TRUE,
effect_at = c(1, 3, 6, 12),
nw_lag = NULL,
post_min = 8L,
time_scale = 1,
include_season = FALSE,
seasonal_harmonics = 1L
),
table = list(
ci_level = .95,
show_model_selection = TRUE,
show_stationarity = TRUE,
show_diagnostics = TRUE,
show_interpretation = TRUE,
show_validation = TRUE
),
interpret = list(
alpha = .05,
include_assumptions = TRUE,
include_model = TRUE,
include_forecast = TRUE,
include_limitations = TRUE,
evidence = TRUE
),
plot = list(
theme = "r4vn",
title = NULL,
subtitle = NULL,
xlab = NULL,
ylab = NULL,
observed = list(
color = "#222222", linewidth = .8, linetype = "solid",
alpha = 1, point = TRUE, point_shape = 16, point_size = 1.8,
point_alpha = 1
),
fitted = list(
color = "#2C6EAA", linewidth = .8, linetype = "solid", alpha = .9
),
forecast = list(
color = "#B23A48", linewidth = .8, linetype = "solid", alpha = 1
),
counterfactual = list(
show = TRUE, color = "#666666", linewidth = .8,
linetype = "dashed", alpha = 1, shade = FALSE, shade_alpha = .10
),
intervention = list(
line = TRUE, label = TRUE, color = "#8B0000",
linewidth = .7, linetype = "dashed", alpha = .9,
label_text = "Intervention", label_angle = 90
),
pi = list(
show = TRUE, alpha = .15, border = FALSE, border_linewidth = .3
),
axis = list(
xlim = NULL, ylim = NULL, x_breaks = NULL, y_breaks = NULL,
date_breaks = NULL, date_labels = NULL, y_log = FALSE,
expand = TRUE
),
legend = list(
show = TRUE, position = "bottom", title = NULL
),
grid = list(
major = TRUE, minor = FALSE
),
reference = list(
xline = NULL, yline = NULL, xline_label = NULL, yline_label = NULL,
color = "#777777", linewidth = .5, linetype = "dotted"
),
facet = list(
show = TRUE, ncol = NULL, nrow = NULL, scales = "fixed"
),
font = list(
family = NULL, base_size = 11, title_size = 13,
subtitle_size = 11, axis_title_size = 11,
axis_text_size = 10, legend_text_size = 10
),
annotation = NULL,
decomposition = list(
color = "#222222", linewidth = .6, free_y = TRUE
),
acf = list(
color = "#2C6EAA", linewidth = .7, ci = TRUE, ci_level = .95
),
pacf = list(
color = "#2C6EAA", linewidth = .7, ci = TRUE, ci_level = .95
),
residual = list(
color = "#222222", linewidth = .6, point = TRUE, point_size = 1.5
)
)
)
}
.r4vn_ts_merge <- function(x, y) {
if (is.null(y) || length(y) == 0L) return(x)
if (!is.list(y)) return(y)
if (!is.list(x)) x <- list()
for (nm in names(y)) {
if (!is.null(names(y)[names(y) == nm])) {
if (is.list(y[[nm]]) && is.list(x[[nm]])) {
x[[nm]] <- .r4vn_ts_merge(x[[nm]], y[[nm]])
} else {
x[[nm]] <- y[[nm]]
}
}
}
x
}
.r4vn_ts_resolve_data <- function(data, env) {
if (!is.null(data)) {
if (!is.data.frame(data)) {
stop("`data` must be a data frame.", call. = FALSE)
}
return(data)
}
# Use the canonical R4VN active-data mechanism first. This keeps tabts()
# consistent with tab(), tabscale(), sum1(), describe(), and other R4VN
# analysis commands, including linked active data created by usedf().
canonical_helpers <- c(".r4vn_resolve_analysis_data", ".r4vn_get_active")
for (nm in canonical_helpers) {
if (exists(nm, mode = "function", inherits = TRUE)) {
fun <- get(nm, mode = "function", inherits = TRUE)
ans <- tryCatch(
{
if (identical(nm, ".r4vn_get_active")) {
fun(required = FALSE)
} else {
fun(NULL)
}
},
error = function(e) NULL
)
if (is.data.frame(ans)) return(ans)
}
}
# Public fallback. This is useful if tabts.R is sourced in a development
# session where usedf() is available but internal helpers are not visible.
if (exists("usedf", mode = "function", inherits = TRUE)) {
ans <- tryCatch(
get("usedf", mode = "function", inherits = TRUE)(quiet = TRUE),
error = function(e) NULL
)
if (is.data.frame(ans)) return(ans)
}
# Backward compatibility with older R4VN active-data implementations.
helper_names <- c(".r4vn_active_data", ".get_active_data")
for (nm in helper_names) {
if (exists(nm, mode = "function", inherits = TRUE)) {
ans <- tryCatch(
get(nm, mode = "function", inherits = TRUE)(),
error = function(e) NULL
)
if (is.data.frame(ans)) return(ans)
}
}
opts <- c("R4VN.active.data", "r4vn.active.data", "R4VN_active_data")
for (op in opts) {
ans <- getOption(op, NULL)
if (is.data.frame(ans)) return(ans)
if (is.character(ans) && length(ans) == 1L &&
exists(ans, envir = .GlobalEnv, inherits = FALSE)) {
z <- get(ans, envir = .GlobalEnv, inherits = FALSE)
if (is.data.frame(z)) return(z)
}
}
stop(
"No active data frame. Use `usedf(data)` or ",
"`opendata(..., active = TRUE)` first, or supply `data =`.",
call. = FALSE
)
}
.r4vn_ts_one_name <- function(expr, data, env, label) {
if (is.symbol(expr)) {
nm <- as.character(expr)
if (nm %in% names(data)) return(nm)
}
val <- tryCatch(eval(expr, envir = env), error = function(e) NULL)
if (is.character(val) && length(val) == 1L && val %in% names(data)) return(val)
stop("Could not resolve `", label, "` to a column in `data`.", call. = FALSE)
}
.r4vn_ts_optional_name <- function(expr, data, env) {
if (identical(expr, quote(NULL))) return(NULL)
if (is.symbol(expr)) {
nm <- as.character(expr)
if (nm %in% names(data)) return(nm)
}
val <- tryCatch(eval(expr, envir = env), error = function(e) NULL)
if (is.null(val)) return(NULL)
if (is.character(val) && length(val) == 1L && val %in% names(data)) return(val)
NULL
}
.r4vn_ts_many_names <- function(expr, data, env) {
if (identical(expr, quote(NULL))) return(character())
extract <- function(e) {
if (is.symbol(e)) return(as.character(e))
if (is.call(e)) {
hd <- as.character(e[[1L]])
if (hd %in% c("vars", "c", "list")) {
return(unlist(lapply(as.list(e)[-1L], extract), use.names = FALSE))
}
}
character()
}
nm <- unique(extract(expr))
nm <- nm[nm %in% names(data)]
if (length(nm)) return(nm)
val <- tryCatch(eval(expr, envir = env), error = function(e) NULL)
if (is.character(val)) return(intersect(val, names(data)))
if (is.list(val) && !is.null(names(val))) return(intersect(names(val), names(data)))
character()
}
.r4vn_ts_resolve_intervention <- function(expr, data, env, time_name) {
if (identical(expr, quote(NULL))) return(NULL)
# Bare column: use first 0 -> 1 transition or all rising edges.
if (is.symbol(expr)) {
nm <- as.character(expr)
if (nm %in% names(data)) {
z <- data[[nm]]
if (is.logical(z) || is.numeric(z) || is.integer(z)) {
zz <- as.numeric(z)
idx <- which(zz == 1 & c(TRUE, head(zz, -1L) != 1))
if (length(idx)) return(unique(data[[time_name]][idx]))
}
}
}
val <- tryCatch(eval(expr, envir = env), error = function(e) NULL)
if (is.null(val)) {
stop("Could not resolve `intervention`.", call. = FALSE)
}
val
}
.r4vn_ts_grouped <- function(
data, group_name, outcome_name, time_name, xreg_names, pop_name,
model, period, order, seasonal, intervention_value, rate, family,
correlation, decompose, stationarity, do_acf, do_pacf, diagnostic,
h, level, future_xreg, test, criterion, missing_time, duplicate_time,
plot, plots, interpret, language, detail, digits, p_digits,
model_options, diagnostic_options, forecast_options, its_options,
table_options, interpret_options, plot_options, export_size) {
g <- data[[group_name]]
lev <- unique(as.character(g[!is.na(g)]))
res <- vector("list", length(lev))
names(res) <- lev
for (i in seq_along(lev)) {
d <- data[as.character(data[[group_name]]) == lev[i], , drop = FALSE]
res[[i]] <- .r4vn_ts_core(
data = d,
outcome_name = outcome_name,
time_name = time_name,
xreg_names = xreg_names,
group_name = NULL,
pop_name = pop_name,
model = model,
period = period,
order = order,
seasonal = seasonal,
intervention_value = intervention_value,
rate = rate,
family = family,
correlation = correlation,
decompose = decompose,
stationarity = stationarity,
do_acf = do_acf,
do_pacf = do_pacf,
diagnostic = diagnostic,
h = h,
level = level,
future_xreg = future_xreg,
test = test,
criterion = criterion,
missing_time = missing_time,
duplicate_time = duplicate_time,
plot = plot,
plots = plots,
interpret = interpret,
language = language,
detail = detail,
digits = digits,
p_digits = p_digits,
model_options = model_options,
diagnostic_options = diagnostic_options,
forecast_options = forecast_options,
its_options = its_options,
table_options = table_options,
interpret_options = interpret_options,
plot_options = plot_options,
export_size = export_size
)
}
structure(
list(
group = group_name,
outcome = outcome_name,
time = time_name,
results = res,
metadata = list(
group = group_name,
groups = lev,
n_groups = length(lev),
forecast_horizon = h,
plots = plots,
revision = .r4vn_tabts_revision
),
call = match.call()
),
class = "r4vn_tabts_grouped"
)
}
.r4vn_ts_core <- function(
data, outcome_name, time_name, xreg_names, group_name, pop_name,
model, period, order, seasonal, intervention_value, rate, family,
correlation, decompose, stationarity, do_acf, do_pacf, diagnostic,
h, level, future_xreg, test, criterion, missing_time, duplicate_time,
plot, plots, interpret, language, detail, digits, p_digits,
model_options, diagnostic_options, forecast_options, its_options,
table_options, interpret_options, plot_options, export_size) {
needed <- unique(c(outcome_name, time_name, xreg_names, group_name, pop_name))
d <- data[, needed, drop = FALSE]
names(d)[names(d) == outcome_name] <- ".y"
names(d)[names(d) == time_name] <- ".time_original"
if (!is.null(group_name)) names(d)[names(d) == group_name] <- ".group"
if (!is.null(pop_name)) names(d)[names(d) == pop_name] <- ".population"
d <- d[!is.na(d$.time_original), , drop = FALSE]
d <- d[order(d$.time_original), , drop = FALSE]
if (!is.null(group_name)) {
# Controlled ITS retains repeated times across groups and regularizes each
# group independently below.
spl <- split(d, as.character(d$.group))
spl <- lapply(spl, function(z) {
gv <- as.character(z$.group[which(!is.na(z$.group))[1L]])
zz <- .r4vn_ts_regularize(
z, ".time_original", ".y", xreg_names,
if (!is.null(pop_name)) ".population" else NULL,
missing_time, duplicate_time
)
zz$.group[is.na(zz$.group)] <- gv
zz
})
d <- do.call(rbind, spl)
rownames(d) <- NULL
d <- d[order(d$.time_original, d$.group), , drop = FALSE]
} else {
d <- .r4vn_ts_regularize(
d, ".time_original", ".y", xreg_names,
if (!is.null(pop_name)) ".population" else NULL,
missing_time, duplicate_time
)
}
step <- .r4vn_ts_infer_step(d$.time_original)
frequency <- .r4vn_ts_period(period, step, d$.time_original)
# Descriptive table.
summary_tab <- .r4vn_ts_summary(d$.y, d$.time_original, frequency)
if (!is.null(pop_name) && all(is.finite(d$.population)) && all(d$.population > 0, na.rm = TRUE)) {
rr <- d$.y / d$.population * rate
summary_tab$Rate_mean <- mean(rr, na.rm = TRUE)
summary_tab$Rate_min <- min(rr, na.rm = TRUE)
summary_tab$Rate_max <- max(rr, na.rm = TRUE)
summary_tab$Rate_multiplier <- rate
}
# For non-controlled analysis, create a ts object.
yts <- NULL
if (is.null(group_name)) {
yts <- stats::ts(d$.y, frequency = frequency)
}
st <- if (stationarity && is.null(group_name)) {
.r4vn_ts_stationarity(yts, frequency)
} else NULL
decomp <- NULL
if (is.null(group_name) && decompose != "none") {
decomp <- .r4vn_ts_decompose(yts, d$.time_original, decompose, frequency)
}
acf_tab <- pacf_tab <- NULL
if (is.null(group_name) && do_acf) {
acf_tab <- .r4vn_ts_acf(yts, diagnostic_options$acf_lag, pacf = FALSE)
}
if (is.null(group_name) && do_pacf) {
pacf_tab <- .r4vn_ts_acf(yts, diagnostic_options$acf_lag, pacf = TRUE)
}
fit <- NULL
fits <- list()
model_selection <- NULL
validation <- NULL
coeff <- NULL
diag_tab <- NULL
forecast_tab <- NULL
counterfactual <- NULL
effect_at <- NULL
fitted_values <- rep(NA_real_, nrow(d))
residual_values <- rep(NA_real_, nrow(d))
chosen_model <- model
its_meta <- NULL
if (model == "its") {
its <- .r4vn_ts_fit_its(
d = d,
xreg_names = xreg_names,
intervention_value = intervention_value,
family = family,
correlation = correlation,
frequency = frequency,
group_name = group_name,
pop_name = pop_name,
rate = rate,
h = h,
level = level,
future_xreg = future_xreg,
its_options = its_options,
diagnostic_options = diagnostic_options,
table_options = table_options,
step = step
)
fit <- its$fit
fits <- list(its = fit)
coeff <- its$coefficients
fitted_values <- its$fitted
residual_values <- its$residuals
diag_tab <- if (diagnostic) its$diagnostics else NULL
forecast_tab <- its$forecast
counterfactual <- its$counterfactual
effect_at <- its$effect_at
its_meta <- its$meta
chosen_model <- "its"
} else {
if (anyNA(d$.y)) {
# ARIMA can sometimes accept NA, but ETS/regression generally cannot.
# Keep behavior explicit and predictable.
if (model %in% c("ets", "compare", "regression", "auto")) {
stop(
"The regularized series contains missing outcome values. ",
"Choose `missing_time = \"interpolate\"` or another explicit strategy ",
"before fitting this model.",
call. = FALSE
)
}
}
xmat <- .r4vn_ts_xreg_matrix(d, xreg_names)
candidate <- switch(
model,
auto = model_options$candidate,
compare = model_options$candidate,
arima = "arima",
ets = "ets",
regression = "regression"
)
candidate <- unique(tolower(candidate))
candidate <- candidate[candidate %in% c("arima", "ets", "regression")]
if (length(xreg_names) && "ets" %in% candidate) {
candidate <- setdiff(candidate, "ets")
}
# ETS is an optional enhancement. A routine tabts() call must remain
# usable on a clean R installation, so automatic/compare analyses simply
# omit ETS when forecast is unavailable; an explicitly requested ETS model
# still gives a clear message.
if ("ets" %in% candidate && !requireNamespace("forecast", quietly = TRUE)) {
if (identical(model, "ets")) {
stop("`model = \"ets\"` requires the optional package `forecast`.", call. = FALSE)
}
candidate <- setdiff(candidate, "ets")
}
if (!length(candidate)) candidate <- "arima"
# Holdout validation, when requested.
split_info <- .r4vn_ts_test_split(length(d$.y), test)
validation_rows <- list()
for (eng in candidate) {
f <- .r4vn_ts_fit_engine(
eng, yts, xmat, frequency, order, seasonal, model_options
)
fits[[eng]] <- f
if (!is.null(split_info)) {
ntr <- split_info$n_train
ytr <- stats::ts(d$.y[seq_len(ntr)], frequency = frequency)
xtr <- if (!is.null(xmat)) xmat[seq_len(ntr), , drop = FALSE] else NULL
xtest <- if (!is.null(xmat)) xmat[(ntr + 1L):nrow(xmat), , drop = FALSE] else NULL
ftr <- .r4vn_ts_fit_engine(
eng, ytr, xtr, frequency, order, seasonal, model_options
)
pred <- .r4vn_ts_forecast_engine(
ftr, eng, split_info$n_test, c(.95), xtest, forecast_options
)
met <- .r4vn_ts_accuracy(
actual = d$.y[(ntr + 1L):length(d$.y)],
predicted = pred$mean,
train = d$.y[seq_len(ntr)],
frequency = frequency
)
validation_rows[[eng]] <- cbind(Model = toupper(eng), met, stringsAsFactors = FALSE)
}
}
if (length(validation_rows)) {
validation <- do.call(rbind, validation_rows)
rownames(validation) <- NULL
}
# Model-selection table.
selection_rows <- lapply(names(fits), function(eng) {
f <- fits[[eng]]
dg <- .r4vn_ts_diagnostics(
f,
diagnostic_options = diagnostic_options,
frequency = frequency,
model_df = .r4vn_ts_model_df(f)
)
data.frame(
Model = toupper(eng),
AIC = .r4vn_ts_aic(f),
AICc = .r4vn_ts_aicc(f),
BIC = .r4vn_ts_bic(f),
Ljung_Box_p = if (!is.null(dg)) dg$p[dg$Diagnostic == "Ljung-Box"][1] else NA_real_,
Residual_OK = if (!is.null(dg)) {
p <- dg$p[dg$Diagnostic == "Ljung-Box"][1]
is.finite(p) && p >= diagnostic_options$alpha
} else NA,
stringsAsFactors = FALSE
)
})
model_selection <- do.call(rbind, selection_rows)
if (!is.null(validation)) {
m <- match(model_selection$Model, validation$Model)
for (nm in c("RMSE", "MAE", "MAPE", "sMAPE", "MASE")) {
model_selection[[nm]] <- validation[[nm]][m]
}
}
chosen_model <- .r4vn_ts_choose_model(
model_selection, criterion, model_options$diagnostic_gate, test
)
fit <- fits[[tolower(chosen_model)]]
coeff <- .r4vn_ts_coefficients(
fit, ci_level = table_options$ci_level, family = "gaussian"
)
fitted_values <- .r4vn_ts_fitted(fit, length(d$.y))
residual_values <- .r4vn_ts_residuals(fit, length(d$.y))
if (diagnostic) {
diag_tab <- .r4vn_ts_diagnostics(
fit,
diagnostic_options = diagnostic_options,
frequency = frequency,
model_df = .r4vn_ts_model_df(fit)
)
}
if (h > 0L) {
fx <- .r4vn_ts_future_xreg(future_xreg, xreg_names, h)
fc <- .r4vn_ts_forecast_engine(
fit, tolower(chosen_model), h, level, fx, forecast_options
)
future_time <- .r4vn_ts_extend_time(d$.time_original, h, step)
forecast_tab <- .r4vn_ts_forecast_table(fc, future_time, level)
}
}
series_data <- d
series_data$.fitted <- fitted_values
series_data$.residual <- residual_values
interpretation <- if (interpret) {
.r4vn_ts_interpret(
summary_tab = summary_tab,
stationarity_tab = st,
model_selection = model_selection,
chosen_model = chosen_model,
coefficients = coeff,
diagnostics = diag_tab,
forecast_tab = forecast_tab,
its_meta = its_meta,
effect_at = effect_at,
language = language,
detail = detail,
options = interpret_options,
digits = digits,
p_digits = p_digits
)
} else NULL
requested_plots <- .r4vn_ts_plot_set(plots, model, h, decomp, do_acf, do_pacf)
plot_list <- if (plot) {
plot_args <- list(
d = series_data,
frequency = frequency,
decomposition = decomp,
acf_tab = acf_tab,
pacf_tab = pacf_tab,
forecast_tab = forecast_tab,
counterfactual = counterfactual,
intervention_value = intervention_value,
model = model,
requested = requested_plots,
options = plot_options,
outcome_label = .r4vn_ts_label(data[[outcome_name]], outcome_name),
time_label = .r4vn_ts_label(data[[time_name]], time_name)
)
if (requireNamespace("ggplot2", quietly = TRUE)) {
do.call(.r4vn_ts_make_plots, plot_args)
} else {
do.call(.r4vn_tabts_make_base_plots, plot_args)
}
} else list()
model_info <- .r4vn_tabts_model_info(
fit = fit, model_name = chosen_model, frequency = frequency,
its_meta = its_meta, n = nrow(d)
)
tables <- list(
descriptive = summary_tab,
stationarity = if (isTRUE(table_options$show_stationarity)) st else NULL,
model_selection = if (isTRUE(table_options$show_model_selection)) model_selection else NULL,
model = model_info,
coefficients = coeff,
diagnostics = if (isTRUE(table_options$show_diagnostics)) diag_tab else NULL,
validation = if (isTRUE(table_options$show_validation)) validation else NULL,
forecast = forecast_tab,
effect_at = effect_at,
interpretation = if (isTRUE(table_options$show_interpretation)) interpretation else NULL
)
tables <- tables[!vapply(tables, is.null, logical(1))]
structure(
list(
data = series_data,
outcome = outcome_name,
time = time_name,
xreg = xreg_names,
group = group_name,
period = frequency,
time_step = step,
summary = summary_tab,
stationarity = st,
decomposition = decomp,
acf = acf_tab,
pacf = pacf_tab,
model_name = chosen_model,
model = fit,
models = fits,
model_selection = model_selection,
model_info = model_info,
coefficients = coeff,
diagnostics = diag_tab,
validation = validation,
forecast = forecast_tab,
counterfactual = counterfactual,
effect_at = effect_at,
its = its_meta,
interpretation = interpretation,
tables = tables,
plots = plot_list,
export_size = export_size,
digits = digits,
p_digits = p_digits,
metadata = list(
outcome = outcome_name,
time = time_name,
xreg = xreg_names,
group = group_name,
model = chosen_model,
period = frequency,
time_step = step,
forecast_horizon = h,
interval_levels = level,
plots = names(plot_list),
plot_engine = if (length(plot_list) && inherits(plot_list[[1L]], "ggplot")) "ggplot2" else if (length(plot_list)) "base" else "none",
optional_packages = c(
forecast = requireNamespace("forecast", quietly = TRUE),
tseries = requireNamespace("tseries", quietly = TRUE),
ggplot2 = requireNamespace("ggplot2", quietly = TRUE),
nlme = requireNamespace("nlme", quietly = TRUE),
sandwich = requireNamespace("sandwich", quietly = TRUE),
MASS = requireNamespace("MASS", quietly = TRUE)
)
),
revision = .r4vn_tabts_revision,
options = list(
model = model_options,
diagnostic = diagnostic_options,
forecast = forecast_options,
its = its_options,
table = table_options,
interpret = interpret_options,
plot = plot_options
)
),
class = "r4vn_tabts"
)
}
.r4vn_ts_regularize <- function(
d, time_col, y_col, xreg_names, pop_col, missing_time, duplicate_time) {
tt <- d[[time_col]]
if (anyDuplicated(tt)) {
if (duplicate_time == "error") {
stop(
"Duplicate time values were found. Use `duplicate_time = \"mean\"`, ",
"`\"sum\"`, or `\"first\"` explicitly.",
call. = FALSE
)
}
d <- .r4vn_ts_collapse_duplicates(d, time_col, duplicate_time)
tt <- d[[time_col]]
}
step <- .r4vn_ts_infer_step(tt)
full_time <- .r4vn_ts_full_time(tt, step)
if (is.null(full_time) || length(full_time) == nrow(d)) {
return(d[order(d[[time_col]]), , drop = FALSE])
}
key <- data.frame(.time_key = full_time)
names(key) <- time_col
z <- merge(key, d, by = time_col, all.x = TRUE, sort = TRUE)
missing_rows <- is.na(z[[y_col]])
if (any(missing_rows)) {
if (missing_time == "error") {
stop(
"Missing time points were detected. Choose an explicit `missing_time` ",
"strategy such as \"NA\", \"zero\", or \"interpolate\".",
call. = FALSE
)
}
if (missing_time == "warn") {
warning(
sum(missing_rows),
" missing time point(s) were inserted as NA. ",
"Choose `missing_time = \"interpolate\"`, `\"zero\"`, or another ",
"explicit strategy before using models that cannot handle NA.",
call. = FALSE
)
}
if (missing_time == "zero") {
z[[y_col]][missing_rows] <- 0
}
if (missing_time == "interpolate") {
cols <- unique(c(y_col, xreg_names, pop_col))
cols <- cols[cols %in% names(z)]
for (nm in cols) {
if (is.numeric(z[[nm]])) {
z[[nm]] <- .r4vn_ts_interp(z[[nm]])
}
}
}
}
z
}
.r4vn_ts_collapse_duplicates <- function(d, time_col, how) {
spl <- split(seq_len(nrow(d)), d[[time_col]])
rows <- lapply(spl, function(ii) {
if (how == "first") return(d[ii[1L], , drop = FALSE])
out <- d[ii[1L], , drop = FALSE]
for (nm in names(d)) {
if (nm == time_col) next
v <- d[[nm]][ii]
if (is.numeric(v)) {
out[[nm]] <- if (how == "sum") sum(v, na.rm = TRUE) else mean(v, na.rm = TRUE)
} else {
out[[nm]] <- v[which(!is.na(v))[1L] %||% 1L]
}
}
out
})
out <- do.call(rbind, rows)
rownames(out) <- NULL
out[order(out[[time_col]]), , drop = FALSE]
}
.r4vn_ts_interp <- function(x) {
if (!is.numeric(x) || all(is.na(x))) return(x)
ok <- which(!is.na(x))
if (length(ok) == 1L) {
x[is.na(x)] <- x[ok]
return(x)
}
stats::approx(ok, x[ok], xout = seq_along(x), rule = 2)$y
}
.r4vn_ts_infer_step <- function(t) {
if (inherits(t, "Date")) {
u <- sort(unique(t))
dd <- as.numeric(diff(u))
med <- if (length(dd)) stats::median(dd, na.rm = TRUE) else NA_real_
if (is.finite(med) && med >= 27 && med <= 32)
return(list(type = "Date", unit = "month", by = "month", delta = med, period = 12L))
if (is.finite(med) && med >= 80 && med <= 100)
return(list(type = "Date", unit = "quarter", by = "3 months", delta = med, period = 4L))
if (is.finite(med) && med >= 6 && med <= 8)
return(list(type = "Date", unit = "week", by = "week", delta = med, period = 52L))
if (is.finite(med) && med >= .5 && med <= 1.5)
return(list(type = "Date", unit = "day", by = "day", delta = med, period = 7L))
if (is.finite(med) && med >= 350 && med <= 380)
return(list(type = "Date", unit = "year", by = "year", delta = med, period = 1L))
return(list(type = "Date", unit = "irregular", by = NULL, delta = med, period = 1L))
}
if (inherits(t, "POSIXt")) {
u <- sort(unique(t))
dd <- as.numeric(diff(u), units = "secs")
med <- if (length(dd)) stats::median(dd, na.rm = TRUE) else NA_real_
day <- 86400
if (is.finite(med) && med >= .8 * day && med <= 1.2 * day)
return(list(type = "POSIXct", unit = "day", by = "day", delta = med, period = 7L))
if (is.finite(med) && med >= 6 * day && med <= 8 * day)
return(list(type = "POSIXct", unit = "week", by = "week", delta = med, period = 52L))
return(list(type = "POSIXct", unit = "irregular", by = NULL, delta = med, period = 1L))
}
if (is.numeric(t) || is.integer(t)) {
u <- sort(unique(as.numeric(t)))
dd <- diff(u)
med <- if (length(dd)) stats::median(dd, na.rm = TRUE) else 1
return(list(type = "numeric", unit = "index", by = med, delta = med, period = 1L))
}
list(type = class(t)[1L], unit = "ordered", by = NULL, delta = NA_real_, period = 1L)
}
.r4vn_ts_full_time <- function(t, step) {
if (length(t) < 2L) return(t)
if (inherits(t, "Date") && !is.null(step$by)) {
return(seq.Date(min(t), max(t), by = step$by))
}
if (inherits(t, "POSIXt") && !is.null(step$by)) {
return(seq.POSIXt(min(t), max(t), by = step$by))
}
if ((is.numeric(t) || is.integer(t)) && is.finite(step$by) && step$by > 0) {
return(seq(min(t), max(t), by = step$by))
}
NULL
}
.r4vn_ts_period <- function(period, step, t) {
if (length(period) == 1L && is.character(period) && period == "auto") {
return(max(1L, as.integer(step$period %||% 1L)))
}
p <- suppressWarnings(as.integer(period))
if (length(p) != 1L || is.na(p) || p < 1L) {
stop("`period` must be \"auto\" or a positive integer.", call. = FALSE)
}
p
}
.r4vn_ts_extend_time <- function(t, h, step) {
if (h <= 0L) return(t[0])
last <- max(t)
if (inherits(t, "Date") && !is.null(step$by)) {
return(seq.Date(last, by = step$by, length.out = h + 1L)[-1L])
}
if (inherits(t, "POSIXt") && !is.null(step$by)) {
return(seq.POSIXt(last, by = step$by, length.out = h + 1L)[-1L])
}
if ((is.numeric(t) || is.integer(t)) && is.finite(step$by)) {
return(last + step$by * seq_len(h))
}
seq_len(h) + length(t)
}
.r4vn_ts_summary <- function(y, time, frequency) {
yy <- y[is.finite(y)]
data.frame(
N_time = length(y),
N_observed = length(yy),
Missing = sum(!is.finite(y)),
Start = as.character(min(time, na.rm = TRUE)),
End = as.character(max(time, na.rm = TRUE)),
Period = frequency,
Mean = if (length(yy)) mean(yy) else NA_real_,
SD = if (length(yy) > 1L) stats::sd(yy) else NA_real_,
Median = if (length(yy)) stats::median(yy) else NA_real_,
Q1 = if (length(yy)) as.numeric(stats::quantile(yy, .25, na.rm = TRUE)) else NA_real_,
Q3 = if (length(yy)) as.numeric(stats::quantile(yy, .75, na.rm = TRUE)) else NA_real_,
Min = if (length(yy)) min(yy) else NA_real_,
Max = if (length(yy)) max(yy) else NA_real_,
Total = if (length(yy)) sum(yy) else NA_real_,
stringsAsFactors = FALSE
)
}
.r4vn_ts_stationarity <- function(yts, frequency) {
y <- as.numeric(yts)
y <- y[is.finite(y)]
if (length(y) < 8L) return(NULL)
rows <- list()
if (requireNamespace("tseries", quietly = TRUE)) {
adf <- tryCatch(
suppressWarnings(tseries::adf.test(y, alternative = "stationary")),
error = function(e) NULL
)
if (!is.null(adf)) {
rows[[length(rows) + 1L]] <- data.frame(
Test = "ADF",
Null = "Unit root",
Statistic = unname(adf$statistic),
Lag = unname(adf$parameter),
p = adf$p.value,
stringsAsFactors = FALSE
)
}
kpss <- tryCatch(
suppressWarnings(tseries::kpss.test(y, null = "Level")),
error = function(e) NULL
)
if (!is.null(kpss)) {
rows[[length(rows) + 1L]] <- data.frame(
Test = "KPSS",
Null = "Level stationary",
Statistic = unname(kpss$statistic),
Lag = unname(kpss$parameter),
p = kpss$p.value,
stringsAsFactors = FALSE
)
}
}
nd <- nsd <- NA_integer_
if (requireNamespace("forecast", quietly = TRUE)) {
nd <- tryCatch(forecast::ndiffs(yts), error = function(e) NA_integer_)
nsd <- if (frequency > 1L && length(y) >= 2L * frequency) {
tryCatch(forecast::nsdiffs(yts), error = function(e) NA_integer_)
} else 0L
}
if (length(rows)) {
out <- do.call(rbind, rows)
} else {
out <- data.frame(
Test = character(), Null = character(), Statistic = numeric(),
Lag = numeric(), p = numeric(), stringsAsFactors = FALSE
)
}
attr(out, "suggested_d") <- nd
attr(out, "suggested_D") <- nsd
out
}
.r4vn_ts_decompose <- function(yts, time, method, frequency) {
if (frequency <= 1L || sum(is.finite(yts)) < 2L * frequency) return(NULL)
if (anyNA(yts)) return(NULL)
if (method == "auto") method <- "stl"
if (method == "stl") {
z <- tryCatch(stats::stl(yts, s.window = "periodic", robust = TRUE),
error = function(e) NULL)
if (is.null(z)) return(NULL)
m <- z$time.series
return(data.frame(
time = time,
observed = as.numeric(yts),
seasonal = as.numeric(m[, "seasonal"]),
trend = as.numeric(m[, "trend"]),
remainder = as.numeric(m[, "remainder"])
))
}
if (method == "classical") {
z <- tryCatch(stats::decompose(yts, type = "additive"), error = function(e) NULL)
if (is.null(z)) return(NULL)
return(data.frame(
time = time,
observed = as.numeric(yts),
seasonal = as.numeric(z$seasonal),
trend = as.numeric(z$trend),
remainder = as.numeric(z$random)
))
}
NULL
}
.r4vn_ts_acf <- function(yts, lag_max = NULL, pacf = FALSE) {
y <- as.numeric(yts)
y <- y[is.finite(y)]
if (length(y) < 3L) return(NULL)
if (is.null(lag_max)) lag_max <- min(floor(length(y) / 3), max(10L, 2L * stats::frequency(yts)))
lag_max <- max(1L, min(as.integer(lag_max), length(y) - 1L))
a <- if (pacf) {
stats::pacf(y, lag.max = lag_max, plot = FALSE, na.action = stats::na.pass)
} else {
stats::acf(y, lag.max = lag_max, plot = FALSE, na.action = stats::na.pass)
}
data.frame(
Lag = as.numeric(a$lag),
Correlation = as.numeric(a$acf),
stringsAsFactors = FALSE
)
}
.r4vn_ts_xreg_matrix <- function(d, xreg_names) {
if (!length(xreg_names)) return(NULL)
x <- as.data.frame(d[, xreg_names, drop = FALSE])
mm <- stats::model.matrix(~ . - 1, data = x)
storage.mode(mm) <- "double"
mm
}
.r4vn_tabts_auto_arima_base <- function(yts, xmat = NULL, frequency = 1L,
seasonal = NULL, options = list()) {
# Dependency-light automatic ARIMA used only when package `forecast` is not
# installed. It searches a deliberately compact, transparent grid of
# stats::arima() models and selects the smallest finite AICc. This is not a
# reimplementation of Hyndman-Khandakar auto.arima(); installing `forecast`
# automatically upgrades tabts() to that engine.
y <- as.numeric(yts)
n <- sum(is.finite(y))
if (n < 8L) stop("Automatic ARIMA requires at least 8 observed time points.", call. = FALSE)
seas_user <- if (is.null(seasonal)) NULL else as.integer(seasonal)
seasonal_ok <- frequency > 1L && n >= 2L * frequency
d_values <- 0:1
D_values <- if (!is.null(seas_user)) seas_user[2L] else if (seasonal_ok) 0:1 else 0L
grid <- expand.grid(
p = 0:2, d = d_values, q = 0:2,
P = if (seasonal_ok || !is.null(seas_user)) 0:1 else 0L,
D = D_values,
Q = if (seasonal_ok || !is.null(seas_user)) 0:1 else 0L,
KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE
)
if (!is.null(seas_user)) {
grid$P <- seas_user[1L]
grid$D <- seas_user[2L]
grid$Q <- seas_user[3L]
grid <- unique(grid)
}
# Keep runtime predictable for routine analyses while retaining common
# parsimonious ARIMA/SARIMA specifications.
grid <- grid[(grid$p + grid$q + grid$P + grid$Q) <= 2L, , drop = FALSE]
grid <- grid[(grid$d + grid$D) <= 2L, , drop = FALSE]
best <- NULL
best_aicc <- Inf
best_spec <- NULL
for (i in seq_len(nrow(grid))) {
g <- grid[i, ]
total_diff <- g$d + g$D
fit <- tryCatch(
suppressWarnings(stats::arima(
yts,
order = c(g$p, g$d, g$q),
seasonal = list(order = c(g$P, g$D, g$Q), period = max(1L, frequency)),
xreg = xmat,
include.mean = isTRUE(options$include_mean) && total_diff == 0L,
method = "ML"
)),
error = function(e) NULL
)
if (is.null(fit)) next
aic <- tryCatch(as.numeric(stats::AIC(fit)), error = function(e) NA_real_)
k <- length(tryCatch(stats::coef(fit), error = function(e) numeric())) + 1L
aicc <- if (is.finite(aic) && n > k + 1L) {
aic + 2 * k * (k + 1) / (n - k - 1)
} else aic
if (is.finite(aicc) && aicc < best_aicc) {
best <- fit
best_aicc <- aicc
best_spec <- c(p = g$p, d = g$d, q = g$q,
P = g$P, D = g$D, Q = g$Q, period = max(1L, frequency))
}
}
if (is.null(best)) {
stop("Automatic base-R ARIMA search could not fit a valid model.", call. = FALSE)
}
if (!is.null(xmat)) attr(best, "r4vn_xreg") <- xmat
attr(best, "r4vn_auto_base") <- TRUE
attr(best, "r4vn_auto_spec") <- best_spec
attr(best, "r4vn_auto_aicc") <- best_aicc
best
}
.r4vn_ts_fit_engine <- function(
engine, yts, xmat, frequency, order, seasonal, options) {
engine <- tolower(engine)
if (engine == "arima") {
if (!requireNamespace("forecast", quietly = TRUE)) {
if (is.null(order)) {
return(.r4vn_tabts_auto_arima_base(
yts = yts, xmat = xmat, frequency = frequency,
seasonal = seasonal, options = options
))
}
seas <- if (is.null(seasonal)) c(0L, 0L, 0L) else as.integer(seasonal)
total_diff <- as.integer(order)[2L] + seas[2L]
fit <- stats::arima(
yts, order = as.integer(order),
seasonal = list(order = seas, period = frequency),
xreg = xmat,
include.mean = isTRUE(options$include_mean) && total_diff == 0L,
method = "ML"
)
# predict.Arima() re-evaluates the original xreg expression. Preserve
# the resolved matrix because the local symbol `xmat` no longer exists
# after this helper returns.
if (!is.null(xmat)) attr(fit, "r4vn_xreg") <- xmat
attr(fit, "r4vn_auto_base") <- FALSE
return(fit)
}
if (is.null(order)) {
args <- c(
list(
y = yts,
xreg = xmat,
seasonal = frequency > 1L,
stepwise = isTRUE(options$stepwise),
approximation = isTRUE(options$approximation)
),
options$auto_arima
)
return(do.call(forecast::auto.arima, args))
}
seas <- if (is.null(seasonal)) c(0L, 0L, 0L) else as.integer(seasonal)
ord <- as.integer(order)
# forecast::Arima() only permits a drift term when the total order of
# differencing is exactly one. Setting include.drift = TRUE with
# d + D >= 2 produces a warning even though the term is then ignored.
# Resolve this here so a valid SARIMA specification runs quietly.
total_diff <- ord[2L] + seas[2L]
use_drift <- isTRUE(options$include_drift) && total_diff == 1L
use_mean <- isTRUE(options$include_mean) && total_diff == 0L
# Keep core structural arguments under tabts() control. This prevents
# nested model_options$arima entries such as include.constant or
# include.drift from silently re-enabling an invalid drift when d + D >= 2.
arima_extra <- options$arima
if (is.null(arima_extra)) arima_extra <- list()
protected <- c(
"y", "order", "seasonal", "xreg", "include.drift",
"include.mean", "include.constant", "method"
)
arima_extra[intersect(names(arima_extra), protected)] <- NULL
args <- c(
list(
y = yts,
order = ord,
seasonal = list(order = seas, period = frequency),
xreg = xmat,
include.drift = use_drift,
include.mean = use_mean,
method = "ML"
),
arima_extra
)
# forecast::Arima() should be silent with the guarded arguments above.
# Muffle only its known benign drift warning as a compatibility safeguard
# across forecast versions; all other warnings remain visible.
fit <- withCallingHandlers(
do.call(forecast::Arima, args),
warning = function(w) {
msg <- conditionMessage(w)
if (grepl(
"No drift term fitted as the order of difference is 2 or more",
msg, fixed = TRUE
)) {
invokeRestart("muffleWarning")
}
}
)
return(fit)
}
if (engine == "ets") {
if (!requireNamespace("forecast", quietly = TRUE)) {
stop("ETS requires package `forecast`.", call. = FALSE)
}
args <- c(list(y = yts), options$ets)
return(do.call(forecast::ets, args))
}
if (engine == "regression") {
tt <- seq_along(yts)
df <- data.frame(.y = as.numeric(yts), .t = tt)
if (!is.null(xmat)) {
xx <- as.data.frame(xmat)
df <- cbind(df, xx)
rhs <- c(".t", names(xx))
} else rhs <- ".t"
if (frequency > 1L) {
df$.season <- factor((tt - 1L) %% frequency + 1L)
rhs <- c(rhs, ".season")
}
fm <- stats::as.formula(paste(".y ~", paste(rhs, collapse = " + ")))
fit <- stats::lm(fm, data = df)
attr(fit, "r4vn_engine") <- "regression"
attr(fit, "r4vn_frequency") <- frequency
attr(fit, "r4vn_xreg_names") <- if (!is.null(xmat)) colnames(xmat) else character()
return(fit)
}
stop("Unsupported engine: ", engine, call. = FALSE)
}
.r4vn_ts_test_split <- function(n, test) {
if (is.null(test)) return(NULL)
if (!is.numeric(test) || length(test) != 1L || is.na(test) || test <= 0) {
stop("`test` must be a positive integer or a proportion between 0 and 1.", call. = FALSE)
}
ntest <- if (test < 1) ceiling(n * test) else as.integer(test)
if (ntest < 1L || ntest >= n - 4L) {
stop("`test` leaves too few observations for model fitting.", call. = FALSE)
}
list(n_train = n - ntest, n_test = ntest)
}
.r4vn_ts_forecast_engine <- function(fit, engine, h, level, xreg, options) {
if (h <= 0L) return(NULL)
if (inherits(fit, "Arima") && requireNamespace("forecast", quietly = TRUE)) {
args <- list(object = fit, h = h, level = level * 100)
if (!is.null(xreg)) args$xreg <- xreg
args$bootstrap <- isTRUE(options$bootstrap)
args$biasadj <- isTRUE(options$biasadj)
z <- do.call(forecast::forecast, args)
return(list(
mean = as.numeric(z$mean),
lower = as.matrix(z$lower),
upper = as.matrix(z$upper),
level = level
))
}
if (inherits(fit, "ets") && requireNamespace("forecast", quietly = TRUE)) {
z <- forecast::forecast(
fit, h = h, level = level * 100,
bootstrap = isTRUE(options$bootstrap)
)
return(list(
mean = as.numeric(z$mean),
lower = as.matrix(z$lower),
upper = as.matrix(z$upper),
level = level
))
}
# stats::arima() and forecast::Arima() both use class "Arima". The latter
# has already been handled above when forecast is available; this branch is
# the dependency-free forecasting path for a fixed stats::arima() model.
if (inherits(fit, "Arima") || inherits(fit, "arima")) {
predict_fit <- fit
fitted_xreg <- attr(fit, "r4vn_xreg", exact = TRUE)
if (identical(predict_fit$call$xreg, quote(xmat))) {
# Assigning NULL removes the stale xreg argument for a model without
# regressors; assigning a matrix makes ARIMAX self-contained.
predict_fit$call$xreg <- fitted_xreg
}
pr <- stats::predict(predict_fit, n.ahead = h, newxreg = xreg)
mu <- as.numeric(pr$pred)
se <- as.numeric(pr$se)
lo <- sapply(level, function(lv) mu + stats::qnorm((1 - lv) / 2) * se)
hi <- sapply(level, function(lv) mu + stats::qnorm(1 - (1 - lv) / 2) * se)
if (is.null(dim(lo))) lo <- matrix(lo, ncol = 1L)
if (is.null(dim(hi))) hi <- matrix(hi, ncol = 1L)
colnames(lo) <- colnames(hi) <- paste0(round(level * 100), "%")
return(list(mean = mu, lower = lo, upper = hi, level = level))
}
if (inherits(fit, "lm") && identical(attr(fit, "r4vn_engine"), "regression")) {
freq <- attr(fit, "r4vn_frequency") %||% 1L
n <- length(stats::fitted(fit))
nd <- data.frame(.t = n + seq_len(h))
xnames <- attr(fit, "r4vn_xreg_names") %||% character()
if (length(xnames)) {
if (is.null(xreg)) stop("Future xreg values are required.", call. = FALSE)
xx <- as.data.frame(xreg)
names(xx) <- xnames
nd <- cbind(nd, xx)
}
if (freq > 1L) nd$.season <- factor((n + seq_len(h) - 1L) %% freq + 1L)
prs <- lapply(level, function(lv) {
stats::predict(fit, newdata = nd, interval = "prediction", level = lv)
})
mu <- prs[[which.max(level)]][, "fit"]
lo <- do.call(cbind, lapply(prs, function(z) z[, "lwr"]))
hi <- do.call(cbind, lapply(prs, function(z) z[, "upr"]))
colnames(lo) <- colnames(hi) <- paste0(round(level * 100), "%")
return(list(mean = mu, lower = lo, upper = hi, level = level))
}
stop("Forecasting is not implemented for this fitted model.", call. = FALSE)
}
.r4vn_ts_accuracy <- function(actual, predicted, train, frequency = 1L) {
ok <- is.finite(actual) & is.finite(predicted)
actual <- actual[ok]
predicted <- predicted[ok]
if (!length(actual)) {
return(data.frame(RMSE = NA, MAE = NA, MAPE = NA, sMAPE = NA, MASE = NA))
}
e <- actual - predicted
rmse <- sqrt(mean(e^2))
mae <- mean(abs(e))
mape <- if (any(actual == 0)) NA_real_ else mean(abs(e / actual)) * 100
denom_smape <- abs(actual) + abs(predicted)
smape <- mean(ifelse(denom_smape == 0, 0, 200 * abs(e) / denom_smape))
lag <- max(1L, as.integer(frequency))
naive <- if (length(train) > lag) abs(train[(lag + 1L):length(train)] - train[seq_len(length(train) - lag)]) else numeric()
scale <- if (length(naive)) mean(naive, na.rm = TRUE) else NA_real_
mase <- if (is.finite(scale) && scale > 0) mae / scale else NA_real_
data.frame(RMSE = rmse, MAE = mae, MAPE = mape, sMAPE = smape, MASE = mase)
}
.r4vn_ts_aic <- function(fit) {
tryCatch(as.numeric(stats::AIC(fit)), error = function(e) NA_real_)
}
.r4vn_ts_bic <- function(fit) {
tryCatch(as.numeric(stats::BIC(fit)), error = function(e) NA_real_)
}
.r4vn_ts_aicc <- function(fit) {
if (!is.null(fit$aicc) && is.finite(fit$aicc)) return(as.numeric(fit$aicc))
aic <- .r4vn_ts_aic(fit)
n <- length(stats::na.omit(.r4vn_ts_residuals(fit)))
k <- .r4vn_ts_model_df(fit) + 1L
if (!is.finite(aic) || n <= k + 1L) return(NA_real_)
aic + 2 * k * (k + 1) / (n - k - 1)
}
.r4vn_tabts_model_info <- function(fit, model_name, frequency = 1L,
its_meta = NULL, n = NA_integer_) {
if (is.null(fit)) return(NULL)
spec <- toupper(model_name)
details <- ""
if (inherits(fit, "Arima") || inherits(fit, "arima")) {
ar <- fit$arma
if (length(ar) >= 7L) {
spec <- sprintf(
"ARIMA(%d,%d,%d)(%d,%d,%d)[%d]",
ar[1L], ar[6L], ar[2L], ar[3L], ar[7L], ar[4L], ar[5L]
)
}
if (isTRUE(attr(fit, "r4vn_auto_base", exact = TRUE))) {
details <- "Automatic compact AICc search using stats::arima (base R)"
} else if (!is.null(fit$call) && grepl("auto.arima", paste(deparse(fit$call), collapse = ""), fixed = TRUE)) {
details <- "Automatic ARIMA selection"
} else {
details <- "ARIMA/SARIMA model"
}
} else if (inherits(fit, "ets")) {
spec <- paste0("ETS ", fit$method %||% "")
details <- "Exponential smoothing state-space model"
} else if (identical(model_name, "its") && !is.null(its_meta)) {
spec <- paste0("ITS (", its_meta$family %||% "model", ")")
details <- paste0(
"Correlation: ", its_meta$correlation %||% "none",
"; pre/post time points: ", its_meta$n_pre %||% NA_integer_,
"/", its_meta$n_post %||% NA_integer_,
if (isTRUE(its_meta$controlled)) paste0(
"; controlled ITS; reference group: ", its_meta$reference_group %||% "first level"
) else ""
)
} else if (inherits(fit, "lm")) {
spec <- if (identical(model_name, "regression")) "Time-series regression" else toupper(model_name)
details <- "Linear regression with time and available seasonal/external regressors"
}
data.frame(
Model = spec,
Details = details,
N = as.integer(n),
Seasonal_period = as.integer(frequency),
AIC = .r4vn_ts_aic(fit),
AICc = .r4vn_ts_aicc(fit),
BIC = .r4vn_ts_bic(fit),
stringsAsFactors = FALSE,
check.names = FALSE
)
}
.r4vn_ts_model_df <- function(fit) {
cf <- tryCatch(stats::coef(fit), error = function(e) numeric())
sum(is.finite(cf))
}
.r4vn_ts_choose_model <- function(tab, criterion, diagnostic_gate, test) {
if (is.null(tab) || nrow(tab) == 1L) return(tolower(tab$Model[1L]))
metric <- switch(
criterion,
aicc = "AICc",
rmse = "RMSE",
mae = "MAE",
mape = "MAPE"
)
if (criterion != "aicc" && is.null(test)) {
warning(
"`criterion = \"", criterion, "\"` requires `test`; using AICc instead.",
call. = FALSE
)
metric <- "AICc"
}
if (!metric %in% names(tab) || all(!is.finite(tab[[metric]]))) metric <- "AICc"
pool <- seq_len(nrow(tab))
if (isTRUE(diagnostic_gate) && "Residual_OK" %in% names(tab)) {
ok <- which(tab$Residual_OK %in% TRUE)
if (length(ok)) pool <- ok
}
vals <- tab[[metric]][pool]
if (all(!is.finite(vals))) vals <- tab$AIC[pool]
idx <- pool[which.min(vals)]
tolower(tab$Model[idx])
}
.r4vn_ts_coefficients <- function(fit, ci_level = .95, family = "gaussian", vcov_override = NULL) {
cf <- tryCatch(stats::coef(fit), error = function(e) NULL)
if (is.null(cf) || !length(cf)) return(NULL)
vc <- vcov_override
if (is.null(vc)) vc <- tryCatch(stats::vcov(fit), error = function(e) NULL)
se <- if (!is.null(vc)) sqrt(pmax(0, diag(vc))) else rep(NA_real_, length(cf))
z <- stats::qnorm(1 - (1 - ci_level) / 2)
stat <- cf / se
p <- 2 * stats::pnorm(abs(stat), lower.tail = FALSE)
if (family %in% c("poisson", "quasipoisson", "negativebinomial")) {
est <- exp(cf)
lo <- exp(cf - z * se)
hi <- exp(cf + z * se)
measure <- "IRR"
} else {
est <- cf
lo <- cf - z * se
hi <- cf + z * se
measure <- "Beta"
}
data.frame(
Term = names(cf),
Measure = measure,
Estimate = unname(est),
SE = unname(se),
CI_lower = unname(lo),
CI_upper = unname(hi),
Statistic = unname(stat),
p = unname(p),
stringsAsFactors = FALSE
)
}
.r4vn_ts_fitted <- function(fit, n = NULL) {
z <- tryCatch(as.numeric(stats::fitted(fit)), error = function(e) numeric())
if (is.null(n)) return(z)
if (length(z) == n) return(z)
out <- rep(NA_real_, n)
if (length(z)) out[(n - length(z) + 1L):n] <- z
out
}
.r4vn_ts_residuals <- function(fit, n = NULL) {
z <- tryCatch(as.numeric(stats::residuals(fit)), error = function(e) numeric())
if (is.null(n)) return(z)
if (length(z) == n) return(z)
out <- rep(NA_real_, n)
if (length(z)) out[(n - length(z) + 1L):n] <- z
out
}
.r4vn_ts_diagnostics <- function(fit, diagnostic_options, frequency, model_df = 0L) {
r <- .r4vn_ts_residuals(fit)
r <- r[is.finite(r)]
if (length(r) < 5L) return(NULL)
lag <- diagnostic_options$ljung_lag
if (is.null(lag)) {
lag <- min(
max(10L, 2L * as.integer(frequency)),
max(model_df + 3L, floor(length(r) / 4))
)
}
lag <- min(as.integer(lag), length(r) - 1L)
lag <- max(lag, min(length(r) - 1L, model_df + 1L))
lb <- tryCatch(
stats::Box.test(r, lag = lag, type = "Ljung-Box", fitdf = min(model_df, lag - 1L)),
error = function(e) NULL
)
rows <- list()
if (!is.null(lb)) {
rows[[1L]] <- data.frame(
Diagnostic = "Ljung-Box",
Statistic = unname(lb$statistic),
df = unname(lb$parameter),
p = lb$p.value,
Value = NA_real_,
stringsAsFactors = FALSE
)
}
rows[[length(rows) + 1L]] <- data.frame(
Diagnostic = "Residual mean",
Statistic = NA_real_, df = NA_real_, p = NA_real_,
Value = mean(r), stringsAsFactors = FALSE
)
rows[[length(rows) + 1L]] <- data.frame(
Diagnostic = "Residual SD",
Statistic = NA_real_, df = NA_real_, p = NA_real_,
Value = stats::sd(r), stringsAsFactors = FALSE
)
if (isTRUE(diagnostic_options$normality) && length(r) >= 3L && length(r) <= 5000L) {
sw <- tryCatch(stats::shapiro.test(r), error = function(e) NULL)
if (!is.null(sw)) {
rows[[length(rows) + 1L]] <- data.frame(
Diagnostic = "Shapiro-Wilk",
Statistic = unname(sw$statistic), df = NA_real_, p = sw$p.value,
Value = NA_real_, stringsAsFactors = FALSE
)
}
}
do.call(rbind, rows)
}
.r4vn_ts_future_xreg <- function(future_xreg, xreg_names, h) {
if (!length(xreg_names)) return(NULL)
if (is.null(future_xreg)) {
stop(
"`future_xreg` is required when forecasting a model with external regressors.",
call. = FALSE
)
}
z <- as.data.frame(future_xreg)
miss <- setdiff(xreg_names, names(z))
if (length(miss)) {
stop("`future_xreg` is missing: ", paste(miss, collapse = ", "), call. = FALSE)
}
if (nrow(z) < h) stop("`future_xreg` must contain at least `forecast` rows.", call. = FALSE)
z <- z[seq_len(h), xreg_names, drop = FALSE]
stats::model.matrix(~ . - 1, data = z)
}
.r4vn_ts_forecast_table <- function(fc, time, level) {
if (is.null(fc)) return(NULL)
out <- data.frame(Time = time, Forecast = fc$mean, stringsAsFactors = FALSE)
lev <- as.numeric(fc$level %||% level)
lo <- as.matrix(fc$lower)
hi <- as.matrix(fc$upper)
if (ncol(lo) != length(lev)) lev <- seq_len(ncol(lo))
for (j in seq_len(ncol(lo))) {
tag <- if (lev[j] <= 1) round(lev[j] * 100) else round(lev[j])
out[[paste0("PI", tag, "_lower")]] <- lo[, j]
out[[paste0("PI", tag, "_upper")]] <- hi[, j]
}
out
}
.r4vn_ts_fit_its <- function(
d, xreg_names, intervention_value, family, correlation, frequency,
group_name, pop_name, rate, h, level, future_xreg, its_options,
diagnostic_options, table_options, step) {
if (is.null(intervention_value) || !length(intervention_value)) {
stop("ITS requires `intervention =`.", call. = FALSE)
}
controlled <- !is.null(group_name)
if (controlled) {
d$.group <- factor(d$.group)
ref_group <- its_options$reference_group
if (!is.null(ref_group)) {
ref_group <- as.character(ref_group)[1L]
if (!ref_group %in% levels(d$.group)) {
stop("`its_options$reference_group` was not found in the grouping variable.", call. = FALSE)
}
d$.group <- stats::relevel(d$.group, ref = ref_group)
}
groups <- levels(d$.group)
if (length(groups) < 2L) stop("Controlled ITS requires at least two groups.", call. = FALSE)
}
# Create time index within the common time grid.
unique_time <- sort(unique(d$.time_original))
d$.time <- match(d$.time_original, unique_time) - 1L
d$.time <- d$.time * (its_options$time_scale %||% 1)
int_times <- .r4vn_ts_match_intervention_times(intervention_value, unique_time)
int_idx <- match(int_times, unique_time) - 1L
if (anyNA(int_idx)) {
stop("One or more intervention times could not be matched to the time series.", call. = FALSE)
}
intervention_cols <- character()
after_cols <- character()
for (j in seq_along(int_idx)) {
ii <- paste0(".int", j)
aa <- paste0(".after", j)
threshold <- int_idx[j] * (its_options$time_scale %||% 1)
d[[ii]] <- as.integer(d$.time >= threshold)
d[[aa]] <- pmax(0, d$.time - threshold)
intervention_cols <- c(intervention_cols, ii)
after_cols <- c(after_cols, aa)
}
if (family == "auto") {
yy <- d$.y[is.finite(d$.y)]
is_count <- length(yy) && all(yy >= 0) && all(abs(yy - round(yy)) < 1e-8)
if (is_count) {
disp <- if (length(yy) > 1L && mean(yy) > 0) stats::var(yy) / mean(yy) else 1
if (is.finite(disp) && disp > 1.5) {
family <- if (requireNamespace("MASS", quietly = TRUE)) "negativebinomial" else "quasipoisson"
} else family <- "poisson"
} else family <- "gaussian"
}
if (identical(correlation, "ar1") && !identical(family, "gaussian")) {
warning(
"`correlation = \"ar1\"` is available only for Gaussian ITS; using `correlation = \"none\"`. ",
"For count ITS, use `correlation = \"nw\"` when Newey-West covariance is desired.",
call. = FALSE
)
correlation <- "none"
}
rhs_base <- c(".time", intervention_cols, after_cols, xreg_names)
if (isTRUE(its_options$include_season) && frequency > 1L) {
kmax <- max(1L, min(
as.integer(its_options$seasonal_harmonics %||% 1L),
max(1L, floor((frequency - 1L) / 2L))
))
tidx <- match(d$.time_original, unique_time) - 1L
seasonal_terms <- character()
for (k in seq_len(kmax)) {
sn <- paste0(".sin", k); cs <- paste0(".cos", k)
d[[sn]] <- sin(2 * pi * k * tidx / frequency)
d[[cs]] <- cos(2 * pi * k * tidx / frequency)
seasonal_terms <- c(seasonal_terms, sn, cs)
}
rhs_base <- c(rhs_base, seasonal_terms)
}
if (controlled) {
rhs <- paste0("(", paste(rhs_base, collapse = " + "), ") * .group")
} else {
rhs <- paste(rhs_base, collapse = " + ")
}
offset_txt <- if (!is.null(pop_name) &&
family %in% c("poisson", "quasipoisson", "negativebinomial")) {
if (any(!is.finite(d$.population) | d$.population <= 0)) {
stop("`population` must be positive and non-missing for count ITS.", call. = FALSE)
}
" + offset(log(.population))"
} else ""
fm <- stats::as.formula(paste(".y ~", rhs, offset_txt))
fit <- NULL
vc_override <- NULL
correlation_used <- correlation
if (correlation == "auto") {
correlation_used <- "none"
if (family == "gaussian" && requireNamespace("nlme", quietly = TRUE)) {
cor_form <- if (controlled) stats::as.formula("~ .time | .group") else stats::as.formula("~ .time")
fit0 <- tryCatch(
nlme::gls(fm, data = d, method = "ML", na.action = stats::na.omit),
error = function(e) NULL
)
fit1 <- tryCatch(
nlme::gls(
fm, data = d,
correlation = nlme::corAR1(form = cor_form),
method = "ML",
na.action = stats::na.omit
),
error = function(e) NULL
)
if (!is.null(fit0) && !is.null(fit1)) {
a0 <- tryCatch(stats::AIC(fit0), error = function(e) Inf)
a1 <- tryCatch(stats::AIC(fit1), error = function(e) Inf)
# Require a small but meaningful improvement before adding AR(1).
if (is.finite(a1) && a1 + 2 < a0) {
fit <- fit1
correlation_used <- "ar1"
} else {
fit <- fit0
correlation_used <- "none"
}
}
}
}
if (family == "gaussian" && correlation_used == "ar1" && is.null(fit)) {
if (!requireNamespace("nlme", quietly = TRUE)) {
warning("Package `nlme` is unavailable; fitting Gaussian ITS without AR(1).", call. = FALSE)
correlation_used <- "none"
} else {
cor_form <- if (controlled) stats::as.formula("~ .time | .group") else stats::as.formula("~ .time")
fit <- nlme::gls(
fm, data = d,
correlation = nlme::corAR1(form = cor_form),
method = "ML",
na.action = stats::na.omit
)
}
}
if (is.null(fit)) {
if (family == "gaussian") {
fit <- stats::lm(fm, data = d)
} else if (family == "poisson") {
fit <- stats::glm(fm, data = d, family = stats::poisson())
} else if (family == "quasipoisson") {
fit <- stats::glm(fm, data = d, family = stats::quasipoisson())
} else if (family == "negativebinomial") {
if (!requireNamespace("MASS", quietly = TRUE)) {
stop("Negative-binomial ITS requires package `MASS`.", call. = FALSE)
}
# MASS::glm.nb() evaluates `link` non-standardly. Inside wrappers,
# passing a symbol or local variable may be reinterpreted as the name
# of that object. Use the literal string "log" so the initial Poisson
# fit and subsequent negative-binomial iterations receive the intended
# canonical log link reliably.
fit <- MASS::glm.nb(fm, data = d)
}
}
if (correlation_used == "nw") {
if (!requireNamespace("sandwich", quietly = TRUE)) {
warning("Package `sandwich` is unavailable; ordinary model covariance is used.", call. = FALSE)
correlation_used <- "none"
} else {
lag <- its_options$nw_lag
if (is.null(lag)) lag <- max(1L, floor(4 * (nrow(d) / 100)^(2 / 9)))
vc_override <- tryCatch(
sandwich::NeweyWest(fit, lag = lag, prewhite = FALSE, adjust = TRUE),
error = function(e) NULL
)
if (is.null(vc_override)) {
warning("Newey-West covariance could not be computed; ordinary model covariance is used.", call. = FALSE)
correlation_used <- "none"
}
}
}
coeff <- .r4vn_ts_coefficients(
fit,
ci_level = table_options$ci_level,
family = family,
vcov_override = vc_override
)
if (!is.null(coeff)) {
coeff$Term_label <- .r4vn_ts_its_term_labels(
coeff$Term, length(int_idx), controlled
)
}
fitted <- tryCatch(as.numeric(stats::fitted(fit)), error = function(e) rep(NA_real_, nrow(d)))
residuals <- tryCatch(as.numeric(stats::residuals(fit)), error = function(e) rep(NA_real_, nrow(d)))
# gls with omitted rows may return shorter vectors.
if (length(fitted) != nrow(d)) {
ff <- rep(NA_real_, nrow(d))
rr <- rep(NA_real_, nrow(d))
ok <- stats::complete.cases(stats::model.frame(fit))
ff[ok] <- fitted
rr[ok] <- residuals
fitted <- ff
residuals <- rr
}
diag_tab <- .r4vn_ts_diagnostics(
fit, diagnostic_options, frequency,
model_df = .r4vn_ts_model_df(fit)
)
cf <- NULL
effect_at <- NULL
if (isTRUE(its_options$counterfactual) && !controlled) {
nd <- d
for (nm in c(intervention_cols, after_cols)) nd[[nm]] <- 0
cf_pred <- .r4vn_ts_predict_mean(fit, nd, family)
cf <- data.frame(
Time = d$.time_original,
Fitted = fitted,
Counterfactual = cf_pred,
Difference = fitted - cf_pred,
Relative_difference_pct = ifelse(cf_pred == 0, NA_real_, (fitted - cf_pred) / cf_pred * 100),
stringsAsFactors = FALSE
)
ea <- unique(as.integer(its_options$effect_at))
ea <- ea[ea >= 0L]
if (length(ea)) {
base_idx <- int_idx[1L] + 1L
ids <- base_idx + ea
ok <- ids >= 1L & ids <= nrow(cf)
ids <- ids[ok]
ea2 <- ea[ok]
if (length(ids)) {
effect_at <- data.frame(
Period_after_intervention = ea2,
Time = cf$Time[ids],
Fitted = cf$Fitted[ids],
Counterfactual = cf$Counterfactual[ids],
Absolute_difference = cf$Difference[ids],
Relative_difference_pct = cf$Relative_difference_pct[ids],
stringsAsFactors = FALSE
)
}
}
}
if (controlled) {
effect_at <- .r4vn_ts_controlled_effect_at(
fit = fit, d = d, int_idx = int_idx,
intervention_cols = intervention_cols, after_cols = after_cols,
its_options = its_options, family = family,
vc_override = vc_override, ci_level = table_options$ci_level
)
}
forecast_tab <- NULL
if (h > 0L && !controlled) {
future_time <- .r4vn_ts_extend_time(unique_time, h, step)
nd <- .r4vn_ts_make_its_future(
d = d, future_time = future_time, xreg_names = xreg_names,
future_xreg = future_xreg, int_idx = int_idx,
frequency = frequency, pop_name = pop_name, its_options = its_options
)
prs <- lapply(level, function(lv) {
.r4vn_ts_predict_interval(fit, nd, family, lv, vc_override)
})
forecast_tab <- data.frame(
Time = future_time,
Forecast = prs[[which.max(level)]]$fit,
stringsAsFactors = FALSE
)
for (j in seq_along(level)) {
tag <- round(level[j] * 100)
forecast_tab[[paste0("CI", tag, "_lower")]] <- prs[[j]]$lower
forecast_tab[[paste0("CI", tag, "_upper")]] <- prs[[j]]$upper
}
}
meta <- list(
family = family,
correlation = correlation_used,
intervention_times = int_times,
intervention_index = int_idx,
controlled = controlled,
groups = if (controlled) levels(d$.group) else NULL,
reference_group = if (controlled) levels(d$.group)[1L] else NULL,
n_pre = int_idx[1L],
n_post = length(unique_time) - int_idx[1L],
post_min = its_options$post_min %||% 8L
)
list(
fit = fit,
coefficients = coeff,
fitted = fitted,
residuals = residuals,
diagnostics = diag_tab,
counterfactual = cf,
effect_at = effect_at,
forecast = forecast_tab,
meta = meta
)
}
.r4vn_ts_match_intervention_times <- function(value, unique_time) {
if (inherits(unique_time, "Date")) {
vv <- if (inherits(value, "Date")) value else as.Date(value)
out <- sapply(vv, function(v) {
if (v %in% unique_time) return(as.character(v))
as.character(unique_time[which.min(abs(as.numeric(unique_time - v)))])
})
return(as.Date(out))
}
if (inherits(unique_time, "POSIXt")) {
vv <- as.POSIXct(value)
idx <- vapply(vv, function(v) which.min(abs(as.numeric(unique_time - v))), integer(1))
return(unique_time[idx])
}
if (is.numeric(unique_time)) {
vv <- as.numeric(value)
return(vapply(vv, function(v) unique_time[which.min(abs(unique_time - v))], numeric(1)))
}
value
}
.r4vn_ts_its_term_labels <- function(terms, n_int, controlled) {
label_piece <- function(piece) {
if (identical(piece, "(Intercept)")) return("Baseline level")
if (identical(piece, ".time")) return("Pre-intervention trend")
for (j in seq_len(n_int)) {
if (identical(piece, paste0(".int", j))) {
return(if (n_int == 1L) "Immediate level change" else
paste0("Immediate level change (intervention ", j, ")"))
}
if (identical(piece, paste0(".after", j))) {
return(if (n_int == 1L) "Trend change" else
paste0("Trend change (intervention ", j, ")"))
}
}
if (grepl("^\\.group", piece)) {
return(paste0("Group: ", sub("^\\.group", "", piece)))
}
if (grepl("^\\.sin[0-9]+$", piece)) {
return(paste0("Seasonal sine harmonic ", sub("^\\.sin", "", piece)))
}
if (grepl("^\\.cos[0-9]+$", piece)) {
return(paste0("Seasonal cosine harmonic ", sub("^\\.cos", "", piece)))
}
if (grepl("^\\.season", piece)) {
return(paste0("Season: ", sub("^\\.season", "", piece)))
}
piece
}
vapply(terms, function(term) {
pieces <- strsplit(term, ":", fixed = TRUE)[[1L]]
paste(vapply(pieces, label_piece, character(1)), collapse = " \u00d7 ")
}, character(1), USE.NAMES = FALSE)
}
.r4vn_ts_controlled_effect_at <- function(
fit, d, int_idx, intervention_cols, after_cols, its_options, family,
vc_override = NULL, ci_level = .95) {
if (!".group" %in% names(d) || !length(int_idx)) return(NULL)
groups <- if (is.factor(d$.group)) levels(d$.group) else unique(as.character(d$.group))
groups <- groups[groups %in% as.character(d$.group)]
if (length(groups) < 2L) return(NULL)
ref <- groups[1L]
unique_time <- sort(unique(d$.time_original))
ea <- unique(as.integer(its_options$effect_at))
ea <- ea[is.finite(ea) & ea >= 0L]
if (!length(ea)) return(NULL)
base_idx <- int_idx[1L] + 1L
ids <- base_idx + ea
ok <- ids >= 1L & ids <= length(unique_time)
ids <- ids[ok]; ea <- ea[ok]
if (!length(ids)) return(NULL)
beta <- tryCatch(stats::coef(fit), error = function(e) NULL)
vc <- vc_override %||% tryCatch(stats::vcov(fit), error = function(e) NULL)
tt <- tryCatch(stats::delete.response(stats::terms(fit)), error = function(e) NULL)
if (is.null(beta) || is.null(vc) || is.null(tt)) return(NULL)
bnames <- names(beta)
if (is.null(bnames) || !length(bnames)) return(NULL)
vc <- tryCatch(vc[bnames, bnames, drop = FALSE], error = function(e) NULL)
if (is.null(vc)) return(NULL)
zcrit <- stats::qnorm(1 - (1 - ci_level) / 2)
rows <- list()
for (k in seq_along(ids)) {
tm <- unique_time[ids[k]]
for (g in groups[-1L]) {
rg <- d[d$.time_original == tm & as.character(d$.group) == g, , drop = FALSE]
rr <- d[d$.time_original == tm & as.character(d$.group) == ref, , drop = FALSE]
if (!nrow(rg) || !nrow(rr)) next
rg <- rg[1L, , drop = FALSE]; rr <- rr[1L, , drop = FALSE]
rg_cf <- rg; rr_cf <- rr
for (nm in c(intervention_cols, after_cols)) {
rg_cf[[nm]] <- 0
rr_cf[[nm]] <- 0
}
nd <- rbind(rg, rg_cf, rr, rr_cf)
X <- tryCatch(stats::model.matrix(tt, data = nd), error = function(e) NULL)
if (is.null(X)) next
miss <- setdiff(bnames, colnames(X))
if (length(miss)) next
X <- X[, bnames, drop = FALSE]
L <- X[1L, ] - X[2L, ] - X[3L, ] + X[4L, ]
est_link <- as.numeric(sum(L * beta))
se_link <- sqrt(max(0, as.numeric(t(L) %*% vc %*% L)))
stat <- if (is.finite(se_link) && se_link > 0) est_link / se_link else NA_real_
pval <- if (is.finite(stat)) 2 * stats::pnorm(abs(stat), lower.tail = FALSE) else NA_real_
lo_link <- est_link - zcrit * se_link
hi_link <- est_link + zcrit * se_link
is_count <- family %in% c("poisson", "quasipoisson", "negativebinomial")
est <- if (is_count) exp(est_link) else est_link
lo <- if (is_count) exp(lo_link) else lo_link
hi <- if (is_count) exp(hi_link) else hi_link
rows[[length(rows) + 1L]] <- data.frame(
Period_after_intervention = ea[k],
Time = tm,
Comparison = paste0(g, " vs ", ref),
Measure = if (is_count) "Ratio of rate ratios" else "Difference in intervention effects",
Estimate = est,
CI_lower = lo,
CI_upper = hi,
p = pval,
Relative_difference_pct = if (is_count) (est - 1) * 100 else NA_real_,
stringsAsFactors = FALSE,
check.names = FALSE
)
}
}
if (!length(rows)) return(NULL)
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_ts_predict_mean <- function(fit, newdata, family) {
if (inherits(fit, "gls")) return(as.numeric(stats::predict(fit, newdata = newdata)))
if (inherits(fit, "glm")) return(as.numeric(stats::predict(fit, newdata = newdata, type = "response")))
as.numeric(stats::predict(fit, newdata = newdata))
}
.r4vn_ts_make_its_future <- function(
d, future_time, xreg_names, future_xreg, int_idx,
frequency, pop_name, its_options) {
h <- length(future_time)
n0 <- length(unique(d$.time_original))
nd <- data.frame(.time_original = future_time)
nd$.time <- (n0 + seq_len(h) - 1L) * (its_options$time_scale %||% 1)
for (j in seq_along(int_idx)) {
threshold <- int_idx[j] * (its_options$time_scale %||% 1)
nd[[paste0(".int", j)]] <- as.integer(nd$.time >= threshold)
nd[[paste0(".after", j)]] <- pmax(0, nd$.time - threshold)
}
if (length(xreg_names)) {
if (is.null(future_xreg)) {
stop("Future xreg values are required for ITS forecasting.", call. = FALSE)
}
fx <- as.data.frame(future_xreg)
miss <- setdiff(xreg_names, names(fx))
if (length(miss)) stop("`future_xreg` is missing: ", paste(miss, collapse = ", "), call. = FALSE)
if (nrow(fx) < h) stop("`future_xreg` has too few rows.", call. = FALSE)
nd[xreg_names] <- fx[seq_len(h), xreg_names, drop = FALSE]
}
if (isTRUE(its_options$include_season) && frequency > 1L) {
kmax <- max(1L, min(
as.integer(its_options$seasonal_harmonics %||% 1L),
max(1L, floor((frequency - 1L) / 2L))
))
tidx <- n0 + seq_len(h) - 1L
for (k in seq_len(kmax)) {
nd[[paste0(".sin", k)]] <- sin(2 * pi * k * tidx / frequency)
nd[[paste0(".cos", k)]] <- cos(2 * pi * k * tidx / frequency)
}
}
if (!is.null(pop_name)) {
if (is.null(future_xreg) || !".population" %in% names(future_xreg)) {
last_pop <- tail(d$.population[is.finite(d$.population)], 1L)
nd$.population <- rep(last_pop, h)
} else {
nd$.population <- future_xreg$.population[seq_len(h)]
}
}
nd
}
.r4vn_ts_predict_interval <- function(fit, newdata, family, level = .95, vc_override = NULL) {
if (inherits(fit, "lm") && !inherits(fit, "glm")) {
p <- stats::predict(fit, newdata = newdata, interval = "confidence", level = level)
return(list(fit = p[, "fit"], lower = p[, "lwr"], upper = p[, "upr"]))
}
if (inherits(fit, "glm")) {
p <- stats::predict(fit, newdata = newdata, type = "link", se.fit = TRUE)
z <- stats::qnorm(1 - (1 - level) / 2)
lo <- p$fit - z * p$se.fit
hi <- p$fit + z * p$se.fit
inv <- fit$family$linkinv
return(list(fit = inv(p$fit), lower = inv(lo), upper = inv(hi)))
}
if (inherits(fit, "gls")) {
tt <- stats::delete.response(stats::terms(fit))
X <- stats::model.matrix(tt, newdata)
beta <- stats::coef(fit)
vc <- vc_override %||% stats::vcov(fit)
mu <- as.numeric(X[, names(beta), drop = FALSE] %*% beta)
se <- sqrt(rowSums((X[, names(beta), drop = FALSE] %*% vc) * X[, names(beta), drop = FALSE]))
z <- stats::qnorm(1 - (1 - level) / 2)
return(list(fit = mu, lower = mu - z * se, upper = mu + z * se))
}
mu <- .r4vn_ts_predict_mean(fit, newdata, family)
list(fit = mu, lower = rep(NA_real_, length(mu)), upper = rep(NA_real_, length(mu)))
}
.r4vn_ts_plot_set <- function(plots, model, h, decomp, do_acf, do_pacf) {
plots <- tolower(as.character(plots))
if (length(plots) == 1L && plots == "none") return(character())
if (length(plots) == 1L && plots %in% c("auto", "all")) {
z <- "series"
if (!is.null(decomp)) z <- c(z, "decomposition")
if (do_acf) z <- c(z, "acf")
if (do_pacf) z <- c(z, "pacf")
z <- c(z, "residual")
if (h > 0L) z <- c(z, "forecast")
if (model == "its") z <- unique(c("its", z, "counterfactual"))
return(unique(z))
}
unique(plots)
}
.r4vn_tabts_plot_spec <- function(type, ...) {
structure(c(list(type = type), list(...)), class = "r4vn_tabts_plot")
}
.r4vn_tabts_make_base_plots <- function(
d, frequency, decomposition, acf_tab, pacf_tab, forecast_tab,
counterfactual, intervention_value, model, requested, options,
outcome_label, time_label) {
out <- list()
common <- list(d = d, counterfactual = NULL,
intervention_value = intervention_value,
options = options, outcome_label = outcome_label,
time_label = time_label)
if ("series" %in% requested) {
out$series <- do.call(.r4vn_tabts_plot_spec,
c(list(type = "series", forecast_tab = NULL), common))
}
if ("forecast" %in% requested && !is.null(forecast_tab)) {
out$forecast <- do.call(.r4vn_tabts_plot_spec,
c(list(type = "series", forecast_tab = forecast_tab), common))
}
if ("its" %in% requested && model == "its") {
out$its <- do.call(.r4vn_tabts_plot_spec,
c(list(type = "series", forecast_tab = NULL), common))
}
if ("decomposition" %in% requested && !is.null(decomposition)) {
out$decomposition <- .r4vn_tabts_plot_spec(
"decomposition", decomposition = decomposition, options = options
)
}
if ("acf" %in% requested && !is.null(acf_tab)) {
out$acf <- .r4vn_tabts_plot_spec(
"acf", acf_tab = acf_tab, n = nrow(d), options = options, kind = "acf"
)
}
if ("pacf" %in% requested && !is.null(pacf_tab)) {
out$pacf <- .r4vn_tabts_plot_spec(
"acf", acf_tab = pacf_tab, n = nrow(d), options = options, kind = "pacf"
)
}
if ("residual" %in% requested && any(is.finite(d$.residual))) {
out$residual <- .r4vn_tabts_plot_spec("residual", d = d, options = options)
}
if ("counterfactual" %in% requested && !is.null(counterfactual)) {
out$counterfactual <- .r4vn_tabts_plot_spec(
"counterfactual", d = d, counterfactual = counterfactual,
intervention_value = intervention_value, options = options,
outcome_label = outcome_label, time_label = time_label
)
}
out
}
.r4vn_tabts_base_x <- function(time) {
if (inherits(time, "Date")) return(list(x = as.numeric(time), axis = "Date", original = time))
if (inherits(time, "POSIXt")) return(list(x = as.numeric(time), axis = "POSIX", original = time))
if (is.numeric(time) || is.integer(time)) return(list(x = as.numeric(time), axis = "numeric", original = time))
list(x = seq_along(time), axis = "label", original = as.character(time))
}
.r4vn_tabts_base_axis <- function(xx, side = 1L) {
if (xx$axis == "Date") {
at <- pretty(xx$x)
graphics::axis(side, at = at, labels = format(as.Date(at, origin = "1970-01-01"), "%Y-%m"), las = 1)
} else if (xx$axis == "POSIX") {
at <- pretty(xx$x)
graphics::axis(side, at = at, labels = format(as.POSIXct(at, origin = "1970-01-01", tz = "UTC"), "%Y-%m"), las = 1)
} else if (xx$axis == "label") {
at <- unique(round(seq(1, length(xx$x), length.out = min(8L, length(xx$x)))))
graphics::axis(side, at = at, labels = xx$original[at], las = 1)
} else {
graphics::axis(side)
}
}
.r4vn_tabts_add_intervention_base <- function(time, intervention_value, options, ylim) {
if (is.null(intervention_value) || !isTRUE(options$intervention$line)) return(invisible(NULL))
ints <- tryCatch(
.r4vn_ts_match_intervention_times(intervention_value, sort(unique(time))),
error = function(e) NULL
)
if (is.null(ints) || !length(ints)) return(invisible(NULL))
xx <- .r4vn_tabts_base_x(time)
for (i in seq_along(ints)) {
iv <- ints[i]
v <- if (inherits(iv, "Date") || inherits(iv, "POSIXt")) as.numeric(iv) else if (xx$axis == "label") match(as.character(iv), xx$original) else as.numeric(iv)
graphics::abline(v = v, lty = 2, lwd = options$intervention$linewidth %||% .8,
col = options$intervention$color %||% "#8B0000")
if (isTRUE(options$intervention$label) && is.finite(v)) {
lab <- options$intervention$label_text %||% "Intervention"
if (length(lab) > 1L) lab <- lab[min(i, length(lab))]
graphics::text(v, ylim[2L], labels = lab, pos = 2, srt = 90,
cex = .75, col = options$intervention$color %||% "#8B0000", xpd = NA)
}
}
invisible(NULL)
}
.r4vn_tabts_draw_base_plot <- function(x) {
stopifnot(inherits(x, "r4vn_tabts_plot"))
tp <- x$type
o <- x$options %||% .r4vn_ts_defaults()$plot
old <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(old), add = TRUE)
graphics::par(family = o$font$family %||% "sans", las = 1,
mar = c(4.2, 4.4, 3.2, 1.2))
if (tp == "series") {
d <- x$d
ft <- x$forecast_tab
grouped <- ".group" %in% names(d) && length(unique(as.character(d$.group[!is.na(d$.group)]))) > 1L
if (grouped) {
groups <- unique(as.character(d$.group[!is.na(d$.group)]))
ng <- length(groups)
nc <- if (ng <= 2L) 1L else 2L
nr <- ceiling(ng / nc)
graphics::par(mfrow = c(nr, nc), mar = c(3.8, 4.2, 2.8, 1.0))
for (gi in seq_along(groups)) {
z <- d[as.character(d$.group) == groups[gi], , drop = FALSE]
xxg <- .r4vn_tabts_base_x(z$.time_original)
yrg <- range(c(z$.y, z$.fitted), finite = TRUE)
if (!all(is.finite(yrg))) yrg <- c(0, 1)
graphics::plot(
xxg$x, z$.y, type = "n", xaxt = "n",
xlab = x$time_label %||% "Time",
ylab = x$outcome_label %||% "Outcome", ylim = yrg,
main = paste0(o$title %||% "Interrupted time series", " - ", groups[gi])
)
.r4vn_tabts_base_axis(xxg)
graphics::lines(xxg$x, z$.y, lwd = o$observed$linewidth %||% .8,
col = o$observed$color %||% "#222222")
if (isTRUE(o$observed$point)) {
graphics::points(xxg$x, z$.y, pch = o$observed$point_shape %||% 16,
cex = (o$observed$point_size %||% 1.8) / 2,
col = o$observed$color %||% "#222222")
}
if (any(is.finite(z$.fitted))) {
graphics::lines(xxg$x, z$.fitted, lwd = o$fitted$linewidth %||% .8,
lty = 2, col = o$fitted$color %||% "#2C6EAA")
}
.r4vn_tabts_add_intervention_base(z$.time_original, x$intervention_value, o, yrg)
if (isTRUE(o$legend$show) && gi == 1L) {
graphics::legend(
"topright", legend = c("Observed", "Fitted"),
col = c(o$observed$color %||% "#222222", o$fitted$color %||% "#2C6EAA"),
lty = c(1, 2), bty = "n", cex = .75
)
}
}
return(invisible(x))
}
xx <- .r4vn_tabts_base_x(d$.time_original)
xf <- if (!is.null(ft)) .r4vn_tabts_base_x(ft$Time) else NULL
all_y <- c(d$.y, d$.fitted)
if (!is.null(ft)) all_y <- c(all_y, ft$Forecast, unlist(ft[grep("_(lower|upper)$", names(ft))], use.names = FALSE))
yr <- range(all_y, finite = TRUE)
if (!all(is.finite(yr))) yr <- c(0, 1)
xr <- range(c(xx$x, if (!is.null(xf)) xf$x else numeric()), finite = TRUE)
graphics::plot(xx$x, d$.y, type = "n", xaxt = "n", xlab = x$time_label %||% "Time",
ylab = x$outcome_label %||% "Outcome", ylim = yr, xlim = xr,
main = o$title %||% if (!is.null(ft)) "Time series and forecast" else "Time series")
.r4vn_tabts_base_axis(xx)
graphics::lines(xx$x, d$.y, lwd = o$observed$linewidth %||% .8,
lty = 1, col = o$observed$color %||% "#222222")
if (isTRUE(o$observed$point)) graphics::points(xx$x, d$.y, pch = o$observed$point_shape %||% 16,
cex = (o$observed$point_size %||% 1.8) / 2,
col = o$observed$color %||% "#222222")
if (any(is.finite(d$.fitted))) graphics::lines(xx$x, d$.fitted, lwd = o$fitted$linewidth %||% .8,
lty = 2, col = o$fitted$color %||% "#2C6EAA")
legend_names <- c("Observed")
legend_col <- c(o$observed$color %||% "#222222")
legend_lty <- 1
if (any(is.finite(d$.fitted))) { legend_names <- c(legend_names, "Fitted"); legend_col <- c(legend_col, o$fitted$color %||% "#2C6EAA"); legend_lty <- c(legend_lty, 2) }
if (!is.null(ft)) {
lower <- grep("^(PI|CI)[0-9]+_lower$", names(ft), value = TRUE)
if (isTRUE(o$pi$show) && length(lower)) {
lc <- tail(lower, 1L); uc <- sub("_lower$", "_upper", lc)
if (uc %in% names(ft)) {
graphics::polygon(c(xf$x, rev(xf$x)), c(ft[[lc]], rev(ft[[uc]])),
border = NA, col = grDevices::adjustcolor(o$forecast$color %||% "#B23A48", alpha.f = o$pi$alpha %||% .15))
}
}
graphics::lines(xf$x, ft$Forecast, lwd = o$forecast$linewidth %||% .9,
col = o$forecast$color %||% "#B23A48")
legend_names <- c(legend_names, "Forecast"); legend_col <- c(legend_col, o$forecast$color %||% "#B23A48"); legend_lty <- c(legend_lty, 1)
}
.r4vn_tabts_add_intervention_base(d$.time_original, x$intervention_value, o, yr)
if (isTRUE(o$legend$show)) graphics::legend("topright", legend = legend_names, col = legend_col,
lty = legend_lty, bty = "n", cex = .8)
} else if (tp == "decomposition") {
z <- x$decomposition
comps <- c("observed", "trend", "seasonal", "remainder")
graphics::par(mfrow = c(4, 1), mar = c(2.1, 4.3, 1.5, 1.0), oma = c(2.0, 0, 1.2, 0))
xx <- .r4vn_tabts_base_x(z$time)
for (nm in comps) {
graphics::plot(xx$x, z[[nm]], type = "l", xaxt = "n", xlab = "", ylab = nm,
lwd = o$decomposition$linewidth %||% .7,
col = o$decomposition$color %||% "#222222")
.r4vn_tabts_base_axis(xx)
}
graphics::mtext("Time-series decomposition", outer = TRUE, side = 3, line = .2, font = 2)
} else if (tp == "acf") {
z <- x$acf_tab; kind <- toupper(x$kind %||% "acf")
ci <- stats::qnorm(.975) / sqrt(x$n %||% 1)
yr <- range(c(z$Correlation, -ci, ci, 0), finite = TRUE)
graphics::plot(z$Lag, z$Correlation, type = "h", lwd = 2,
col = o[[tolower(kind)]]$color %||% "#2C6EAA",
xlab = "Lag", ylab = "Correlation", ylim = yr, main = kind)
graphics::abline(h = 0, lwd = .6)
if (isTRUE(o[[tolower(kind)]]$ci)) graphics::abline(h = c(-ci, ci), lty = 2, col = "#777777")
} else if (tp == "residual") {
d <- x$d
grouped <- ".group" %in% names(d) && length(unique(as.character(d$.group[!is.na(d$.group)]))) > 1L
groups <- if (grouped) unique(as.character(d$.group[!is.na(d$.group)])) else NA_character_
if (grouped) {
ng <- length(groups); nc <- if (ng <= 2L) 1L else 2L; nr <- ceiling(ng / nc)
graphics::par(mfrow = c(nr, nc), mar = c(3.8, 4.2, 2.8, 1.0))
for (gi in seq_along(groups)) {
z <- d[as.character(d$.group) == groups[gi], , drop = FALSE]
xxg <- .r4vn_tabts_base_x(z$.time_original)
graphics::plot(xxg$x, z$.residual, type = "n", xaxt = "n",
xlab = o$xlab %||% "Time", ylab = "Residual",
main = paste0("Residuals - ", groups[gi]))
.r4vn_tabts_base_axis(xxg); graphics::abline(h = 0, lty = 2, col = "#777777")
graphics::lines(xxg$x, z$.residual, col = o$residual$color %||% "#222222",
lwd = o$residual$linewidth %||% .7)
if (isTRUE(o$residual$point)) {
graphics::points(xxg$x, z$.residual, pch = 16, cex = .65,
col = o$residual$color %||% "#222222")
}
}
return(invisible(x))
}
xx <- .r4vn_tabts_base_x(d$.time_original)
graphics::plot(xx$x, d$.residual, type = "n", xaxt = "n", xlab = o$xlab %||% "Time",
ylab = "Residual", main = "Residual time series")
.r4vn_tabts_base_axis(xx); graphics::abline(h = 0, lty = 2, col = "#777777")
graphics::lines(xx$x, d$.residual, col = o$residual$color %||% "#222222",
lwd = o$residual$linewidth %||% .7)
if (isTRUE(o$residual$point)) graphics::points(xx$x, d$.residual, pch = 16, cex = .65,
col = o$residual$color %||% "#222222")
} else if (tp == "counterfactual") {
d <- x$d; cf <- x$counterfactual; xx <- .r4vn_tabts_base_x(cf$Time)
yr <- range(c(d$.y, cf$Fitted, cf$Counterfactual), finite = TRUE)
graphics::plot(xx$x, d$.y, type = "n", xaxt = "n", ylim = yr,
xlab = x$time_label %||% "Time", ylab = x$outcome_label %||% "Outcome",
main = o$title %||% "Interrupted time-series analysis")
.r4vn_tabts_base_axis(xx)
graphics::lines(xx$x, d$.y, col = o$observed$color %||% "#222222", lwd = o$observed$linewidth %||% .8)
graphics::points(xx$x, d$.y, pch = 16, cex = .65, col = o$observed$color %||% "#222222")
graphics::lines(xx$x, cf$Fitted, col = o$fitted$color %||% "#2C6EAA", lwd = o$fitted$linewidth %||% .9)
graphics::lines(xx$x, cf$Counterfactual, col = o$counterfactual$color %||% "#666666",
lwd = o$counterfactual$linewidth %||% .9, lty = 2)
.r4vn_tabts_add_intervention_base(cf$Time, x$intervention_value, o, yr)
if (isTRUE(o$legend$show)) graphics::legend("topright", c("Observed", "Fitted", "Counterfactual"),
col = c(o$observed$color %||% "#222222", o$fitted$color %||% "#2C6EAA", o$counterfactual$color %||% "#666666"),
lty = c(1, 1, 2), bty = "n", cex = .8)
} else {
graphics::plot.new(); graphics::title(main = "tabts plot")
}
invisible(x)
}
.r4vn_ts_make_plots <- function(
d, frequency, decomposition, acf_tab, pacf_tab, forecast_tab,
counterfactual, intervention_value, model, requested, options,
outcome_label, time_label) {
out <- list()
if ("series" %in% requested) {
out$series <- .r4vn_ts_plot_series(
d, NULL, NULL, intervention_value,
options, outcome_label, time_label
)
}
if ("forecast" %in% requested && !is.null(forecast_tab)) {
out$forecast <- .r4vn_ts_plot_series(
d, forecast_tab, NULL, intervention_value,
options, outcome_label, time_label
)
}
if ("its" %in% requested && model == "its") {
out$its <- .r4vn_ts_plot_series(
d, NULL, NULL, intervention_value,
options, outcome_label, time_label
)
}
if ("decomposition" %in% requested && !is.null(decomposition)) {
out$decomposition <- .r4vn_ts_plot_decomposition(decomposition, options)
}
if ("acf" %in% requested && !is.null(acf_tab)) {
out$acf <- .r4vn_ts_plot_acf(acf_tab, nrow(d), options, "ACF", "acf")
}
if ("pacf" %in% requested && !is.null(pacf_tab)) {
out$pacf <- .r4vn_ts_plot_acf(pacf_tab, nrow(d), options, "PACF", "pacf")
}
if ("residual" %in% requested && any(is.finite(d$.residual))) {
out$residual <- .r4vn_ts_plot_residual(d, options)
}
if ("counterfactual" %in% requested && !is.null(counterfactual)) {
out$counterfactual <- .r4vn_ts_plot_counterfactual(
d, counterfactual, intervention_value, options, outcome_label, time_label
)
}
out
}
.r4vn_ts_plot_series <- function(
d, forecast_tab, counterfactual, intervention_value,
o, outcome_label, time_label) {
grouped <- ".group" %in% names(d) && length(unique(as.character(d$.group[!is.na(d$.group)]))) > 1L
df <- data.frame(
Time = d$.time_original,
Observed = d$.y,
Fitted = d$.fitted,
stringsAsFactors = FALSE
)
if (grouped) df$Group <- factor(d$.group)
mapping <- if (grouped) {
ggplot2::aes(x = Time, y = Observed, group = Group)
} else {
ggplot2::aes(x = Time, y = Observed)
}
p <- ggplot2::ggplot(df, mapping)
if (isTRUE(o$observed$point)) {
p <- p + ggplot2::geom_point(
ggplot2::aes(color = "Observed"),
shape = o$observed$point_shape,
size = o$observed$point_size, alpha = o$observed$point_alpha,
na.rm = TRUE
)
}
p <- p + ggplot2::geom_line(
ggplot2::aes(color = "Observed", linetype = "Observed"),
linewidth = o$observed$linewidth, alpha = o$observed$alpha,
na.rm = TRUE
)
has_fitted <- any(is.finite(df$Fitted))
if (has_fitted) {
p <- p + ggplot2::geom_line(
ggplot2::aes(y = Fitted, color = "Fitted", linetype = "Fitted"),
linewidth = o$fitted$linewidth, alpha = o$fitted$alpha,
na.rm = TRUE
)
}
has_forecast <- !is.null(forecast_tab)
if (has_forecast) {
ft <- forecast_tab
lower_cols <- grep("^(PI|CI)[0-9]+_lower$", names(ft), value = TRUE)
upper_cols <- grep("^(PI|CI)[0-9]+_upper$", names(ft), value = TRUE)
if (isTRUE(o$pi$show) && length(lower_cols) && length(upper_cols)) {
tags <- as.numeric(gsub("\\D", "", lower_cols))
ord <- order(tags, decreasing = TRUE)
for (j in ord) {
lc <- lower_cols[j]
uc <- sub("_lower$", "_upper", lc)
if (!uc %in% names(ft)) next
rib <- data.frame(Time = ft$Time, ymin = ft[[lc]], ymax = ft[[uc]])
p <- p + ggplot2::geom_ribbon(
data = rib,
ggplot2::aes(x = Time, ymin = ymin, ymax = ymax),
inherit.aes = FALSE,
fill = o$forecast$color,
alpha = o$pi$alpha,
color = if (isTRUE(o$pi$border)) o$forecast$color else NA,
linewidth = o$pi$border_linewidth
)
}
}
p <- p + ggplot2::geom_line(
data = data.frame(Time = ft$Time, Forecast = ft$Forecast),
ggplot2::aes(x = Time, y = Forecast, color = "Forecast", linetype = "Forecast"),
inherit.aes = FALSE,
linewidth = o$forecast$linewidth, alpha = o$forecast$alpha
)
}
vals_col <- c(Observed = o$observed$color)
vals_lty <- c(Observed = o$observed$linetype)
if (has_fitted) {
vals_col <- c(vals_col, Fitted = o$fitted$color)
vals_lty <- c(vals_lty, Fitted = o$fitted$linetype)
}
if (has_forecast) {
vals_col <- c(vals_col, Forecast = o$forecast$color)
vals_lty <- c(vals_lty, Forecast = o$forecast$linetype)
}
p <- p +
ggplot2::scale_color_manual(values = vals_col, name = o$legend$title) +
ggplot2::scale_linetype_manual(values = vals_lty, name = o$legend$title)
p <- .r4vn_ts_add_interventions(p, intervention_value, o, d$.time_original)
if (grouped && isTRUE(o$facet$show)) {
p <- p + ggplot2::facet_wrap(
~ Group, ncol = o$facet$ncol, nrow = o$facet$nrow,
scales = o$facet$scales %||% "fixed"
)
}
.r4vn_ts_apply_plot_options(p, o, outcome_label, time_label)
}
.r4vn_ts_plot_decomposition <- function(x, o) {
nm <- c("observed", "trend", "seasonal", "remainder")
pieces <- lapply(nm, function(v) {
data.frame(Time = x$time, Component = v, Value = x[[v]])
})
z <- do.call(rbind, pieces)
z$Component <- factor(z$Component, levels = nm)
p <- ggplot2::ggplot(z, ggplot2::aes(x = Time, y = Value)) +
ggplot2::geom_line(
color = o$decomposition$color,
linewidth = o$decomposition$linewidth,
na.rm = TRUE
) +
ggplot2::facet_wrap(
~ Component, ncol = 1,
scales = if (isTRUE(o$decomposition$free_y)) "free_y" else "fixed"
) +
ggplot2::labs(x = o$xlab %||% "Time", y = NULL)
.r4vn_ts_apply_theme(p, o)
}
.r4vn_ts_plot_acf <- function(x, n, o, title, kind = c("acf", "pacf")) {
kind <- match.arg(kind)
s <- o[[kind]]
ci <- stats::qnorm(1 - (1 - s$ci_level) / 2) / sqrt(n)
p <- ggplot2::ggplot(x, ggplot2::aes(x = Lag, y = Correlation)) +
ggplot2::geom_hline(yintercept = 0, linewidth = .4) +
ggplot2::geom_segment(
ggplot2::aes(xend = Lag, y = 0, yend = Correlation),
color = s$color, linewidth = s$linewidth
) +
ggplot2::labs(title = title, x = "Lag", y = "Correlation")
if (isTRUE(s$ci)) {
p <- p + ggplot2::geom_hline(
yintercept = c(-ci, ci), linetype = "dashed", linewidth = .4
)
}
.r4vn_ts_apply_theme(p, o)
}
.r4vn_ts_plot_residual <- function(d, o) {
grouped <- ".group" %in% names(d) && length(unique(as.character(d$.group[!is.na(d$.group)]))) > 1L
z <- data.frame(Time = d$.time_original, Residual = d$.residual, stringsAsFactors = FALSE)
if (grouped) z$Group <- factor(d$.group)
mapping <- if (grouped) {
ggplot2::aes(x = Time, y = Residual, group = Group)
} else {
ggplot2::aes(x = Time, y = Residual)
}
p <- ggplot2::ggplot(z, mapping) +
ggplot2::geom_hline(yintercept = 0, linetype = "dashed", linewidth = .4) +
ggplot2::geom_line(
color = o$residual$color, linewidth = o$residual$linewidth, na.rm = TRUE
)
if (isTRUE(o$residual$point)) {
p <- p + ggplot2::geom_point(
color = o$residual$color, size = o$residual$point_size, na.rm = TRUE
)
}
if (grouped && isTRUE(o$facet$show)) {
p <- p + ggplot2::facet_wrap(
~ Group, ncol = o$facet$ncol, nrow = o$facet$nrow,
scales = o$facet$scales %||% "fixed"
)
}
p <- p + ggplot2::labs(x = o$xlab %||% "Time", y = "Residual")
.r4vn_ts_apply_theme(p, o)
}
.r4vn_ts_plot_counterfactual <- function(
d, cf, intervention_value, o, outcome_label, time_label) {
z <- data.frame(
Time = cf$Time,
Observed = d$.y,
Fitted = cf$Fitted,
Counterfactual = cf$Counterfactual
)
p <- ggplot2::ggplot(z, ggplot2::aes(x = Time, y = Observed)) +
ggplot2::geom_point(
ggplot2::aes(color = "Observed"),
size = o$observed$point_size, shape = o$observed$point_shape, na.rm = TRUE
) +
ggplot2::geom_line(
ggplot2::aes(color = "Observed", linetype = "Observed"),
linewidth = o$observed$linewidth, na.rm = TRUE
) +
ggplot2::geom_line(
ggplot2::aes(y = Fitted, color = "Fitted", linetype = "Fitted"),
linewidth = o$fitted$linewidth, na.rm = TRUE
) +
ggplot2::geom_line(
ggplot2::aes(y = Counterfactual, color = "Counterfactual",
linetype = "Counterfactual"),
linewidth = o$counterfactual$linewidth,
alpha = o$counterfactual$alpha, na.rm = TRUE
) +
ggplot2::scale_color_manual(
values = c(
Observed = o$observed$color,
Fitted = o$fitted$color,
Counterfactual = o$counterfactual$color
),
name = o$legend$title
) +
ggplot2::scale_linetype_manual(
values = c(
Observed = o$observed$linetype,
Fitted = o$fitted$linetype,
Counterfactual = o$counterfactual$linetype
),
name = o$legend$title
)
if (isTRUE(o$counterfactual$shade)) {
p <- p + ggplot2::geom_ribbon(
ggplot2::aes(ymin = pmin(Fitted, Counterfactual),
ymax = pmax(Fitted, Counterfactual)),
inherit.aes = TRUE,
fill = o$counterfactual$color, alpha = o$counterfactual$shade_alpha,
color = NA
)
}
p <- .r4vn_ts_add_interventions(p, intervention_value, o, d$.time_original)
.r4vn_ts_apply_plot_options(p, o, outcome_label, time_label)
}
.r4vn_ts_add_interventions <- function(p, intervention_value, o, time) {
if (is.null(intervention_value) || !isTRUE(o$intervention$line)) return(p)
ints <- tryCatch(
.r4vn_ts_match_intervention_times(intervention_value, sort(unique(time))),
error = function(e) NULL
)
if (is.null(ints)) return(p)
for (i in seq_along(ints)) {
v <- ints[i]
p <- p + ggplot2::geom_vline(
xintercept = v,
color = o$intervention$color,
linewidth = o$intervention$linewidth,
linetype = o$intervention$linetype,
alpha = o$intervention$alpha
)
if (isTRUE(o$intervention$label)) {
lab <- o$intervention$label_text
if (length(lab) > 1L) lab <- lab[min(i, length(lab))]
p <- p + ggplot2::annotate(
"text", x = v, y = Inf, label = lab,
angle = o$intervention$label_angle,
vjust = 1.2, hjust = 1.05,
color = o$intervention$color
)
}
}
p
}
.r4vn_ts_apply_plot_options <- function(p, o, outcome_label, time_label) {
p <- p + ggplot2::labs(
title = o$title,
subtitle = o$subtitle,
x = o$xlab %||% time_label,
y = o$ylab %||% outcome_label
)
ax <- o$axis
is_date <- !is.null(p$data$Time) && inherits(p$data$Time, "Date")
if (is_date && (!is.null(ax$date_breaks) || !is.null(ax$date_labels) ||
!is.null(ax$x_breaks))) {
args <- list()
if (!is.null(ax$date_breaks)) args$date_breaks <- ax$date_breaks
if (!is.null(ax$date_labels)) args$date_labels <- ax$date_labels
if (!is.null(ax$x_breaks)) args$breaks <- ax$x_breaks
p <- p + do.call(ggplot2::scale_x_date, args)
} else if (!is_date && !is.null(ax$x_breaks)) {
p <- p + ggplot2::scale_x_continuous(breaks = ax$x_breaks)
}
if (isTRUE(ax$y_log)) {
p <- p + ggplot2::scale_y_log10(
breaks = ax$y_breaks %||% ggplot2::waiver()
)
} else if (!is.null(ax$y_breaks)) {
p <- p + ggplot2::scale_y_continuous(breaks = ax$y_breaks)
}
if (!is.null(ax$xlim) || !is.null(ax$ylim)) {
p <- p + ggplot2::coord_cartesian(xlim = ax$xlim, ylim = ax$ylim)
}
if (!is.null(o$reference$yline)) {
p <- p + ggplot2::geom_hline(
yintercept = o$reference$yline,
color = o$reference$color,
linewidth = o$reference$linewidth,
linetype = o$reference$linetype
)
}
if (!is.null(o$reference$xline)) {
p <- p + ggplot2::geom_vline(
xintercept = o$reference$xline,
color = o$reference$color,
linewidth = o$reference$linewidth,
linetype = o$reference$linetype
)
}
if (!is.null(o$annotation) && is.data.frame(o$annotation) &&
all(c("x", "y", "label") %in% names(o$annotation))) {
p <- p + ggplot2::geom_text(
data = o$annotation,
ggplot2::aes(x = x, y = y, label = label),
inherit.aes = FALSE
)
}
.r4vn_ts_apply_theme(p, o)
}
.r4vn_ts_apply_theme <- function(p, o) {
family <- o$font$family %||% ""
th <- switch(
o$theme,
minimal = ggplot2::theme_minimal(base_size = o$font$base_size, base_family = family),
classic = ggplot2::theme_classic(base_size = o$font$base_size, base_family = family),
bw = ggplot2::theme_bw(base_size = o$font$base_size, base_family = family),
r4vn = ggplot2::theme_minimal(base_size = o$font$base_size, base_family = family)
)
p <- p + th + ggplot2::theme(
plot.title = ggplot2::element_text(size = o$font$title_size),
plot.subtitle = ggplot2::element_text(size = o$font$subtitle_size),
axis.title = ggplot2::element_text(size = o$font$axis_title_size),
axis.text = ggplot2::element_text(size = o$font$axis_text_size),
legend.text = ggplot2::element_text(size = o$font$legend_text_size),
panel.grid.minor = if (isTRUE(o$grid$minor)) ggplot2::element_line() else ggplot2::element_blank(),
panel.grid.major = if (isTRUE(o$grid$major)) ggplot2::element_line() else ggplot2::element_blank()
)
if (o$theme == "r4vn") {
p <- p + ggplot2::theme(plot.title.position = "plot")
}
if (!isTRUE(o$legend$show) || identical(o$legend$position, "none")) {
p <- p + ggplot2::theme(legend.position = "none")
} else {
p <- p + ggplot2::theme(legend.position = o$legend$position)
}
p
}
.r4vn_ts_label <- function(x, fallback) {
lab <- attr(x, "label", exact = TRUE)
if (is.null(lab) || !nzchar(as.character(lab)[1L])) fallback else as.character(lab)[1L]
}
.r4vn_ts_interpret <- function(
summary_tab, stationarity_tab, model_selection, chosen_model,
coefficients, diagnostics, forecast_tab, its_meta, effect_at,
language, detail, options, digits, p_digits) {
alpha <- options$alpha %||% .05
rows <- list()
add <- function(section, finding, evidence, status, source) {
rows[[length(rows) + 1L]] <<- data.frame(
Section = section,
Finding = finding,
Evidence = evidence,
Status = status,
Source = source,
stringsAsFactors = FALSE
)
}
fmt <- function(x) formatC(x, digits = digits, format = "f")
fmtp <- function(p) {
if (!is.finite(p)) return(NA_character_)
cut <- 10^-p_digits
if (p < cut) paste0("<", formatC(cut, digits = p_digits, format = "f"))
else formatC(p, digits = p_digits, format = "f")
}
# Assumptions/stationarity.
if (isTRUE(options$include_assumptions) && !is.null(stationarity_tab) && nrow(stationarity_tab)) {
adf <- stationarity_tab[stationarity_tab$Test == "ADF", , drop = FALSE]
kpss <- stationarity_tab[stationarity_tab$Test == "KPSS", , drop = FALSE]
if (nrow(adf) && nrow(kpss)) {
evidence <- paste0(
"ADF=", fmt(adf$Statistic[1]), ", p=", fmtp(adf$p[1]),
"; KPSS=", fmt(kpss$Statistic[1]), ", p=", fmtp(kpss$p[1])
)
if (adf$p[1] < alpha && kpss$p[1] >= alpha) {
finding <- "The tests provide evidence consistent with stationarity in the analysed series."
add("Stationarity", finding, evidence, "Supported", "stationarity")
} else if (adf$p[1] >= alpha && kpss$p[1] < alpha) {
finding <- "The tests provide evidence against stationarity; differencing or explicit trend modelling should be considered."
add("Stationarity", finding, evidence, "Not supported", "stationarity")
} else {
finding <- "ADF and KPSS do not provide fully concordant evidence about stationarity; interpretation should be cautious."
add("Stationarity", finding, evidence, "Mixed", "stationarity")
}
}
}
# Model selection.
if (isTRUE(options$include_model) && !is.null(model_selection) && nrow(model_selection)) {
row <- model_selection[tolower(model_selection$Model) == tolower(chosen_model), , drop = FALSE]
if (nrow(row)) {
evidence <- paste0(
"Model=", row$Model[1],
", AICc=", fmt(row$AICc[1]),
if ("RMSE" %in% names(row) && is.finite(row$RMSE[1])) paste0(", RMSE=", fmt(row$RMSE[1])) else "",
if (is.finite(row$Ljung_Box_p[1])) paste0(", Ljung-Box p=", fmtp(row$Ljung_Box_p[1])) else ""
)
finding <- paste0(
"The ", row$Model[1],
" model was selected using the requested selection criterion and available diagnostics."
)
add("Model selection", finding, evidence, "Selected", "model_selection")
}
}
# Residual autocorrelation.
if (!is.null(diagnostics) && nrow(diagnostics)) {
lb <- diagnostics[diagnostics$Diagnostic == "Ljung-Box", , drop = FALSE]
if (nrow(lb) && is.finite(lb$p[1])) {
evidence <- paste0(
"Ljung-Box Q=", fmt(lb$Statistic[1]),
", df=", fmt(lb$df[1]), ", p=", fmtp(lb$p[1])
)
if (lb$p[1] >= alpha) {
finding <- "There is no detected evidence of substantial remaining residual autocorrelation."
add("Residual diagnostics", finding, evidence, "Supported", "diagnostics")
} else {
finding <- "Residual autocorrelation remains detectable; the model may not fully capture temporal dependence."
add("Residual diagnostics", finding, evidence, "Warning", "diagnostics")
}
}
}
# ITS effects.
if (!is.null(its_meta) && !is.null(coefficients) && nrow(coefficients)) {
is_count <- its_meta$family %in% c("poisson", "quasipoisson", "negativebinomial")
int_rows <- grep("^\\.int[0-9]+$", coefficients$Term)
after_rows <- grep("^\\.after[0-9]+$", coefficients$Term)
for (ii in int_rows) {
r <- coefficients[ii, ]
evidence <- paste0(
r$Measure, "=", fmt(r$Estimate),
", 95% CI ", fmt(r$CI_lower), " to ", fmt(r$CI_upper),
", p=", fmtp(r$p)
)
if (is_count) {
pct <- (r$Estimate - 1) * 100
finding <- paste0(
if (isTRUE(its_meta$controlled)) paste0(
"In the reference group (", its_meta$reference_group %||% "first group", "), "
) else "",
"immediately after the intervention, the estimated outcome rate ",
if (pct < 0) "decreased by " else "increased by ",
fmt(abs(pct)), "% relative to the model-implied pre-intervention expectation."
)
} else {
finding <- paste0(
if (isTRUE(its_meta$controlled)) paste0(
"In the reference group (", its_meta$reference_group %||% "first group", "), "
) else "The ITS model ",
"estimates an immediate level change of ",
fmt(r$Estimate), " units."
)
}
add("ITS effect", finding, evidence,
if (is.finite(r$p) && r$p < alpha) "Evidence" else "Uncertain",
"coefficients")
}
for (ii in after_rows) {
r <- coefficients[ii, ]
evidence <- paste0(
r$Measure, "=", fmt(r$Estimate),
", 95% CI ", fmt(r$CI_lower), " to ", fmt(r$CI_upper),
", p=", fmtp(r$p)
)
if (is_count) {
pct <- (r$Estimate - 1) * 100
finding <- paste0(
if (isTRUE(its_meta$controlled)) paste0(
"In the reference group (", its_meta$reference_group %||% "first group", "), the "
) else "The ",
"post-intervention trend changed by about ", fmt(abs(pct)),
"% per time unit ", if (pct < 0) "downward" else "upward",
" relative to the pre-intervention trend."
)
} else {
finding <- paste0(
if (isTRUE(its_meta$controlled)) paste0(
"In the reference group (", its_meta$reference_group %||% "first group", "), the "
) else "The ",
"post-intervention trend change is estimated at ",
fmt(r$Estimate), " units per time unit."
)
}
add("ITS effect", finding, evidence,
if (is.finite(r$p) && r$p < alpha) "Evidence" else "Uncertain",
"coefficients")
}
if (isTRUE(its_meta$controlled) && !is.null(effect_at) && nrow(effect_at)) {
for (jj in seq_len(nrow(effect_at))) {
r <- effect_at[jj, , drop = FALSE]
evidence <- paste0(
r$Measure[1], "=", fmt(r$Estimate[1]),
", 95% CI ", fmt(r$CI_lower[1]), " to ", fmt(r$CI_upper[1]),
", p=", fmtp(r$p[1])
)
if (is_count) {
pct <- r$Relative_difference_pct[1]
finding <- paste0(
"At ", r$Period_after_intervention[1],
" period(s) after intervention, the controlled contrast ",
r$Comparison[1], " corresponds to an estimated ",
if (is.finite(pct) && pct < 0) "lower" else "higher",
" intervention-associated rate by ", fmt(abs(pct)), "% relative to the reference group."
)
} else {
finding <- paste0(
"At ", r$Period_after_intervention[1],
" period(s) after intervention, the controlled contrast ",
r$Comparison[1], " is ", fmt(r$Estimate[1]), " outcome units."
)
}
add("Controlled ITS effect", finding, evidence,
if (is.finite(r$p[1]) && r$p[1] < alpha) "Evidence" else "Uncertain",
"effect_at")
}
}
}
# Forecast.
if (isTRUE(options$include_forecast) && !is.null(forecast_tab) && nrow(forecast_tab)) {
last <- forecast_tab[nrow(forecast_tab), , drop = FALSE]
pi_cols <- grep("^PI[0-9]+_(lower|upper)$", names(last), value = TRUE)
ci_cols <- grep("^CI[0-9]+_(lower|upper)$", names(last), value = TRUE)
interval_cols <- if (length(pi_cols)) pi_cols else ci_cols
evidence <- paste0(
"Horizon=", nrow(forecast_tab),
", forecast=", fmt(last$Forecast)
)
if (length(interval_cols) >= 2L) {
low <- interval_cols[grepl("_lower$", interval_cols)]
up <- interval_cols[grepl("_upper$", interval_cols)]
if (length(low) && length(up)) {
evidence <- paste0(
evidence, ", interval=", fmt(last[[low[length(low)]]]),
" to ", fmt(last[[up[length(up)]]])
)
}
}
finding <- "The longest-horizon forecast should be interpreted together with its uncertainty interval rather than as a certain point value."
add("Forecast", finding, evidence, "Forecast", "forecast")
}
# Limitations based on actual numerical conditions.
if (isTRUE(options$include_limitations)) {
n <- summary_tab$N_observed[1]
if (is.finite(n) && n < 24L) {
finding <- "The series contains relatively few observed time points, so seasonality estimates and forecasts may be unstable."
add("Limitation", finding, paste0("N observed=", n), "Caution", "descriptive")
}
if (!is.null(its_meta) && is.finite(its_meta$n_post) &&
its_meta$n_post < (its_meta$post_min %||% 8L)) {
finding <- "Few post-intervention time points are available; the post-intervention trend estimate should be interpreted cautiously."
add("Limitation", finding, paste0("Post-intervention time points=", its_meta$n_post),
"Caution", "its")
}
}
if (!length(rows)) return(NULL)
out <- do.call(rbind, rows)
rownames(out) <- NULL
if (detail == "brief" && nrow(out) > 4L) out <- out[seq_len(4L), , drop = FALSE]
out
}
.r4vn_ts_ai <- function(x, ai) {
if (identical(ai, FALSE) || is.null(ai)) return(NULL)
if (!exists("aiask", mode = "function", inherits = TRUE)) {
warning("`ai` was requested but R4VN `aiask()` is not available.", call. = FALSE)
return(NULL)
}
fun <- get("aiask", mode = "function", inherits = TRUE)
fml <- names(formals(fun))
args <- list(x)
if (is.character(ai) && length(ai) == 1L && "api" %in% fml) args$api <- ai
tryCatch(do.call(fun, args), error = function(e) {
warning("AI interpretation failed: ", conditionMessage(e), call. = FALSE)
NULL
})
}
.r4vn_ts_label_df <- function(df, digits = 2, p_digits = 3) {
if (is.null(df) || !is.data.frame(df)) return(df)
out <- df
for (nm in names(out)) {
if (is.numeric(out[[nm]])) {
if (tolower(nm) == "p" || grepl("_p$", tolower(nm))) {
out[[nm]] <- vapply(out[[nm]], function(v) {
if (!is.finite(v)) return("")
cut <- 10^-p_digits
if (v < cut) paste0("<", formatC(cut, digits = p_digits, format = "f"))
else formatC(v, digits = p_digits, format = "f")
}, character(1))
} else {
out[[nm]] <- ifelse(
is.finite(out[[nm]]),
formatC(out[[nm]], digits = digits, format = "f"),
""
)
}
}
}
out
}
.r4vn_ts_html_escape <- function(x) {
x <- as.character(x)
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub('"', """, x, fixed = TRUE)
x
}
.r4vn_tabts_transpose_table <- function(df) {
if (!is.data.frame(df) || !nrow(df)) return(df)
if (!(nrow(df) <= 5L && ncol(df) >= 9L && ncol(df) > nrow(df) + 4L)) return(df)
row_names <- if ("Statistic" %in% names(df)) as.character(df$Statistic) else {
rn <- rownames(df)
if (is.null(rn) || identical(rn, as.character(seq_len(nrow(df))))) paste0("Result ", seq_len(nrow(df))) else rn
}
keep <- if ("Statistic" %in% names(df)) setdiff(names(df), "Statistic") else names(df)
m <- t(as.matrix(df[, keep, drop = FALSE]))
out <- data.frame(Statistic = rownames(m), m, check.names = FALSE, stringsAsFactors = FALSE)
names(out)[-1L] <- make.unique(row_names)
rownames(out) <- NULL
out
}
.r4vn_ts_df_html <- function(df, title = NULL, digits = 2, p_digits = 3,
smart_transpose = TRUE) {
if (is.null(df) || !nrow(df)) return("")
if (isTRUE(smart_transpose)) df <- .r4vn_tabts_transpose_table(df)
z <- .r4vn_ts_label_df(df, digits, p_digits)
th <- paste0("<th>", .r4vn_ts_html_escape(names(z)), "</th>", collapse = "")
trs <- apply(z, 1, function(r) {
paste0("<tr>", paste0("<td>", .r4vn_ts_html_escape(r), "</td>", collapse = ""), "</tr>")
})
paste0(
if (!is.null(title)) paste0("<h3>", .r4vn_ts_html_escape(title), "</h3>") else "",
"<div class='table-wrap'><table><thead><tr>", th, "</tr></thead><tbody>",
paste(trs, collapse = ""), "</tbody></table></div>"
)
}
.r4vn_tabts_base64 <- function(x) {
bytes <- as.integer(x)
if (!length(bytes)) return("")
alphabet <- strsplit(
"ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/",
"", fixed = TRUE
)[[1L]]
padding <- (3L - length(bytes) %% 3L) %% 3L
if (padding) bytes <- c(bytes, rep.int(0L, padding))
z <- matrix(bytes, ncol = 3L, byrow = TRUE)
code <- cbind(
bitwShiftR(z[, 1L], 2L),
bitwOr(bitwShiftL(bitwAnd(z[, 1L], 3L), 4L), bitwShiftR(z[, 2L], 4L)),
bitwOr(bitwShiftL(bitwAnd(z[, 2L], 15L), 2L), bitwShiftR(z[, 3L], 6L)),
bitwAnd(z[, 3L], 63L)
)
encoded <- as.vector(t(matrix(alphabet[code + 1L], ncol = 4L)))
if (padding) encoded[(length(encoded) - padding + 1L):length(encoded)] <- "="
paste0(encoded, collapse = "")
}
.r4vn_tabts_draw_plot_object <- function(p) {
if (inherits(p, "r4vn_tabts_plot")) {
.r4vn_tabts_draw_base_plot(p)
} else {
print(p)
}
invisible(p)
}
.r4vn_tabts_plot_png_html <- function(p, alt = "tabts plot",
width = 1200L, height = 760L,
res = 144L) {
path <- tempfile(fileext = ".png")
on.exit(unlink(path), add = TRUE)
args <- list(filename = path, width = width, height = height, units = "px",
res = res, bg = "white")
if (isTRUE(capabilities("cairo"))) args$type <- "cairo-png"
do.call(grDevices::png, args)
ok <- tryCatch({
.r4vn_tabts_draw_plot_object(p)
TRUE
}, error = function(e) FALSE, finally = grDevices::dev.off())
if (!ok || !file.exists(path) || !is.finite(file.info(path)$size) || file.info(path)$size <= 0) return("")
bytes <- readBin(path, what = "raw", n = file.info(path)$size)
paste0(
"<img class='ts-plot-image' alt='", .r4vn_ts_html_escape(alt),
"' src='data:image/png;base64,", .r4vn_tabts_base64(bytes), "'>"
)
}
.r4vn_tabts_plot_title <- function(name) {
c(
series = "Observed and fitted time series",
decomposition = "Time-series decomposition",
acf = "Autocorrelation function (ACF)",
pacf = "Partial autocorrelation function (PACF)",
residual = "Residual time series",
forecast = "Forecast",
its = "Interrupted time-series model",
counterfactual = "Observed, fitted and counterfactual series"
)[[name]] %||% gsub("_", " ", name, fixed = TRUE)
}
.r4vn_tabts_result_body <- function(x, heading = NULL) {
digits <- x$digits %||% 2L
p_digits <- x$p_digits %||% 3L
tabs <- x$tables
titles <- c(
descriptive = "Descriptive statistics",
stationarity = "Stationarity assessment",
model_selection = "Model selection",
model = "Final model specification",
coefficients = "Model estimates",
diagnostics = "Residual diagnostics",
validation = "Out-of-sample validation",
forecast = "Forecast",
effect_at = if (!is.null(x$its) && isTRUE(x$its$controlled))
"Controlled intervention effects at requested times" else
"Counterfactual effects at requested times",
interpretation = "Suggested interpretation"
)
body <- character()
if (!is.null(heading)) body <- c(body, paste0("<h2>", .r4vn_ts_html_escape(heading), "</h2>"))
body <- c(body, paste0(
"<p class='sub'>Model: <b>", .r4vn_ts_html_escape(toupper(x$model_name)),
"</b> • Seasonal period: ", x$period,
" • N: ", if (!is.null(x$summary$N_observed)) x$summary$N_observed[1] else nrow(x$data),
"</p>"
))
for (nm in names(tabs)) {
if (is.data.frame(tabs[[nm]]) && nrow(tabs[[nm]])) {
body <- c(body, .r4vn_ts_df_html(
tabs[[nm]], titles[[nm]] %||% nm,
digits = digits, p_digits = p_digits
))
}
}
if (length(x$plots)) {
body <- c(body, "<h2 class='plots-heading'>Figures</h2>")
for (nm in names(x$plots)) {
img <- tryCatch(
.r4vn_tabts_plot_png_html(x$plots[[nm]], .r4vn_tabts_plot_title(nm)),
error = function(e) ""
)
if (nzchar(img)) {
body <- c(body, paste0(
"<section class='chart'><h3>", .r4vn_ts_html_escape(.r4vn_tabts_plot_title(nm)),
"</h3>", img, "</section>"
))
}
}
}
body
}
.r4vn_tabts_wrap_html <- function(body, subtitle = NULL) {
paste0(
"<!doctype html><html><head><meta charset='utf-8'>",
"<meta name='viewport' content='width=device-width,initial-scale=1'>",
"<style>",
"body{font-family:Arial,Helvetica,sans-serif;margin:24px;color:#222;line-height:1.4;max-width:1180px}",
"h1{font-size:23px;margin-bottom:3px}h2{font-size:18px;margin:28px 0 8px}h3{margin:22px 0 8px;font-size:15px}",
".sub{color:#666;margin-top:0}.table-wrap{overflow-x:auto;max-width:100%}",
"table{border-collapse:collapse;width:100%;font-size:12.5px;min-width:560px}",
"th{background:#f3f5f7;text-align:right;border-top:2px solid #333;border-bottom:1px solid #777;padding:7px 8px;white-space:nowrap}",
"td{border-bottom:1px solid #ddd;padding:7px 8px;vertical-align:top;text-align:right;white-space:nowrap}",
"th:first-child,td:first-child{text-align:left}tr:nth-child(even){background:#fafafa}",
".chart{max-width:960px;margin:18px 0 30px;page-break-inside:avoid}.chart img{display:block;width:100%;height:auto;max-width:960px;background:white}",
".plots-heading{border-top:1px solid #ddd;padding-top:18px}.group-block{border-top:3px solid #222;margin-top:34px;padding-top:4px}",
"</style></head><body><h1>R4VN tabts</h1>",
if (!is.null(subtitle)) paste0("<p class='sub'>", .r4vn_ts_html_escape(subtitle), "</p>") else "",
paste(body, collapse = "\n"), "</body></html>"
)
}
.r4vn_tabts_build_html <- function(x) {
if (inherits(x, "r4vn_tabts_grouped")) {
body <- character()
for (nm in names(x$results)) {
block <- .r4vn_tabts_result_body(x$results[[nm]], paste0("Group: ", nm))
body <- c(body, "<div class='group-block'>", block, "</div>")
}
return(.r4vn_tabts_wrap_html(body, paste0("Grouped analysis by ", x$group)))
}
.r4vn_tabts_wrap_html(.r4vn_tabts_result_body(x), "Time-series analysis")
}
.r4vn_tabts_open_viewer <- function(html) {
f <- tempfile(pattern = "r4vn-tabts-", fileext = ".html")
writeLines(enc2utf8(html), f, useBytes = TRUE)
viewer <- getOption("viewer")
if (is.function(viewer)) {
try(viewer(f), silent = TRUE)
} else if (requireNamespace("rstudioapi", quietly = TRUE) &&
isTRUE(tryCatch(rstudioapi::isAvailable(), error = function(e) FALSE))) {
try(rstudioapi::viewer(f), silent = TRUE)
} else {
try(utils::browseURL(f), silent = TRUE)
}
invisible(f)
}
.r4vn_ts_primary_plot <- function(x) {
if (inherits(x, "r4vn_tabts_grouped")) {
if (!length(x$results)) return(NULL)
return(.r4vn_ts_primary_plot(x$results[[1L]]))
}
x$plots$forecast %||% x$plots$counterfactual %||% x$plots$its %||% x$plots$series %||%
(if (length(x$plots)) x$plots[[1L]] else NULL)
}
.r4vn_ts_show <- function(x) {
if (!inherits(x, "r4vn_tabts")) return(invisible(x))
p <- .r4vn_ts_primary_plot(x)
if (!is.null(p)) try(.r4vn_tabts_draw_plot_object(p), silent = TRUE)
.r4vn_tabts_open_viewer(.r4vn_tabts_build_html(x))
invisible(x)
}
.r4vn_ts_show_grouped <- function(x) {
if (!inherits(x, "r4vn_tabts_grouped")) return(invisible(x))
p <- .r4vn_ts_primary_plot(x)
if (!is.null(p)) try(.r4vn_tabts_draw_plot_object(p), silent = TRUE)
.r4vn_tabts_open_viewer(.r4vn_tabts_build_html(x))
invisible(x)
}
#' @export
print.r4vn_tabts <- function(x, ..., digits = 3) {
cat("\nR4VN time-series analysis\n")
cat("Outcome:", x$outcome, "\n")
cat("Time:", x$time, "\n")
cat("Seasonal period:", x$period, "\n")
cat("Model:", toupper(x$model_name), "\n\n")
for (nm in names(x$tables)) {
cat("== ", gsub("_", " ", nm, fixed = TRUE), " ==\n", sep = "")
print(.r4vn_ts_label_df(x$tables[[nm]], digits = digits, p_digits = 3), row.names = FALSE)
cat("\n")
}
invisible(x)
}
#' @export
print.r4vn_tabts_grouped <- function(x, ...) {
cat("\nR4VN grouped time-series analysis\n")
cat("Grouping variable:", x$group, "\n")
cat("Groups:", paste(names(x$results), collapse = ", "), "\n\n")
for (nm in names(x$results)) {
cat("---- ", nm, " ----\n", sep = "")
print(x$results[[nm]], ...)
}
invisible(x)
}
#' @export
summary.r4vn_tabts <- function(object, ...) {
list(
descriptive = object$summary,
model = object$model_name,
model_info = object$model_info,
model_selection = object$model_selection,
coefficients = object$coefficients,
diagnostics = object$diagnostics,
validation = object$validation,
forecast = object$forecast,
interpretation = object$interpretation
)
}
#' Plot an R4VN tabts result
#'
#' @param x An `r4vn_tabts` object.
#' @param which Plot name such as `"series"`, `"forecast"`, `"acf"`,
#' `"pacf"`, `"decomposition"`, `"residual"`, `"its"`, or
#' `"counterfactual"`.
#' @param type Optional alias for `which`.
#' @param file Optional output path. Supported extensions are `.png`, `.jpg`,
#' `.jpeg`, `.tif`, `.tiff`, `.pdf`, and `.svg`.
#' @param width,height Figure width and height in inches when saving.
#' @param dpi Raster resolution when saving PNG/JPEG/TIFF files.
#' @param ... Reserved for future use.
#' @return Invisibly returns the plot object that is displayed or saved.
#' @examples
#' \donttest{
#' data(dengue_ts)
#' x <- tabts(cases, time = month, data = dengue_ts,
#' model = "arima", order = c(1, 0, 1), forecast = 6,
#' show = FALSE)
#' plot(x, "series")
#' plot(x, type = "forecast")
#' plot(x, "forecast", file = file.path(tempdir(), "forecast.png"),
#' width = 8, height = 5, dpi = 300)
#' }
#' @export
plot.r4vn_tabts <- function(
x,
which = c("series", "forecast", "its", "counterfactual",
"decomposition", "acf", "pacf", "residual"),
type = NULL,
file = NULL,
width = NULL,
height = NULL,
dpi = NULL,
...) {
if (!is.null(type)) which <- type
which <- match.arg(which, c("series", "forecast", "its", "counterfactual",
"decomposition", "acf", "pacf", "residual"))
p <- x$plots[[which]]
if (is.null(p)) {
stop("Plot `", which, "` is not available in this result.", call. = FALSE)
}
if (!is.null(file)) {
.r4vn_tabts_save_plot_object(
p, file = file,
width = width %||% x$export_size$width %||% 8,
height = height %||% x$export_size$height %||% 5,
dpi = dpi %||% x$export_size$dpi %||% 300
)
} else {
.r4vn_tabts_draw_plot_object(p)
}
invisible(p)
}
.r4vn_tabts_save_plot_object <- function(p, file, width = 8, height = 5, dpi = 300) {
if (!is.character(file) || length(file) != 1L || !nzchar(file)) {
stop("`file` must be one valid output path.", call. = FALSE)
}
ext <- tolower(tools::file_ext(file))
if (!ext %in% c("png", "jpg", "jpeg", "tif", "tiff", "pdf", "svg")) {
stop("Unsupported plot format. Use PNG, JPEG, TIFF, PDF, or SVG.", call. = FALSE)
}
dir.create(dirname(file), recursive = TRUE, showWarnings = FALSE)
if (ext == "png") {
args <- list(filename = file, width = width, height = height,
units = "in", res = dpi, bg = "white")
if (isTRUE(capabilities("cairo"))) args$type <- "cairo-png"
do.call(grDevices::png, args)
} else if (ext %in% c("jpg", "jpeg")) {
grDevices::jpeg(file, width = width, height = height, units = "in", res = dpi,
quality = 95, bg = "white")
} else if (ext %in% c("tif", "tiff")) {
grDevices::tiff(file, width = width, height = height, units = "in", res = dpi,
compression = "lzw", bg = "white")
} else if (ext == "pdf") {
grDevices::pdf(file, width = width, height = height, onefile = TRUE,
useDingbats = FALSE, bg = "white")
} else {
grDevices::svg(file, width = width, height = height, onefile = TRUE, bg = "white")
}
tryCatch(.r4vn_tabts_draw_plot_object(p), finally = grDevices::dev.off())
invisible(normalizePath(file, winslash = "/", mustWork = FALSE))
}
#' @export
print.r4vn_tabts_plot <- function(x, ...) {
.r4vn_tabts_draw_base_plot(x)
invisible(x)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.