quick and ugly answer, but it does the job until someone delivers a better one ...
ggplot(df, aes(x, y, fill = color)) +
geom_point(size = 4, pch = 21) +
guides(fill = guide_legend(
title = expression(atop(Median~Nitrate-Nitrogen~(NO[3]^{textstyle("-")}-N), "Concentration"~(mg~L^{textstyle("-")})~phantom (1000000)~phantom (1000000)))))
source
share