varin sacha
2016-Jul-03 11:47 UTC
[R] BCa Bootstrapped regression coefficients from lmrob function not working
Dear R-experts,
I am trying to calculate the bootstrapped (BCa) regression coefficients for a
robust regression using MM-type estimator (lmrob function from robustbase
package).
My R code here below is showing a warning message ([1] "All values of t are
equal to
22.2073014256803\n Can not calculate confidence intervals" NULL), I was
wondering if it was because I am trying to fit a robust regression with lmrob
function rather than a simple lm ? I mean maybe the boot.ci function does not
work with lmrob function ? If not, I was wondering what was going on ?
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))
library("robustbase")
newdata=na.omit(Dataset)
a=Dataset$PIBparHab
b=Dataset$QUALITESANSREDONDANCE
c=Dataset$competitivite
d=Dataset$innovation
fm.lmrob=lmrob(a~b+c+d,data=newdata)
fm.lmrob
boot.Lmrob=function(formula,data,indices) {
d=data[indices,]
fit=lmrob(formula,data=d)
return(coef(fit))
}
library(boot)
results=boot(data=newdata, statistic=boot.Lmrob, R=1000,formula=a~b+c+d)
boot.ci(results, type= "bca",index=2)
Any help would be highly appreciated,
S
peter dalgaard
2016-Jul-03 16:19 UTC
[R] BCa Bootstrapped regression coefficients from lmrob function not working
> On 03 Jul 2016, at 13:47 , varin sacha via R-help <r-help at r-project.org> wrote: > > Dear R-experts, > > I am trying to calculate the bootstrapped (BCa) regression coefficients for a robust regression using MM-type estimator (lmrob function from robustbase package). > > My R code here below is showing a warning message ([1] "All values of t are equal to > 22.2073014256803\n Can not calculate confidence intervals" NULL), I was wondering if it was because I am trying to fit a robust regression with lmrob function rather than a simple lm ? I mean maybe the boot.ci function does not work with lmrob function ? If not, I was wondering what was going on ?You need to review your code. You calculate a,b,c,d in the global environment and create newdata as a subset of Dataset, then use a,b,c,d in the formula, but no such variables are in newdata. AFAICT, all your bootstrap fits use the _same_ global values for a,b,c,d hence give the same result 1000 times... -pd> > 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)) > > library("robustbase") > newdata=na.omit(Dataset) > a=Dataset$PIBparHab > b=Dataset$QUALITESANSREDONDANCE > c=Dataset$competitivite > d=Dataset$innovation > > fm.lmrob=lmrob(a~b+c+d,data=newdata) > fm.lmrob > > boot.Lmrob=function(formula,data,indices) { > d=data[indices,] > fit=lmrob(formula,data=d) > return(coef(fit)) > } > > library(boot) > results=boot(data=newdata, statistic=boot.Lmrob, R=1000,formula=a~b+c+d) > boot.ci(results, type= "bca",index=2) > > > Any help would be highly appreciated, > 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.-- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Office: A 4.23 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
varin sacha
2016-Jul-05 21:51 UTC
[R] BCa Bootstrapped regression coefficients from lmrob function not working
Dear Professor Dalgaard,
I really thank you lots for your response. I have solved my problem. Now, I have
tried to do the same (calculate the BCa bootstrapped CIs) for the MARS
regression, and I get an error message. If somebody has a hint to solve my
problem, 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)
summary(model)
plot(model)
plotmo(model)
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)
Best,
S
________________________________
De : peter dalgaard <pdalgd at gmail.com>
Cc : R-help Mailing List <r-help at r-project.org>
Envoy? le : Dimanche 3 juillet 2016 18h19
Objet : Re: [R] BCa Bootstrapped regression coefficients from lmrob function not
working
> On 03 Jul 2016, at 13:47 , varin sacha via R-help <r-help at
r-project.org> wrote:
>
> Dear R-experts,
>
> I am trying to calculate the bootstrapped (BCa) regression coefficients for
a robust regression using MM-type estimator (lmrob function from robustbase
package).
>
> My R code here below is showing a warning message ([1] "All values of
t are equal to
> 22.2073014256803\n Can not calculate confidence intervals" NULL), I
was wondering if it was because I am trying to fit a robust regression with
lmrob function rather than a simple lm ? I mean maybe the boot.ci function does
not work with lmrob function ? If not, I was wondering what was going on ?
You need to review your code. You calculate a,b,c,d in the global environment
and create newdata as a subset of Dataset, then use a,b,c,d in the formula, but
no such variables are in newdata. AFAICT, all your bootstrap fits use the _same_
global values for a,b,c,d hence give the same result 1000 times...
-pd
>
> 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))
>
> library("robustbase")
> newdata=na.omit(Dataset)
> a=Dataset$PIBparHab
> b=Dataset$QUALITESANSREDONDANCE
> c=Dataset$competitivite
> d=Dataset$innovation
>
> fm.lmrob=lmrob(a~b+c+d,data=newdata)
> fm.lmrob
>
> boot.Lmrob=function(formula,data,indices) {
> d=data[indices,]
> fit=lmrob(formula,data=d)
> return(coef(fit))
> }
>
> library(boot)
> results=boot(data=newdata, statistic=boot.Lmrob, R=1000,formula=a~b+c+d)
> boot.ci(results, type= "bca",index=2)
>
>
> Any help would be highly appreciated,
> 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.
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com