在单个页面中按变量绘制分位数回归

Plotting quantile regression by variables in a single page

我是运行分别对几个自变量(同因变量)进行分位数回归。我只想在单个图中绘制每个变量的几个分位数的斜率估计值。

这是一个玩具数据:

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

我只想绘制分位数的斜率估计值。因此,我在 plot.

中给出 parm=2
 plot(fit1, parm=2)
 plot(fit2, parm=2)

现在,我想将这两个图合并到一个页面中。

到目前为止我尝试了什么;

  1. 我尝试设置 par(mfrow=c(2,2)) 并绘制它们。但它正在生成一个空白页。
  2. 我曾尝试使用 gridExtra 和 gridGraphics,但没有成功。尝试将基本图转换为 Grob 对象,如所述
  3. 尝试使用函数 layout 中的函数 this 文档
  4. 我正在尝试查看 plot.rqs 的源代码。但我无法理解它是如何绘制置信带的(我只能绘制分位数上的系数)或在那里更改 mfrow 参数。

谁能指出我哪里错了?我应该查看 plot.rqs 的源代码并更改其中的任何参数吗?

quantreg 包使用的函数 plot 有它自己的 mfrow 参数。如果你不指定它,它会强制执行它自己选择的一些选项(因此会覆盖你的 par(mfrow = c(2,2)).

plot.rqs 中使用 mfrow 参数:

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

虽然 quantreg::plot.summary.rqs 有一个 mfrow 参数,但它使用它来覆盖 par('mfrow') 以便分面 parm 值,这不是你想要做的.

一种替代方法是解析对象并手动绘制。您可以从 fit1fit2 中提取 tau 值和系数矩阵,它们只是每个 tau 的值列表,所以在 tidyverse 语法中,

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`

如果你喜欢,可以从对象中拉出更多,添加原始的水平线,否则就疯狂。


另一种选择是使用magick捕捉原始图像(或用任何设备保存并重新读取)并手动组合它们:

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