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

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.

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.

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