Uso de error estándar de distribución bootstrap

19

(ignore el código R si es necesario, ya que mi pregunta principal es independiente del idioma)

Si quiero ver la variabilidad de una estadística simple (ej .: media), sé que puedo hacerlo a través de una teoría como:

x = rnorm(50)

# Estimate standard error from theory
summary(lm(x~1))
# same as...
sd(x) / sqrt(length(x))

o con el bootstrap como:

library(boot)

# Estimate standard error from bootstrap
(x.bs = boot(x, function(x, inds) mean(x[inds]), 1000))
# which is simply the standard *deviation* of the bootstrap distribution...
sd(x.bs$t)

Sin embargo, lo que me pregunto es, ¿ puede ser útil / válido (?) Mirar el error estándar de una distribución de arranque en ciertas situaciones? La situación con la que estoy lidiando es una función no lineal relativamente ruidosa, como:

# Simulate dataset
set.seed(12345)
n   = 100
x   = runif(n, 0, 20)
y   = SSasymp(x, 5, 1, -1) + rnorm(n, sd=2)
dat = data.frame(x, y)

Aquí el modelo ni siquiera converge usando el conjunto de datos original,

> (fit = nls(y ~ SSasymp(x, Asym, R0, lrc), dat))
Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

por lo tanto, las estadísticas que me interesan son estimaciones más estabilizadas de estos parámetros nls, tal vez sus medias en varias réplicas de arranque.

# Obtain mean bootstrap nls parameter estimates
fit.bs = boot(dat, function(dat, inds)
              tryCatch(coef(nls(y ~ SSasymp(x, Asym, R0, lrc), dat[inds, ])),
                       error=function(e) c(NA, NA, NA)), 100)
pars = colMeans(fit.bs$t, na.rm=T)

Aquí están, de hecho, en el campo de juego de lo que solía simular los datos originales:

> pars
[1]  5.606190  1.859591 -1.390816

Una versión trazada se ve así:

# Plot
with(dat, plot(x, y))

newx = seq(min(x), max(x), len=100)
lines(newx, SSasymp(newx, pars[1], pars[2], pars[3]))

lines(newx, SSasymp(newx, 5, 1, -1), col='red')
legend('bottomright', c('Actual', 'Predicted'), bty='n', lty=1, col=2:1)

ingrese la descripción de la imagen aquí

Ahora, si quiero la variabilidad de estas estimaciones de parámetros estabilizados , creo que puedo, suponiendo la normalidad de esta distribución de arranque, simplemente calcular sus errores estándar:

> apply(fit.bs$t, 2, function(x) sd(x, na.rm=T) / sqrt(length(na.omit(x))))
[1] 0.08369921 0.17230957 0.08386824

¿Es este un enfoque sensato? ¿Existe un mejor enfoque general para la inferencia sobre los parámetros de modelos no lineales inestables como este? (Supongo que podría hacer una segunda capa de remuestreo aquí, en lugar de depender de la teoría para el último bit, pero eso podría llevar mucho tiempo dependiendo del modelo. Aún así, no estoy seguro de si estos errores estándar podrían sería útil para cualquier cosa, ya que se acercarían a 0 si solo aumentara el número de réplicas de arranque).

Muchas gracias y, por cierto, soy ingeniero, así que perdóname por ser un novato relativo por aquí.

John Colby
fuente

Respuestas:

13

Hay varios problemas en esta pregunta. Primero, está la cuestión de si los promedios bootstrapped serán estimadores razonables incluso cuando algunos de los estimadores bootstrapped individuales no sean computables (falta de convergencia, inexistencia de soluciones). En segundo lugar, dado que los estimadores bootstrap son razonables, hay una cuestión de cómo obtener intervalos de confianza o quizás solo errores estándar para estas estimaciones.

La idea de promediar estimaciones bootstrap está estrechamente relacionada con, si no es realmente la misma, que la agregación bootstrap o el embolsado, utilizado en el aprendizaje automático para mejorar el rendimiento de predicción de predictores débiles. Ver ESL , Sección 8.7. En ciertos casos también para estimar parámetros el promedio de las estimaciones de bootstrap puede reducir la varianza del estimador resultante en comparación con solo usar el estimador en el conjunto de datos original.---

Sin embargo, el propósito de la pregunta es producir estimaciones incluso en los casos en que el algoritmo para calcular las estimaciones puede fallar ocasionalmente o en ocasiones el estimador no está definido. Como enfoque general hay un problema:

  • Hacer un promedio de las estimaciones de bootstrap y tirar a ciegas las muestras de bootstrap para las cuales las estimaciones no son computables en general dará resultados sesgados.

La gravedad del problema general depende de varias cosas. Por ejemplo, con qué frecuencia la estimación no es computable y si la distribución condicional de la muestra dado que la estimación no es computable difiere de la distribución condicional de la muestra dado que la estimación es computable. No recomendaría usar el método.

Xθ^θ^(X)Y

θ~(X)=mi(θ^(Y)X,UN(X))
UN(X)Xθ^(Y)N / A-XUN(X)θ~(X)

θ^(Y)XUN(X)θ~(X)

UN(X)θ~(X)

Editar :

El muy buen artículo Estimación y precisión después de la selección del modelo por Efron ofrece un método general para estimar el error estándar de un estimador en bolsas sin usar una segunda capa de arranque. El documento no trata explícitamente con estimadores que ocasionalmente no son computables.

NRH
fuente
Gracias por la excelente respuesta. El punto de sesgo es especialmente bien tomado. Puede imaginar un caso extremo donde la nube de puntos es totalmente uniforme, salvo por un solo conjunto de puntos lejanos que se ajustan muy bien al modelo. La gran mayoría de los nlsajustes pueden fallar, pero, de los que convergen, el sesgo será enorme y los errores estándar / ICs predichos espuriosamente pequeños. nlsBootutiliza un requisito ad hoc del 50% de ajustes exitosos, pero estoy de acuerdo con usted en que la (des) similitud de las distribuciones condicionales es igualmente preocupante.
John Colby
(Trataré de darte una bonificación mañana si este sitio me permite SO)
John Colby