How do I add multiple subplots into a multirow figure in R?

Question

i need to overlay multiple subplots onto a single plot which is already contained inside a multirow figure (see image)
the reason why i need subplots instead of screen layout is because the figure will be possibly multicolumn, also (a 5 by 3 plot, in total)
there are packages which assist in doing subplots, but they break when you use multirow figures, and sequential subplots, except the first one, are rendered next to the overall figure border, not relative to the current row/column plot borders
i understand large packages such as ggplot2 allow this relatively easily, but base R plots are highly preferable

UPD: the minimum reproducible example depicting the problem is here:

require(Hmisc)

COL.1 <- c('red','orange','yellow'); COL.2 <- c('blue','green','turquoise')

SUBPLOT.FUN <- function(COL) {plot(rnorm(100), type='l', col=COL)}
PLOT.FUN <- function(i) {
  plot(rnorm(100),ylim=c(-1,1))
  subplot(SUBPLOT.FUN(COL.1[i]), 100,1, vadj=1,hadj=1,pars=list(mfg=c(1,i)))
  subplot(SUBPLOT.FUN(COL.2[i]), 100,-1,vadj=0,hadj=1,pars=list(mfg=c(1,i)))
}
plot.new(); par(mfrow=c(1,3))
for (i in 1:3) {
  PLOT.FUN(i)
}

which looks like that:
while what is required is shown on the first image (meaning, EACH of the three plots must contain 3 subplots in their respective locations (along the right border, arranged vertically))
N.B. either the figure is multirow or multicolumn (as depicted) does not matter

Answer 1

Something like this? Inspired in this R-bloggers post.

# reproducible test data
set.seed(2022)
x <- rnorm(1000)
y <- rbinom(1000, 1, 0.5)
z <- rbinom(1000, 4, 0.5)

# save default values and prepare 
# to create custom plot areas
old_par <- par(fig = c(0,1,0,1))

# set x axis limits based on data
h <- hist(x, plot = FALSE)
xlim <- c(h$breaks[1] - 0.5, h$breaks[length(h$breaks)] + 2)

hist(x, xlim = xlim)

# x = c(0.6, 1)  right part of plot
# y = c(0.5, 1)  top part of plot
par(fig = c(0.6, 1, 0.5, 1), new = TRUE)  
boxplot(x ~ y)  

# x = c(0.6, 1)   right part of plot
# y = c(0.1, 0.6) bottom part of plot
par(fig = c(0.6, 1, 0.1, 0.6), new = TRUE)  
boxplot(x ~ z)

# put default values back
par(old_par)

^{Created on 2022-08-18 by the reprex package (v2.0.1)}