Cómo: Intervalos de predicción para regresión lineal a través de bootstrapping

14

Tengo problemas para entender cómo usar bootstrapping para calcular los intervalos de predicción para un modelo de regresión lineal. ¿Alguien puede describir un procedimiento paso a paso? Busqué a través de google pero nada realmente tiene sentido para mí.

Entiendo cómo usar bootstrapping para calcular intervalos de confianza para los parámetros del modelo.

Max
fuente
1
Esto se discute en detalle en el libro de Davison y Hinkley, Bootstrap Methods and These Application , junto con un algoritmo explícito (Algoritmo 6.4). Explican conceptos, dificultades y detalles con mayor extensión de lo que es posible en una respuesta razonable aquí.
Glen_b -Reinstate Monica
@Glen_b Gracias por la referencia. Lamentablemente, no estoy en una universidad o empresa, así que no tengo los recursos para adquirir el libro.
Max
Se puede pedir a Amazon; Una explicación completa del algoritmo y todas las advertencias y problemas asociados no es realmente el tipo de cosa que puede cubrir en unos cientos de palabras o incluso una respuesta de una página.
Glen_b -Reinstale a Monica el
1
@Glen_b Escribí el algoritmo de Davison y HInkely --- hay varias preguntas en CV sobre esto, así que pensé que valía la pena el esfuerzo. Cualquier comentario que tenga será apreciado. stats.stackexchange.com/questions/226565/…
Bill
Este hilo parece responder a su pregunta: stats.stackexchange.com/questions/226565/…
user2683832

Respuestas:

7

Los intervalos de confianza tienen en cuenta la incertidumbre de la estimación. Los intervalos de predicción añaden a esto la incertidumbre fundamental. R predict.lmle dará el intervalo de predicción para un modelo lineal. A partir de ahí, todo lo que tiene que hacer es ejecutarlo repetidamente en muestras de arranque.

n <- 100
n.bs <- 30

dat <- data.frame( x<-runif(n), y=x+runif(n) )
plot(y~x,data=dat)


regressAndPredict <- function( dat ) {
  model <- lm( y~x, data=dat )
  predict( model, interval="prediction" )
}

regressAndPredict(dat)

replicate( n.bs, regressAndPredict(dat[ sample(seq(n),replace=TRUE) ,]) )

El resultado de replicatees una matriz tridimensional ( nx 3x n.bs). La dimensión de longitud 3 consiste en el valor ajustado para cada elemento de datos y los límites inferior / superior del intervalo de predicción del 95%.

Método de Gary King

Dependiendo de lo que quieras, hay un método genial de King, Tomz y Wittenberg . Es relativamente fácil de implementar y evita los problemas de arranque para ciertas estimaciones (por ejemplo max(Y)).

Citaré su definición de incertidumbre fundamental aquí, ya que es razonablemente agradable:

Una segunda forma de variabilidad, la incertidumbre fundamental representada por el componente estocástico (la distribución f) en la Ecuación 1, resulta de innumerables eventos fortuitos como el clima o la enfermedad que pueden influir en Y pero no están incluidos en X. Incluso si conocíamos los valores exactos de los parámetros (eliminando así la incertidumbre de la estimación), la incertidumbre fundamental nos impediría predecir Y sin error.

Ari B. Friedman
fuente
3
No estoy seguro de cómo construir un intervalo de confianza a partir de esta matriz de intervalos de predicción n.bs.
B_Miner
1

Bootstrapping no asume ningún conocimiento de la forma de la distribución principal subyacente de la cual surgió la muestra. Las estimaciones de parámetros estadísticos clásicos tradicionales se basan en el supuesto de normalidad. Bootstrap trata la no normalidad y es más preciso en la práctica que los métodos clásicos.

Bootstrapping sustituye la potencia informática sin procesar de las computadoras por un análisis teórico riguroso. Es una estimación de la distribución de muestreo de un término de error del conjunto de datos. Bootstrapping incluye: volver a muestrear el conjunto de datos un número específico de veces, calcular la media de cada muestra y encontrar el error estándar de la media.

El siguiente código "R" demuestra el concepto:

Este ejemplo práctico demuestra la utilidad de bootstrapping y estima el error estándar. Se requiere el error estándar para calcular el intervalo de confianza.

Supongamos que tiene un conjunto de datos sesgado "a":

a<-rexp(395, rate=0.1)          # Create skewed data

visualización del conjunto de datos sesgados

plot(a,type="l")                # Scatter plot of the skewed data
boxplot(a,type="l")             # Box plot of the skewed data
hist(a)                         # Histogram plot of the skewed data

Realice el procedimiento de arranque:

n <- length(a)                  # the number of bootstrap samples should equal the original data set
    xbarstar <- c()                 # Declare the empty set “xbarstar” variable which will be holding the mean of every bootstrap iteration
    for (i in 1:1000) {             # Perform 1000 bootstrap iteration
        boot.samp <- sample(a, n, replace=TRUE) #”Sample” generates the same number of elements as the original data set
    xbarstar[i] <- mean(boot.samp)} # “xbarstar” variable  collects 1000 averages of the original data set
    ## 
    plot(xbarstar)                  # Scatter plot of the bootstrapped data
    boxplot(xbarstar)               # Box plot of the bootstrapped data
    hist(xbarstar)                  # Histogram plot of the bootstrapped data

    meanOfMeans <- mean(xbarstar)
    standardError <- sd(xbarstar)    # the standard error is the standard deviation of the mean of means
    confidenceIntervalAboveTheMean <- meanOfMeans + 1.96 * standardError # for 2 standard deviation above the mean 
    confidenceIntervalBelowTheMean <- meanOfMeans - 1.96 * standardError # for 2 standard deviation above the mean 
    confidenceInterval <- confidenceIntervalAboveTheMean + confidenceIntervalBelowTheMean
    confidenceInterval
Isaac Ragy
fuente
1
Gracias Ragy por el ejemplo. Sin embargo, por lo que puedo ver, su respuesta no cubrió el cálculo de los intervalos de predicción utilizando bootstrap. Como dije en mi respuesta, ya entiendo cómo usar bootstrapping para calcular los intervalos de confianza, lo que parece hacer su código.
Max