[R] BCa Bootstrapped regression coefficients from lmrob function not working

Thu Jul 7 00:37:57 CEST 2016

Dear Duncan,

Many thanks for your reply. What you propose is a very good idea yes. However, I am interested in the coefficients...

Best,
S

________________________________
De : Duncan Murdoch <murdoch.duncan at gmail.com>

Cc : R-help Mailing List <r-help at r-project.org>
Envoyé le : Mercredi 6 juillet 2016 14h13
Objet : Re: [R] BCa Bootstrapped regression coefficients from lmrob function not working

On 06/07/2016 7:56 AM, varin sacha via R-help wrote:
> Dear Professor Dalgaard,
>
> Okay, this is what I was afraid of.
> Many thanks for your response. I know that my next question is off-topic here (on this website) but is it nevertheless possible to calculate confidence intervals for MARS regression or CIs for MARS will just be impossible to calculate ?

If you think of the fitted curve or surface as the output rather than 
the coefficients that define it, it makes sense to get pointwise 
confidence intervals, or joint confidence regions for multiple 
locations.  Just do predictions at a set of fixed locations.

Duncan Murdoch

> Best
> S
>
>
>
> ----- Mail original -----
> De : peter dalgaard <pdalgd at gmail.com>

> Cc : Bert Gunter <bgunter.4567 at gmail.com>; R-help Mailing List <r-help at r-project.org>
> Envoyé le : Mercredi 6 juillet 2016 11h05
> Objet : Re: [R] BCa Bootstrapped regression coefficients from lmrob function not working
>
> Offhand, I would suspect that the cause is that earth() does not return a coefficient vector of the same length in every bootstrap iteration (_adaptive_ regression splines). This makes the bootstrap rather tricky to even define.
>
> -pd
>

>>
>> Dear Bert,
>>
>> You are right.
>>
>>> results=boot(data=newdata, statistic=boot.MARS, R=1000,formula=PIBparHab ~ QUALITESANSREDONDANCE + competitivite + innovation)
>> Erreur dans boot(data = newdata, statistic = boot.MARS, R = 1000, formula = PIBparHab ~  :
>> le nombre d'objets à remplacer n'est pas multiple de la taille du remplacement
>>
>> In English it would be something like : number of items to replace is not a multiple of replacement length
>>
>>
>> Best
>> S
>>
>> ________________________________
>> De : Bert Gunter <bgunter.4567 at gmail.com>

>> Cc : peter dalgaard <pdalgd at gmail.com>; R-help Mailing List <r-help at r-project.org>
>> Envoyé le : Mercredi 6 juillet 2016 1h19
>> Objet : Re: [R] BCa Bootstrapped regression coefficients from lmrob function not working
>>
>>
>> It would help to show your  error message, n'est-ce pas?
>>
>> 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 Tue, Jul 5, 2016 at 2:51 PM, varin sacha via R-help
>> <r-help at r-project.org> wrote:
>>> 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
>
>>
>>>
>>> ______________________________________________
>>> 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.
>