[R] BCa confidence bands around fitted curves MARS regression

[varin sacha]

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.4567 at gmail.com>

Cc : R-help Mailing List <r-help at 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 at 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 at 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.