Quando Breiman introduziu o algoritmo Random Forest (RF) em 2001, ele sabia o tremendo efeito que teria? Atualmente, a RF é uma ferramenta muito utilizada em muitas partes da ciência de dados. É fácil entender por quê — o RF é fácil de usar e conhecido por seu alto desempenho em uma gama extremamente ampla de tarefas. Isso por si só é impressionante, mas o que o torna ainda mais interessante é que nenhum ajuste geralmente é necessário para obter esses resultados.

Neste artigo, discutimos uma extensão bastante massiva do RF original, a Distributional Random Forest (DRF) que desenvolvemos recentemente neste artigo:

Florestas aleatórias distributivas: ajuste de heterogeneidade e regressão distributiva multivariada (jmlr.org)

trabalho conjunto com Loris Michel , que também contribuiu para este artigo. Código adicional e uma implementação python podem ser encontrados aqui:https://github.com/lorismichel/drf

Começaremos com uma pequena introdução ao tópico Random Forests e depois seguiremos as etapas para desenvolver o DRF. Uma pequena lista de exemplos aguarda no final, que estendemos alegremente no futuro.

Antes de prosseguirmos, primeiro precisamos tirar algumas notações do caminho: Neste artigo, nos referiremos vagamente a X , Y como vetores aleatórios a partir dos quais supomos observar uma amostra iid (yi, x i), i= 1,… n. Isso é bastante comum em trabalhos de pesquisa. Se, em vez disso, lidarmos com variáveis ​​aleatórias (ou seja, univariadas), apenas escrevemos Y ou X. O objetivo principal em muitas aplicações da ciência de dados é prever Y a partir de alguns recursos X , digamos que queremos prever um rótulo “cachorro” (Y) de digamos um vetor de intensidades de pixel ( X ). Formulada com menos glamour, uma previsão é, na verdade, quase sempre uma função da distribuição condicional P(Y| X = x). Então perguntamos “se o vetor X é fixo em um valor x , então o que é um certo aspecto distributivo de Y?”. O exemplo mais proeminente é a média ou expectativa condicional, ou seja, o valor esperado sob a distribuição P(Y| X = x ). Isso é o que as pessoas geralmente querem quando pedem uma previsão. Por exemplo, se prevermos com regressão linear, inserimos um valor x na fórmula linear e obtemos um valor esperado de Y como a previsão. O exemplo de classificação da imagem acima já é diferente — aqui podemos perguntar pelo valor mais provável, ou seja, para qual rótulo l é P(Y= l | X = x) maximizado. Resumindo: a maioria das previsões são aspectos de uma distribuição condicional, geralmente expectativa condicional.

Curiosamente, as pessoas descobriram mais tarde que o RF pode ser usado não apenas para expectativa condicional ou previsões médias, mas também para prever outros aspectos de uma distribuição (condicional). Por exemplo, pode-se prever quantis e, assim, obter intervalos de previsão medindo a incerteza. O DRF leva essa abordagem a um novo nível: você pode alimentá-lo com uma resposta multivariada Y junto com o vetor típico de recursos X e ele estimará toda a distribuição condicional do Y multivariado dado qualquer (realização) x . Essa estimativa é dada na forma de pesos simples, facilitando simular a partir da distribuição ou calcular quantidades de interesse, ou seja, “previsões”.

Uma explicação completa de RF é dada aqui, por exemplo, então apenas tocamos nisso brevemente: Em um RF, várias árvores de decisão independentes são ajustadas. Cada árvore obtém apenas uma parte da variável e, em seguida, divide o espaço resultante de acordo com os recursos em X .

Assim, quando "descemos" um novo ponto x , ele terminará em uma folha para cada árvore. Uma folha é um conjunto com observações i e tirando a média de todos os yi naquela folha dá a previsão para uma árvore. Essas previsões são então calculadas para dar o resultado final. Assim, para um dado x , se você quiser prever a média condicional de Y dado que x , você:

