如何衡量一些回归模型对异常值的稳健性

如何解决如何衡量一些回归模型对异常值的稳健性

我将模型系数的方差相加，然后返回总和的均值。

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

但是，我的代码告诉我，普通最小二乘法是最好的，但这不是预期的结果，因为 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 ~ .,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,mu)
#####################################

更新由于第一个提供的答案，我改变了衡量稳健性的想法。

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

############ R CODE ##############
rmn <- function(n,seed = TRUE) {
 if (seed) set.seed(12345)
 p <- length(mu)
 matrix( rnorm(n * p,ncol = p)
}
##################################
out.cv <- function(y,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 ~.,-1])) ) - be) )
 b2[i] <- mean( abs( coef( rlm(z ~ .,data = data.frame(x[,-1]),method = "MM") ) - be ) )
 b3[i] <- mean( abs( coef( rlm(z ~ .,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),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[,method = "MM") ) - be ) )
 b3[i] <- mean( abs( coef( rlm(y ~ .,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,mu)
###########################
mu <-colMeans(x)+15 ## outliers in x
out.cv(y,mu)
###################

解决方法

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

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