【統計学】外れ値があってもいい感じで回帰直線を引いてみる

最小刈込み2乗法というものについてのメモです。
通常、回帰直線は誤差に正規分布を想定することが多く、通常の最小二乗法を適用すると外れ値の影響を受けやすいですが、外れ値を除外（刈り込む）ことでその影響を除いて直線を引いてみる手法です。

パッケージのインストールと読み込み。

install.packages("galts")
install.packages("ggplot2")
require(galts)
require(ggplot2)

通常の最小二乗法と最小刈込み2乗法の両方で回帰直線を引いてみる関数を定義

lts_test <- function(x, y) {
  df <- data.frame(x, y)
  
  result <- lm(y ~ x, data=df)
  summary(result)
  
  #Estimating LTS with ga (Default optimization method)
  lts <- ga.lts(y ~ x, popsize=40, iters=2, lower=-20, upper=20)
  print(lts)
  
  x2 <- seq(-4, 10, 0.1)
  y2 <- lts$coefficients[1] + x2*lts$coefficients[2]
  df2 <- data.frame(x2, y2)
  
  y1 = result[[1]][1] + result[[1]][2]*x2
  df1 <- data.frame(x2, y1)
  
  
  # 赤い線：通常の回帰直線、青い線：刈込み平均を適用した回帰直線
  ggplot() +
    layer(data=df, mapping=aes(x=x, y=y), geom="point",
          size=3, na.rm=TRUE) +
    layer(data=df1, mapping=aes(x=x2, y=y1), geom="line", color="Red") +
    layer(data=df2, mapping=aes(x=x2, y=y2), geom="line", color="Blue")
}

x軸方向に外れ値があった場合

################### x方向の外れ値
x <- rnorm(100)
e  <- rnorm(100)
y <- 5 + 5*x + e
outlyings <- sample(1:100,48)
x[outlyings] <- 10
lts_test(x, y)

y軸方向に外れ値があった場合

################### y方向の外れ値
x <- rnorm(100)
e  <- rnorm(100)
y <- 5 + 5*x + e
outlyings <- sample(1:100,48)
y[outlyings] <- 30
lts_test(x, y)

参考

経済データの回帰分析と頑健性 (美添泰人)
　　http://www.econ.aoyama.ac.jp/~yasuto_yoshizoe/econstat/stat200809.pdf

LTS (Least Trimmed Squares)
　　https://cran.r-project.org/web/packages/galts/galts.pdf