Je ne suis pas sûr que les intervalles de confiance pour smooth.spline ont « belles » intervalles de confiance comme ceux sous forme lowess font. Mais j'ai trouvé un échantillon de code à partir d'un CMU Data Analysis course pour créer des intervalles de confiance bayésiens.
Voici les fonctions utilisées et un exemple. La fonction principale est spline.cis où le premier paramètre est une trame de données où la première colonne sont les valeurs x et la deuxième colonne les valeurs y. L'autre paramètre important est B qui indique le nombre de réplications bootstrap à effectuer. (Voir le PDF lien ci-dessus pour tous les détails.)
# Helper functions
resampler <- function(data) {
n <- nrow(data)
resample.rows <- sample(1:n,size=n,replace=TRUE)
return(data[resample.rows,])
}
spline.estimator <- function(data,m=300) {
fit <- smooth.spline(x=data[,1],y=data[,2],cv=TRUE)
eval.grid <- seq(from=min(data[,1]),to=max(data[,1]),length.out=m)
return(predict(fit,x=eval.grid)$y) # We only want the predicted values
}
spline.cis <- function(data,B,alpha=0.05,m=300) {
spline.main <- spline.estimator(data,m=m)
spline.boots <- replicate(B,spline.estimator(resampler(data),m=m))
cis.lower <- 2*spline.main - apply(spline.boots,1,quantile,probs=1-alpha/2)
cis.upper <- 2*spline.main - apply(spline.boots,1,quantile,probs=alpha/2)
return(list(main.curve=spline.main,lower.ci=cis.lower,upper.ci=cis.upper,
x=seq(from=min(data[,1]),to=max(data[,1]),length.out=m)))
}
#sample data
data<-data.frame(x=rnorm(100), y=rnorm(100))
#run and plot
sp.cis <- spline.cis(data, B=1000,alpha=0.05)
plot(data[,1],data[,2])
lines(x=sp.cis$x,y=sp.cis$main.curve)
lines(x=sp.cis$x,y=sp.cis$lower.ci, lty=2)
lines(x=sp.cis$x,y=sp.cis$upper.ci, lty=2)
Et cela donne quelque chose comme
En fait, il semble qu'il y ait peut-être une façon plus paramétrique pour calculer les intervalles de confiance en utilisant la les résidus de jackknife. Ce code provient de la S+ help page for smooth.spline
fit <- smooth.spline(data$x, data$y) # smooth.spline fit
res <- (fit$yin - fit$y)/(1-fit$lev) # jackknife residuals
sigma <- sqrt(var(res)) # estimate sd
upper <- fit$y + 2.0*sigma*sqrt(fit$lev) # upper 95% conf. band
lower <- fit$y - 2.0*sigma*sqrt(fit$lev) # lower 95% conf. band
matplot(fit$x, cbind(upper, fit$y, lower), type="plp", pch=".")
et que les résultats dans
Et pour autant que les intervalles de confiance gam vont, si vous lisez le fichier d'aide print.gam, il y a un paramètre se= avec par défaut TRUE et les docs disent
lorsque TRUE (de faille) les lignes supérieure et inférieure sont ajoutées aux tracés 1-d à 2 erreurs-types au-dessus et en-dessous de l'estimation de l'être lisse tracée tandis que pour les tracés 2-d, les surfaces aux erreurs -1 et -1 sont contournées et superposées sur le tracé de contour pour l'estimation. Si un nombre positif est fourni, ce nombre est multiplié par les erreurs standard lors du calcul des courbes d'erreur standard ou des surfaces. Voir aussi l'ombre, ci-dessous.
Vous pouvez donc ajuster l'intervalle de confiance en ajustant ce paramètre. (Ce serait dans l'appel print().
La première méthode en utilisant la méthode bootstrap de sens, mais les résultats donnent des motifs complètement différents comparant à celui que je reçois de 'gam()'. Celui qui utilise les résidus jackknife donne aussi un schéma très différent et les CI dans cette parcelle sont très très cahoteux. –
@YuDeng Eh bien, différentes méthodes donnent des résultats différents. Je suppose que vous devez décider si vous êtes un fréquentiste ou un bayésien. Vous voudrez peut-être consulter un statisticien pour voir quelle méthode est la plus appropriée pour vos données. – MrFlick