Same results but different functions ?

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Same results but different functions ?

R help mailing list-2
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" )
# # # # # # # # # # # # # # # # # # # # # # # # #

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Same results but different functions ?

Michael Dewey-3
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

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" )
> # # # # # # # # # # # # # # # # # # # # # # # # #
>
> ______________________________________________
> [hidden email] 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Same results but different functions ?

Martin Maechler
>>>>> 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" )
    >> # # # # # # # # # # # # # # # # # # # # # # # # #
    >>
    >> ______________________________________________
    >> [hidden email] 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

    > ______________________________________________
    > [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: Same results but different functions ?

R help mailing list-2
Dear Michael,
Dear Martin,

Many thanks for your suggestions.

Best,







Le lundi 23 mars 2020 à 22:34:41 UTC+1, Martin Maechler <[hidden email]> 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" )
    >> # # # # # # # # # # # # # # # # # # # # # # # # #
    >>
    >> ______________________________________________
    >> [hidden email] 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


    > ______________________________________________
    > [hidden email] 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.

______________________________________________
[hidden email] 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.