1. “Drop down” x cada árvore (isso é indicado em vermelho na figura acima). Como as regras de divisão foram feitas em X , seu novo ponto x ficará seguro em algum lugar de um nó folha.

2. Para cada árvore, você calcula a média das respostas yi naquela folha para obter uma estimativa da média condicional de cada árvore.

3. Você calcula a média de cada média condicional sobre as árvores para obter a previsão final.

A média importante da previsão de todas as árvores leva a uma redução acentuada na variância. Os erros ficam "em média", por assim dizer. É por isso que é importante tornar as árvores o mais independentes possível.

Agora trabalhamos para uma extensão dessa ideia para não apenas fazer uma previsão média, mas, em vez disso, prever toda a distribuição (condicional). Ou seja, queremos prever P( Y | X = x ) para um ponto de teste x . Como já mencionado, muitas coisas que geralmente se deseja prever são funções diretas dessa distribuição condicional: a média (condicional), quantis condicionais, variâncias condicionais, covariâncias e assim por diante. Acontece que a forma do DRF também o torna adequado para análise causal, além do que é possível até mesmo para métodos de ponta. Damos um pequeno exemplo no final deste documento e potencialmente muito mais em postagens posteriores.

Por enquanto, derivamos passo a passo como chegar à Floresta Aleatória Distributiva. O primeiro passo é obter os pesos.

A teoria

Passo 1: Ganhar algum peso

A chave para isso é obter os pesos que o RF oferece. Isso geralmente é esquecido em uma primeira introdução ao RF, mas é vital aqui:

Em vez de calcular diretamente a média em um nó folha, calculam-se os pesos que são usados ​​implicitamente ao fazer o cálculo da média. O peso wi( x ) é uma função de (1) o ponto de teste xe (2) uma observação i . Ou seja, se diminuirmos xdescemos uma árvore, observamos em que folha ela termina. Todas as observações que estão nessa folha recebem 1, todas as outras 0. Portanto, se terminarmos em uma folha com observações (1,3,10), o peso das observações 1,3,10 para essa árvore é 1, enquanto todas as outras observações recebem 0. Em seguida, dividimos esse peso pelo número de elementos no nó folha. No exemplo anterior, tínhamos 3 observações no nó folha, então os pesos para as observações 1,3, 10 são 1/3 cada, enquanto todas as outras observações ainda recebem peso 0 nesta árvore. A média desses pesos sobre todas as árvores dá o peso final wi( x ).

Por que esses pesos são úteis? Primeiro, pode-se mostrar que calcular a média, como fazemos em uma Random Forest tradicional, é o mesmo que somar wi( x )*yi. Isso parece um pouco feio, mas basicamente se trata apenas de trocar duas somas. Portanto, se essa estimativa for sensata, deve fazer sentido que, se você tomar outras funções de yi, digamos yi², obtenha a média de Y² dado X . Ou você, para algum número t, poderia usar I(yi<t) que é 1 se yi < t e 0 caso contrário e ver a soma sobre wi( x)I(yi<t) como uma estimativa da probabilidade condicional de que Y < t. Assim, você pode obter diretamente a função de distribuição cumulativa condicional (cdf), que caracteriza totalmente a distribuição condicional em teoria. Em segundo lugar, esses pesos realmente fornecem uma estimativa do vizinho mais próximo sobre observações que são semelhantes em certo sentido à consulta x . Como veremos a seguir, essa similaridade é medida de forma que, idealmente, seja possível ver os yi's associados como uma amostra iid da distribuição condicional.

Isso pode ser levado ainda mais longe. Em particular, se você amostrar de Y com esses pesos, obterá uma amostra aproximada da distribuição condicional. Antes de continuarmos, vamos fazer um rápido exemplo usando o famoso pacote ranger em R. O pacote não fornece esses pesos diretamente, então precisamos implementá-los manualmente, o que é um bom exercício aqui:

