--- title: "Getting Started with psvr" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Getting Started with psvr} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4 ) ``` ## Introduction Classical SVR minimises absolute-error losses (MAE, MSE), which are misaligned with the scale-free accuracy criteria standard in forecasting. An error of 1 unit is negligible when the target is 1 000 but large when it is 2. **psvr** implements four SVR variants derived from percentage-error loss functions (Benavides-Herrera et al., 2026), accessed through one fitter per model family: | Model | Function | `sym_type` | Solver | |-------|----------|------------|------------------| | 1 | `psvr_mape()` | `"none"` | quadratic program | | 2 | `psvr_mape()` | `"even"` / `"odd"` | quadratic program | | 3 | `psvr_rmspe()` | `"none"` | linear system | | 4 | `psvr_rmspe()` | `"even"` / `"odd"` | linear system | `sym_type = "even"` enforces even symmetry `f(x) = f(-x)`; `"odd"` enforces odd symmetry. Use the symmetric variants only with kernels that satisfy Assumption 3 of the paper (RBF and even-degree polynomial kernels do). These are the two direct fitters. They replaced the unified `psvr()` entry point in 0.0.2.9012, which had in turn replaced four separate wrappers in 0.0.2.9010: seven of `psvr()`'s eleven arguments were conditional on which family you were fitting, so the families are now separate functions. Both return the same classes the tidymodels engine returns, so a fit obtained either way behaves identically. All models require **strictly positive targets** (`y > 0`), which is the condition under which percentage residuals are well-defined. ## Fuel-economy data This vignette demonstrates the package on a small public dataset. We use `ggplot2::mpg`: 234 records of US passenger vehicles. The target is `hwy`, highway fuel economy in miles per gallon. It is strictly positive on every observation, so the percentage residuals — and therefore MAPE and RMSPE — are well-defined throughout. ```{r data} library(psvr) library(ggplot2) # Target: highway fuel economy (all > 0) y_all <- mpg$hwy # Predictors: engine displacement, model year, cylinder count. # # `cty` is deliberately EXCLUDED. City and highway fuel economy are two # measurements of the same property of the same vehicle (they correlate at # 0.96), so predicting one from the other is leakage rather than modelling. X_raw <- as.matrix(mpg[, c("displ", "year", "cyl")]) stopifnot(all(y_all > 0)) cat("N =", nrow(X_raw), " p =", ncol(X_raw), " y range: [", min(y_all), ",", max(y_all), "]\n") ``` ### 70 / 30 train–test split Features are standardised using training-set statistics so that the RBF kernel operates on a comparable scale across all three predictors. ```{r split} set.seed(42) n <- nrow(X_raw) tr_idx <- sample(n, floor(0.7 * n)) X_raw_tr <- X_raw[tr_idx, ]; y_tr <- y_all[tr_idx] X_raw_te <- X_raw[-tr_idx, ]; y_te <- y_all[-tr_idx] # Standardise: centre and scale by training mean/sd col_mean <- colMeans(X_raw_tr) col_sd <- apply(X_raw_tr, 2, sd) X_tr <- scale(X_raw_tr, center = col_mean, scale = col_sd) X_te <- scale(X_raw_te, center = col_mean, scale = col_sd) ``` ### Helper metrics ```{r metrics} mape <- function(y, yhat) mean(abs(y - yhat) / y) * 100 rmspe <- function(y, yhat) sqrt(mean(((y - yhat) / y)^2)) * 100 r2 <- function(y, yhat) 1 - sum((y - yhat)^2) / sum((y - mean(y))^2) ``` ### Baseline: linear regression ```{r baseline} lm_df_tr <- as.data.frame(X_tr) lm_df_te <- as.data.frame(X_te) lm_fit <- lm(y_tr ~ ., data = lm_df_tr) lm_pred <- predict(lm_fit, newdata = lm_df_te) cat(sprintf("Linear regression — MAPE: %.2f%% RMSPE: %.2f%% R²: %.4f\n", mape(y_te, lm_pred), rmspe(y_te, lm_pred), r2(y_te, lm_pred))) ``` ## Model 1: ε-SVR with MAPE The ε-SVR formulation optimises a QP with **sample-dependent box constraints** $|\beta_k| \le 100C/y_k$: tighter bounds for small targets, concentrating model capacity on low-magnitude observations. ```{r mape-svr} # make_kernel() returns a closure K(xi, xj) = exp(-||xi - xj||^2 / (2 sigma^2)). # sigma is a LENGTH in the units of the preprocessed feature space, so it has # to be set on that scale -- sigma_heuristic() reads it off the data instead of # guessing. See "Hyperparameter search ranges" below. K <- make_kernel("rbf", sigma = sigma_heuristic(X_tr)) # C = 10: per-sample box bound |beta_k| <= 100*C/y_k; eps = 1: tube width (% of y_k) fit_ep <- psvr_mape(X_tr, y_tr, kernel = K, C = 10, eps = 1) pred_ep <- predict(fit_ep, X_te) cat(sprintf("ε-SVR MAPE — MAPE: %.2f%% RMSPE: %.2f%% R²: %.4f\n", mape(y_te, pred_ep), rmspe(y_te, pred_ep), r2(y_te, pred_ep))) cat(sprintf("Support vectors: %d / %d\n", length(fit_ep$beta), fit_ep$n_train)) print(fit_ep) ``` ```{r mape-plot, echo = FALSE} lim <- range(c(y_te, pred_ep)) data.frame(actual = y_te, predicted = pred_ep) |> ggplot(aes(actual, predicted)) + geom_point(colour = "#d7191c", alpha = 0.65, size = 1.8) + geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey40") + coord_equal(xlim = lim, ylim = lim) + labs(x = "Actual hwy (mpg)", y = "Predicted hwy (mpg)", title = "Model 1: \u03b5-SVR with MAPE") ``` ```{r mape-coef} cf_ep <- coef(fit_ep) # alpha, alpha_star: length-N dual variables (paired); the pre-pruning # solution. Useful as a warm-start across CV folds; for # prediction use `beta` instead. # beta: beta_k = alpha_k - alpha_k* for each SUPPORT VECTOR only # (non-zero only for training points outside the # percentage-error ε-tube — sparse) # b: bias / intercept term # support_data: training rows corresponding to support vectors only cat(sprintf("b = %.4f | beta range: [%.4f, %.4f]\n", cf_ep$b, min(cf_ep$beta), max(cf_ep$beta))) ``` ## Model 3: LS-SVR with RMSPE The LS-SVR formulation replaces the QP with a linear system by using a quadratic penalty on percentage residuals. The dual reduces to: $$ \begin{bmatrix} 0 & \mathbf{1}^\top \\ \mathbf{1} & \Omega + Y_\Gamma \end{bmatrix} \begin{bmatrix} b \\ \boldsymbol{\alpha} \end{bmatrix} = \begin{bmatrix} 0 \\ \mathbf{y} \end{bmatrix} $$ where $Y_\Gamma = \operatorname{diag}(y_1^2/\Gamma, \ldots, y_N^2/\Gamma)$. ```{r rmspe-lssvr} # gamma = 5000: regularisation; larger gamma -> smaller Y_Gamma diagonal -> tighter fit. # This is roughly var(y_tr) * N, the scale cost_psvr_ls_data() computes -- and # already five times the ceiling of the registered `cost` default. See below. fit_ls <- psvr_rmspe(X_tr, y_tr, kernel = K, gamma = 5000) pred_ls <- predict(fit_ls, X_te) cat(sprintf("LS-SVR RMSPE — MAPE: %.2f%% RMSPE: %.2f%% R²: %.4f\n", mape(y_te, pred_ls), rmspe(y_te, pred_ls), r2(y_te, pred_ls))) print(fit_ls) ``` ```{r rmspe-plot, echo = FALSE} lim <- range(c(y_te, pred_ls)) data.frame(actual = y_te, predicted = pred_ls) |> ggplot(aes(actual, predicted)) + geom_point(colour = "#2c7bb6", alpha = 0.65, size = 1.8) + geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey40") + coord_equal(xlim = lim, ylim = lim) + labs(x = "Actual hwy (mpg)", y = "Predicted hwy (mpg)", title = "Model 3: LS-SVR with RMSPE") ``` ```{r rmspe-coef} cf_ls <- coef(fit_ls) # alpha: N dual variables; weight each training point's kernel # contribution in f(x) = sum_k alpha_k K(x_k, x) + b # (all N points, no sparsity) # b: bias / intercept term # support_data: all N training inputs stored for prediction cat(sprintf("b = %.4f | alpha range: [%.4f, %.4f]\n", cf_ls$b, min(cf_ls$alpha), max(cf_ls$alpha))) ``` ## Comparing objectives ```{r comparison} results <- data.frame( Model = c("Linear regression", "\u03b5-SVR MAPE (Model 1)", "LS-SVR RMSPE (Model 3)"), MAPE = c(mape(y_te, lm_pred), mape(y_te, pred_ep), mape(y_te, pred_ls)), RMSPE = c(rmspe(y_te, lm_pred), rmspe(y_te, pred_ep), rmspe(y_te, pred_ls)), R2 = c(r2(y_te, lm_pred), r2(y_te, pred_ep), r2(y_te, pred_ls)) ) results[, 2:4] <- round(results[, 2:4], 2) knitr::kable(results, col.names = c("Model", "MAPE (%)", "RMSPE (%)", "R²"), align = "lrrr", caption = paste("Test-set performance on ggplot2::mpg", "(70/30 split, RBF kernel, single run,", "untuned hyperparameters).")) ``` Both psvr models improve on the linear baseline under their respective percentage-error objectives. The ε-SVR formulation (Model 1) targets MAPE through sample-dependent box constraints on the dual variables; the LS-SVR formulation (Model 3) minimises RMSPE directly. Read the table for the ordering, not for the magnitudes. It is one split of one small dataset at hyperparameters nobody tuned, so the margins are not evidence about how these models perform in general — that is what the resampled comparison in the [When to Use Percentage-Error SVR](https://pbenavidesh.github.io/psvr/articles/when-to-use-psvr.html) article is for. ## Using psvr with tidymodels All four models are registered as parsnip engines and integrate seamlessly with the tidymodels ecosystem. This enables hyperparameter tuning via `tune_grid()`, resampling via `rsample`, and unified model comparison via `workflow_set()`. See the [tidymodels workflow](tidymodels-workflow.html) article for a complete example using `psvr_rmspe_rbf()` with `tune_grid()` and data-driven hyperparameter ranges via `rbf_sigma_psvr_data()`. The [When to Use Percentage-Error SVR](https://pbenavidesh.github.io/psvr/articles/when-to-use-psvr.html) article, on the package website rather than in the installed package, shows a full `workflow_set()` comparison of all four models against standard baselines. ## Hyperparameter search ranges The parsnip specs register **static** default ranges — a `dials` parameter object has to exist before any data does. Three of them are deliberately conservative, one of them is actively wrong for LS-SVR, and one searches a level you probably do not want, so tuning a psvr model means supplying your own ranges through `param_info`. Nothing below fits a model. It builds the parameter set and prints it, before and after. ```{r ranges-defaults, message = FALSE} library(parsnip) library(tune) show_ranges <- function(ps) { for (i in seq_len(nrow(ps))) { ob <- ps$object[[i]] cat(sprintf(" %-10s %s\n", ps$id[i], if (inherits(ob, "quant_param")) sprintf("[%s] on the %s scale", paste(signif(unlist(ob$range), 4), collapse = ", "), if (is.null(ob$trans)) "identity" else ob$trans$name) else sprintf("{%s}", paste(ob$values, collapse = ", ")))) } invisible(ps) } spec_mape <- psvr_mape_rbf(cost = tune(), margin = tune(), rbf_sigma = tune(), sym_type = tune()) |> set_engine("psvr") extract_parameter_set_dials(spec_mape) |> show_ranges() ``` For an ε-SVR (Models 1–2) only `rbf_sigma` really needs replacing; the `cost` and `margin` defaults are usable as they stand. ```{r ranges-mape} extract_parameter_set_dials(spec_mape) |> update( cost = cost_psvr(), # [-2, 10] log2 — fine for C margin = margin_percentage(), # 1-20% of each target rbf_sigma = rbf_sigma_psvr_data(X_tr), # data-driven; see below sym_type = sym_type_param(c("even", "odd")) # drops "none"; see below ) |> show_ranges() ``` For an LS-SVR (Models 3–4) `cost` is $\Gamma$, and there the registered default is **not** usable. ```{r ranges-rmspe} spec_ls <- psvr_rmspe_rbf(cost = tune(), rbf_sigma = tune(), sym_type = tune()) |> set_engine("psvr") extract_parameter_set_dials(spec_ls) |> update( cost = cost_psvr_ls_data(y_tr), rbf_sigma = rbf_sigma_psvr_data(X_tr), sym_type = sym_type_param(c("even", "odd")) ) |> show_ranges() ``` No single spec carries all four helpers: `margin` exists only on the MAPE specs, and `cost_psvr_ls_data()` applies only to the RMSPE ones. ### `rbf_sigma` is a length scale, not a precision Of everything on this page, this is the one most likely to cost you a day. psvr's RBF kernel is $$K(\mathbf{x}_i, \mathbf{x}_j) = \exp\left(-\frac{\|\mathbf{x}_i - \mathbf{x}_j\|^2}{2\sigma^2}\right)$$ so `rbf_sigma` is a **length**, measured in the units of the preprocessed feature space, and *larger* means a *wider* kernel. `kernlab::rbfdot(sigma =)`, `e1071::svm(gamma =)` and `parsnip::svm_rbf(rbf_sigma =)` all use the reciprocal convention, $\exp(-\sigma\|\cdot\|^2)$, where larger means *narrower*. Carrying a tuned value across from one of those gives a silently wrong kernel width — no error, no warning, just a worse model. ### `rbf_sigma_psvr_data()` needs preprocessed predictors `rbf_sigma_psvr_data(X)` centres its range on the median pairwise Euclidean distance between the rows of `X`, spanning one decade either side of it on the log10 scale. Being a distance-based heuristic, it means nothing except on the scale the model is actually fitted on: pass the standardised predictors (`X_tr` here, or the baked output of a recipe), never the raw ones. ```{r ranges-sigma} sigma_heuristic(X_tr) # the geometric centre of the range printed above, # and the value the fits at the top of this page used ``` Above `sample_size` rows (default 500) the median is taken on a random subsample, so the centre becomes an estimate rather than the exact median — and a seed-dependent one. Pass `seed` if you need it reproducible. ### LS-SVR needs `cost_psvr_ls_data()`, and it cannot be automated `cost` maps to `C` on the ε-SVR models and to $\Gamma$ on the LS-SVR ones, but both register the same default of $[-2, 10]$ on the log2 scale, i.e. $\Gamma \le 1024$. That is the ε-SVR range. $\Gamma$ enters the LS-SVR system only through the $y_k^2/\Gamma$ diagonal, so the value that balances that term against the kernel scales with `var(y) * n` — it is not a fixed magnitude, and it grows with both the spread of the outcome and the size of the training set. The fit at the top of this page is already past the default ceiling: it used `gamma = 5000`, and $[-2, 10]$ stops at 1024. A grid over the default would be **boundary-trapped** — every candidate is legal, the search reports a plausible number, and the optimum was never inside the range. The [tidymodels workflow](tidymodels-workflow.html) article shows the size of the effect. `cost_psvr_ls_data(y)` sets the ceiling from `var(y) * n`. It has to be passed by hand: `tune` finalises parameters from the molded **predictors** alone and never passes the outcome to `dials::finalize()`, so no machinery could compute a `var(y)`-based bound on your behalf. ### `cost` is not comparable across datasets either The same trap has a milder form on the ε-SVR side, and `cost_psvr()` does not protect you from it. The dual box constraint is $$|\beta_k| \le \frac{100\,C}{y_k}$$ so the bound a given `C` imposes depends on the **magnitude of the outcome**. On data where `y` is of order $10^4$ rather than order $10$, the same `C` yields a box three decades tighter, every multiplier saturates against it, and the fit degenerates towards a constant. `cost_psvr()` returns a static $[-2, 10]$ calibrated for outcomes of order 10, so on large-`y` data `C` has to be raised by hand to compensate — there is currently no `cost_psvr_data()` counterpart to `cost_psvr_ls_data()` for the ε-SVR family. ### Tuning `sym_type` cannot answer "does symmetry help" `sym_type` registers all three levels, `{none, even, odd}`, and `"none"` *is* the asymmetric model. So `sym_type = tune()` puts Models 1 and 2 — or 3 and 4 — inside a single search space, and `select_best()` returns whichever candidate won. That is a selection, not a comparison. To contrast the two families, fix `sym_type` on two specs and compare their resampled metrics. Restricting to `c("even", "odd")`, as above, tunes *within* the symmetric family, which is a third question again. ### The outcome stays in original units `recipes::step_normalize()` applies to predictors. Leave the outcome alone: the percentage-error losses divide by `y`, so they need `y > 0` on the original scale, and centring the outcome would destroy positivity and change what "percentage error" even refers to. This is a genuine departure from classical SVR, where rescaling the outcome is routine and harmless. ## References Benavides-Herrera, P., Álvarez, G., Ruiz-Cruz, R., & Sánchez-Torres, J. D. (2026). A unified family of percentage-error support vector regression models with symmetric kernel extensions. *Mathematics*, 14(10), 1679.