In zmjones/fanova: Functional ANOVA Decomposition of Prediction Functions

library(fanova)
library(data.table)
library(ggplot2)
library(randomForest)

data(Boston, package = "MASS")
boston.features = Boston[, -ncol(Boston)]
boston.target = Boston$medv

model = randomForest(boston.features, boston.target)

fa = functionalANOVA(boston.features, c("rm", "lstat"), c(10, 2), model)
print(fa)

plt = melt(fa[fa$effect %in% c("rm", "lstat"), ],
  id.vars = c("f", "effect"), na.rm = TRUE)
ggplot(plt, aes(value, f)) +
  geom_point() + geom_line() + facet_wrap(~ variable, scales = "free_x")

ggplot(fa[fa$effect == "rm:lstat", ], aes(rm, lstat, z = f, fill = f)) + geom_raster()

## from section 9 of the paper
data = data.frame(
  x = runif(100),
  y = runif(100),
  z = runif(100)
)

F = function(x, y, z) x * y + x * z + y * z
G = function(x, y, z) 2 * ifelse(x > .5 & y > .5 & z > .5, 1, 0)
predict.fun = function(object, newdata)
  F(newdata$x, newdata$y, newdata$z) + G(newdata$x, newdata$y, newdata$z)

no.weight = function(design, data) rep(1, nrow(design))
corner.weight = function(design, data)
  ifelse(design$x > .5 & design$y > .5 & design$z > .5, 0, 1)

vars = c("x", "y")
n = c(10, 2)
plt = list(
  "unweighted" = functionalANOVA(data, vars, n, NULL, predict.fun, no.weight),
  "corner-deleted" = functionalANOVA(data, vars, n, NULL, predict.fun, corner.weight)
)
plt = rbindlist(plt, idcol = "method")
ggplot(plt[plt$effect == "x:y", ], aes(x, y, fill = f, z = f)) +
  geom_raster() + 
  facet_wrap(~ method)