Hi, So, I can draw/plot the usual 95% least squares confidence bands around the 3 fitted curves for MARS regression, but I don't know how to get the 95% BCa bootstrapped confidence bands.
Any help would be highly appreciated. Reproducible example : ################## Dataset = data.frame(PIBparHab=c(43931,67524,48348,44827,52409,15245,24453,57636,28992,17102,51495,47243,40908,22494,12784,48391,44221,32514,35132,46679,106022,9817,99635,38678,49128,12876,20732,17151,19670,41053,22488,57134,83295,10660), QUALITESANSREDONDANCE=c(1082.5,1066.6,1079.3,1079.9,1074.9,1008.6,1007.5,1111.3,1108.2,1109.7,1059.6,1165.1,1026.7,1035.1,997.8,1044.8,1073.6,1085.7,1083.8,1021.6,1036.2,1075.3,1069.3,1101.4,1086.9,1072.1,1166.7,983.9,1004.5,1082.5,1123.5,1094.9,1105.1,1010.8), competitivite=c(89,83,78,73,90,71,77,85,61,67,98,82,70,43,57,78,72,79,61,71,86,63,90,75,87,64,60,56,66,80,53,91,97,62), innovation=c(56,52,53,54,57,43,54,60,47,55,58,62,52,35,47,59,56,56,45,52,58,33,57,57,61,40,45,41,50,61,50,65,68,34)) install.packages("earth") library(earth) newdata=na.omit(Dataset) model=earth(PIBparHab ~ QUALITESANSREDONDANCE + competitivite + innovation,data=newdata, penalty=-1) summary(model) plot(model) plotmo(model) model2 <- earth(PIBparHab ~ QUALITESANSREDONDANCE + competitivite + innovation,data=newdata, penalty=-1, varmod.method="lm",nfold=10,ncross=3) plotmo(model2, pt.col=1, level=.95) boot.MARS=function(formula,data,indices) { d=data[indices,] fit=earth(formula,data=d) return(coef(fit)) } library(boot) results=boot(data=newdata, statistic=boot.MARS, R=1000,formula=PIBparHab ~ QUALITESANSREDONDANCE + competitivite + innovation) boot.ci(results, type= "bca", index=2) ################## ________________________________ De : Bert Gunter <bgunter.4...@gmail.com> Cc : R-help Mailing List <r-help@r-project.org> Envoyé le : Dimanche 25 septembre 2016 22h39 Objet : Re: [R] BCa confidence bands around fitted curves MARS regression Presumably the "earth" package lacks this functionality ...? So, obvious query: did you try using the boot package? If not, why not? If so, show us the code that failed. Or am I missing the point? Cheers, Bert Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip ) On Sun, Sep 25, 2016 at 12:19 PM, varin sacha via R-help <r-help@r-project.org> wrote: > Dear R-experts, > > I have fitted a MARS regression and am trying now to plot/draw the BCa > confidence bands around the 3 fitted curves (QUALITESANSREDONDANC, > competitivite and innovation). > > > Here is the reproducible example. > > > ############################ > > Dataset = > data.frame(PIBparHab=c(43931,67524,48348,44827,52409,15245,24453,57636,28992,17102,51495,47243,40908,22494,12784,48391,44221,32514,35132,46679,106022,9817,99635,38678,49128,12876,20732,17151,19670,41053,22488,57134,83295,10660), > > QUALITESANSREDONDANCE=c(1082.5,1066.6,1079.3,1079.9,1074.9,1008.6,1007.5,1111.3,1108.2,1109.7,1059.6,1165.1,1026.7,1035.1,997.8,1044.8,1073.6,1085.7,1083.8,1021.6,1036.2,1075.3,1069.3,1101.4,1086.9,1072.1,1166.7,983.9,1004.5,1082.5,1123.5,1094.9,1105.1,1010.8), > > competitivite=c(89,83,78,73,90,71,77,85,61,67,98,82,70,43,57,78,72,79,61,71,86,63,90,75,87,64,60,56,66,80,53,91,97,62), > > innovation=c(56,52,53,54,57,43,54,60,47,55,58,62,52,35,47,59,56,56,45,52,58,33,57,57,61,40,45,41,50,61,50,65,68,34)) > > install.packages("earth") > > library(earth) > > newdata=na.omit(Dataset) > > model=earth(PIBparHab ~ QUALITESANSREDONDANCE + competitivite + > innovation,data=newdata, penalty=-1) > > summary(model) > > plot(model) > > plotmo(model) > > > ############################ > > Best Regards, > S > > ______________________________________________ > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.