library(ranger)
library(ggplot2)
n<-1000
## Simulate some data:
X<-matrix(rnorm(n*2),ncol=2)
Y=2*X[,1] + 0.5*X[,2]^2 + rnorm(n)
Y<-as.matrix(Y)
## Get a test point
x<-matrix(c(1,-0.5),ncol=2)
# Fit Random forest
RF<-ranger(Y~., data=data.frame(Y=Y, X=X),min.node.size = 10, num.trees=500)
## Get the leaf/terminal node indices of each observation in the training samples
terminalX <- predict(RF, data=data.frame(X=X), type="terminalNodes")$predictions
## Get the leaf/terminal node of the test point
terminalxnew <- predict(RF, data=data.frame(X=x), type="terminalNodes")$predictions
## For the leafs in which x ends up, what are the number of observations (we need to normalize by that)
divid<-sapply(1:ncol(terminalxnew[1,, drop=F]), function(j) {sum(terminalX[,j]==terminalxnew[1,j])} )
  
# Average the weights over the trees:
weights<-rowSums(t(sapply(1:nrow(terminalX), function(j) as.numeric(terminalX[j,]==terminalxnew[1,])/divid )), na.rm=T)/length(terminalxnew[1,])
## We can now sample according to those weights from Y
Ysample<-Y[sample(1:dim(Y)[1], prob=weights, replace = T),]
## True conditional distribution
Ytrue<-2*x[,1] + 0.5*x[,2]^2 + rnorm(n) 
p1 <- hist(Ytrue, freq=F)                     
p2 <- hist(Ysample, freq=F)                     
plot( p2, col=rgb(0,0,1,1/4), freq=F)  
plot( p1, col=rgb(1,0,0,1/4), add=T, freq=F)
## Plot the weights by looking at the point x and its nearest neighbors according to the weights
ggplot(data = data.frame(X=X), aes(x = X.1, y = X.2)) +
  geom_point(aes(alpha = exp(weights ))) + geom_point(data=data.frame(X=x), color = "red")
## Get the estimated conditional mean:
mean(Ytrue)
mean(Ysample)
sum(weights*Y)

      
                

                    
                    
                
            

Há outra maneira de olhar para o algoritmo Random Forest: é uma máquina de homogeneidade. Em cada divisão em uma árvore, a divisão em X é escolhida de forma que as duas amostras de Y nos nós resultantes sejam tão ''diferentes'' quanto possível. A figura abaixo mostra um pequeno exemplo para X e Y univariados.

Nesse exemplo, uma determinada árvore provavelmente dividirá X em S como na figura. Então todo Yi com Xi < S será lançado no nó 1 e todo Xi >= S no nó 2, identificando assim grupos de dados. Se fizermos isso por tempo suficiente, cada nó folha terá uma amostra muito homogênea de Y, ou em outras palavras, todos os yi em um nó folha serão semelhantes de alguma forma, porque todas as observações ''dissimilares'' foram cortadas. Assim, para cada árvore, você agrupa seus dados em baldes de coisas semelhantes.

Por que isso faz sentido? Porque, em essência, RF é um algoritmo de vizinho mais próximo. Se você der um x , ele cairá das árvores até cair em um balde ou folha com observações “semelhantes”. A partir dessas observações, a média condicional é então calculada. Ou seja, em cada árvore são consideradas apenas as observações da folha, nada mais. Então é como um k-NN onde a distância não é medida pela distância euclidiana, mas decidida pela floresta. A floresta, por sua vez, decide rotular esses x i's como ''semelhantes'' que possuem yi's ''semelhantes''. então a semelhança de xi é decidido com base em seu yi associado, o que faz muito sentido se seu objetivo é inferir coisas sobre Y. Na verdade, até mesmo os métodos k-NN assumem algo como “a distribuição condicional P( Y | X = x i) para x i próximo (em distância euclidiana) de x é aproximadamente o mesmo”. A Figura abaixo mostra uma ilustração disso: Você pode ver para cada valor x i na amostra a distribuição condicional verdadeira associada, da qual yi é amostrado. Uma versão perfeita do DRF reconheceria que as distribuições condicionais para ( x 1, x 4, x 6) e ( x 3, x 5, x7) são semelhantes (não importa qual seja realmente a distância euclidiana) e tratam o grupo correspondente de yi, (y1, y4, y6) e (y3, y5, y7), cada um como uma amostra iid da distribuição condicional.

