Same results but different functions ?

4 messages
Open this post in threaded view
|

Same results but different functions ?

 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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

Re: Same results but different functions ?

 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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.