The par() command before making the first plot is to reset the graphics device. The inset will need tiny or no margins, and once you create the smaller graphics device all plotting will be done in it until you remake everything or create a new device.

Not a map, but:

par(fig=c(0,1,0,1))
hist(rnorm(100))
par(mar=c(0,0,0,0),
fig=c(0.8,0.99,0.8,0.99), new=TRUE)
hist(runif(10))

par()$fig will create a new graphics window inside your current graphics window, with its extent defined where 0 is the left or bottom edge, and 1 is the right or top edge

That's probably what I'll end up doing, or just use an older installation of R that still has it. I was just wondering if there was a newer implementation somewhere for the convenience.

Hard way was lm, then anova on the result, then a separate call to aov with the same formula to feed to TukeyHSD. Turns out you can use the output of aov in summary AND TukeyHSD and get everything you need. I feel silly for not realizing this sooner. This will save me dozens of minutes from now on.

Nice post! Even as a base R user I've never used subset or transform much but might have to try. People don't realize how powerful base aggregate can be.