Idealmente, isso significaria que, na prática, as amostras yi homogêneas que acabamos dentro de uma folha são, na verdade, aproximadamente uma amostra iid da distribuição condicional Y| X=x . Isso é o que justifica tomar a média (ponderada).

Infelizmente, no RF original, essa abordagem não funciona como pretendido fora da previsão de média condicional. Novamente, o que queremos é um critério de divisão que torne a distribuição de Y nas duas divisões o mais diferente possível. O que obtemos no RF original é simplesmente uma divisão que torna a diferença de médias entre as duas amostras tão grande quanto possível. Na figura acima, tal método pode agrupar todos menos x 2 em um grupo, porque x 1, x 3, x 4, x 6, x7 todos têm meios muito semelhantes. Mas claro, como já é visível na foto acima, uma distribuição não se define pela sua média. Uma distribuição normal pode ter a mesma média, mas variância ou outros momentos muito diferentes. De um modo geral, você pode ter muitas distribuições com a mesma média, mas sendo muito diferentes.

A chave é que cada divisão deve depender de uma medida de diferença na distribuição entre os dois nós resultantes. Portanto, não apenas diferenças na média ou variância são detectadas, mas qualquer diferença na distribuição. O DRF resolve esse problema ajustando o critério de divisão normalmente empregado em um RF para aproveitar o poder teórico e prático dos kernelse o chamado critério MMD. O MMD pode ser aproximado com muita eficiência e é, em princípio, capaz de detectar qualquer diferença na distribuição. Teoricamente falando, nós assim enviamos cada ponto yi para um espaço de dimensão infinita, o Reproduzindo Kernel Hilbert Space, e realmente comparamos as médias naquele espaço. Através da mágica dos métodos do kernel, essa comparação entre médias é, na verdade, uma comparação de distribuições! Acontece que neste espaço especial os meios são distribuições. O que isso significa na prática é o seguinte: um nó folha conterá x i’s semelhantes, no sentido de que a distribuição de yi naquele balde é semelhante. Assim, se a distribuição condicional de Y dado x i e xj são semelhantes, eles serão agrupados no mesmo balde. Isso pode, em princípio, ser verdade, mesmo se x i e x j estiverem distantes no espaço euclidiano (ou seja, se não forem vizinhos mais próximos no sentido k-NN). Portanto, se usarmos esses pesos para calcular coisas condicionais nas quais estamos interessados, usaremos um método de vizinho mais próximo que considera x i e x j semelhantes, quando as distribuições de seus associados yi, yj são semelhantes. E, em particular, sob algumas suposições de suavidade, a amostra no nó folha x termina em aproximadamente uma amostra iid da distribuição P( Y | X = x ).

Passo 3: Use uma resposta multivariada

Esta etapa é realmente fácil, pois o MMD também permite comparar distribuições multivariadas. É importante ressaltar que distinguir mais do que apenas a média torna-se ainda mais importante para respostas multivariadas, pois as diferenças nas distribuições podem ser ainda mais complicadas. Por exemplo, duas distribuições multivariadas podem ter a mesma média e variância para todos os marginais, mas diferentes covariâncias entre os elementos.

Exemplos

Vamos fazer alguns pequenos exemplos. Aqui, o objetivo é apenas fornecer exemplos simulados muito simples para ter uma ideia do método. Primeiro, refazemos o que fizemos manualmente acima:

library(drf)
# Fit DRF
DRF<-drf(X,Y)
weights<-predict(DRF,x)$weights[1,]
### We can now sample according to those weights from Y
Ysample<-Y[sample(1:dim(Y)[1], prob=weights, replace = T),]
## True conditional distribution
Ytrue<-2*x[,1] + 0.5*x[,2]^2 + rnorm(n) 
p1 <- hist(Ytrue, freq=F)                     
p2 <- hist(Ysample, freq=F)                     
plot( p2, col=rgb(0,0,1,1/4), freq=F)  
plot( p1, col=rgb(1,0,0,1/4), add=T, freq=F)
ggplot(data = data.frame(X=X), aes(x = X.1, y = X.2)) +
  geom_point(aes(alpha = exp(weights ))) + geom_point(data=data.frame(X=x), color = "red")
