[ ] display single legend for all subplots: could be in an empty subplot (e.g. if grid is (nrow, ncol) and there are fewer than nrox*ncol variables, the empty plot in the bottom right could display only the legend) or could be outside the grid of plots.
[ ] support optional keyword arguments nrow and ncol. If both are given, that'd be the same as layout = (nrow, ncol). We need the case of only one of them given and the other computed from the number of variables, with some subplots potentially left empty.