>>>>> Michael Dewey
>>>>> on Mon, 23 Mar 2020 13:45:44 +0000 writes:
> The documentation suggests that the rlm method for a formula does not
> have psi as a parameter. Perhaps try using the method for a matrix x and
> a vector y.
> Michael
or use lmrob() from pkg robustbase which is at least one
generation more recent and also with many more options than
rlm().
rlm() has been fantastic when it was introduced (into S /
S-plus, before R existed [in a publicly visible way]) but it had
been based of what was available back then, end of the 80's, beginning 90's.
Martin
> On 23/03/2020 12:39, varin sacha via R-help wrote:
>> Dear R-experts,
>>
>> The rlm command in the MASS package command implements several versions
of robust regression, for example the Huber and the Tukey (bisquare weighting
function) estimators.
>> In my R code here below I try to get the Tukey (bisquare weighting
function) estimation, R gives me an error message : Error in statistic(data,
original, ...) : unused argument (psi = psi.bisquare)
>> If I cancel psi=psi.bisquare my code is working but IMHO I will get the
Huber estimation and not the Tukey. So how can I get the Tukey ? Many thanks
for your help.
>>
>>
>> # # # # # # # # # # # # # # # # # # # # # # # #
>> install.packages( "boot",dependencies=TRUE )
>> install.packages( "MASS",dependencies=TRUE )
>> library(boot)
>> library(MASS)
>>
>> n<-50
>> b<-runif(n, 0, 5)
>> z <- rnorm(n, 2, 3)
>> a <- runif(n, 0, 5)
>>
>> y_model<- 0.1*b - 0.5 * z - a + 10
>> y_obs <- y_model +c( rnorm(n*0.9, 0, 0.1), rnorm(n*0.1, 0, 0.5) )
>> df<-data.frame(b,z,a,y_obs)
>>
>> # function to obtain MSE
>> MSE <- function(data, indices, formula) {
>> d <- data[indices, ] # allows boot to select sample
>> fit <- rlm(formula, data = d)
>> ypred <- predict(fit)
>> d[["y_obs "]] <-y_obs
>> mean((d[["y_obs"]]-ypred)^2)
>> }
>>
>> # Make the results reproducible
>> set.seed(1234)
>>
>> # bootstrapping with 600 replications
>> results <- boot(data = df, statistic = MSE,
>> R = 600, formula = y_obs ~ b+z+a, psi = psi.bisquare)
>>
>> str(results)
>>
>> boot.ci(results, type="bca" )
>> # # # # # # # # # # # # # # # # # # # # # # # # #
>>
>> ______________________________________________
>> 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.
>>
> --
> Michael
> http://www.dewey.myzen.co.uk/home.html
> ______________________________________________
> 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.