## Get the estimated conditional mean:
mean(Ytrue)
mean(Ysample)
sum(weights*Y)

Também podemos prever quantis condicionais. Fazer isso, por exemplo, fornece um intervalo de previsão para os valores de Y| x , de modo que os valores extraídos dessa distribuição devem ser aproximadamente 95% do tempo no intervalo:

# Predict quantiles for a test point.
quants<-predict(DRF,x, functional = “quantile”, quantiles=c(0.025, 0.975))$quantile 
q1<-quants[1]
q2<-quants[2]
mean(Ytrue >=q1 & Ytrue <=q2)

Resposta Bidimensional

Aqui construímos um exemplo louco e difícil com uma resposta bidimensional e apenas calculamos uma variedade de previsões. Primeiro simulamos os dados:

n<-5000
d<-2
p<-3
## Simulate X and Y
X<-matrix( cbind(runif(n)*2, rnorm(n, mean = 1), rt(n, df=8)),nrow=n)
Y<-matrix(NA, ncol=2, nrow=n)
Y[,1]<-X[,1]^2 + X[,2]*X[,1] + X[,3] + rnorm(n)
Y[,2] <- 2*Y[,1]*X[,2] + 0.1*X[,3] + rnorm(n)

Agora pegamos dois pontos de teste e fazemos algumas previsões malucas, só porque podemos:

library(drf)
# Fit DRF
DRF<- drf(X=X, Y=Y, num.features=200)
# Choose a few test point:
x= matrix(c(0.093, -0.5,  1.37, 0.093, 0.5,  1.37) , ncol=p, nrow=2, byrow=T)
# mean prediction
(predict(DRF, newdata=x, functional="mean")$mean)
# correlation prediction
matrix(predict(DRF, newdata=x[1,], functional="cor")$cor, nrow=d,ncol=d)
matrix(predict(DRF, newdata=x[2,], functional="cor")$cor, nrow=d,ncol=d)
# Estimated probability that Y is smaller than 0:
weightstotal<-predict(DRF, newdata=x)$weights
p1<-sum(weightstotal[1,]* (Y[,1] + Y[,2]<= 0) )
p2<-sum(weightstotal[2,]* (Y[,1] + Y[,2] <= 0) )
# Bootstrapping the estimated probability of the sum being <=0 for both points:
B<-100
pb1<-matrix(NA, ncol=1,nrow=B)
pb2<-matrix(NA, ncol=1,nrow=B)
for (b in 1:B){
  Ybx1<-Y[sample(1:n, size=n, replace=T, prob=weightstotal[1,]), ]
  Ybx2<-Y[sample(1:n, size=n, replace=T, prob=weightstotal[2,]), ]
  pb1[b] <- mean(Ybx1[,1] + Ybx1[,2] <= 0)
  pb2[b] <- mean(Ybx2[,1] + Ybx2[,2] <= 0)
}
ggplot(data.frame(x=1:2, y=c(p1, p2)), aes(x = x, y = y)) +
  geom_point(size = 4) +
  geom_errorbar(aes(ymax = c(quantile(pb1,1- 0.025), quantile(pb2, 1-0.025)), ymin = c(quantile(pb1, 0.025),quantile(pb2, 0.025) )))

Da previsão ao efeito causal

Para tornar as coisas mais interessantes, estudamos um exemplo médico onde gostaríamos de obter um efeito causal (é inteiramente feito com números completamente irrealistas, embora seja inspirado em problemas reais - a diferença de reação que homens e mulheres podem ter em medicamentos) .

