## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>",
  fig.width = 7, fig.height = 5,
  warning = FALSE, message = FALSE
)
library(SQIpro)

## ----install, eval=FALSE------------------------------------------------------
# # From CRAN (once published)
# install.packages("SQIpro")
# 
# # Development version from GitHub
# # install.packages("remotes")
# remotes::install_github("yourname/SQIpro")

## ----data---------------------------------------------------------------------
data(soil_data)
dim(soil_data)
head(soil_data)

## ----data-summary-------------------------------------------------------------
table(soil_data$LandUse, soil_data$Depth)

## ----validate-----------------------------------------------------------------
result <- validate_data(
  soil_data,
  group_cols = c("LandUse", "Depth")
)
result$n_per_group

## ----config-------------------------------------------------------------------
cfg <- make_config(
  variable    = c("pH",   "EC",   "BD",   "OC",   "MBC",
                  "PMN",  "Clay", "WHC",  "DEH",  "AP",  "TN"),
  type        = c("opt",  "less", "less", "more", "more",
                  "more", "opt",  "more", "more", "more","more"),
  opt_low     = c(6.0,   NA,     NA,     NA,     NA,
                  NA,     20,     NA,     NA,     NA,    NA),
  opt_high    = c(7.0,   NA,     NA,     NA,     NA,
                  NA,     35,     NA,     NA,     NA,    NA),
  description = c(
    "Soil pH (H2O 1:2.5)",
    "Electrical Conductivity (dS/m)",
    "Bulk Density (g/cm3)",
    "Organic Carbon (%)",
    "Microbial Biomass Carbon (mg/kg)",
    "Potentially Mineralizable N (mg/kg)",
    "Clay content (%)",
    "Water Holding Capacity (%)",
    "Dehydrogenase Activity (ug TPF/g/day)",
    "Available Phosphorus (mg/kg)",
    "Total Nitrogen (%)"
  )
)
print(cfg)

## ----scoring-curves, fig.height=7---------------------------------------------
plot_scoring_curves(soil_data, cfg,
                    group_cols = c("LandUse", "Depth"))

## ----score--------------------------------------------------------------------
scored <- score_all(soil_data, cfg,
                    group_cols = c("LandUse", "Depth"))
head(scored[, c("LandUse", "Depth", "pH", "OC", "MBC")])

## ----mds----------------------------------------------------------------------
mds <- select_mds(scored,
                  group_cols     = c("LandUse", "Depth"),
                  load_threshold = 0.6)
mds$mds_vars

## ----linear-------------------------------------------------------------------
res_lin <- sqi_linear(scored, cfg,
                      group_cols = c("LandUse", "Depth"),
                      mds_vars   = mds$mds_vars)
print(res_lin)

## ----regression---------------------------------------------------------------
res_reg <- sqi_regression(scored, cfg,
                           dep_var    = "OC",
                           group_cols = c("LandUse", "Depth"),
                           mds_vars   = mds$mds_vars)
print(res_reg)

## ----pca----------------------------------------------------------------------
res_pca <- sqi_pca(scored, cfg,
                   group_cols = c("LandUse", "Depth"),
                   mds        = mds)
print(res_pca)

## ----fuzzy--------------------------------------------------------------------
res_fuz <- sqi_fuzzy(scored, cfg,
                     group_cols = c("LandUse", "Depth"),
                     mds_vars   = mds$mds_vars,
                     operator   = "mean")
print(res_fuz)

## ----entropy------------------------------------------------------------------
res_ent <- sqi_entropy(scored, cfg,
                       group_cols = c("LandUse", "Depth"),
                       mds_vars   = mds$mds_vars)
print(res_ent)
cat("\nEntropy weights:\n")
print(attr(res_ent, "entropy_weights"))

## ----topsis-------------------------------------------------------------------
res_top <- sqi_topsis(scored, cfg,
                      group_cols = c("LandUse", "Depth"),
                      mds_vars   = mds$mds_vars)
print(res_top)

## ----compare------------------------------------------------------------------
comparison <- sqi_compare(scored, cfg,
                           group_cols = c("LandUse", "Depth"),
                           dep_var    = "OC",
                           mds        = mds)
print(comparison)

## ----heatmap------------------------------------------------------------------
plot_scores(scored, cfg,
            group_cols = c("LandUse", "Depth"),
            group_by   = "LandUse",
            facet_by   = "Depth")

## ----bar----------------------------------------------------------------------
plot_sqi(res_lin,
         sqi_col   = "SQI_linear",
         group_col = "LandUse",
         fill_col  = "Depth")

## ----radar, fig.height=6, eval=requireNamespace("fmsb", quietly=TRUE)---------
# # Requires the 'fmsb' package: install.packages("fmsb")
# plot_radar(scored, cfg,
#            group_col  = "LandUse",
#            group_cols = c("LandUse", "Depth"))

## ----biplot-------------------------------------------------------------------
plot_pca_biplot(mds, scored, group_col = "LandUse")

## ----anova--------------------------------------------------------------------
# Compute per-observation index
scored$SQI_obs <- rowMeans(scored[, mds$mds_vars], na.rm = TRUE)

aov_res <- sqi_anova(scored,
                     sqi_col   = "SQI_obs",
                     group_col = "LandUse")
cat("ANOVA significant:", aov_res$significant, "\n")
print(aov_res$compact_letters)

## ----sensitivity--------------------------------------------------------------
sa <- sqi_sensitivity(scored, cfg,
                       group_cols = c("LandUse", "Depth"),
                       method     = "linear",
                       mds_vars   = mds$mds_vars)
print(sa)

## ----tornado------------------------------------------------------------------
plot_sensitivity(sa)

## ----bib, echo=FALSE----------------------------------------------------------