## ----setup, include = FALSE---------------------------------------------------
knitr::knit_hooks$set(pngquant = knitr::hook_pngquant)
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE
)
shapr_eval <- requireNamespace("shapr", quietly = TRUE)

## -----------------------------------------------------------------------------
library(nestedcv)
library(mlbench)  # Boston housing dataset

data(BostonHousing2)
dat <- BostonHousing2
y <- dat$cmedv
x <- subset(dat, select = -c(cmedv, medv, town, chas))

# Fit a glmnet model using nested CV
set.seed(1, "L'Ecuyer-CMRG")
fit <- nestcv.glmnet(y, x, family = "gaussian",
                     min_1se = 1, alphaSet = 1, cv.cores = 2)
vs <- var_stability(fit)
vs

## ----fig.dim = c(10, 5)-------------------------------------------------------
p1 <- plot_var_stability(fit)

# overlay directionality using colour
p2 <- plot_var_stability(fit, final = FALSE, direction = 1)

# or show directionality with the sign of the variable importance
# plot_var_stability(fit, final = FALSE, percent = F)

ggpubr::ggarrange(p1, p2, ncol=2)

## ----eval=FALSE---------------------------------------------------------------
# # change bubble colour scheme
# p1 + scale_fill_manual(values=c("orange", "green3"))

## ----fig.dim = c(10, 5)-------------------------------------------------------
# beeswarm plot of variable ranks 
p1 <- plot_var_ranks(fit)

# histogram of variable ranks
p2 <- hist_var_ranks(fit)

ggpubr::ggarrange(p1, p2, ncol=2)

## ----eval = shapr_eval, message=FALSE-----------------------------------------
library(shapr)

# Generate SHAP values using shapr via nestcv.explain()
# For regression, alternative `approach` could be "empirical" or "gaussian".
sh <- nestcv.explain(fit, pred_nestcv_glmnet)

# Plot overall variable importance
plot_shap_bar(sh, x)

## ----eval = shapr_eval--------------------------------------------------------
# Plot beeswarm plot
plot_shap_beeswarm(sh, x, size = 1)

## ----eval=FALSE---------------------------------------------------------------
# # Only 3 outer folds to speed up process
# fit_gbm <- nestcv.train(y, x,
#                         method = "gbm",
#                         n_outer_folds = 3, cv.cores = 2)
# 
# # phi0 and approach are both auto-computed/auto-selected
# sh <- nestcv.explain(fit_gbm, pred_train)
# plot_shap_beeswarm(sh, x, size = 1)

## ----eval = shapr_eval, message=FALSE, fig.width = 9, fig.height = 3.5--------
library(ggplot2)
data("iris")
dat <- iris
y <- dat$Species
x <- dat[, 1:4]

# Only 3 outer folds to speed up process
fit_iris <- nestcv.glmnet(y, x, family = "multinomial", n_outer_folds = 3,
                          alphaSet = 0.6)

# SHAP values for each of the 3 classes
# phi0 is auto-computed per class as mean predicted probability
sh1 <- nestcv.explain(fit_iris,
                      predict_model = pred_nestcv_glmnet_class(1))
sh2 <- nestcv.explain(fit_iris,
                      predict_model = pred_nestcv_glmnet_class(2))
sh3 <- nestcv.explain(fit_iris,
                      predict_model = pred_nestcv_glmnet_class(3))

s1 <- plot_shap_bar(sh1, x, sort = FALSE) +
  ggtitle("Setosa")
s2 <- plot_shap_bar(sh2, x, sort = FALSE) +
  ggtitle("Versicolor")
s3 <- plot_shap_bar(sh3, x, sort = FALSE) +
  ggtitle("Virginica")

ggpubr::ggarrange(s1, s2, s3, ncol=3, legend = "bottom", common.legend = TRUE)

## ----eval = shapr_eval, fig.width = 9.5, fig.height = 3.5---------------------
s1 <- plot_shap_beeswarm(sh1, x, sort = FALSE, cex = 0.7) +
  ggtitle("Setosa")
s2 <- plot_shap_beeswarm(sh2, x, sort = FALSE, cex = 0.7) +
  ggtitle("Versicolor")
s3 <- plot_shap_beeswarm(sh3, x, sort = FALSE, cex = 0.7) +
  ggtitle("Virginica")

ggpubr::ggarrange(s1, s2, s3, ncol=3, legend = "right", common.legend = TRUE)

## ----eval=FALSE---------------------------------------------------------------
# sh <- nestcv.explain(fit, x_explain = fit$xsub,
#                      predict_model = pred_nestcv_glmnet)
# plot_shap_bar(sh, fit$xsub)
# 
# sh <- nestcv.explain(fit, x_explain = x[, fit$final_vars],
#                      predict_model = pred_nestcv_glmnet)
# plot_shap_bar(sh, x[, fit$final_vars])

