## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4
)
set.seed(20260307)
library(ReproStat)

## ----run-diag-----------------------------------------------------------------
mt_diag <- run_diagnostics(
  mpg ~ wt + hp + disp,
  data = mtcars,
  B = 200,
  method = "bootstrap"
)
mt_diag

## ----metrics------------------------------------------------------------------
coef_stability(mt_diag)
pvalue_stability(mt_diag)
selection_stability(mt_diag)
prediction_stability(mt_diag)$mean_variance

## ----ri-----------------------------------------------------------------------
reproducibility_index(mt_diag)

## ----plots--------------------------------------------------------------------
oldpar <- par(mfrow = c(1, 2))
plot_stability(mt_diag, "coefficient")
plot_stability(mt_diag, "selection")
par(oldpar)

## ----methods, eval = FALSE----------------------------------------------------
# # Bootstrap resampling (default)
# run_diagnostics(mpg ~ wt + hp, mtcars, method = "bootstrap")
# 
# # Subsampling (80 % of rows without replacement)
# run_diagnostics(mpg ~ wt + hp, mtcars, method = "subsample", frac = 0.8)
# 
# # Gaussian noise injection (5 % of each column's SD)
# run_diagnostics(mpg ~ wt + hp, mtcars, method = "noise", noise_sd = 0.05)

## ----cv-----------------------------------------------------------------------
models <- list(
  baseline = mpg ~ wt + hp + disp,
  compact = mpg ~ wt + hp,
  expanded = mpg ~ wt + hp + disp + qsec
)

cv_result <- cv_ranking_stability(models, mtcars, v = 5, R = 40)
cv_result$summary

## ----cv-plot------------------------------------------------------------------
plot_cv_stability(cv_result, metric = "top1_frequency")

## ----iris-example-------------------------------------------------------------
iris_diag <- run_diagnostics(
  Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width,
  data = iris,
  B = 150,
  method = "noise",
  noise_sd = 0.03
)
reproducibility_index(iris_diag)

## ----airquality-example-------------------------------------------------------
aq <- na.omit(airquality[, c("Ozone", "Solar.R", "Wind", "Temp")])
aq_diag <- run_diagnostics(
  Ozone ~ Solar.R + Wind + Temp,
  data = aq,
  B = 150,
  method = "bootstrap"
)
reproducibility_index(aq_diag)

## ----ci-----------------------------------------------------------------------
set.seed(20260307)
d <- run_diagnostics(mpg ~ wt + hp + disp, mtcars, B = 150)
ri_confidence_interval(d, level = 0.95, R = 500)

## ----backends, eval = FALSE---------------------------------------------------
# # Generalized linear model (logistic)
# run_diagnostics(
#   am ~ wt + hp,
#   mtcars,
#   B = 100,
#   family = stats::binomial()
# )
# 
# # Robust regression (requires MASS)
# if (requireNamespace("MASS", quietly = TRUE)) {
#   run_diagnostics(mpg ~ wt + hp, mtcars, B = 100, backend = "rlm")
# }
# 
# # LASSO (requires glmnet)
# if (requireNamespace("glmnet", quietly = TRUE)) {
#   run_diagnostics(
#     mpg ~ wt + hp + disp + qsec,
#     mtcars,
#     B = 100,
#     backend = "glmnet",
#     en_alpha = 1
#   )
# }

## ----ggplot-helpers, eval = requireNamespace("ggplot2", quietly = TRUE)-------
if (requireNamespace("ggplot2", quietly = TRUE)) {
  set.seed(1)
  d_gg <- run_diagnostics(mpg ~ wt + hp + disp, mtcars, B = 50)
  print(plot_stability_gg(d_gg, "coefficient"))
  print(plot_stability_gg(d_gg, "selection"))
}