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

## ----setup--------------------------------------------------------------------
library(CCI)

## -----------------------------------------------------------------------------
normal_data <- function(n) {
  Z1 <- rnorm(n)
  Z2 <- rnorm(n)
  X <- Z1 + Z2 + rnorm(n)
  Y <- Z1 + Z2 + rnorm(n)
  data.frame(Z1, Z2, X, Y)
}
set.seed(1)
dat <- normal_data(500)

## ----eval = FALSE-------------------------------------------------------------
# my_metric <- function(actual, predictions) {
#   # compute and return one number
# }

## -----------------------------------------------------------------------------
r_squared <- function(actual, predictions) {
  1 - sum((actual - predictions)^2) / sum((actual - mean(actual))^2)
}
res_r2 <- CCI.test(Y ~ X | Z1, data = dat, metricfunc = r_squared, tail = "right",
                   seed = 1, progress = FALSE)
summary(res_r2)

## -----------------------------------------------------------------------------
mae <- function(actual, predictions) mean(abs(actual - predictions))
summary(CCI.test(Y ~ X | Z1 + Z2, data = dat, method = "KNN", metricfunc = mae, tail = "left",
                 seed = 1, progress = FALSE))

## -----------------------------------------------------------------------------
set.seed(2)
cat_data <- normal_data(500)
cat_data$Y <- factor(ifelse(cat_data$Y > 1, "high", "low"))   # unequal class sizes
table(cat_data$Y)

balanced_accuracy <- function(actual, predictions) {
  mean(tapply(as.character(predictions) == as.character(actual), actual, mean))
}
summary(CCI.test(Y ~ X | Z1, data = cat_data, metricfunc = balanced_accuracy, tail = "right",
                 seed = 1, progress = FALSE))

## -----------------------------------------------------------------------------
brier <- function(actual, predictions) {
  mean((as.numeric(actual == levels(actual)[2]) - predictions)^2)
}
summary(CCI.test(Y ~ X | Z1, data = cat_data, method = "xgboost", nrounds = 100, eta = 0.1,
                 metricfunc = brier, tail = "left", seed = 1, progress = FALSE))

## ----eval = FALSE-------------------------------------------------------------
# my_wrapper <- function(formula, data, train_indices, test_indices, ...) {
#   model <- train_model(formula, data = data[train_indices, ], ...)
#   predictions <- predict(model, data[test_indices, ])
#   actual <- data[test_indices, all.vars(formula)[1]]
#   compute_metric(actual, predictions)
# }

## -----------------------------------------------------------------------------
lm_wrapper <- function(formula, data, train_indices, test_indices, ...) {
  model <- lm(formula, data = data[train_indices, ])
  predictions <- predict(model, newdata = data[test_indices, ])
  actual <- data[test_indices, all.vars(formula)[1]]
  sqrt(mean((actual - predictions)^2))   # RMSE: lower is better
}
summary(CCI.test(Y ~ X | Z1 + Z2, data = dat, mlfunc = lm_wrapper, tail = "left",
                 seed = 1, progress = FALSE))
summary(CCI.test(Y ~ X | Z1, data = dat, mlfunc = lm_wrapper, tail = "left",
                 seed = 1, progress = FALSE))

## -----------------------------------------------------------------------------
logistic_wrapper <- function(formula, data, train_indices, test_indices, clip = 1e-6, ...) {
  model <- glm(formula, data = data[train_indices, ], family = binomial)
  prob <- predict(model, newdata = data[test_indices, ], type = "response")
  prob <- pmin(pmax(prob, clip), 1 - clip)
  actual <- data[test_indices, all.vars(formula)[1]]
  is_second <- actual == levels(actual)[2]    # glm models the probability of the second level
  -mean(ifelse(is_second, log(prob), log(1 - prob)))   # log loss: lower is better
}
summary(CCI.test(Y ~ X | Z1, data = cat_data, mlfunc = logistic_wrapper, tail = "left",
                 clip = 1e-4, poly = FALSE, interaction = FALSE, seed = 1, progress = FALSE))

## -----------------------------------------------------------------------------
caret_wrapper <- function(formula, data, train_indices, test_indices, caret_method, ...) {
  model <- caret::train(formula, data = data[train_indices, ], method = caret_method,
                        trControl = caret::trainControl(method = "none"), ...)
  predictions <- predict(model, newdata = data[test_indices, ])
  actual <- data[test_indices, all.vars(formula)[1]]
  sqrt(mean((actual - predictions)^2))
}
summary(CCI.test(Y ~ X | Z1, data = dat, mlfunc = caret_wrapper, tail = "left",
                 caret_method = "knn", tuneGrid = data.frame(k = 15),
                 seed = 1, progress = FALSE))

## -----------------------------------------------------------------------------
QQplot(res_r2, nperm = 40, progress = FALSE)