1

I am building a plot with ggplot. I have data where y is mostly independent of X, but I randomly have a few extreme values of Y at low values of X. Like this:

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)

I want to make the point that the MEDIAN y-value is still constant over X values. I can see that this is basically true here:

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

If I make a geom_smooth() plot, it does mean, and is very affected by outliers:

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

I have a few potential fixes. For example, I could first use group_by/summarize to make a dataset of binned medians and then plot that. I would rather NOT do this because in my real data I have a lot of facetting and grouping variables, and it would be a lot to keep track of (non-ideal). A lot plot definitely looks better, but log does not have nice interpretation in my application (median does have nice interpretation)

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

Finally, I know about geom_quantile but I think I'm using it wrong. Is there a way to add an error bar? Also- this geom_quantile plot looks way too smooth, and I don't understand why it is sloping down. Am I using it wrong?

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

I realize that this problem probably has a LOT of workarounds, but if possible I would love to use geom_smooth and just provide an argument that tells it to use a median. I want geom_smooth for a side-by-side comparison with consistency. I want to put the mean and median geom_smooths side-by-side to show "hey look, super strong pattern between Y and X is driven by a few large outliers, if we look only at median the pattern disappears".

Thanks!!

2
  • Under the hood, geom_smooth fits a model. Depending on the size of data, loess or gam is used, which both fit to the mean. When you use geom_quantile, it's using quantile regression from the quantreg package, which fits a linear approximation of the quantiles Dec 4, 2020 at 16:58
  • You can see the built-in smoothing methods options at the ?geom_smooth help page. None of them (to my knowledge) are designed to be robust to outliers. You may have some luck using method.args() to adjust the ?loess or ?mgcv::gam smoothers, but you may be better off asking a methods question on stats.stackexchange for recommendations on how to get a smoothed estimate of a median. Dec 4, 2020 at 17:03

2 Answers 2

Reset to default

Trending sort

Trending sort is based off of the default sorting method — by highest score — but it boosts votes that have happened recently, helping to surface more up-to-date answers.

It falls back to sorting by highest score if no posts are trending.

6

You can create your own method to use in geom_smooth. As long as you have a function that produces an object on which the predict generic works to take a data frame with a column called x and translate into appropriate values of y.

As an example, let's create a simple model that interpolates along a running median. We wrap it in its own class and give it its own predict method:

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)
}

Now we can use our method in geom_smooth:

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

enter image description here

Now of course, this doesn't look very "flat", but it is way flatter than the line calculated by the loess method of the standard 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)

enter image description here

Now, I understand that this is not the same thing as "regressing on the median", so you may wish to explore different methods, but if you want to get geom_smooth to plot them, this is how you can go about it. Note that if you want standard errors, you will need to have your predict function return a list with members called fit and se.fit

1
  • 5
    Very nice demonstration and explanation! Dec 4, 2020 at 21:22
0

Here's a modification of @Allan's answer that uses a fixed x window rather than a fixed number of points. This is useful for irregular time series and series with multiple observations at the same time (x value). It uses a loop so it's not very efficient and will be slow for larger data sets.

# 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 & 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))

Created on 2022-01-05 by the reprex package (v2.0.1)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.