Now the challenge: I'd like to get fitted outputs using newdata, where I've incremented wt by cyl/4:

newdata <-
mtcars2 %>%
mutate(
wt = wt + cyl/4)

I expect that this will produce a data frame of the same size as lm1 %>% augment(fit): one row for each row in newdata, because broom will match up models and newdata by the grouping variables cyl and am.

From reading elsewhere, I gather that group_by and rowwise data frames aren't compatible, so lm1 is ungrouped, and augment can't associate models and newdata. Is there another design pattern that lets me do this? It would be nice if it were as simple and transparent as the above attempt, but it's more important that it work.

@aosmith: I have been exploring your second option, and I like it. When I try it on my real data, though, I have a problem in the mutate command: it returns "Error: augment doesn't know how to deal with data of class list".

Where I say it looks like yours, I mean I have the following columns (renamed here for consistency): ID (chr), attr1 (dbl), cyl (dbl), am (chr), fit (list), and data (list). You have cyl, am (dbl), fit, and data. I changed my am to dbl, but that didn't help.

I think the difference is that I have 3 (ID ... similar to the rownames in mtcars) x 2 (cyl) x 2 (am) units in this sample (with each sample having 12 measurements), while the mtcars example has 3 (cyl) x 2 (am) cells x a random number of car types per cell. In my analysis, I need to see the ID values, but newdata applies equally to all units. If it helps, think of it as the speed of a headwind applied to each car in the test. Does that suggest a cause for augment's complaint it can't deal with data of class list?

EDIT: Merging the ID with the newdata (using full=TRUE) solved the last problem. I'm currently using your first proposed solution.

I've used map2 from package purrr for this sort of situation. map2 loops through the elements of two lists simultaneously. The lists must be the same length and be in the same order.

The elements of the lists are used as arguments for some function you want to apply (augment, in your case). Here your two lists would be a list of models and a list of datasets (one list for each cyl/am combination).

Using map2_df returns the results as a data.frame instead of a list.

library(purrr)

I made the list of data.frames to predict with using split. The order of the factors to split on determined the list order, so I made sure it was in the same order as lm1.