R: O que vejo nos gráficos de dependência parcial de gbm e RandomForest?

Na verdade, eu pensei que tinha entendido o que se pode mostrar com um gráfico de dependência parcial, mas usando um exemplo hipotético muito simples, fiquei bastante confuso. No seguinte pedaço de código que geram três variáveis independentes ( a , b , c ) e uma variável dependente ( y ) com c que mostra uma relação linear próxima com y , ao passo que um e b não estão correlacionados com y . Eu faço uma análise de regressão com uma árvore de regressão impulsionada usando o pacote R gbm:

a <- runif(100, 1, 100)
b <- runif(100, 1, 100)
c <- 1:100 + rnorm(100, mean = 0, sd = 5)
y <- 1:100 + rnorm(100, mean = 0, sd = 5)
par(mfrow = c(2,2))
plot(y ~ a); plot(y ~ b); plot(y ~ c)
Data <- data.frame(matrix(c(y, a, b, c), ncol = 4))
names(Data) <- c("y", "a", "b", "c")
library(gbm)
gbm.gaus <- gbm(y ~ a + b + c, data = Data, distribution = "gaussian")
par(mfrow = c(2,2))
plot(gbm.gaus, i.var = 1)
plot(gbm.gaus, i.var = 2)
plot(gbm.gaus, i.var = 3)

Não surpreendentemente, para as variáveis de um e b as parcelas dependência parciais produzir linhas horizontais em torno da média de um . O que eu acho é o enredo para a variável c . Eu recebo linhas horizontais para os intervalos c <40 ec > 60 e o eixo y é restrito a valores próximos à média de y . Uma vez que um e b são completamente alheios ao y (e, portanto, não importância variável em que o modelo é 0), I esperado que cmostraria dependência parcial ao longo de toda a faixa, em vez da forma sigmóide, para uma faixa muito restrita de seus valores. Tentei encontrar informações em Friedman (2001) "Aproximação da função gananciosa: uma máquina de aumento de gradiente" e em Hastie et al. (2011) "Elements of Statistical Learning", mas minhas habilidades matemáticas são muito baixas para entender todas as equações e fórmulas contidas nela. Assim, minha pergunta: o que determina a forma do gráfico de dependência parcial para a variável c ? (Por favor, explique em palavras compreensíveis para um não matemático!)

ADICIONADO em 17 de abril de 2014:

Enquanto esperava por uma resposta, usei o mesmo exemplo de dados para uma análise com o R-package randomForest. As parcelas dependência parciais de Floresta aleatória assemelham-se muito mais para o que o esperado a partir das parcelas gbm: a dependência parcial de variáveis explanatórias a e b variam de forma aleatória e estreitamente em torno de 50, enquanto explicativas variável c dependência mostra parcial ao longo de toda a sua gama (e ao longo de quase todo o intervalo de y ). Quais poderiam ser as razões para essas diferentes formas das parcelas de dependência parcial em gbme randomForest?

Aqui, o código modificado que compara os gráficos:

a <- runif(100, 1, 100)
b <- runif(100, 1, 100)
c <- 1:100 + rnorm(100, mean = 0, sd = 5)
y <- 1:100 + rnorm(100, mean = 0, sd = 5)
par(mfrow = c(2,2))
plot(y ~ a); plot(y ~ b); plot(y ~ c)
Data <- data.frame(matrix(c(y, a, b, c), ncol = 4))
names(Data) <- c("y", "a", "b", "c")

library(gbm)
gbm.gaus <- gbm(y ~ a + b + c, data = Data, distribution = "gaussian")

library(randomForest)
rf.model <- randomForest(y ~ a + b + c, data = Data)

x11(height = 8, width = 5)
par(mfrow = c(3,2))
par(oma = c(1,1,4,1))
plot(gbm.gaus, i.var = 1)
partialPlot(rf.model, Data[,2:4], x.var = "a")
plot(gbm.gaus, i.var = 2)
partialPlot(rf.model, Data[,2:4], x.var = "b")
plot(gbm.gaus, i.var = 3)
partialPlot(rf.model, Data[,2:4], x.var = "c")
title(main = "Boosted regression tree", outer = TRUE, adj = 0.15)
title(main = "Random forest", outer = TRUE, adj = 0.85)

Passei algum tempo escrevendo meu próprio "plotador de função parcial." Antes de perceber que ele já estava incluído na biblioteca R randomForest.

[EDITAR ... mas depois passei um ano fabricando o pacote CRAN forestFloor , o que, na minha opinião, é significativamente melhor do que as parcelas clássicas de dependência parcial]

O gráfico Partial.function é excelente em exemplos como este exemplo de simulação que você mostra aqui, onde a variável explicativa não interage com outras variáveis. Se cada variável explicativa contribuir de forma aditiva para o alvo Y por alguma função desconhecida, esse método é ótimo para mostrar essa função oculta estimada. Muitas vezes vejo esse achatamento nas fronteiras das funções parciais.

