(ignorez le code R si nécessaire, car ma question principale est indépendante de la langue)
Si je veux regarder la variabilité d'une statistique simple (ex: moyenne), je sais que je peux le faire via la théorie comme:
x = rnorm(50)
# Estimate standard error from theory
summary(lm(x~1))
# same as...
sd(x) / sqrt(length(x))
ou avec le bootstrap comme:
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)
Cependant, ce que je me demande, peut-il être utile / valide (?) De rechercher l' erreur standard d'une distribution bootstrap dans certaines situations? La situation à laquelle je fais face est une fonction non linéaire relativement bruyante, telle que:
# 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)
Ici, le modèle ne converge même pas en utilisant l'ensemble de données d'origine,
> (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
donc les statistiques qui m'intéressent plutôt sont des estimations plus stabilisées de ces paramètres nls - peut-être leurs moyens à travers un certain nombre de réplications bootstrap.
# 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)
Les voici, en effet, dans le parc à billes de ce que j'ai utilisé pour simuler les données originales:
> pars
[1] 5.606190 1.859591 -1.390816
Une version tracée ressemble à:
# 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)
Maintenant, si je veux la variabilité de ces estimations de paramètres stabilisés , je pense que je peux, en supposant la normalité de cette distribution bootstrap, simplement calculer leurs erreurs standard:
> apply(fit.bs$t, 2, function(x) sd(x, na.rm=T) / sqrt(length(na.omit(x))))
[1] 0.08369921 0.17230957 0.08386824
Est-ce une approche sensée? Existe-t-il une meilleure approche générale de l'inférence sur les paramètres de modèles non linéaires instables comme celui-ci? (Je suppose que je pourrais plutôt faire une deuxième couche de rééchantillonnage ici, au lieu de compter sur la théorie pour le dernier bit, mais cela pourrait prendre beaucoup de temps selon le modèle. Même encore, je ne suis pas sûr si ces erreurs standard être utile pour n'importe quoi, car ils approcheraient de 0 si j'augmentais simplement le nombre de réplications de bootstrap.)
Merci beaucoup, et en passant, je suis ingénieur alors s'il vous plaît pardonnez-moi d'être un novice relatif ici.
la source
nls
ajustements peuvent échouer, mais, parmi ceux qui convergent, le biais sera énorme et les erreurs / IC standard prédits sont extrêmement faibles.nlsBoot
utilise une exigence ad hoc de 50% d'ajustements réussis, mais je suis d'accord avec vous que la (dis) similitude des distributions conditionnelles est également une préoccupation.