## ----include = FALSE----------------------------------------------------------
ok <- requireNamespace("rpart", quietly = TRUE) && requireNamespace("MASS", quietly = TRUE)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4, dpi = 72, out.width = "100%", eval = ok)

## -----------------------------------------------------------------------------
library(SynthesizerPlus)
library(rpart)
library(ggplot2)
data(Boston, package = "MASS")

rmse <- function(pred, obs) sqrt(mean((pred - obs)^2))
col_real <- "#2a78d6"
col_aug <- "#eb6834"

## -----------------------------------------------------------------------------
# four learners with very different capacity
deep_tree <- function(d, te) {
  predict(rpart(medv ~ ., d, control = rpart.control(cp = 0, minsplit = 2,
                                                     minbucket = 1, xval = 0)), te)
}
pruned_tree <- function(d, te) {     # complexity chosen by 10-fold cross-validation
  f <- rpart(medv ~ ., d, control = rpart.control(cp = 0, minsplit = 2,
                                                  minbucket = 1, xval = 10))
  cp <- f$cptable[which.min(f$cptable[, "xerror"]), "CP"]
  predict(prune(f, cp = cp), te)
}
knn <- function(d, te, k) {
  X <- scale(as.matrix(d[, names(d) != "medv"]))
  Xt <- scale(as.matrix(te[, names(d) != "medv"]),
              attr(X, "scaled:center"), attr(X, "scaled:scale"))
  D <- outer(rowSums(Xt^2), rowSums(X^2), "+") - 2 * Xt %*% t(X)
  apply(D, 1, function(dd) mean(d$medv[order(dd)[seq_len(k)]]))
}
learners <- list(
  "fully grown tree" = deep_tree,
  "CV-pruned tree" = pruned_tree,
  "1-nearest neighbour" = function(d, te) knn(d, te, 1),
  "5-nearest neighbours" = function(d, te) knn(d, te, 5)
)

## -----------------------------------------------------------------------------
sizes <- c(0, 100, 300, 1000, 2000)
curve <- do.call(rbind, lapply(1:8, function(r) {
  set.seed(r)
  idx <- sample(nrow(Boston), 60)
  train <- Boston[idx, ]
  test <- Boston[-idx, ]
  syn <- generate(fit_synthesizer(train), n = max(sizes), seed = r)
  do.call(rbind, lapply(sizes, function(s) {
    d <- rbind(train, syn[seq_len(s), ])
    data.frame(rep = r, synthetic = s, learner = names(learners),
               rmse = vapply(learners, function(f) rmse(f(d, test), test$medv), numeric(1)))
  }))
}))
curve$learner <- factor(curve$learner, levels = names(learners))
avg <- aggregate(rmse ~ synthetic + learner, curve, mean)
avg$rmse <- round(avg$rmse, 2)
reshape(avg, idvar = "learner", timevar = "synthetic", direction = "wide")

## ----fig.height = 4.5---------------------------------------------------------
ggplot(avg, aes(synthetic, rmse)) +
  geom_line(colour = col_aug, linewidth = 0.8) +
  geom_point(colour = col_aug, size = 2.2) +
  facet_wrap(~learner, nrow = 1) +
  labs(x = "synthetic records added to 60 real ones", y = "test RMSE",
       title = "Augmentation helps learners that control their complexity") +
  theme_minimal(base_size = 10) +
  theme(panel.grid.minor = element_blank())

## -----------------------------------------------------------------------------
set.seed(42)
test_idx <- sample(nrow(Boston), 200)
test <- Boston[test_idx, ]
pool <- Boston[-test_idx, ]

R <- 20
pred_real <- pred_aug <- matrix(NA_real_, R, nrow(test))
for (r in seq_len(R)) {
  train <- pool[sample(nrow(pool), 60), ]
  pred_real[r, ] <- pruned_tree(train, test)
  pred_aug[r, ] <- pruned_tree(augment_data(train, n = 2000, seed = r), test)
}
decompose <- function(P) {
  c(bias2 = mean((colMeans(P) - test$medv)^2),
    variance = mean(apply(P, 2, var)),
    mse = mean(sweep(P, 2, test$medv)^2))
}
bv <- rbind(real = decompose(pred_real), augmented = decompose(pred_aug))
round(bv, 1)

