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

## ----setup, results = "hide", message = FALSE---------------------------------
library(mrIML)
library(tidymodels)
library(flashlight)
library(ggplot2)
library(patchwork)

set.seed(7007)

## ----load-data----------------------------------------------------------------
data <- MRFcov::Bird.parasites
head(data)

## ----prep-data----------------------------------------------------------------
# Responses
Y <- select(data, "Hzosteropis", "Hkillangoi", "Plas", "Microfilaria")
# Covariates
X <- select(data, "scale.prop.zos")

## ----setup-cluster, eval=FALSE------------------------------------------------
# future::plan(multisession, workers = 2)

## ----define-tidymodels--------------------------------------------------------
model_rf <- rand_forest(
  trees = 100, # 100 trees are set for brevity. Aim to start with 1000.
  mode = "classification",
  mtry = tune(),
  min_n = tune()
) %>%
  set_engine("randomForest")

model_lm <- logistic_reg() %>%
  set_engine("glm")

## ----fit-mrIML, results = "hide"----------------------------------------------
mrIML_rf <- mrIMLpredicts(
  X = X,
  Y = Y,
  X1 = Y,
  Model = model_rf,
  prop = 0.7,
  k = 5,
  racing = TRUE
)

mrIML_lm <- mrIMLpredicts(
  X = X,
  Y = Y,
  X1 = Y,
  Model = model_lm , 
  balance_data = 'no',
  prop = 0.6,
  k = 5,
  racing = FALSE
)

## ----performance-mrIML--------------------------------------------------------
perf_rf <- mrIML_rf %>%
  mrIMLperformance()
perf_rf$model_performance

perf_lm <- mrIML_lm %>%
  mrIMLperformance()
perf_rf$model_performance

## ----compare-mrIML------------------------------------------------------------
perf_comp <- mrPerformancePlot(perf_rf, perf_lm)
perf_comp$performance_plot + perf_comp$performance_diff_plot

## -----------------------------------------------------------------------------
fl_rf <- mrIML_rf %>%
  mrFlashlight()

fl_rf$Microfilaria %>%
  light_profile(data = data, v = "scale.prop.zos") %>%
  plot() +
  ggtitle("Effect of scale.prop.zos on Microfilaria") +
  theme_bw()

## ----fig.height = 6-----------------------------------------------------------
PD_scale.prop.zos_rf <- mrIML_rf %>%
  mrCovar(var = "scale.prop.zos", sdthresh = 0)

PD_scale.prop.zos_rf[[1]] /
  PD_scale.prop.zos_rf[[2]] /
  PD_scale.prop.zos_rf[[3]] +
  plot_layout(axis = "collect")

## ----fig.height = 6-----------------------------------------------------------
PD_scale.prop.zos_lm <- mrIML_lm %>%
  mrCovar(var = "scale.prop.zos", sdthresh = 0)

PD_scale.prop.zos_lm[[1]] /
  PD_scale.prop.zos_lm[[2]] /
  PD_scale.prop.zos_lm[[3]] +
  plot_layout(axis = "collect")

## -----------------------------------------------------------------------------
vip_rf <- mrIML_rf %>%
  mrVip()

vip_rf[[3]]

## -----------------------------------------------------------------------------
mrIML_boot_rf <- mrIML_rf %>%
  mrBootstrap()

mrPdPlotBootstrap(
  mrIML_rf,
  mrBootstrap_obj = mrIML_boot_rf,
  target = "Plas",
  global_top_var = 4
)[[2]]

## -----------------------------------------------------------------------------
vip_boot_rf <- mrVip(
  mrIMLobj = mrIML_rf,
  mrBootstrap_obj = mrIML_boot_rf
)

vip_boot_rf[[3]]

## ----fig.height = 10----------------------------------------------------------
interactions_rf <- mrIML_rf %>%
  mrInteractions(feature = "Plas", num_bootstrap = 10)

interactions_rf[[1]] /
  interactions_rf[[2]] /
  interactions_rf[[3]]