How to make one qq-plot with two variables

First, we create a data set.

x <- c(33.5, 33.8, 31.4, 30.2, 30, 27.1, 24.1, 21.5, 9.2)
y <- c(51.7, 48.1, 47.2, 46.9, 46.7, 40.9, 39.9, 38.2, 30.9)

Next, we will create a qq-plot for x. We supply arguments for xlim and ylim based on the min and max of the output from the qqplot function for each variable, otherwise the second qq-plot might not fit using the automatic axis scales for the first variable.

qqnorm(x, 
  xlim = c(
    min(qqnorm(x, plot.it = FALSE)$x, qqnorm(y, plot.it = FALSE)$x), 
    max(qqnorm(x, plot.it = FALSE)$x, qqnorm(y, plot.it = FALSE)$x)
    ),
  ylim = c(
    min(qqnorm(x, plot.it = FALSE)$y, qqnorm(y, plot.it = FALSE)$y), 
    max(qqnorm(x, plot.it = FALSE)$y, qqnorm(y, plot.it = FALSE)$y)
    ),
  pch = 16,
  col = 'brown',
  bty = 'l'
)

The function qqline adds the theoretically expected line for comparison.

qqnorm(x, 
  xlim = c(
    min(qqnorm(x, plot.it = FALSE)$x, qqnorm(y, plot.it = FALSE)$x), 
    max(qqnorm(x, plot.it = FALSE)$x, qqnorm(y, plot.it = FALSE)$x)
    ),
  ylim = c(
    min(qqnorm(x, plot.it = FALSE)$y, qqnorm(y, plot.it = FALSE)$y), 
    max(qqnorm(x, plot.it = FALSE)$y, qqnorm(y, plot.it = FALSE)$y)
    ),
  pch = 16,
  col = 'brown',
  bty = 'l'
)
qqline(x, col = 'brown')

The qq-plot for the second variable is added using the points function, with the result of qqnorm as the first argument. HERE IS THE TRICK IN USING points: you must specify plot.it = FALSE in the qqnorm function so that qqnorm just supplies points instead of it attempting to create a new plot.

qqnorm(x, 
  xlim = c(
    min(qqnorm(x, plot.it = FALSE)$x, qqnorm(y, plot.it = FALSE)$x), 
    max(qqnorm(x, plot.it = FALSE)$x, qqnorm(y, plot.it = FALSE)$x)
    ),
  ylim = c(
    min(qqnorm(x, plot.it = FALSE)$y, qqnorm(y, plot.it = FALSE)$y), 
    max(qqnorm(x, plot.it = FALSE)$y, qqnorm(y, plot.it = FALSE)$y)
    ),
  pch = 16,
  col = 'brown',
  bty = 'l'
)
qqline(x, col = 'brown')

points(qqnorm(y, plot.it = FALSE), pch = 16, col = 'navy')
qqline(y, col = 'navy')

legend('bottomright', legend = c('x', 'y'), fill = c('brown', 'navy'), bty = 'n')