时间戳

首页 > 解决方案 > 如何衡量一些回归模型对异常值的鲁棒性

r - 如何衡量一些回归模型对异常值的鲁棒性

问题描述

我正在添加模型系数的方差，然后返回总和的平均值。

我只是想检查哪种回归方法对异常值更稳健。我将研究许多场景。

但是，我的代码告诉我普通的最小二乘是最好的，但这不是预期的结果，因为 MM 估计和 Huber 被称为鲁棒回归方法。

我的代码做错了吗？

#####################################
rmn <- function(n, mu) {
  p <- length(mu)
  matrix(rnorm(n*p, mean = mu), ncol = p)
}
#####################################
RI<-function(y,x,a,mu,R=30,t=1000){
  x <- as.matrix(x)
  dm <- dim(x)
  n <- dm[1]  
  bias1 <- bias2 <- bias3 <- numeric(t)
  b1 <- b2<- b3 <- numeric(R) 
  ### Outliers in X ######
  for (j in 1:t) {
    for (i in 1:R) {
      id <- sample(n, a * n)
      z <- x
      z[id, ] <- rmn(length(id), mu)  
      b1[i] <- var(coef(lm(y ~., data = data.frame(z))))
      b2[i] <- var(coef(rlm(y ~ ., data = data.frame(z), maxit = 2000, method = "MM")))
      b3[i] <- var(coef(rlm(y ~ ., data = data.frame(z), psi = psi.huber,maxit = 300)))
    }
    bias1[j] <- sum(b1); bias2[j] <- sum(b2); bias3[j] <- sum(b3)
  }
  bias <- cbind("lm" = bias1,"MM-rlm" = bias2, "H-rlm" = bias3)
  colMeans(bias)
}
#####################################
p <- 5
n <- 300

x<- matrix(rnorm(n * p), ncol = p)
y<-rnorm(n)
a=0.2
mu <-colMeans(x)+10
#####################################
RI(y,x,a,mu)
#####################################

更新 由于第一个提供的答案，我改变了测量鲁棒性的想法。

我通过计算数据未受污染和数据受污染时的系数之间的平均绝对差来衡量稳健性。我首先在 y 中引入异常值，然后在 x 中引入异常值。我还是有问题。

############ R CODE ##############
rmn <- function(n, mu, seed = TRUE) {
 if (seed) set.seed(12345)
 p <- length(mu)
 matrix( rnorm(n * p, mean = mu), ncol = p)
}
##################################
out.cv <- function(y, x, a, mu, R = 500, seed = TRUE) {
 ## y: response variable
 ## x: independent variables
 ## a: percent of outliers
 ## mu: how far should the outliers be. A vector if outliers in x,
 ## or a single number if outliers in y
 ## R: how many times to repeat this process
 x <- as.matrix(x)
 dm <- dim(x)
 n <- dm[1] ; d <- dm[2] + 1
 b1 <- b2<- b3 <- numeric(R)
 be <- coef( lm(y ~., data = as.data.frame(x[,-1]) ) )
####################################
 ### Outliers in Y ######
 if ( length(mu) == 1 ) {
 for (i in 1:R) {
 if (seed) set.seed(12345)
 id <- sample(n, a * n)
 z <- y
 if (seed) set.seed(12345)
 z[id] <- rnorm(id, mu) ## mu has to be a single number here
 ## mean absolute difference between coefficients of clean data
 ## and coefficients with contaminated data
 b1[i] <- mean( abs( coef( lm(z ~., data = as.data.frame(x[,-1])) ) - be) )
 b2[i] <- mean( abs( coef( rlm(z ~ ., data = data.frame(x[,-1]), maxit = 2000, method = "MM") ) - be ) )
 b3[i] <- mean( abs( coef( rlm(z ~ ., data = data.frame(x[,-1]), psi = psi.huber,maxit = 300) ) - be ) )
 }
########################
##### Outliers in X #########
 } else {
 for (i in 1:R) {
 if (seed) set.seed(12345)
 id <- sample(n, a * n)
 z <- x
 z[id, ] <- rmn( length(id), mu, seed ) ## mu must be a vector
 b1[i] <- mean( abs( coef( lm(y ~., data = as.data.frame(z[,-1])) )- be) )
 b2[i] <- mean( abs( coef( rlm(y ~ ., data = data.frame(z[,-1]), maxit = 2000, method = "MM") ) - be ) )
 b3[i] <- mean( abs( coef( rlm(y ~ ., data = data.frame(z[,-1]), psi = psi.huber,maxit = 300) ) - be ) )
 }
 }
 bias1 <- mean(b1) ; bias2 <- mean(b2); bias3 <- mean(b3)
 bias <- c(bias1, bias2, bias3)
 names(bias) <- c("lm", "MM-rlm","Huber-rlm")
 bias
}
################################
p <- 5
n <- 200
##############################
# Independent X and Y ####
#set.seed(12345)
#x<- matrix( rnorm(n * p), ncol = p)
#y<-rnorm(n)

## Related X and Y ####
set.seed(12345)
x <- rmn(n, numeric(p))
ber <- rnorm(p)
m <- x %*% ber
y <- rnorm(n, m, 1)

############################
a <- 0.2 #outliers 10%
mu <- 15 ## outliers in y
out.cv(y, x, a, mu)
###########################
mu <-colMeans(x)+15 ## outliers in x
out.cv(y, x, a, mu)
###################

标签： rregressionlinear-regression

解决方案

首先，我没有看到您从长尾分布中生成样本。请使用rt(n, 3)具有非常小的 df 的 t 学生来获得这样的分布或与其他像对数正态分布一样玩。因此不要rnorm确定使用。我看到你使用了一些看起来过于复杂的注射产品。另一件事是 MASS::rlm 的规范并不是那么简单。在我看来，从quantreg::rq分位数回归开始，并将其视为稳健的基准方法。

此外，您的抽样程序看起来不是一个有效的程序。您每次迭代都会生成一个新的观察结果，这些观察结果是先验未知的。我希望在训练集或测试集上进行引导。

推荐阅读