[R] Same results but different functions ?

[varin sacha]

Dear Michael, 
Dear Martin,

Many thanks for your suggestions.

Best,







Le lundi 23 mars 2020 à 22:34:41 UTC+1, Martin Maechler <maechler using stat.math.ethz.ch> a écrit : 





>>>>> 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 using 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 using 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.