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"?)
# # # # # # # # # # # # # # # # # # # # # # # # #
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"?) > # # # # # # # # # # # # # # # # # # # # # # # # # > > ______________________________________________ > 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. >-- Michael http://www.dewey.myzen.co.uk/home.html
>>>>> 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 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. >> > -- > Michael > http://www.dewey.myzen.co.uk/home.html > ______________________________________________ > 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.