Theil-sen comparison

TheilSen.R

library(zyp)
library(mblm)

# From https://library.virginia.edu/data/articles/theil-sen-regression-programming-and-understanding-an-outlier-resistant-alternative-to-least-squares
theil_sen <- function(x,y){
  n <- length(x)
  max_n_slopes <- (n * (n - 1)) / 2
  slopes <- vector(mode = 'list', length = max_n_slopes) # list of NULL elements
  add_at <- 1
  # Calculate point-to-point slopes if x_i != x_j
  for (i in 1:(n - 1)) {
    slopes_i <- lapply((i + 1):n, function(j) if (x[j] != x[i]) { (y[j] - y[i]) / (x[j] - x[i]) })
    slopes[add_at:(add_at + length(slopes_i) - 1)] <- slopes_i
    add_at <- add_at + length(slopes_i)
  }
  # Calculate Theil-Sen slope
  slopes <- unlist(slopes) # drop NULL elements (indefinite slopes/unfilled list entries)
  theil_sen_slope <- median(slopes, na.rm = TRUE)
  # Calculate Theil-Sen intercept
  intercepts <- y - (theil_sen_slope * x)
  theil_sen_intercept <- median(intercepts, na.rm = TRUE)
  # Return
  c('Theil-Sen intercept' = theil_sen_intercept, 'Theil-Sen slope' = theil_sen_slope)
}


set.seed(123)

# Generate linear data
x <- 1:100
y <- 2*x + rnorm(100, 0, 10)

# Add some outliers
y[c(10, 30, 50, 70, 90)] <- y[c(10, 30, 50, 70, 90)] + 50

plot(y~x)

zyp.mod <- zyp.sen(y~x)
zyp.mod

mblm.mod <- mblm(y~x,repeated = F)
summary(mblm.mod)

UV.theilsen <- theil_sen(x,y)
UV.theilsen