Neste exemplo, simulamos um resultado, digamos, um efeito de afinamento do sangue (B) que deve ser regulado por algum medicamento. Conhecemos também a idade e o sexo dos pacientes e simulamos a seguinte relação: Para um paciente do sexo masculino, independente da idade, a medicação aumenta linearmente o efeito anticoagulante. Para pacientes do sexo feminino, a medicação também aumenta B, mas em maior proporção do que para os homens, se tiverem mais de 50 anos de idade. Se estiverem abaixo de 50, no entanto, o efeito se reverte completamente e a medicação leva a um menor efeito de afinamento do sangue. O processo exato de geração de dados está aqui:

n<-1000
# We randomly sample from different ages…
Age<-sample(20:70, n, replace=T)
# and sexes
Sex <- rbinom(n, size=1,prob=0.6) # 1=woman, 0=man
# W denotes the dosis of the medication, the causal effect we are directly interested in
W<-sample( c(5,10,50,100),n, replace=T)
# B is the thing we want to understand
B<- 60+ (0.5*W)*(Sex==0) + (-0.5*W)*(Sex==1)*(Age<50) + (0.8*W)*(Sex==1)*(Age>=50) + rnorm(n)

Uma maneira de abordar isso com DRF é (1) tomar Y =(B,W) (então, novamente, nosso Y é bidimensional) e X =(Idade,Sexo), (2) obter os pesos para um determinado x e então (3) estime uma regressão linear ponderada com esses pesos. O que isso dá é um estimador do efeito, dado que X é fixo em x :

get_CATE_DRF = function(fit, newdata){
 out = predict(fit, newdata)
 ret = c()
 for(i in 1:nrow(newdata)){
 ret = c(ret, lm(out$y[,1]~out$y[,2], weights=out$weights[i,])$coefficients[2])
 }
 return(ret)
}

library(drf)
# Construct the data matrices
X<-matrix(cbind(Age, Sex), ncol=2)
Y<-matrix(cbind(B,W), ncol=2)
# Fit
DRF<-drf(Y=Y, X=X)
# Get a test point (changing this test point gives different queries)
x<-matrix(c(40, 1), nrow=1)
# Get the weights
weights <- predict(DRF, newdata=x)$weights[1,]
# Study the weights
head(X[order(weights, decreasing=T),])
# Result for given test point
(CATE_DRF = get_CATE_DRF(DRF, newdata=x))

De certa forma, este ainda é um exemplo muito simples e até mesmo uma regressão linear com idade e sexo como preditores pode funcionar bem. O importante é que DRF não tem suposições anteriores aqui (como linearidade) e aprende as relações por si só, mesmo quando o efeito de X é não linear. Fica muito mais difícil estimar esse efeito para tamanhos de amostra menores, mas a direção geral geralmente não é tão ruim.

Conclusão

Este artigo explicou o método Distributional Random Forest (espero que de uma forma compreensível). O método é uma Floresta Aleatória, onde cada árvore divide a resposta Y de acordo com X de forma que observações com distribuições semelhantes terminem em um nó folha. Se então um novo ponto x cair em uma árvore, ele terminará em um nó folha com outros x i's que tenham uma distribuição condicional semelhante de Y . Isso resulta em pesos que, calculados sobre todas as árvores, fornecem uma estimativa da distribuição condicional de forma simples. Isso dá uma estimativa puramente não paramétrica de P( Y | X = x) a partir do qual muitas quantidades interessantes podem ser estimadas.

No final deste artigo, queremos apenas alertar que estimar uma distribuição condicional multivariada não parametricamente é uma tarefa assustadora. Faz sentido principalmente quando há muitas observações e quando se suspeita de um relacionamento complicado. No entanto, às vezes, apenas assumir um modelo linear com uma distribuição Gaussiana também funciona. O que torna o DRF tão versátil é que mesmo nos casos em que um modelo paramétrico é mais apropriado, os pesos ainda podem ser úteis para uma abordagem semi-paramétrica,

Um grande número de exemplos adicionais pode ser encontrado no artigo original ou em futuros artigos potenciais do Medium. Esperamos que o DRF ajude em muitas tarefas dependentes de dados!