Plotting quantile regression by variables in a single page

Question

I am running quantile regressions for several independent variables separately (same dependent). I want to plot only the slope estimates over several quantiles of each variable in a single plot.

Here's a toy data:

set.seed(1988)

y <- rnorm(50, 5, 3)
x1 <- rnorm(50, 3, 1)
x2 <- rnorm(50, 1, 0.5)

# Running Quantile Regression 
require(quantreg)
fit1 <- summary(rq(y~x1, tau=1:9/10), se="boot")
fit2 <- summary(rq(y~x2, tau=1:9/10), se="boot")

I want to plot only the slope estimates over quantiles. Hence, I am giving parm=2 in plot.

 plot(fit1, parm=2)
 plot(fit2, parm=2)

Now, I want to combine both these plots in a single page.

What I have tried so far;

1. I tried setting par(mfrow=c(2,2)) and plotting them. But it's producing a blank page.
2. I have tried using gridExtra and gridGraphics without success. Tried to convert base graphs into Grob objects as stated here
3. Tried using function layout function as in this document
4. I am trying to look into the source code of plot.rqs. But I am unable to understand how it's plotting confidence bands (I'm able to plot only the coefficients over quantiles) or to change mfrow parameter there.

Can anybody point out where am I going wrong? Should I look into the source code of plot.rqs and change any parameters there?

Answer 1

While quantreg::plot.summary.rqs has an mfrow parameter, it uses it to override par('mfrow') so as to facet over parm values, which is not what you want to do.

One alternative is to parse the objects and plot manually. You can pull the tau values and coefficient matrix out of fit1 and fit2, which are just lists of values for each tau, so in tidyverse grammar,

library(tidyverse)

c(fit1, fit2) %>%    # concatenate lists, flattening to one level
    # iterate over list and rbind to data.frame
    map_dfr(~cbind(tau = .x[['tau']],    # from each list element, cbind the tau...
                   coef(.x) %>%    # ...and the coefficient matrix, 
                       data.frame(check.names = TRUE) %>%    # cleaned a little
                       rownames_to_column('term'))) %>% 
    filter(term != '(Intercept)') %>%    # drop intercept rows
    # initialize plot and map variables to aesthetics (positions)
    ggplot(aes(x = tau, y = Value, 
               ymin = Value - Std..Error, 
               ymax = Value + Std..Error)) + 
    geom_ribbon(alpha = 0.5) + 
    geom_line(color = 'blue') + 
    facet_wrap(~term, nrow = 2)    # make a plot for each value of `term`

Pull more out of the objects if you like, add the horizontal lines of the original, and otherwise go wild.

---

Another option is to use magick to capture the original images (or save them with any device and reread them) and manually combine them:

library(magick)

plots <- image_graph(height = 300)    # graphics device to capture plots in image stack
plot(fit1, parm = 2)
plot(fit2, parm = 2)
dev.off()

im1 <- image_append(plots, stack = TRUE)    # attach images in stack top to bottom

image_write(im1, 'rq.png')

Answer 2

The function plot used by quantreg package has it's own mfrow parameter. If you do not specify it, it enforces some option which it chooses on it's own (and thus overrides your par(mfrow = c(2,2)).

Using the mfrow parameter within plot.rqs:

# make one plot, change the layout
plot(fit1, parm = 2, mfrow = c(2,1))
# add a new plot
par(new = TRUE)
# create a second plot
plot(fit2, parm = 2, mfrow = c(2,1))