Si ejecuto un randomForest
modelo, puedo entonces hacer predicciones basadas en el modelo. ¿Existe una manera de obtener un intervalo de predicción de cada una de las predicciones de manera que sepa lo "seguro" que está el modelo de su respuesta? Si esto es posible, ¿se basa simplemente en la variabilidad de la variable dependiente para todo el modelo o tendrá intervalos más amplios y más estrechos dependiendo del árbol de decisión particular que se siguió para una predicción particular?
Respuestas
¿Demasiados anuncios?Esto es en parte una respuesta a @Sashikanth Dareddy (ya que no cabe en un comentario) y en parte una respuesta al post original.
Recuerda lo que es un intervalo de predicción, es un intervalo o conjunto de valores en los que predecimos que se encontrarán las observaciones futuras. Por lo general, el intervalo de predicción tiene dos partes principales que determinan su anchura: una parte que representa la incertidumbre sobre la media predicha (u otro parámetro), que es la parte del intervalo de confianza, y una parte que representa la variabilidad de las observaciones individuales en torno a esa media. El intervalo de confianza es bastante robusto debido al Teorema Central del Límite y, en el caso de un bosque aleatorio, el bootstrapping también ayuda. Pero el intervalo de predicción depende completamente de las suposiciones sobre cómo se distribuyen los datos dadas las variables predictoras, el CLT y el bootstrapping no tienen efecto en esa parte.
El intervalo de predicción debe ser más amplio cuando el correspondiente intervalo de confianza también lo sea. Otras cosas que afectarían a la anchura del intervalo de predicción son las suposiciones sobre la igualdad de varianza o no, esto tiene que venir del conocimiento del investigador, no del modelo de bosque aleatorio.
Un intervalo de predicción no tiene sentido para un resultado categórico (se podría hacer un conjunto de predicciones en lugar de un intervalo, pero la mayoría de las veces probablemente no sería muy informativo).
Podemos ver algunos de los problemas en torno a los intervalos de predicción simulando datos en los que conocemos la verdad exacta. Consideremos los siguientes datos:
set.seed(1)
x1 <- rep(0:1, each=500)
x2 <- rep(0:1, each=250, length=1000)
y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000)
Estos datos en particular siguen los supuestos de una regresión lineal y son bastante sencillos para un ajuste de bosque aleatorio. Sabemos por el modelo "verdadero" que cuando ambos predictores son 0 la media es 10, también sabemos que los puntos individuales siguen una distribución normal con una desviación estándar de 1. Esto significa que el intervalo de predicción del 95% basado en un conocimiento perfecto para estos puntos sería de 8 a 12 (bueno, en realidad de 8,04 a 11,96, pero el redondeo lo hace más sencillo). Cualquier intervalo de predicción estimado debería ser más amplio que esto (no tener información perfecta añade amplitud para compensar) e incluir este rango.
Veamos los intervalos de la regresión:
fit1 <- lm(y ~ x1 * x2)
newdat <- expand.grid(x1=0:1, x2=0:1)
(pred.lm.ci <- predict(fit1, newdat, interval='confidence'))
# fit lwr upr
# 1 10.02217 9.893664 10.15067
# 2 14.90927 14.780765 15.03778
# 3 20.02312 19.894613 20.15162
# 4 21.99885 21.870343 22.12735
(pred.lm.pi <- predict(fit1, newdat, interval='prediction'))
# fit lwr upr
# 1 10.02217 7.98626 12.05808
# 2 14.90927 12.87336 16.94518
# 3 20.02312 17.98721 22.05903
# 4 21.99885 19.96294 24.03476
Podemos ver que hay cierta incertidumbre en las medias estimadas (intervalo de confianza) y eso nos da un intervalo de predicción que es más amplio (pero incluye) el rango de 8 a 12.
Ahora veamos el intervalo basado en las predicciones individuales de los árboles individuales (deberíamos esperar que fueran más amplias, ya que el bosque aleatorio no se beneficia de las suposiciones (que sabemos que son ciertas para estos datos) que hace la regresión lineal):
library(randomForest)
fit2 <- randomForest(y ~ x1 + x2, ntree=1001)
pred.rf <- predict(fit2, newdat, predict.all=TRUE)
pred.rf.int <- apply(pred.rf$individual, 1, function(x) {
c(mean(x) + c(-1, 1) * sd(x),
quantile(x, c(0.025, 0.975)))
})
t(pred.rf.int)
# 2.5% 97.5%
# 1 9.785533 13.88629 9.920507 15.28662
# 2 13.017484 17.22297 12.330821 18.65796
# 3 16.764298 21.40525 14.749296 21.09071
# 4 19.494116 22.33632 18.245580 22.09904
Los intervalos son más amplios que los intervalos de predicción de la regresión, pero no cubren todo el rango. Incluyen los valores verdaderos y, por tanto, pueden ser legítimos como intervalos de confianza, pero sólo predicen dónde está la media (valor predicho), sin la pieza añadida para la distribución alrededor de esa media. En el primer caso, en el que x1 y x2 son ambos 0, los intervalos no bajan de 9,7, lo que es muy diferente del verdadero intervalo de predicción, que baja hasta 8. Si generamos nuevos puntos de datos, habrá varios puntos (mucho más del 5%) que están en los intervalos verdadero y de regresión, pero no caen en los intervalos del bosque aleatorio.
Para generar un intervalo de predicción tendrá que hacer algunas suposiciones fuertes sobre la distribución de los puntos individuales alrededor de las medias predichas, entonces podría tomar las predicciones de los árboles individuales (la pieza del intervalo de confianza bootstrapped) y luego generar un valor aleatorio de la distribución asumida con ese centro. Los cuantiles de esos trozos generados pueden formar el intervalo de predicción (pero aún así lo probaría, puede que tenga que repetir el proceso varias veces más y combinarlo).
He aquí un ejemplo de cómo hacerlo añadiendo desviaciones normales (ya que sabemos que los datos originales utilizaban una normal) a las predicciones con la desviación estándar basada en el MSE estimado de ese árbol:
pred.rf.int2 <- sapply(1:4, function(i) {
tmp <- pred.rf$individual[i, ] + rnorm(1001, 0, sqrt(fit2$mse))
quantile(tmp, c(0.025, 0.975))
})
t(pred.rf.int2)
# 2.5% 97.5%
# [1,] 7.351609 17.31065
# [2,] 10.386273 20.23700
# [3,] 13.004428 23.55154
# [4,] 16.344504 24.35970
Estos intervalos contienen los basados en el conocimiento perfecto, por lo que parecen razonables. Pero, dependerán en gran medida de las suposiciones que se hagan (las suposiciones son válidas aquí porque utilizamos el conocimiento de cómo se simularon los datos, pueden no ser tan válidas en casos de datos reales). Aún así, yo repetiría las simulaciones varias veces para datos que se parezcan más a los reales (pero simulados para saber la verdad) varias veces antes de confiar plenamente en este método.
Si usas R puedes producir fácilmente intervalos de predicción para las predicciones de una regresión de bosques aleatorios: Sólo tiene que utilizar el paquete quantregForest
(disponible en CRAN ) y leer el papel de N. Meinshausen sobre cómo se pueden inferir los cuantiles condicionales con los bosques de regresión cuantílica y cómo se pueden utilizar para construir intervalos de predicción. ¡Muy informativo incluso si no trabajas con R!
Me doy cuenta de que este es un post antiguo, pero he estado haciendo algunas simulaciones sobre esto y pensé en compartir mis hallazgos.
@GregSnow hizo un post muy detallado sobre esto, pero creo que al calcular el intervalo utilizando las predicciones de los árboles individuales estaba mirando $[ \mu + \sigma, \mu - \sigma]$ que es sólo un intervalo de predicción del 70%. Tenemos que ver $[\mu+1.96*\sigma,\mu-1.96*\sigma]$ para obtener el intervalo de predicción del 95%.
Haciendo este cambio en el código de @GregSnow, obtenemos los siguientes resultados
set.seed(1)
x1 <- rep( 0:1, each=500 )
x2 <- rep( 0:1, each=250, length=1000 )
y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000)
library(randomForest)
fit2 <- randomForest(y~x1+x2)
pred.rf <- predict(fit2, newdat, predict.all=TRUE)
pred.rf.int <- t(apply( pred.rf$individual, 1, function(x){
c( mean(x) + c(-1.96,1.96)*sd(x), quantile(x, c(0.025,0.975)) )}))
pred.rf.int
2.5% 97.5%
1 7.826896 16.05521 9.915482 15.31431
2 11.010662 19.35793 12.298995 18.64296
3 14.296697 23.61657 14.749248 21.11239
4 18.000229 23.73539 18.237448 22.10331
Ahora, comparando estos con los intervalos generados al añadir la desviación normal a las predicciones con la desviación estándar como MSE como sugirió @GregSnow obtenemos,
pred.rf.int2 <- sapply(1:4, function(i) {
tmp <- pred.rf$individual[i,] + rnorm(1000, 0, sqrt(fit2$mse))
quantile(tmp, c(0.025, 0.975))
})
t(pred.rf.int2)
2.5% 97.5%
[1,] 7.486895 17.21144
[2,] 10.551811 20.50633
[3,] 12.959318 23.46027
[4,] 16.444967 24.57601
El intervalo de estos dos enfoques está ahora muy cerca. El gráfico del intervalo de predicción de los tres enfoques frente a la distribución de errores en este caso es el siguiente
- Líneas negras = intervalos de predicción de la regresión lineal,
- Líneas rojas = intervalos del bosque aleatorio calculados sobre las predicciones individuales,
- Líneas azules = Intervalos de bosque aleatorio calculados añadiendo la desviación normal a las predicciones
Ahora, volvamos a realizar la simulación pero esta vez aumentando la varianza del término de error. Si nuestros cálculos del intervalo de predicción son buenos, deberíamos terminar con intervalos más amplios que los que obtuvimos anteriormente.
set.seed(1)
x1 <- rep( 0:1, each=500 )
x2 <- rep( 0:1, each=250, length=1000 )
y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000,mean=0,sd=5)
fit1 <- lm(y~x1+x2)
newdat <- expand.grid(x1=0:1,x2=0:1)
predict(fit1,newdata=newdat,interval = "prediction")
fit lwr upr
1 10.75006 0.503170 20.99695
2 13.90714 3.660248 24.15403
3 19.47638 9.229490 29.72327
4 22.63346 12.386568 32.88035
set.seed(1)
fit2 <- randomForest(y~x1+x2,localImp=T)
pred.rf.int <- t(apply( pred.rf$individual, 1, function(x){
c( mean(x) + c(-1.96,1.96)*sd(x), quantile(x, c(0.025,0.975)) )}))
pred.rf.int
2.5% 97.5%
1 7.889934 15.53642 9.564565 15.47893
2 10.616744 18.78837 11.965325 18.51922
3 15.024598 23.67563 14.724964 21.43195
4 17.967246 23.88760 17.858866 22.54337
pred.rf.int2 <- sapply(1:4, function(i) {
tmp <- pred.rf$individual[i,] + rnorm(1000, 0, sqrt(fit2$mse))
quantile(tmp, c(0.025, 0.975))
})
t(pred.rf.int2)
2.5% 97.5%
[1,] 1.291450 22.89231
[2,] 4.193414 25.93963
[3,] 7.428309 30.07291
[4,] 9.938158 31.63777
Ahora bien, esto deja claro que el cálculo de los intervalos de predicción mediante el segundo enfoque es mucho más exacto y está dando resultados bastante cercanos al intervalo de predicción de la regresión lineal.
Tomando el supuesto de normalidad, hay otra forma más sencilla de calcular los intervalos de predicción del bosque aleatorio. De cada uno de los árboles individuales tenemos el valor predicho ( $\mu_i$ ) así como el error medio cuadrático ( $MSE_i$ ). Así, la predicción de cada árbol individual puede considerarse como $ N(\mu_i,RMSE_i)$ . Utilizando las propiedades de la distribución normal, nuestra predicción del bosque aleatorio tendría la distribución $N(\sum \mu_i/n, \sum RMSE_i/n)$ . Aplicando esto al ejemplo que comentamos anteriormente, obtenemos los siguientes resultados
mean.rf <- pred.rf$aggregate
sd.rf <- mean(sqrt(fit2$mse))
pred.rf.int3 <- cbind(mean.rf - 1.96*sd.rf, mean.rf + 1.96*sd.rf)
pred.rf.int3
1 1.332711 22.09364
2 4.322090 25.08302
3 8.969650 29.73058
4 10.546957 31.30789
Esto coincide muy bien con los intervalos del modelo lineal y también con el enfoque sugerido por @GregSnow. Pero tenga en cuenta que el supuesto subyacente en todos los métodos que hemos discutido es que los errores siguen una distribución normal.
Esto es fácil de resolver con randomForest.
En primer lugar, permítame abordar la tarea de regresión (suponiendo que su bosque tiene 1000 árboles). En el predict
tiene la opción de devolver resultados de árboles individuales. Esto significa que recibirá una salida de 1000 columnas. Podemos tomar la media de las 1000 columnas para cada fila - esta es la salida regular que RF habría producido de cualquier manera. Ahora, para obtener el intervalo de predicción, digamos +/- 2 desviaciones estándar, todo lo que tiene que hacer es, para cada fila, a partir de los 1000 valores, calcular +/-2 desviaciones estándar y convertirlas en los límites superior e inferior de su predicción.
En segundo lugar, en el caso de la clasificación, recuerde que cada árbol da como resultado 1 o 0 (por defecto) y la suma de todos los 1000 árboles dividida por 1000 da la probabilidad de la clase (en el caso de la clasificación binaria). Para poner un intervalo de predicción en la probabilidad necesitas modificar la opción min. nodesize (ver la documentación de randomForest para el nombre exacto de esa opción) una vez que establezcas un valor >>1 entonces los árboles individuales darán salida a números entre 1 y 0. Ahora, a partir de aquí puedes repetir el mismo proceso descrito anteriormente para la tarea de regresión.
Espero que tenga sentido.
El problema de la construcción de intervalos de predicción para las predicciones de bosques aleatorios se ha abordado en el siguiente documento:
Zhang, Haozhe, Joshua Zimmerman, Dan Nettleton y Daniel J. Nordman. "Intervalos de predicción de bosques aleatorios". The American Statistician, 2019.
El paquete R "rfinterval" es su implementación disponible en CRAN.
Instalación
Para instalar el paquete R rfinterval :
#install.packages("devtools")
#devtools::install_github(repo="haozhestat/rfinterval")
install.packages("rfinterval")
library(rfinterval)
?rfinterval
Uso
Inicio rápido:
train_data <- sim_data(n = 1000, p = 10)
test_data <- sim_data(n = 1000, p = 10)
output <- rfinterval(y~., train_data = train_data, test_data = test_data,
method = c("oob", "split-conformal", "quantreg"),
symmetry = TRUE,alpha = 0.1)
### print the marginal coverage of OOB prediction interval
mean(output$oob_interval$lo < test_data$y & output$oob_interval$up > test_data$y)
### print the marginal coverage of Split-conformal prediction interval
mean(output$sc_interval$lo < test_data$y & output$sc_interval$up > test_data$y)
### print the marginal coverage of Quantile regression forest prediction interval
mean(output$quantreg_interval$lo < test_data$y & output$quantreg_interval$up > test_data$y)
Ejemplo de datos:
oob_interval <- rfinterval(pm2.5 ~ .,
train_data = BeijingPM25[1:1000, ],
test_data = BeijingPM25[1001:2000, ],
method = "oob",
symmetry = TRUE,
alpha = 0.1)
str(oob_interval)
- Ver respuestas anteriores
- Ver más respuestas