op <- par()

par(mfrow = c(2, 2))
plot(rnorm(10),
     rnorm(10),
     col = gl(2, 5),
     pch = as.integer(gl(2, 5)))
plot(rnorm(10),
     rnorm(10),
     col = gl(2, 5),
     pch = as.integer(gl(2, 5)))
plot(rnorm(10),
     rnorm(10),
     col = gl(2, 5),
     pch = as.integer(gl(2, 5)))

par(xaxt = "n", yaxt = "n", bty = "n")
plot(1, xlab = "", ylab = "")
legend(
  "left",
  c("Group A", "Group B"),
  pch = 1:2,
  col = 1:2,
  ncol = 2
)

## At end of plotting, reset to previous settings:
par(op)