Блог

Я строю график с помощью ggplot. У меня есть данные, где y в основном не зависит от X, но у меня случайно есть несколько экстремальных значений Y при низких значениях X. Вот так:

 set.seed(1)
X <- rnorm(500, mean=5)
y <- rnorm(500)
y[X < 3] <- sample(c(0, 1000), size=length(y[X < 3]),prob=c(0.9, 0.1),  
    replace=TRUE)
 

Я хочу подчеркнуть, что МЕДИАННОЕ значение y по-прежнему остается постоянным при значениях X. Я вижу, что здесь это в основном верно:

 mean(y[X < 3])
median(y[X < 3])
 

Если я создам график geom_smooth (), он действительно означает и сильно зависит от выбросов:

 ggplot(data=NULL, aes(x=X, y=y))   geom_smooth()
 

У меня есть несколько потенциальных исправлений. Например, я мог бы сначала использовать group_by/summary для создания набора данных из привязанных медиан, а затем построить его. Я бы предпочел ЭТОГО НЕ делать, потому что в моих реальных данных у меня много переменных фасетирования и группировки, и было бы сложно отслеживать (не идеально). Большой график определенно выглядит лучше, но в моем приложении log не имеет хорошей интерпретации (у медианы есть хорошая интерпретация)

 ggplot(data=NULL, aes(x=X, y=y))   geom_smooth()   
  scale_y_log10()
 

Наконец, я знаю о geom_quantile, но я думаю, что использую его неправильно. Есть ли способ добавить панель ошибок? Кроме того, этот график geom_quantile выглядит слишком гладким, и я не понимаю, почему он наклонен вниз. Я использую это неправильно?

 ggplot(data=NULL, aes(x=X, y=y))   
  geom_quantile(quantiles=c(0.5))
 

Я понимаю, что у этой проблемы, вероятно, МНОГО обходных путей, но, если возможно, я бы хотел использовать geom_smooth и просто предоставить аргумент, который указывает ему использовать медиану. Я хочу geom_smooth для параллельного сравнения с согласованностью. Я хочу поместить среднее и среднее значения geom_smooths рядом друг с другом, чтобы показать «эй, смотри, супер сильная закономерность между Y и X обусловлена несколькими большими выбросами, если мы посмотрим только на медиану, шаблон исчезнет».

Спасибо!!

Вы можете создать свой собственный метод для использования geom_smooth . До тех пор, пока у вас есть функция, которая создает объект, над которым predict работает generic, чтобы взять фрейм данных с вызываемым столбцом x и преобразовать в соответствующие значения y .

В качестве примера давайте создадим простую модель, которая интерполирует вдоль текущей медианы. Мы оборачиваем его в его собственный класс и даем ему его собственный predict метод:

 rolling_median <- function(formula, data, n_roll = 11, ...) {
  x <- data$x[order(data$x)]
  y <- data$y[order(data$x)]
  y <- zoo::rollmedian(y, n_roll, na.pad = TRUE)
  structure(list(x = x, y = y, f = approxfun(x, y)), class = "rollmed")
}

predict.rollmed <- function(mod, newdata, ...) {
  setNames(mod$f(newdata$x), newdata$x)
}
 

Теперь мы можем использовать наш метод в geom_smooth :

 ggplot(data = NULL, aes(x = X, y = y))   
  geom_smooth(formula = y ~ x, method = "rolling_median", se = FALSE) 
 

Теперь, конечно, это выглядит не очень «плоским», но оно намного более плоское, чем линия, рассчитанная стандартным методом Лесса geom_smooth() :

 ggplot(data = NULL, aes(x = X, y = y))   
  geom_smooth(formula = y ~ x, color = "red", se = FALSE)  
  geom_smooth(formula = y ~ x, method = "rolling_median", se = FALSE)
 

Теперь я понимаю, что это не то же самое, что «регрессия по медиане», поэтому вы можете изучить другие методы, но если вы хотите geom_smooth их построить, вот как вы можете это сделать. Обратите внимание, что если вам нужны стандартные ошибки, вам нужно, чтобы ваша predict функция возвращала список с вызываемыми элементами fit и se.fit

Вот модификация ответа @Allan, которая использует фиксированное окно x, а не фиксированное количество точек. Это полезно для нерегулярных временных рядов и рядов с несколькими наблюдениями одновременно (значение x). Он использует цикл, поэтому он не очень эффективен и будет медленным для больших наборов данных.

 # running median with time window

library(dplyr)
library(ggplot2)
library(zoo)

# some irregular and skewed data
set.seed(1)
x <- seq(2000, 2020, length.out = 400) # normal time series, gives same result for both methods
x <- sort(rep(runif(40, min = 2000, max = 2020), 10)) # irregular and repeated time series
y <- exp(runif(length(x), min = -1, max = 3))
data <- data.frame(x = x, y = y)
# ggplot(data)   geom_point(aes(x = x, y = y))

# 2 year window
xwindow <- 2
nwindow <- xwindow * length(x) / 20 - 1

# rolling median
rolling_median <- function(formula, data, n_roll = 11, ...) {
  x <- data$x[order(data$x)]
  y <- data$y[order(data$x)]
  y <- zoo::rollmedian(y, n_roll, na.pad = TRUE)
  structure(list(x = x, y = y, f = approxfun(x, y)), class = "rollmed")
}

predict.rollmed <- function(mod, newdata, ...) {
  setNames(mod$f(newdata$x), newdata$x)
}

# rolling time window median
rolling_median2 <- function(formula, data, xwindow = 2, ...) {
  x <- data$x[order(data$x)]
  y <- data$y[order(data$x)]
  ys <- rep(NA, length(x)) # for the smoothed y values
  xs <- setdiff(unique(x), NA) # the unique x values
  i <- 1 # for testing
  for (i in seq_along(xs)){
    j <- xs[i] - xwindow/2 < x amp; x < xs[i]   xwindow/2 # x points in this window
    ys[x == xs[i]] <- median(y[j], na.rm = TRUE) # y median over this window
  }
  y <- ys
  structure(list(x = x, y = y, f = approxfun(x, y)), class = "rollmed2")
}

predict.rollmed2 <- function(mod, newdata, ...) {
  setNames(mod$f(newdata$x), newdata$x)
}

# plot smooth
ggplot(data)  
  geom_point(aes(x = x, y = y))  
  geom_smooth(aes(x = x, y = y, colour = "nwindow"), formula = y ~ x, method = "rolling_median", se = FALSE, method.args = list(n_roll = nwindow))  
  geom_smooth(aes(x = x, y = y, colour = "xwindow"), formula = y ~ x, method = "rolling_median2", se = FALSE, method.args = list(xwindow = xwindow))
 

^{Создано 2022-01-05 пакетом reprex (v2.0.1)}