R: ¿Qué veo en las parcelas de dependencia parcial de gbm y RandomForest?

En realidad, pensé que había entendido lo que se puede mostrar con un diagrama de dependencia parcial, pero usando un ejemplo hipotético muy simple, me quedé bastante perplejo. En el siguiente fragmento de código, genero tres variables independientes ( a , b , c ) y una variable dependiente ( y ) con c que muestra una relación lineal estrecha con y , mientras que a y b son correlacionados con y . Hago un análisis de regresión con un árbol de regresión potenciado usando el paquete 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)

No es sorprendente que, para las variables a y b las parcelas de dependencia parciales dió líneas horizontales alrededor de la media de una . Lo que me intriga es la trama de la variable c . Obtengo líneas horizontales para los rangos c <40 yc > 60 y el eje y está restringido a valores cercanos a la media de y . Desde un y b son completamente no relacionada con y (e importancia por lo tanto no variable en el modelo es 0), que esperaba que cmostraría una dependencia parcial a lo largo de todo su rango en lugar de esa forma sigmoidea para un rango muy restringido de sus valores. Traté de encontrar información en Friedman (2001) "Aproximación de la función codiciosa: una máquina de aumento de gradiente" y en Hastie et al. (2011) "Elementos de aprendizaje estadístico", pero mis habilidades matemáticas son demasiado bajas para comprender todas las ecuaciones y fórmulas que contiene. Por lo tanto, mi pregunta: ¿qué determina la forma del gráfico de dependencia parcial para la variable c ? (¡Por favor explique en palabras comprensibles para un no matemático!)

AGREGADO el 17 de abril de 2014:

Mientras esperaba una respuesta, utilicé los mismos datos de ejemplo para un análisis con R-package randomForest. Las parcelas de dependencia parciales de randomForest se parecen mucho más a lo que esperaba de las parcelas de GBM: la dependencia parcial de variables explicativas a y b varían aleatoriamente y estrechamente alrededor de 50, mientras variable explicativa c dependencia shows parcial en toda su gama (y sobre casi rango completo de y ). ¿Cuáles podrían ser las razones de estas diferentes formas de las parcelas de dependencia parcial en gbmy randomForest?

Aquí el código modificado que compara las parcelas:

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)

Pasé algún tiempo escribiendo mi propio "plotter de funciones parciales" antes de darme cuenta de que ya estaba incluido en la biblioteca R randomForest.

[EDITAR ... pero luego pasé un año haciendo el paquete CRAN forestFloor , que en mi opinión es significativamente mejor que las parcelas clásicas de dependencia parcial]

El gráfico de función parcial es excelente en casos como este ejemplo de simulación que muestra aquí, donde la variable explicativa no interactúa con otras variables. Si cada variable explicativa contribuye de manera aditiva al objetivo Y mediante alguna función desconocida, este método es excelente para mostrar esa función oculta estimada. A menudo veo tal aplanamiento en los bordes de las funciones parciales.

Algunas razones: randomForsest tiene un argumento llamado 'nodeize = 5' que significa que ningún árbol subdividirá un grupo de 5 miembros o menos. Por lo tanto, cada árbol no puede distinguir con mayor precisión. La capa de ensacado / arranque de conjunto suaviza votando las muchas funciones de paso de los árboles individuales, pero solo en el medio de la región de datos. Al acercarse a los bordes de los datos representados espacio, la 'amplitud' de la función parcial caerá. Establecer nodosize = 3 y / u obtener más observaciones en comparación con el ruido puede reducir este efecto de aplanamiento de borde ... Cuando la relación señal / ruido cae en general en un bosque aleatorio, la escala de predicciones se condensa. Por lo tanto, las predicciones no son términos absolutamente precisos, sino que solo se correlacionan linealmente con el objetivo. Puede ver los valores ayb como ejemplos de una relación señal / ruido extremadamente baja, y por lo tanto estas funciones parciales son muy planas. Es una buena característica del bosque aleatorio que usted, desde el rango de predicciones del conjunto de entrenamiento, puede adivinar qué tan bien se está desempeñando el modelo. OOB.predictions es genial también ...

El aplanamiento de parcelas parciales en regiones sin datos es razonable: como el bosque aleatorio y el CART son modelos basados ​​en datos, personalmente me gusta el concepto de que estos modelos no se extrapolan. Por lo tanto, la predicción de c = 500 o c = 1100 es exactamente igual a c = 100 o, en la mayoría de los casos, también c = 98.

Aquí hay un ejemplo de código con el aplanamiento del borde reducido:

No he probado el paquete gbm ...

Aquí hay un código ilustrativo basado en su ejemplo ...

#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 se mencionó en los comentarios anteriores, el modelo gbm sería mejor con algunos ajustes de parámetros. Una manera fácil de detectar problemas en el modelo y la necesidad de tales parámetros es generar algunos gráficos de diagnóstico. Por ejemplo, para el modelo gbm anterior con los parámetros predeterminados (y usando el paquete plotmo para crear los gráficos) tenemos

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

lo que da

En la gráfica de la izquierda vemos que la curva de error no ha tocado fondo. Y en la gráfica de la derecha, los residuos no son lo que queremos.

Si reconstruimos el modelo con una mayor cantidad de árboles

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

obtenemos

Vemos la curva de error al fondo con una gran cantidad de árboles, y la gráfica de residuos es más saludable. También podemos trazar las parcelas de dependencia parcial para el nuevo modelo gbm y el modelo de bosque aleatorio

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

lo que da

Las parcelas de modelo de bosque aleatorio y gbm ahora son similares, como se esperaba.

Debe actualizar su interaction.depthparámetro cuando construya su modelo mejorado. El valor predeterminado es 1 y eso hará que todos los árboles que gbmconstruye el algoritmo se dividan solo una vez cada uno. Esto significaría que cada árbol se está dividiendo en variables cy, dependiendo de la muestra de observaciones que use, se dividirá en algún lugar alrededor de 40-60.

Aquí están las parcelas parciales con interaction.depth = 3