He estado buscando en el arranque del paquete en R y aunque he encontrado un buen número de cartillas sobre cómo usarlo, todavía tengo que encontrar algo que describe exactamente lo que sucede "detrás de las escenas". Por ejemplo, en este ejemplo, la guía muestra cómo utilizar el estándar de los coeficientes de regresión como un punto de partida para un bootstrap de regresión, pero no explica cuál es el procedimiento de arranque está haciendo derivar el bootstrap de los coeficientes de regresión. Parece que hay algún tipo de proceso iterativo que está pasando, pero me parece que no puede averiguar exactamente lo que está pasando.
Respuestas
¿Demasiados anuncios?
Hay varios "sabores" o formas de bootstrap no paramétrico paramétrico, remuestreo residual y muchos más). El bootstrap en el ejemplo se llama un no-paramétrico de bootstrap, o en caso de remuestreo (ver aquí, aquí, aquí y aquí para aplicaciones en la regresión). La idea básica es que el tratamiento de la muestra de la población, y en varias ocasiones dibujar nuevas muestras de que con la sustitución. Todas las observaciones originales tienen la misma probabilidad de ser atraído hacia el nuevo modelo. A continuación, calcular y almacenar la estadística(s) de interés, esto puede ser la media, la mediana o coeficientes de regresión utilizando el recién extraídas de la muestra. Esto se repite a lo $n$ veces. En cada iteración, algunas observaciones de su muestra original se dibujan varias veces mientras algunas observaciones pueden no ser dibujado. Después de $n$ iteraciones, ha $n$ almacenado bootstrap estimaciones de la estadística(s) de interés (por ejemplo, si $n=1000$ y el estadístico de interés es la media, tiene 1000 bootstrap estimaciones de la media). Por último, resumen de los datos estadísticos como la media, la mediana y la desviación estándar de la $n$ bootstrap-las estimaciones se calculan.
De arranque se utiliza a menudo para:
- Cálculo de intervalos de confianza (y la estimación de los errores estándar)
- La estimación del sesgo de las estimaciones puntuales
Hay varios métodos para calcular los intervalos de confianza basados en el proceso de arranque de muestras (este documento proporciona una explicación y orientación). Un método muy sencillo para el cálculo de un 95%-intervalo de confianza es simplemente el cálculo de la empíricos 2.5 y 97.5 th percentiles de las muestras bootstrap (este intervalo se llama bootstrap percentil intervalo; véase el código de abajo). El simple percentil intervalo método se utiliza raramente en la práctica, ya que hay mejores métodos, tales como la corrección de sesgo y acelerado de bootstrap (BCa). BCa intervalos de ajustar tanto el sesgo y la asimetría en la distribución bootstrap.
El sesgo es simplemente calcula como la diferencia entre la media de los $n$ almacenan muestras bootstrap y la estimación original(s).
Vamos a replicar el ejemplo de la página web, pero utilizando nuestro propio bucle de la incorporación de las ideas que he descrito arriba (dibujo repetidamente con reemplazo):
#-----------------------------------------------------------------------------
# Load packages
#-----------------------------------------------------------------------------
require(ggplot2)
require(pscl)
require(MASS)
require(boot)
#-----------------------------------------------------------------------------
# Load data
#-----------------------------------------------------------------------------
zinb <- read.csv("http://www.ats.ucla.edu/stat/data/fish.csv")
zinb <- within(zinb, {
nofish <- factor(nofish)
livebait <- factor(livebait)
camper <- factor(camper)
})
#-----------------------------------------------------------------------------
# Calculate zero-inflated regression
#-----------------------------------------------------------------------------
m1 <- zeroinfl(count ~ child + camper | persons, data = zinb,
dist = "negbin", EM = TRUE)
#-----------------------------------------------------------------------------
# Store the original regression coefficients
#-----------------------------------------------------------------------------
original.estimates <- as.vector(t(do.call(rbind, coef(summary(m1)))[, 1:2]))
#-----------------------------------------------------------------------------
# Set the number of replications
#-----------------------------------------------------------------------------
n.sim <- 2000
#-----------------------------------------------------------------------------
# Set up a matrix to store the results
#-----------------------------------------------------------------------------
store.matrix <- matrix(NA, nrow=n.sim, ncol=12)
#-----------------------------------------------------------------------------
# The loop
#-----------------------------------------------------------------------------
set.seed(123)
for(i in 1:n.sim) {
#-----------------------------------------------------------------------------
# Draw the observations WITH replacement
#-----------------------------------------------------------------------------
data.new <- zinb[sample(1:dim(zinb)[1], dim(zinb)[1], replace=TRUE),]
#-----------------------------------------------------------------------------
# Calculate the model with this "new" data
#-----------------------------------------------------------------------------
m <- zeroinfl(count ~ child + camper | persons,
data = data.new, dist = "negbin",
start = list(count = c(1.3711, -1.5152, 0.879),
zero = c(1.6028, -1.6663)))
#-----------------------------------------------------------------------------
# Store the results
#-----------------------------------------------------------------------------
store.matrix[i, ] <- as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
}
#-----------------------------------------------------------------------------
# Save the means, medians and SDs of the bootstrapped statistics
#-----------------------------------------------------------------------------
boot.means <- colMeans(store.matrix, na.rm=T)
boot.medians <- apply(store.matrix,2,median, na.rm=T)
boot.sds <- apply(store.matrix,2,sd, na.rm=T)
#-----------------------------------------------------------------------------
# The bootstrap bias is the difference between the mean bootstrap estimates
# and the original estimates
#-----------------------------------------------------------------------------
boot.bias <- colMeans(store.matrix, na.rm=T) - original.estimates
#-----------------------------------------------------------------------------
# Basic bootstrap CIs based on the empirical quantiles
#-----------------------------------------------------------------------------
conf.mat <- matrix(apply(store.matrix, 2 ,quantile, c(0.025, 0.975), na.rm=T),
ncol=2, byrow=TRUE)
colnames(conf.mat) <- c("95%-CI Lower", "95%-CI Upper")
Y aquí está nuestro resumen de la tabla:
#-----------------------------------------------------------------------------
# Set up summary data frame
#-----------------------------------------------------------------------------
summary.frame <- data.frame(mean=boot.means, median=boot.medians,
sd=boot.sds, bias=boot.bias, "CI_lower"=conf.mat[,1], "CI_upper"=conf.mat[,2])
summary.frame
mean median sd bias CI_lower CI_upper
1 1.2998 1.3013 0.39674 -0.0712912 0.51960 2.0605
2 0.2527 0.2486 0.03208 -0.0034461 0.19898 0.3229
3 -1.5662 -1.5572 0.26220 -0.0509239 -2.12900 -1.0920
4 0.2005 0.1986 0.01949 0.0049019 0.16744 0.2418
5 0.9544 0.9252 0.48915 0.0753405 0.03493 1.9025
6 0.2702 0.2688 0.02043 0.0009583 0.23272 0.3137
7 -0.8997 -0.9082 0.22174 0.0856793 -1.30664 -0.4380
8 0.1789 0.1781 0.01667 0.0029513 0.14494 0.2140
9 2.0683 1.7719 1.59102 0.4654898 0.44150 8.0471
10 4.0209 0.8270 13.23434 3.1845710 0.58114 57.6417
11 -2.0969 -1.6717 1.56311 -0.4306844 -8.43440 -1.1156
12 3.8660 0.6435 13.27525 3.1870642 0.33631 57.6062
Algunas explicaciones
- La diferencia entre la media de las estimaciones bootstrap y las estimaciones originales es lo que se llama "sesgo" en la salida de
boot - Lo que la salida de
bootllamadas "std. error" es la desviación estándar de las estimaciones bootstrap
Comparar con el resultado de boot:
#-----------------------------------------------------------------------------
# Compare with boot output and confidence intervals
#-----------------------------------------------------------------------------
set.seed(10)
res <- boot(zinb, f, R = 2000, parallel = "snow", ncpus = 4)
res
Bootstrap Statistics :
original bias std. error
t1* 1.3710504 -0.076735010 0.39842905
t2* 0.2561136 -0.003127401 0.03172301
t3* -1.5152609 -0.064110745 0.26554358
t4* 0.1955916 0.005819378 0.01933571
t5* 0.8790522 0.083866901 0.49476780
t6* 0.2692734 0.001475496 0.01957823
t7* -0.9853566 0.083186595 0.22384444
t8* 0.1759504 0.002507872 0.01648298
t9* 1.6031354 0.482973831 1.58603356
t10* 0.8365225 3.240981223 13.86307093
t11* -1.6665917 -0.453059768 1.55143344
t12* 0.6793077 3.247826469 13.90167954
perc.cis <- matrix(NA, nrow=dim(res$t)[2], ncol=2)
for( i in 1:dim(res$t)[2] ) {
perc.cis[i,] <- boot.ci(res, conf=0.95, type="perc", index=i)$percent[4:5]
}
colnames(perc.cis) <- c("95%-CI Lower", "95%-CI Upper")
perc.cis
95%-CI Lower 95%-CI Upper
[1,] 0.52240 2.1035
[2,] 0.19984 0.3220
[3,] -2.12820 -1.1012
[4,] 0.16754 0.2430
[5,] 0.04817 1.9084
[6,] 0.23401 0.3124
[7,] -1.29964 -0.4314
[8,] 0.14517 0.2149
[9,] 0.29993 8.0463
[10,] 0.57248 56.6710
[11,] -8.64798 -1.1088
[12,] 0.33048 56.6702
#-----------------------------------------------------------------------------
# Our summary table
#-----------------------------------------------------------------------------
summary.frame
mean median sd bias CI_lower CI_upper
1 1.2998 1.3013 0.39674 -0.0712912 0.51960 2.0605
2 0.2527 0.2486 0.03208 -0.0034461 0.19898 0.3229
3 -1.5662 -1.5572 0.26220 -0.0509239 -2.12900 -1.0920
4 0.2005 0.1986 0.01949 0.0049019 0.16744 0.2418
5 0.9544 0.9252 0.48915 0.0753405 0.03493 1.9025
6 0.2702 0.2688 0.02043 0.0009583 0.23272 0.3137
7 -0.8997 -0.9082 0.22174 0.0856793 -1.30664 -0.4380
8 0.1789 0.1781 0.01667 0.0029513 0.14494 0.2140
9 2.0683 1.7719 1.59102 0.4654898 0.44150 8.0471
10 4.0209 0.8270 13.23434 3.1845710 0.58114 57.6417
11 -2.0969 -1.6717 1.56311 -0.4306844 -8.43440 -1.1156
12 3.8660 0.6435 13.27525 3.1870642 0.33631 57.6062
Comparar el "sesgo" columnas " y la "enfermedad de transmisión sexual. error" con la "sd" en la columna de nuestra propia tabla de resumen. Nuestra 95%-los intervalos de confianza son muy similares a los intervalos de confianza calculados por boot.ci utilizando el método del percentil (no todos, sin embargo: mira el límite inferior del parámetro con índice 9).
aron
Puntos
174
Usted debe centrarse en la función que se pasa a boot como la "estadística" de los parámetros y observe cómo se construye.
f <- function(data, i) {
require(pscl)
m <- zeroinfl(count ~ child + camper | persons,
data = data[i, ], dist = "negbin",
start = list(count = c(1.3711, -1.5152, 0.879), zero = c(1.6028, -1.6663)))
as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
}
Los "datos" argumento va a recibir toda una trama de datos, pero el "yo" argumento va a recibir una muestra de la fila de los índices generados por el "arranque" de 1:NROW(datos). Como se puede ver en el código "i" se utiliza a continuación para crear un neo-muestra que se pasa a zeroinl y, a continuación, sólo partes seleccionadas de la que se devuelven los resultados.
Vamos a imaginar que el "yo" es {1,2,3,3,3,6,7,7,10}. "[" La función devolverá sólo aquellas filas con 3 copias de la fila 3 y 2 copias de la fila 7. Que sería la base para una sola zeroinl() de cálculo y, a continuación, los coeficientes serán devueltos a boot como resultado de que replicar el proceso. El número de repeticiones es controlado por la "R" de los parámetros.
Dado que sólo los coeficientes de regresión son devueltos desde statistic en este caso, el boot función devolverá estos acumulado coeficientes como el valor de "t". Más que las comparaciones pueden ser realizadas por otras de arranque-paquete de funciones.