## ----fig.height = 3.5---------------------------------------------------------
bv_long <- data.frame(
  training = factor(rep(c("real only", "real + synthetic"), 3),
                    levels = c("real only", "real + synthetic")),
  component = factor(rep(c("squared bias", "variance", "mean squared error"), each = 2),
                     levels = c("squared bias", "variance", "mean squared error")),
  value = as.vector(bv)
)
ggplot(bv_long, aes(component, value, fill = training)) +
  geom_col(position = position_dodge2(padding = 0.15), width = 0.7) +
  scale_fill_manual(values = c(col_real, col_aug), name = NULL) +
  labs(x = NULL, y = NULL, title = "CV-pruned tree trained on 60 records") +
  theme_minimal(base_size = 10) +
  theme(legend.position = "top", panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank())

## ----fig.height = 4-----------------------------------------------------------
grid <- data.frame(lstat = seq(2, 35, length.out = 200))
fits <- do.call(rbind, lapply(1:15, function(r) {
  set.seed(100 + r)
  train <- pool[sample(nrow(pool), 60), c("lstat", "medv")]
  aug <- augment_data(train, n = 2000, seed = r)
  rbind(data.frame(grid, medv = pruned_tree(train, grid), run = r, training = "real only"),
        data.frame(grid, medv = pruned_tree(aug, grid), run = r, training = "real + synthetic"))
}))
fits$training <- factor(fits$training, levels = c("real only", "real + synthetic"))
ggplot(fits, aes(lstat, medv, group = run, colour = training)) +
  geom_point(data = Boston, aes(lstat, medv), inherit.aes = FALSE,
             colour = "grey80", size = 0.6) +
  geom_step(alpha = 0.6, linewidth = 0.5) +
  facet_wrap(~training) +
  scale_colour_manual(values = c(col_real, col_aug), guide = "none") +
  labs(title = "15 trees, each trained on a different sample of 60 records") +
  theme_minimal(base_size = 10) +
  theme(panel.grid.minor = element_blank())

## -----------------------------------------------------------------------------
bagged_trees <- function(d, te, B = 25) {
  rowMeans(sapply(seq_len(B), function(b) {
    db <- d[sample(nrow(d), replace = TRUE), ]
    predict(rpart(medv ~ ., db, control = rpart.control(cp = 0, minsplit = 5, xval = 0)), te)
  }))
}
cmp <- t(sapply(1:8, function(r) {
  set.seed(r)
  idx <- sample(nrow(Boston), 60)
  train <- Boston[idx, ]
  test <- Boston[-idx, ]
  aug <- augment_data(train, n = 2000, seed = r)
  c(`pruned tree` = rmse(pruned_tree(train, test), test$medv),
    `pruned tree + synthetic` = rmse(pruned_tree(aug, test), test$medv),
    `bagged trees` = rmse(bagged_trees(train, test), test$medv),
    `bagged trees + synthetic` = rmse(bagged_trees(aug, test), test$medv))
}))
round(colMeans(cmp), 2)

## -----------------------------------------------------------------------------
cv_augmentation <- function(train, sizes, learner, folds = 5, seed = 1) {
  set.seed(seed)
  fold <- sample(rep_len(seq_len(folds), nrow(train)))
  err <- sapply(sizes, function(s) {
    mean(sapply(seq_len(folds), function(k) {
      fit_part <- train[fold != k, ]
      valid <- train[fold == k, ]
      d <- if (s == 0) fit_part else augment_data(fit_part, n = s, seed = k)
      rmse(learner(d, valid), valid$medv)
    }))
  })
  data.frame(synthetic = sizes, cv_rmse = err)
}
set.seed(7)
train <- Boston[sample(nrow(Boston), 60), ]
cv_augmentation(train, sizes = c(0, 500, 2000), learner = pruned_tree)