En términos de líneas de tiempo de aprendizaje automático, los bosques aleatorios (RF), introducidos en el artículo fundamental de Breimann ([1]), son antiguos. A pesar de su edad, siguen impresionando por su rendimiento y son objeto de investigación activa. El objetivo de este artículo es resaltar en qué caja de herramientas versátil se han convertido los métodos de Random Forest, centrándose en Bosque aleatorio generalizado (GRF) y Bosque aleatorio distributivo (DRF).
En resumen, la idea principal que subyace a ambos métodos es que las ponderaciones implícitamente producidas por RF pueden usarse para estimar objetivos distintos de la expectativa condicional. La idea de GRF es utilizar un bosque aleatorio con un criterio de división que se adapta al objetivo que se tiene en mente (por ejemplo, media condicional, cuantiles condicionales o efecto de tratamiento condicional). La idea de DRF es adaptar el criterio de división de manera que se pueda estimar toda la distribución condicional. A partir de este objeto se pueden derivar en un segundo paso muchos objetivos diferentes. De hecho, en este artículo hablo principalmente de DRF, ya que estoy más familiarizado con este método y es algo más sencillo (solo se debe adaptar un bosque para una amplia gama de objetivos). Sin embargo, todas las ventajas indicadas en la figura anterior también se aplican a GRF y, de hecho, el paquete DRF en R se basa en el implementación profesional de GRF. Además, el hecho de que el criterio de división de los bosques GRF se adapte al objetivo significa que puede tener un mejor rendimiento que el DRF. Esto es particularmente cierto para el sistema binario. Y, donde se debe utilizar probabilidad_forests(). Entonces, aunque hablo principalmente de DRF, se debe tener en cuenta GRF a lo largo de este artículo.
El objetivo de este artículo es proporcionar una descripción general con enlaces a una lectura más profunda en las secciones correspondientes. Revisaremos cada uno de los puntos de la figura anterior en el sentido de las agujas del reloj, haremos referencia a los artículos correspondientes y los resaltaremos con un pequeño ejemplo. Primero, resumo rápidamente los enlaces más importantes para leer más a continuación:
Versatilidad/Rendimiento: Artículo mediano y artículos originales (DRF/GRF)
Valores faltantes incorporados: Artículo mediano
Medidas de incertidumbre: Artículo mediano
Importancia variable: Artículo mediano
El ejemplo
Nosotros tomamos X_1, X_2, X_4,…, X_10 uniforme independientemente entre (-1,1) y crea dependencia entre X_1 y X_3 tomando X_3=X_1 + error uniforme. Luego simulamos Y como
## Load packages and functions needed
library(drf)
library(mice)
source("drfnew_v2.R")
## The function drfnew_v2.R is available below, or on
## https://github.com/JeffNaef/drfupdate## Set parameters
set.seed(10)
n<-1000
##Simulate Data that experiences both a mean as well as sd shift
# Simulate from X
x1 <- runif(n,-1,1)
x2 <- runif(n,-1,1)
x3 <- x1+ runif(n,-1,1)
X0 <- matrix(runif(7*n,-1,1), nrow=n, ncol=7)
Xfull <- cbind(x1,x2, x3, X0)
colnames(Xfull)<-paste0("X", 1:10)
# Simulate dependent variable Y
Y <- as.matrix(rnorm(n,mean = 0.8*(x1 > 0), sd = 1 + 1*(x2 > 0)))
##Also add MAR missing values using ampute from the mice package
X<-ampute(Xfull)$amp
head(cbind(Y,X))
Y X1 X2 X3 X4 X5
1 -3.0327466 -0.4689827 0.06161759 0.27462737 NA -0.624463079
2 1.2582824 -0.2557522 0.36972181 NA -0.04100963 0.009518047
3 -0.8781940 0.1457067 -0.23343321 NA -0.64519687 -0.945426305
4 3.1595623 0.8164156 0.90997600 0.69184618 -0.20573331 -0.007404298
5 1.1176545 -0.5966361 NA -1.21276055 0.62845399 0.894703422
6 -0.4377359 0.7967794 -0.92179989 -0.03863182 0.88271415 -0.237635732
X6 X7 X8 X9 X10
1 -0.9290009 0.5401628 0.39735433 -0.7434697 0.8807558
2 -0.2885927 0.3805251 -0.09051334 -0.7446170 0.9935311
3 -0.5022541 0.3009541 0.29424395 0.5554647 -0.5341800
4 0.7583608 -0.8506881 0.22758566 -0.1596993 -0.7161976
5 -0.3640422 0.8051613 -0.46714833 0.4318039 -0.8674060
6 -0.3577590 -0.7341207 0.85504668 -0.6933918 0.4656891
Observe que con la función amputar desde el paquete de ratonesnosotros ponemos Desaparecido no al azar (MAR) valores faltantes en X resaltar la capacidad de GRF/DRF para abordar los valores faltantes. Además, en el proceso anterior sólo X_1 y X_2 son relevantes para predecir Y, todas las demás variables son variables de “ruido”. Una configuración tan “escasa” podría en realidad ser común en conjuntos de datos de la vida real.
Ahora elegimos un punto de prueba para este ejemplo que usaremos en todo momento:
x<-matrix(c(0.2, 0.4, runif(8,-1,1)), nrow=1, ncol=10)
print(x)[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.2 0.4 0.7061058 0.8364877 0.2284314 0.7971179 0.78581 0.5310279
[,9] [,10]
[1,] -0.5067102 0.6918785
Versatilidad
DRF estima la distribución condicional P_ {Y|X=X} en forma de pesos simples:
A partir de estos pesos, se puede calcular una amplia gama de objetivos o se pueden utilizar para simular a partir de la distribución condicional. Una buena referencia por su versatilidad es el original. Artículo de investigacióndonde se utilizaron muchos ejemplos, así como los correspondientes artículo mediano.
En el ejemplo, primero simulamos a partir de esta distribución:
DRF<-drfCI(X=X, Y=Y, B=50,num.trees=1000, min.node.size = 5)
DRFpred<-predictdrf(DRF, newdata=x)## Sample from P_{Y| X=x}
Yxs<-Y[sample(1:n, size=n, replace = T, DRFpred$weights[1,])]
hist(Yxs, prob=T)
z<-seq(-6,7,by=0.01)
d<-dnorm(z, mean=0.8 * (x[1] > 0), sd=(1+(x[2] > 0)))
lines(z,d, col="darkred" )
El gráfico muestra los sorteos aproximados de la distribución condicional superpuestos con la densidad real en rojo. Ahora usamos esto para estimar el expectativa condicional y el cuantiles condicionales (0,05, 0,95) en X:
# Calculate quantile prediction as weighted quantiles from Y
qx <- quantile(Yxs, probs = c(0.05,0.95))# Calculate conditional mean prediction
mux <- mean(Yxs)
# True quantiles
q1<-qnorm(0.05, mean=0.8 * (x[1] > 0), sd=(1+(x[2] > 0)))
q2<-qnorm(0.95, mean=0.8 * (x[1] > 0), sd=(1+(x[2] > 0)))
mu<-0.8 * (x[1] > 0)
hist(Yxs, prob=T)
z<-seq(-6,7,by=0.01)
d<-dnorm(z, mean=0.8 * (x[1] > 0), sd=(1+(x[2] > 0)))
lines(z,d, col="darkred" )
abline(v=q1,col="darkred" )
abline(v=q2, col="darkred" )
abline(v=qx[1], col="darkblue")
abline(v=qx[2], col="darkblue")
abline(v=mu, col="darkred")
abline(v=mux, col="darkblue")
Del mismo modo, se pueden calcular muchos objetivos con GRF, sólo que en este caso para cada uno de esos dos objetivos sería necesario ajustar un bosque diferente. En particular, bosque_regresión() para la expectativa condicional y bosque_cuantil() para los cuantiles.
Lo que GRF no puede hacer es abordar objetivos multivariados, lo que también es posible con DRF.
Actuación
A pesar de todo el trabajo realizado en métodos nuevos y potentes (no paramétricos), como las redes neuronales, los métodos basados en árboles son consistentemente capaces de vencer a los competidores en datos tabulares. Ver por ejemplo, esto papel fascinanteo este artículo más antiguo sobre el fuerza de RF en la clasificación.
Para ser justos, con el ajuste de parámetros, métodos de árbol impulsados, como XGboost, a menudo toman la delantera, al menos cuando se trata de predicción clásica (que corresponde a la estimación de expectativas condicionales). No obstante, es notable el rendimiento robusto que tienden a tener los métodos de RF sin ningún ajuste. Además, también se ha trabajado para mejorar el rendimiento de los bosques aleatorios, por ejemplo, el enfoque de bosque aleatorio cubierto.
Valores faltantes incorporados
“Criterio faltante incorporado en atributos” (MIA) de este papel es una idea muy simple pero muy poderosa que permite que los métodos basados en árboles manejen los datos faltantes. Esto se implementó en el paquete GRF R y por eso también está disponible en DRF. Los detalles también se explican en este artículo mediano. Por simple que sea el concepto, funciona notablemente bien en la práctica: en el ejemplo anterior, DRF no tuvo problemas para manejar una falta sustancial de MAR en los datos de entrenamiento. X (!)
Medidas de incertidumbre
Como estadístico no quiero sólo estimaciones puntuales (incluso de una distribución), sino también una medida de incertidumbre en la estimación de mis parámetros (incluso si el “parámetro” es toda mi distribución). Resulta que un submuestreo adicional simple integrado en DRF/GRF permite una cuantificación de la incertidumbre basada en principios para tamaños de muestra grandes. La teoría detrás de esto en el caso de DRF se deriva de este Artículo de investigaciónpero también lo explico en este artículo mediano. GRF tiene toda la teoría en el papel original.
Adaptamos esto para el ejemplo anterior:
# Calculate uncertainty
alpha<-0.05
B<-length(DRFpred$weightsb)
qxb<-matrix(NaN, nrow=B, ncol=2)
muxb<-matrix(NaN, nrow=B, ncol=1)
for (b in 1:B){
Yxsb<-Y[sample(1:n, size=n, replace = T, DRFpred$weightsb[[b]][1,])]
qxb[b,] <- quantile(Yxsb, probs = c(0.05,0.95))
muxb[b] <- mean(Yxsb)
}
CI.lower.q1 <- qx[1] - qnorm(1-alpha/2)*sqrt(var(qxb[,1]))
CI.upper.q1 <- qx[1] + qnorm(1-alpha/2)*sqrt(var(qxb[,1]))CI.lower.q2 <- qx[2] - qnorm(1-alpha/2)*sqrt(var(qxb[,2]))
CI.upper.q2 <- qx[2] + qnorm(1-alpha/2)*sqrt(var(qxb[,2]))
CI.lower.mu <- mux - qnorm(1-alpha/2)*sqrt(var(muxb))
CI.upper.mu <- mux + qnorm(1-alpha/2)*sqrt(var(muxb))
hist(Yxs, prob=T)
z<-seq(-6,7,by=0.01)
d<-dnorm(z, mean=0.8 * (x[1] > 0), sd=(1+(x[2] > 0)))
lines(z,d, col="darkred" )
abline(v=q1,col="darkred" )
abline(v=q2, col="darkred" )
abline(v=qx[1], col="darkblue")
abline(v=qx[2], col="darkblue")
abline(v=mu, col="darkred")
abline(v=mux, col="darkblue")
abline(v=CI.lower.q1, col="darkblue", lty=2)
abline(v=CI.upper.q1, col="darkblue", lty=2)
abline(v=CI.lower.q2, col="darkblue", lty=2)
abline(v=CI.upper.q2, col="darkblue", lty=2)
abline(v=CI.lower.mu, col="darkblue", lty=2)
abline(v=CI.upper.mu, col="darkblue", lty=2)
Como se puede ver en el código anterior, esencialmente tenemos B subárboles que se pueden utilizar para calcular la medida cada vez. A partir de estos B muestras de media y cuantiles, luego podemos calcular las varianzas y usar una aproximación normal para obtener los intervalos de confianza (asintóticos) que se ven en la línea discontinua de la Figura. Nuevamente, todo esto se puede hacer a pesar de los valores faltantes en X(!)
Importancia variable
Un último aspecto importante de Random Forests son las medidas de importancia variable calculadas de manera eficiente. Si bien las medidas tradicionales son algo ad hoc, para el FR tradicional y el DRF, ahora hay medidas basadas en principios disponibles, como se explica en este artículo mediano. Para RF, el método Sobol-MDA identifica de manera confiable las variables más importantes para la estimación de expectativas condicionales, mientras que para DRF, el MMD-MDA identifica las variables más importantes para la estimación de la distribución general. Como se analiza en el artículo, utilizando la idea de Bosques aleatorios proyectados, estas medidas pueden implementarse de manera muy eficiente. Demostramos esto en el ejemplo con una implementación menos eficiente de la medida de importancia de la variable MMD:
## Variable importance for conditional Quantile Estimation## For the conditional quantiles we use a measure that considers the whole distribution,
## i.e. the MMD based measure of DRF.
MMDVimp <- compute_drf_vimp(X=X,Y=Y, print=F)
sort(MMDVimp, decreasing = T)
X2 X1 X8 X6 X3 X10
0.852954299 0.124110913 0.012194176 0.009578300 0.008191663 0.007517931
X9 X7 X5 X4
0.006861688 0.006632175 0.005257195 0.002401974
aqui ambos X_1 y X_2 se identifican correctamente como la variable más relevante al intentar estimar la distribución. Sorprendentemente, a pesar de la dependencia de X_3 y X_1la medida cuantifica correctamente que X_3 no es importante para la predicción de la distribución de Y. Esto es algo que la medida MDA original de Random Forests tiende a hacer mal, como se demuestra en el artículo mediano. Además, observe nuevamente que los valores faltantes en X No hay problema aquí.
Conclusión
GRF/DRF y también el tradicional Random Forest no deberían faltar en la caja de herramientas de ningún científico de datos. Si bien métodos como XGboost pueden tener un mejor rendimiento en la predicción tradicional, las muchas fortalezas de los enfoques modernos basados en RF los convierten en una herramienta increíblemente versátil.
Por supuesto, hay que tener en cuenta que estos métodos todavía son totalmente no paramétricos y que se necesitan muchos puntos de datos para que el ajuste tenga sentido. Esto es particularmente cierto para la cuantificación de la incertidumbre, que sólo es válida asintóticamente, es decir, para muestras “grandes”.
Literatura
[1] Breiman, L. (2001). Bosques aleatorios. Aprendizaje automático, 45(1):5–32.
Apéndice: Código
require(drf)
require(Matrix)
require(kernlab)drfCI <- function(X, Y, B, sampling = "binomial", ...) {
n <- dim(X)[1]
# compute point estimator and DRF per halfsample
# weightsb: B times n matrix of weights
DRFlist <- lapply(seq_len(B), function(b) {
# half-sample index
indexb <- if (sampling == "binomial") {
seq_len(n)[as.logical(rbinom(n, size = 1, prob = 0.5))]
} else {
sample(seq_len(n), floor(n / 2), replace = FALSE)
}
## Using normal Bootstrap on the data and refitting DRF
DRFb <-
drfown(X = X[indexb, , drop = F], Y = Y[indexb, , drop = F], ...)
return(list(DRF = DRFb, indices = indexb))
})
return(list(DRFlist = DRFlist, X = X, Y = Y))
}
predictdrf <- function(DRF, newdata, functional = NULL, ...) {
##Predict for testpoints in newdata, with B weights for each point, representing
##uncertainty
ntest <- nrow(x)
n <- nrow(DRF$Y)
weightsb <- lapply(DRF$DRFlist, function(l) {
weightsbfinal <- Matrix(0, nrow = ntest, ncol = n, sparse = TRUE)
weightsbfinal[, l$indices] <- predict(l$DRF, x)$weights
return(weightsbfinal)
})
weightsall <- Reduce("+", weightsb) / length(weightsb)
if (!is.null(functional)) {
stopifnot("Not yet implemented for several x" = ntest == 1)
thetahatb <-
lapply(weightsb, function(w)
functional(weights = w, X = DRF$X, Y = DRF$Y, x = x))
thetahatbforvar <- do.call(rbind, thetahatb)
thetahat <- functional(weights = weightsall, X = DRF$X, Y = DRF$Y, x = x)
thetahat <- matrix(thetahat, nrow = dim(x)[1])
var_est <- if (dim(thetahat)[2] > 1) {
a <- sweep(thetahatbforvar, 2, thetahat, FUN = "-")
crossprod(a, a) / B
} else {
mean((c(thetahatbforvar) - c(thetahat)) ^ 2)
}
return(list(weights = weightsall, thetahat = thetahat, weightsb = weightsb,
var_est = var_est))
} else {
return(list(weights = weightsall, weightsb = weightsb))
}
}
drfown <- function(X, Y,
num.trees = 500,
splitting.rule = "FourierMMD",
num.features = 10,
bandwidth = NULL,
response.scaling = TRUE,
node.scaling = FALSE,
sample.weights = NULL,
sample.fraction = 0.5,
mtry = min(ceiling(sqrt(ncol(X)) + 20), ncol(X)),
min.node.size = 15,
honesty = TRUE,
honesty.fraction = 0.5,
honesty.prune.leaves = TRUE,
alpha = 0.05,
imbalance.penalty = 0,
compute.oob.predictions = TRUE,
num.threads = NULL,
seed = stats::runif(1, 0, .Machine$integer.max),
compute.variable.importance = FALSE) {
# initial checks for X and Y
if (is.data.frame(X)) {
if (is.null(names(X))) {
stop("the regressor should be named if provided under data.frame format.")
}
if (any(apply(X, 2, class) %in% c("factor", "character"))) {
any.factor.or.character <- TRUE
X.mat <- as.matrix(fastDummies::dummy_cols(X, remove_selected_columns = TRUE))
} else {
any.factor.or.character <- FALSE
X.mat <- as.matrix(X)
}
mat.col.names.df <- names(X)
mat.col.names <- colnames(X.mat)
} else {
X.mat <- X
mat.col.names <- NULL
mat.col.names.df <- NULL
any.factor.or.character <- FALSE
}
if (is.data.frame(Y)) {
if (any(apply(Y, 2, class) %in% c("factor", "character"))) {
stop("Y should only contain numeric variables.")
}
Y <- as.matrix(Y)
}
if (is.vector(Y)) {
Y <- matrix(Y,ncol=1)
}
#validate_X(X.mat)
if (inherits(X, "Matrix") && !(inherits(X, "dgCMatrix"))) {
stop("Currently only sparse data of class 'dgCMatrix' is supported.")
}
drf:::validate_sample_weights(sample.weights, X.mat)
#Y <- validate_observations(Y, X)
# set legacy GRF parameters
clusters <- vector(mode = "numeric", length = 0)
samples.per.cluster <- 0
equalize.cluster.weights <- FALSE
ci.group.size <- 1
num.threads <- drf:::validate_num_threads(num.threads)
all.tunable.params <- c("sample.fraction", "mtry", "min.node.size", "honesty.fraction",
"honesty.prune.leaves", "alpha", "imbalance.penalty")
# should we scale or not the data
if (response.scaling) {
Y.transformed <- scale(Y)
} else {
Y.transformed <- Y
}
data <- drf:::create_data_matrices(X.mat, outcome = Y.transformed, sample.weights = sample.weights)
# bandwidth using median heuristic by default
if (is.null(bandwidth)) {
bandwidth <- drf:::medianHeuristic(Y.transformed)
}
args <- list(num.trees = num.trees,
clusters = clusters,
samples.per.cluster = samples.per.cluster,
sample.fraction = sample.fraction,
mtry = mtry,
min.node.size = min.node.size,
honesty = honesty,
honesty.fraction = honesty.fraction,
honesty.prune.leaves = honesty.prune.leaves,
alpha = alpha,
imbalance.penalty = imbalance.penalty,
ci.group.size = ci.group.size,
compute.oob.predictions = compute.oob.predictions,
num.threads = num.threads,
seed = seed,
num_features = num.features,
bandwidth = bandwidth,
node_scaling = ifelse(node.scaling, 1, 0))
if (splitting.rule == "CART") {
##forest <- do.call(gini_train, c(data, args))
forest <- drf:::do.call.rcpp(drf:::gini_train, c(data, args))
##forest <- do.call(gini_train, c(data, args))
} else if (splitting.rule == "FourierMMD") {
forest <- drf:::do.call.rcpp(drf:::fourier_train, c(data, args))
} else {
stop("splitting rule not available.")
}
class(forest) <- c("drf")
forest[["ci.group.size"]] <- ci.group.size
forest[["X.orig"]] <- X.mat
forest[["is.df.X"]] <- is.data.frame(X)
forest[["Y.orig"]] <- Y
forest[["sample.weights"]] <- sample.weights
forest[["clusters"]] <- clusters
forest[["equalize.cluster.weights"]] <- equalize.cluster.weights
forest[["tunable.params"]] <- args[all.tunable.params]
forest[["mat.col.names"]] <- mat.col.names
forest[["mat.col.names.df"]] <- mat.col.names.df
forest[["any.factor.or.character"]] <- any.factor.or.character
if (compute.variable.importance) {
forest[['variable.importance']] <- variableImportance(forest, h = bandwidth)
}
forest
}
#' Variable importance for Distributional Random Forests
#'
#' @param X Matrix with input training data.
#' @param Y Matrix with output training data.
#' @param X_test Matrix with input testing data. If NULL, out-of-bag estimates are used.
#' @param num.trees Number of trees to fit DRF. Default value is 500 trees.
#' @param silent If FALSE, print variable iteration number, otherwise nothing is print. Default is FALSE.
#'
#' @return The list of importance values for all input variables.
#' @export
#'
#' @examples
compute_drf_vimp <- function(X, Y, X_test = NULL, num.trees = 500, silent = FALSE){
# fit initial DRF
bandwidth_Y <- drf:::medianHeuristic(Y)
k_Y <- rbfdot(sigma = bandwidth_Y)
K <- kernelMatrix(k_Y, Y, Y)
DRF <- drfown(X, Y, num.trees = num.trees)
wall <- predict(DRF, X_test)$weights
# compute normalization constant
wbar <- colMeans(wall)
wall_wbar <- sweep(wall, 2, wbar, "-")
I0 <- as.numeric(sum(diag(wall_wbar %*% K %*% t(wall_wbar))))
# compute drf importance dropping variables one by one
I <- sapply(1:ncol(X), function(j) {
if (!silent){print(paste0('Running importance for variable X', j, '...'))}
DRFj <- drfown(X = X[, -j, drop=F], Y = Y, num.trees = num.trees)
DRFpredj <- predict(DRFj, X_test[, -j])
wj <- DRFpredj$weights
Ij <- sum(diag((wj - wall) %*% K %*% t(wj - wall)))/I0
return(Ij)
})
# compute retraining bias
DRF0 <- drfown(X = X, Y = Y, num.trees = num.trees)
DRFpred0 = predict(DRF0, X_test)
w0 <- DRFpred0$weights
vimp0 <- sum(diag((w0 - wall) %*% K %*% t(w0 - wall)))/I0
# compute final importance (remove bias & truncate negative values)
vimp <- sapply(I - vimp0, function(x){max(0,x)})
names(vimp)<-colnames(X)
return(vimp)
}