Multiple plots on one device using stl

If you 

?plot.stl

you will see that that the second argument, set.pars, is a list of argument settings for par(), including a (variable) default setting for mfrow.  I.e. plot.stl overrides your external setting (which will also override any layout() setting).  

It looks like to override it back, you may need to take charge yourself.  Here is the gist of a partial way round it, perhaps, sort of. ...

################
par(mar = c(0, 6, 0, 6), oma = c(6, 0, 4, 0),
  tck = -0.01, mfcol = c(4, 2))

plot(stl.1, set.pars = list())
plot(stl.2, set.pars = list())
################

Things might be a bit more flexible if you were to use xyplot.stl in the latticeExtra package.  With lattice the operations of making the plot and displaying it are more cleanly separated.

I am looking at monthly reports, and have three series of monthly data from 2007 to 2009. I would like to show the season decomposition of these two series side by side on the one device, however using plot doesn't seem to respect any use of layout(matrix(1:3, ncol=3)) or par(mfcol=c(1,3)).

I'm guessing that this means that the plot(stl) perhaps uses them, but I can't find anywhere the / a plot.stl() - ie, I can't work out where the plot() call is going? 

I've attached a small example of data and some probably overly verbose code. There are only two sets of example data, not three as mentioned above.

# load the data
vals.1 <-
structure(list(year = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("2007", 
"2008", "2009"), class = c("ordered", "factor")), month = structure(c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 
8L, 9L, 10L, 11L, 12L), .Label = c("Jan", "Feb", "Mar", "Apr", 
"May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"), class = c("ordered", 
"factor")), count = c(105, 100, 64, 44, 49, 65, 88, 90, 99, 92, 
93, 88, 212, 146, 96, 131, 220, 143, 137, 138, 395, 362, 349, 
222, 294, 268, 298, 310, 426, 348, 287, 101, 66, 112, 105, 4)), .Names = c("year", 
"month", "count"), row.names = c("1", "2", "3", "4", "5", "6", 
"7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", 
"18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", 
"29", "30", "31", "32", "33", "34", "35", "36"), class = "data.frame")
vals.2 <-
structure(list(year = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("2007", 
"2008", "2009"), class = c("ordered", "factor")), month = structure(c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 
5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 
8L, 9L, 10L, 11L), .Label = c("Jan", "Feb", "Mar", "Apr", "May", 
"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"), class = c("ordered", 
"factor")), count = c(45, 34, 17, 6, 7, 16, 12, 11, 14, 17, 11, 
20, 27, 12, 10, 14, 22, 23, 92, 144, 385, 274, 320, 252, 240, 
146, 222, 142, 122, 117, 163, 89, 51, 89, 108)), .Names = c("year", 
"month", "count"), row.names = c("1", "2", "3", "4", "5", "6", 
"7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", 
"18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", 
"29", "30", "31", "32", "33", "34", "35"), class = "data.frame")

# create the time series
series.1 <- ts(vals.1$count, start=c(2007,1), freq=12)
series.2 <- ts(vals.2$count, start=c(2007,1), freq=12)

# apply the seasonal decomposition
stl.1 <- stl(series.1, "per", robust=TRUE)
stl.2 <- stl(series.2, "per", robust=TRUE)

# set up the device for side by side display
layout(matrix(1:2, ncol=2))

# little check
layout.show(2)

# plot the first series and second series
plot(stl.1, labels = c('Count by month','Seasonal Component','Trend','Remainder'), main='Series.1 Decomposition')
plot(stl.2, labels = c('Count by month','Seasonal Component','Trend','Remainder'), main='Series.2 Decomposition')

# hmm, used whole device twice, try again

par(mfcol=c(1,2))

# and now the second

plot(stl.1, labels = c('Count by month','Seasonal Component','Trend','Remainder'), main='Series.1 Decomposition')
plot(stl.2, labels = c('Count by month','Seasonal Component','Trend','Remainder'), main='Series.2 Decomposition')

# oh what about split.screen()

split.screen(c(1,2))
screen(1)
# now plot
plot(stl.1, labels = c('Count by month','Seasonal Component','Trend','Remainder'), main='Series.1 Decomposition') 

# something wrong with the plot, not seeing original series at the top.
screen(2)

# Error in par(split.screens[[n]]) : parameter "j" in "mfg" is out of range