My goal is to show multiple plots of residual vs. fitted values for a dataframe of 13 stations each having a separate lm() model. Working from an example in Chapter 25 (Many Models) in the R for Data Science (1st ed) textbook.