How to put multiple plots with inner plots in a M x N layout without messing up the par()s?

I make a plot within a plot.

(I use curve() here for simplicity, but this also refers to plot().)

curve(exp(x), 0, 1, col=4)
.op <- par(
  fig=c(grconvertX(c(.05, .4), to='ndc'),
        grconvertY(c(2, 2.75), to='ndc')),
  mar=c(1, 1, 1, 1),
  new=TRUE
)
curve(sin(x), 0, 2*pi)
par(.op)

Repeating this in a for loop basically works well. However, when I try use a 2x2 layout, the figures are spread over four different plots instead of being shown on one.

layout(matrix(1:4, 2, 2))
for (i in 1:4) {
  curve(exp(x), 0, 1, col=4)
  .op <- par(
    fig=c(grconvertX(c(.05, .3), to='ndc'),
          grconvertY(c(2, 2.75), to='ndc')),
    mar=c(1, 1, 1, 1),
    new=TRUE
  )
  curve(sin(x), 0, 2*pi)
  par(.op)
}

For some reason the inner par messes up the outer one. I also tried to op <- par(mfrow=); ... par(op) instead of layout, but with same result.

How can I fix that?

Try using screen instead.

split.screen(c(2,2))
for (i in 1:4) {
  screen(i)
  curve(exp(x), 0, 1, col=4)
  title(paste("i=",i))

  .op <- par(
    fig=c(grconvertX(c(.05, .3), to='ndc'),
          grconvertY(c(2, 2.75), to='ndc')),
    mar=c(1, 1, 1, 1),
    new=TRUE
  )
  curve(sin(x), 0, 2*pi)
  par(.op)
}

If you want to stick to using layout (sort of), you can first create a series of empty plots to find the coordinates for the layout regions. Then, not relying on layout anymore, you can define them again with par(fig=...) which will allow you to draw the inner plots without things getting messed up.

layout(matrix(1:4, 2, 2))
# store coordinates of layout regions
figs <- sapply(1:4, \(i) {
  frame()
  par('fig')
})
# to overwrite the empty plots
par(new=T)

for (i in 1:4) {
  
  # set figure region to the current layout window
  par(fig = figs[, i])
  
  curve(exp(x), 0, 1, col=4)
  .op <- par(
    fig=c(grconvertX(c(.05, .3), to='ndc'),
          grconvertY(c(2, 2.75), to='ndc')),
    mar=c(1, 1, 1, 1),
    new=TRUE
  )
  curve(sin(x), 0, 2*pi)
  par(.op)
  par(new=T)
}

Here's how I combined the insights from the two answers in a function plot_fun(). From @Edward I take split.screen, and from @Robert Hacken the fact that curve (or plot respectively) produces margins relative to the font size that we may manipulate using cex.

plot_fun <- \(spl, .cex=3, iadj=0) {
  scr <- split.screen(spl)
  for (i in scr) {
    screen(i)
    par(mar=c(5, 6, 4, 3) + .1)
    curve(exp(x), 0, 1, col=4, cex.axis=.cex, cex.lab=.cex)
    .op <- par(
      fig=c(grconvertX(c(.1, .4), to='ndc'),
            grconvertY(c(2 + iadj, 2.75), to='ndc')),
      mar=c(1, 1, 1, 1),
      cex=.cex*.1,
      new=TRUE
    )
    curve(sin(x), 0, 2*pi, cex.axis=.cex*3)
    par(.op)
  }
  close.screen(all.screens=TRUE)
}
    

f <- 4  ## appears to work fine with the two versions below

pdf("tmp/p1.pdf", 16*f, 9*f)
plot_fun(spl=c(6, 7))
dev.off()

pdf("tmp/p1.pdf", 16*f, 9*f)
plot_fun(spl=c(2, 3), iadj=.2)
dev.off()

If we wish cleaner tick labels on the axes for publishing or so, we may still use curve(., axes=FALSE) and mess around with axis() and mtext(), but this already gives an acceptable result for analyses without great effort.

If your device gets messed up during setup and gives weird results, use graphics.off() and close.screen(all.screens=TRUE) to reset the device.