(ignore o código R, se necessário, pois minha pergunta principal é independente do idioma)
Se eu quiser olhar para a variabilidade de uma estatística simples (ex: média), sei que posso fazê-lo via teoria como:
x = rnorm(50)
# Estimate standard error from theory
summary(lm(x~1))
# same as...
sd(x) / sqrt(length(x))
ou com o 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)
No entanto, o que eu quero saber é: pode ser útil / válido (?) Procurar o erro padrão de uma distribuição de inicialização em determinadas situações? A situação com a qual estou lidando é uma função não linear relativamente barulhenta, 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)
Aqui, o modelo nem converge usando o conjunto de dados 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
portanto, as estatísticas nas quais estou interessado são estimativas mais estabilizadas desses parâmetros nls - talvez seus meios em várias replicações de autoinicialização.
# 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)
Aqui estão, de fato, no parque de bolas do que eu costumava simular os dados originais:
> pars
[1] 5.606190 1.859591 -1.390816
Uma versão plotada se parece com:
# 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)
Agora, se eu quiser a variabilidade dessas estimativas de parâmetros estabilizados , acho que posso, assumindo a normalidade dessa distribuição de bootstrap, apenas calcular seus erros padrão:
> apply(fit.bs$t, 2, function(x) sd(x, na.rm=T) / sqrt(length(na.omit(x))))
[1] 0.08369921 0.17230957 0.08386824
Essa é uma abordagem sensata? Existe uma abordagem geral melhor para inferência sobre os parâmetros de modelos não-lineares instáveis como este? (Suponho que eu poderia fazer uma segunda camada de reamostragem aqui, em vez de confiar na teoria do último bit, mas isso pode levar muito tempo dependendo do modelo. Mesmo assim, não tenho certeza se esses erros padrão ser útil para qualquer coisa, pois eles se aproximariam de 0 se eu apenas aumentasse o número de replicações de inicialização.)
Muito obrigado e, a propósito, sou engenheiro, por favor, perdoe-me sendo um novato por aqui.
fonte
nls
ajustes pode falhar, mas, dentre os que convergirem, o viés será enorme e os erros / ICs padrão previstos espuriosamente pequenos.nlsBoot
usa um requisito ad hoc de 50% de ajustes bem-sucedidos, mas concordo com você que a (des) semelhança das distribuições condicionais é igualmente uma preocupação.