varin sacha
2016-Sep-25 19:19 UTC
[R] BCa confidence bands around fitted curves MARS regression
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
Bert Gunter
2016-Sep-25 20:39 UTC
[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.
I have a data frame that contains data for multiple (seven) subjects. Each subject is represented by a new value of PID. I would like to plot the data for all seven subjects. For each subject I want to plot a line showing CT as a function of Nit, with the dots for each subject joined. I have tried to accomplish this using the by function. I get an error message, Error in match.fun(panel) : 'xx[, "CT"]' is not a function, character or symbol I have no idea what is causing the error, nor how to correct the error, nor how to get the dots for each point be connected by a line. Any help would be appreciated! PID <- c( 1 , 1 , 1 , 1 , 2, 2, 2, 2, 3 , 3 , 3 , 3 , 4 , 4, 4, 4 , 5, 5, 5, 5, 6, 6, 6, 6, 7 , 7 , 7 , 7) Nit <- c(NA , -9.23,-11.61,-7.88,NA,NA,NA,NA,-5.59, 0.73,-10.55, -9.13, 3.67, NA, NA,-13.26,NA,NA,NA,NA,NA,NA,NA,NA,-9.36, 5.08, -5.73, 2.02) CT <- c(544,459 ,432 ,NA ,NA,NA,NA,NA,1398 ,1287 ,1049 , NA ,543 ,474,507,NA ,NA,NA,NA,NA,NA,NA,NA,NA,992 ,992 ,1078 ,NA) xx <- data.frame(PID=PID,Nit=Nit,CT=CT) xx by(xx,as.factor(xx[,"PID"]),plot,xx[,"Nit"],xx[,"CT"]) John David Sorkin M.D., Ph.D. Professor of Medicine Chief, Biostatistics and Informatics University of Maryland School of Medicine Division of Gerontology and Geriatric Medicine Baltimore VA Medical Center 10 North Greene Street GRECC (BT/18/GR) Baltimore, MD 21201-1524 (Phone) 410-605-7119 (Fax) 410-605-7913 (Please call phone number above prior to faxing) Confidentiality Statement: This email message, including any attachments, is for the sole use of the intended recipient(s) and may contain confidential and privileged information. Any unauthorized use, disclosure or distribution is prohibited. If you are not the intended recipient, please contact the sender by reply email and destroy all copies of the original message.
varin sacha
2016-Oct-02 22:11 UTC
[R] BCa confidence bands around fitted curves MARS regression
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.