R/tabts.R

Defines functions print.r4vn_tabts_plot .r4vn_tabts_save_plot_object plot.r4vn_tabts summary.r4vn_tabts print.r4vn_tabts_grouped print.r4vn_tabts .r4vn_ts_show_grouped .r4vn_ts_show .r4vn_ts_primary_plot .r4vn_tabts_open_viewer .r4vn_tabts_build_html .r4vn_tabts_wrap_html .r4vn_tabts_result_body .r4vn_tabts_plot_title .r4vn_tabts_plot_png_html .r4vn_tabts_draw_plot_object .r4vn_tabts_base64 .r4vn_ts_df_html .r4vn_tabts_transpose_table .r4vn_ts_html_escape .r4vn_ts_label_df .r4vn_ts_ai .r4vn_ts_interpret .r4vn_ts_label .r4vn_ts_apply_theme .r4vn_ts_apply_plot_options .r4vn_ts_add_interventions .r4vn_ts_plot_counterfactual .r4vn_ts_plot_residual .r4vn_ts_plot_acf .r4vn_ts_plot_decomposition .r4vn_ts_plot_series .r4vn_ts_make_plots .r4vn_tabts_draw_base_plot .r4vn_tabts_add_intervention_base .r4vn_tabts_base_axis .r4vn_tabts_base_x .r4vn_tabts_make_base_plots .r4vn_tabts_plot_spec .r4vn_ts_plot_set .r4vn_ts_predict_interval .r4vn_ts_make_its_future .r4vn_ts_predict_mean .r4vn_ts_controlled_effect_at .r4vn_ts_its_term_labels .r4vn_ts_match_intervention_times .r4vn_ts_fit_its .r4vn_ts_forecast_table .r4vn_ts_future_xreg .r4vn_ts_diagnostics .r4vn_ts_residuals .r4vn_ts_fitted .r4vn_ts_coefficients .r4vn_ts_choose_model .r4vn_ts_model_df .r4vn_tabts_model_info .r4vn_ts_aicc .r4vn_ts_bic .r4vn_ts_aic .r4vn_ts_accuracy .r4vn_ts_forecast_engine .r4vn_ts_test_split .r4vn_ts_fit_engine .r4vn_tabts_auto_arima_base .r4vn_ts_xreg_matrix .r4vn_ts_acf .r4vn_ts_decompose .r4vn_ts_stationarity .r4vn_ts_summary .r4vn_ts_extend_time .r4vn_ts_period .r4vn_ts_full_time .r4vn_ts_infer_step .r4vn_ts_interp .r4vn_ts_collapse_duplicates .r4vn_ts_regularize .r4vn_ts_core .r4vn_ts_grouped .r4vn_ts_resolve_intervention .r4vn_ts_many_names .r4vn_ts_optional_name .r4vn_ts_one_name .r4vn_ts_resolve_data .r4vn_ts_merge .r4vn_ts_defaults `%||%` tabts

Documented in plot.r4vn_tabts tabts

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("&", "&amp;", x, fixed = TRUE)
  x <- gsub("<", "&lt;", x, fixed = TRUE)
  x <- gsub(">", "&gt;", x, fixed = TRUE)
  x <- gsub('"', "&quot;", 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> &nbsp;&bull;&nbsp; Seasonal period: ", x$period,
    " &nbsp;&bull;&nbsp; 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)
}

Try the R4VN package in your browser

Any scripts or data that you put into this service are public.

R4VN documentation built on Sept. 30, 2026, 5:13 p.m.