(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)
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í.
fuente
nls
ajustes pueden fallar, pero, de los que convergen, el sesgo será enorme y los errores estándar / ICs predichos espuriosamente pequeños.nlsBoot
utiliza 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.