Quantile-based bias correction

# simulate data
x <- rnorm(1000)
y <- x + rnorm(1000)

# split into training/validation/test
x.train <- x[1:600]
y.train <- y[1:600]
x.valid <- x[601:800]
y.valid <- y[601:800]
x.test <- x[801:1000]
y.test <- y[801:1000]

# generate linear model and predict on validation data
mod <- lm(y.train ~ x.train)
pred.valid <- predict(mod, list(x.train = x.valid))

# set plot parameters
par(mar = c(5,5,1,1), las = 1)

# plot predicted vs. observed
ax.min <- min(c(y.valid, pred.valid))
ax.max <- max(c(y.valid, pred.valid))
plot(pred.valid ~ y.valid, xlim = c(ax.min, ax.max), ylim = c(ax.min, ax.max),
     ylab = "Predicted", xlab = "Observed")
grid()
lines(c(-100,100), c(-100,100))
abline(lm(pred.valid ~ y.valid), lwd = 2, col = "red")

# get quantiles
qs <- seq(0,1,0.001)
q.obs <- quantile(y.valid, probs = qs)
q.pred <- quantile(pred.valid, probs = qs)

# plot them out
plot(q.obs ~ qs, type = "l", col = "blue", ylim = c(ax.min, ax.max),
     xlab = "Quantiles", ylab = "Predicted/Observed Values")
lines(q.pred ~ qs, col = "red")
legend("topleft", legend = c("Observed", "Predicted"),
       col = c("blue", "red"), lty = 1)

# get quantile differences
q.diff <- q.obs - q.pred

# use differences to correct for bias in test data
pred.valid.corr <- pred.valid + approx(q.pred, q.diff, pred.valid)$y

# plot the data out
par(mfrow = c(1,2))
plot(pred.valid ~ y.valid, xlim = c(ax.min, ax.max), ylim = c(ax.min, ax.max))
grid()
lines(c(-100,100), c(-100,100))
abline(lm(pred.valid ~ y.valid), lwd = 2, col = "red")
plot(pred.valid.corr ~ y.valid, xlim = c(ax.min, ax.max), ylim = c(ax.min, ax.max))
grid()
lines(c(-100,100), c(-100,100))
abline(lm(pred.valid.corr ~ y.valid), lwd = 2, col = "red")

# apply bias correction to test data
pred.test <- predict(mod, list(x.train = x.test))
pred.test.corr <- pred.test + approx(q.pred, q.diff, pred.test)$y
par(mfrow = c(1,2))
plot(pred.test ~ y.test, xlim = c(ax.min, ax.max), ylim = c(ax.min, ax.max))
grid()
lines(c(-100,100), c(-100,100))
abline(lm(pred.test ~ y.test), lwd = 2, col = "red")
plot(pred.test.corr ~ y.test, xlim = c(ax.min, ax.max), ylim = c(ax.min, ax.max))
grid()
lines(c(-100,100), c(-100,100))
abline(lm(pred.test.corr ~ y.test), lwd = 2, col = "red")