Determine how variables are grouped in subplots

[itsdfish]

As described in this discourse question, when working with tabular data, it is common to plot a dependent variable according to groups and organize subplots so that each subplot has a fixed value of one of the factors. Currently, it does not appear to be possible to specify how to group factors in the subplots. The code below demonstrates this limitation:

using StatsPlots, DataFrames

linear(x, β0, β1) = β0 + β1 * x

exponential(x, β, λ) = β * exp(λ * x)

xs = 1:10

np1 = [(p1,p2,x,y = linear(x, p1, p2)) for p1 in [0.0,1.0] for p2 in [.2,.6] for x in xs]

df = DataFrame(np1)
df[!,:function] .= "linear"

np2 = [(p1,p2,x,y = exponential(x, p1, p2)) for p1 in [.0002,.0003] for p2 in [1.0,1.1] for x in xs]
temp =  DataFrame(np2)
temp[!,:function] .= "exponential"

append!(df, temp)

@df df plot(:x, :y, group=(:p1,:p2,:function), layout=(2,1), ylims=(0,20))

A proposed workaround is to group the dataframe according to the factor that varies across the subplots:

gdf = groupby(df, :function)
fig = plot(layout=(2,1), legend=:topleft)
@df gdf[1] plot!(fig[1],:x, :y, group=(:function, :p1, :p2), ylims=(0,20))
@df gdf[2] plot!(fig[2],:x, :y, group=(:function, :p1, :p2), ylims=(0,20))

Although this works, it places a lot of demands on the user when generalizing to more factors or different layout schemes. I was wondering if it would be possible to elegantly handle this within StatsPlots? For example grouped_layout = (:function, 1) or grouped_layout = (:factor1, :factor2).