Algumas razões: randomForsest possui um argumento chamado 'nodesize = 5', o que significa que nenhuma árvore subdividirá um grupo de 5 membros ou menos. Portanto, cada árvore não pode distinguir com mais precisão. A camada de ensacamento / bootstrapping de ensemple suaviza votando as muitas funções de etapas das árvores individuais - mas apenas no meio da região de dados. Aproximando as bordas dos dados representados no espaço, a 'amplitude' da função parcial. Definir o tamanho do nó = 3 e / ou obter mais observações em comparação ao ruído pode reduzir esse efeito de nivelamento da borda ... Quando a relação sinal / ruído cai em geral em florestas aleatórias, a escala de previsões condensa. Assim, as previsões não são absolutamente termos precisos, mas apenas linearmente correlacionados com o alvo. Você pode ver os valores aeb como exemplos de uma relação sinal / ruído extremamente baixa, e, portanto, essas funções parciais são muito planas. É um recurso interessante da floresta aleatória que você, do conjunto de previsões do conjunto de treinamento, pode adivinhar o desempenho do modelo. OOB.predictions também é ótimo ..

o achatamento de parcelas parciais em regiões sem dados é razoável: como a floresta aleatória e o CART são modelagem orientada a dados, eu pessoalmente gosto do conceito de que esses modelos não extrapolam. Assim, a previsão de c = 500 ou c = 1100 é exatamente a mesma que c = 100 ou, na maioria dos casos, também c = 98.

Aqui está um exemplo de código com o nivelamento da borda reduzido:

Eu não tentei o pacote gbm ...

aqui está um código ilustrativo baseado em seu exemplo ...

#more observations are created...
a <- runif(5000, 1, 100)
b <- runif(5000, 1, 100)
c <- (1:5000)/50 + rnorm(100, mean = 0, sd = 0.1)
y <- (1:5000)/50 + rnorm(100, mean = 0, sd = 0.1)
par(mfrow = c(1,3))
plot(y ~ a); plot(y ~ b); plot(y ~ c)
Data <- data.frame(matrix(c(y, a, b, c), ncol = 4))
names(Data) <- c("y", "a", "b", "c")
library(randomForest)
#smaller nodesize "not as important" when there number of observartion is increased
#more tress can smooth flattening so boundery regions have best possible signal to             noise, data specific how many needed

plot.partial = function() {
partialPlot(rf.model, Data[,2:4], x.var = "a",xlim=c(1,100),ylim=c(1,100))
partialPlot(rf.model, Data[,2:4], x.var = "b",xlim=c(1,100),ylim=c(1,100))
partialPlot(rf.model, Data[,2:4], x.var = "c",xlim=c(1,100),ylim=c(1,100))
}

#worst case! : with 100 samples from Data and nodesize=30
rf.model <- randomForest(y ~ a + b + c, data = Data[sample(5000,100),],nodesize=30)
plot.partial()

#reasonble settings for least partial flattening by few observations: 100 samples and nodesize=3 and ntrees=2000
#more tress can smooth flattening so boundery regions have best possiblefidelity
rf.model <- randomForest(y ~ a + b + c, data = Data[sample(5000,100),],nodesize=5,ntress=2000)
plot.partial()

#more observations is great!
rf.model <- randomForest(y ~ a + b + c,
 data = Data[sample(5000,5000),],
 nodesize=5,ntress=2000)
plot.partial()

Como mencionado nos comentários acima, o modelo gbm seria melhor com algum ajuste de parâmetro. Uma maneira fácil de detectar problemas no modelo e a necessidade de tais parâmetros é gerar alguns gráficos de diagnóstico. Por exemplo, para o modelo gbm acima com os parâmetros padrão (e usando o pacote plotmo para criar os gráficos), temos

gbm.gaus <- gbm(y~., data = Data, dist = "gaussian")
library(plotmo)   # for the plotres function
plotres(gbm.gaus) # plot the error per ntrees and the residuals

que dá

No gráfico à esquerda, vemos que a curva de erro não chegou ao fundo do poço. E no gráfico à direita, os resíduos não são o que gostaríamos.

Se reconstruirmos o modelo com um número maior de árvores

gbm.gaus1 <- gbm(y~., data = Data, dist = "gaussian",
                 n.trees=5000, interact=3)
plotres(gbm.gaus1)

Nós temos

Vemos a curva de erro com um grande número de árvores e o gráfico de resíduos é mais saudável. Também podemos traçar os gráficos de dependência parcial para o novo modelo gbm e o modelo de floresta aleatória

library(plotmo)
plotmo(gbm.gaus1, pmethod="partdep", all1=TRUE, all2=TRUE)
plotmo(rf.model,  pmethod="partdep", all1=TRUE, all2=TRUE)

que dá

Os gráficos de modelo de floresta aleatória e gbm agora são semelhantes, como esperado.

Você precisa atualizar seu interaction.depthparâmetro ao criar seu modelo otimizado. O padrão é 1 e isso fará com que todas as árvores que o gbmalgoritmo construa se dividam apenas uma vez cada. Isso significa que todas as árvores estão apenas dividindo as variáveis ce, dependendo da amostra de observações que ela usa, elas se dividem em torno de 40 - 60.

Aqui estão as parcelas parciais com interaction.depth = 3