diff --git a/R/by_aesthetics.R b/R/by_aesthetics.R
index dd76f2ab..6562c4a2 100755
--- a/R/by_aesthetics.R
+++ b/R/by_aesthetics.R
@@ -1,14 +1,18 @@
-by_col = function(ngrps = 1L, col = NULL, palette = NULL, gradient = NULL, ordered = NULL, alpha = NULL) {
+by_col = function(ngrps = 1L, col = NULL, palette = NULL, gradient = NULL, ordered = NULL, alpha = NULL, nested = FALSE) {
+  gradient = isTRUE(gradient)
+  col_index = FALSE # flag for subsetting/index over colors
   if (is.null(ordered)) ordered = FALSE
   if (is.null(alpha)) alpha = 1
   if (is.null(gradient)) gradient = FALSE
-  if (isTRUE(gradient)) {
-    ngrps = 100L
-  }
+  if (gradient) ngrps = 100L
+  # if (gradient && nested) col = round(rescale_num(col, to = c(1,100)))
+  if (gradient && nested) col = NULL
+  
+  ncomp = if (gradient || !nested) ngrps else length(unique(col)) 
 
   if (is.null(palette)) {
     pal_qual = get_tpar("palette.qualitative", default = NULL)
-    if (ngrps <= max(c(length(pal_qual), 8))) {
+    if (ncomp <= max(c(length(pal_qual), 8))) {
       palette = pal_qual
     } else {
       palette = get_tpar("palette.sequential", default = NULL)
@@ -37,10 +41,15 @@ by_col = function(ngrps = 1L, col = NULL, palette = NULL, gradient = NULL, order
     }
     if (isTRUE(gradient)) {
       col = rev(col)
-    } else if (!ordered && is.numeric(col)) {
-      col = palette()[col]
+    # } else if (!ordered && is.numeric(col)) {
+    #   col = palette()[col]
     }
-    if (anyNA(col) || is.character(col)) {
+    # if (anyNA(col) || is.character(col)) {
+    if (!gradient && nested && !ordered && is.numeric(col)) { ## double check this
+      if (is.null(palette)) palette = grDevices::palette
+      cidx = col
+      col_index = TRUE
+    } else if (anyNA(col) || is.character(col)) {
       if (alpha) col = adjustcolor(col, alpha.f = alpha)
       return(col)
     }
@@ -157,12 +166,18 @@ by_col = function(ngrps = 1L, col = NULL, palette = NULL, gradient = NULL, order
       )
     }
   }
+  
+  # if (col_index && !is.null(args$n)) args$n = max(args$n, cidx)
+  if (col_index && !is.null(args$n)) args$n = max(cidx)
+  
 
   cols = tryCatch(
     do.call(palette_fun, args),
     error = function(e) do.call(eval(palette), args) # catch for bespoke palette generating funcs
   )
 
+  if (col_index) cols = cols[cidx]
+  
   if (length(cols) > ngrps) cols = cols[1:ngrps]
 
   # For gradient and ordered colors, we'll run high to low
@@ -307,7 +322,8 @@ by_bg = function(
     ngrps,
     palette,
     ribbon.alpha,
-    type) {
+    type,
+    nested = FALSE) {
   if (is.null(bg) && !is.null(fill)) bg = fill
   if (!is.null(bg) && length(bg) == 1 && is.numeric(bg) && bg >= 0 && bg <= 1) {
     alpha = bg
@@ -320,7 +336,8 @@ by_bg = function(
       palette = palette,
       gradient = by_continuous,
       ordered = by_ordered,
-      alpha = alpha
+      alpha = alpha,
+      nested = nested
     )
   } else if (length(bg) != ngrps) {
     bg = rep(bg, ngrps)
diff --git a/R/tinyformula.R b/R/tinyformula.R
index 56fd546e..4475aed7 100644
--- a/R/tinyformula.R
+++ b/R/tinyformula.R
@@ -1,9 +1,10 @@
 ## auxiliary functions for formula/facet parsing
 
-tinyformula = function(formula, facet = NULL) {
+tinyformula = function(formula, facet = NULL, col_fml = NULL) {
   ## input
   ## - formula:       y ~ x or y ~ x | z or ~ x or ~ x | z
   ## - facet:         ~ a or ~ a + b or b ~ a
+  ## - col:           ~ c
   ##
   ## output:
   ## - x:      ~ x
@@ -11,7 +12,8 @@ tinyformula = function(formula, facet = NULL) {
   ## - by:     NULL or ~ z or ~ z1 + z2 + ... (use interaction of all)
   ## - xfacet: NULL or ~ a or ~ a + b etc.
   ## - yfacet: NULL or ~ b
-  ## - full:   e.g. ~ x + y + z + a + b
+  ## - col:    NULL or ~ c
+  ## - full:   e.g. ~ x + y + z + a + b + c
 
   ## preliminaries
   if (!inherits(formula, "formula")) formula = as.formula(formula)
@@ -28,6 +30,11 @@ tinyformula = function(formula, facet = NULL) {
     xfacet = ~ a
     yfacet = if (length(facet) == 2L) NULL else ~ b
   }
+  if (is.null(col_fml) || !inherits(col_fml, "formula")) {
+    col = NULL
+  } else {
+    col = ~ c
+  }
 
   ## fill with actual terms
   environment(x) = environment(formula)
@@ -50,6 +57,10 @@ tinyformula = function(formula, facet = NULL) {
     environment(yfacet) = environment(formula)
     yfacet[[2L]] = facet[[2L]]
   }
+  if (!is.null(col)) {
+    environment(col) = environment(formula)
+    col[[2L]] = col_fml[[2L]]
+  }
 
   ## combine everything
   full = x
@@ -57,6 +68,7 @@ tinyformula = function(formula, facet = NULL) {
   if (!is.null(by))     full[[2L]] = call("+", full[[2L]], by[[2L]])
   if (!is.null(xfacet)) full[[2L]] = call("+", full[[2L]], xfacet[[2L]])
   if (!is.null(yfacet)) full[[2L]] = call("+", full[[2L]], yfacet[[2L]])
+  if (!is.null(col))    full[[2L]] = call("+", full[[2L]], col[[2L]])
 
   ## return list of all formulas
   return(list(
@@ -65,6 +77,7 @@ tinyformula = function(formula, facet = NULL) {
     by = by,
     xfacet = xfacet,
     yfacet = yfacet,
+    col = col,
     full = full
   ))
 }
diff --git a/R/tinyplot.R b/R/tinyplot.R
index 2ac3727e..789f0c2b 100644
--- a/R/tinyplot.R
+++ b/R/tinyplot.R
@@ -847,9 +847,30 @@ tinyplot.default = function(
   # split data
   by_ordered = FALSE
   by_continuous = !is.null(by) && inherits(datapoints$by, c("numeric", "integer"))
+  
+  nested = isTRUE(dots$nested)
+  if (nested) {
+    by_nested = col
+    if (inherits(by_nested, c("numeric", "integer"))) {
+      ngrps_nested = 100
+      by_continuous = TRUE
+      datapoints$col = col
+      
+    } else {
+      col = factor(tapply(col, by, FUN = `[[`, 1))  # grab group colours
+      nested_labs = levels(col)
+      col = as.integer(col)
+      ngrps_nested = length(col)
+    }
+    # need to solve by_continuous flag below...
+  }
   if (isTRUE(by_continuous) && type %in% c("l", "b", "o", "ribbon", "polygon", "polypath", "boxplot")) {
-    warning("\nContinuous legends not supported for this plot type. Reverting to discrete legend.")
-    by_continuous = FALSE
+    # warning("\nContinuous legends not supported for this plot type. Reverting to discrete legend.")
+    # by_continuous = FALSE
+    if (!nested) {  #TEMPORARY!!!
+      warning("\nContinuous legends not supported for this plot type. Reverting to discrete legend.")
+      by_continuous = FALSE
+    }
   } else if (!is.null(by)) {
     by_ordered = is.ordered(by)
   }
@@ -866,16 +887,17 @@ tinyplot.default = function(
 
   # aesthetics by group: col, bg, etc.
   ngrps = if (is.null(by)) 1L else if (is.factor(by)) length(levels(by)) else if (by_continuous) 100L else length(unique(by))
+  ngrps_col = if (!nested) ngrps else ngrps_nested
   pch = by_pch(ngrps = ngrps, type = type, pch = pch)
   lty = by_lty(ngrps = ngrps, type = type, lty = lty)
   lwd = by_lwd(ngrps = ngrps, type = type, lwd = lwd)
   col = by_col(
-    ngrps = ngrps, col = col, palette = palette,
-    gradient = by_continuous, ordered = by_ordered, alpha = alpha)
+    ngrps = ngrps_col, col = col, palette = palette,
+    gradient = by_continuous, ordered = by_ordered, alpha = alpha, nested = nested)
   bg = by_bg(
     adjustcolor = adjustcolor, alpha = alpha, bg = bg, by = by, by_continuous = by_continuous,
     by_ordered = by_ordered, col = col, fill = fill, palette = substitute(palette),
-    ribbon.alpha = ribbon.alpha, ngrps = ngrps, type = type)
+    ribbon.alpha = ribbon.alpha, ngrps = ngrps_col, type = type, nested = nested)
   
   ncolors = length(col)
   lgnd_labs = rep(NA, times = ncolors)
@@ -883,7 +905,7 @@ tinyplot.default = function(
     ## Identify the pretty break points for our labels
     nlabs = 5
     ncolors = length(col)
-    ubyvar = unique(by)
+    ubyvar = if (!nested) unique(by) else(unique(by_nested))
     byvar_range = range(ubyvar)
     pbyvar = pretty(byvar_range, n = nlabs)
     pbyvar = pbyvar[pbyvar >= byvar_range[1] & pbyvar <= byvar_range[2]]
@@ -929,7 +951,7 @@ tinyplot.default = function(
     legend_args[["x"]] = "none"
   }
 
-  if (is.null(by)) {
+  if (is.null(by) && !nested) {
     if (is.null(legend)) {
       legend = "none"
       legend_args[["x"]] = "none"
@@ -938,7 +960,9 @@ tinyplot.default = function(
 
   if ((is.null(legend) || legend != "none") && isFALSE(add)) {
     if (isFALSE(by_continuous)) {
-      if (ngrps > 1) {
+      if (nested) {
+        lgnd_labs = nested_labs
+      } else if (ngrps > 1) {
         lgnd_labs = if (is.factor(datapoints$by)) levels(datapoints$by) else unique(datapoints$by)
       } else {
         lgnd_labs = ylab
@@ -961,7 +985,7 @@ tinyplot.default = function(
       pch = pch,
       lty = lty,
       lwd = lwd,
-      col = col,
+      col = if (nested && !by_continuous) unique(col) else col,
       bg = bg,
       gradient = by_continuous,
       cex = cex * cex_fct_adj,
@@ -1153,9 +1177,17 @@ tinyplot.default = function(
     iby = idata[["by"]]
     if (!is.null(by)) { ## maybe all(iby=="")
       if (isTRUE(by_continuous)) {
-        idata[["col"]] = col[round(rescale_num(idata$by, from = range(datapoints$by), to = c(1, 100)))]
-        idata[["bg"]] = bg[round(rescale_num(idata$by, from = range(datapoints$by), to = c(1, 100)))]
-        idata = list(idata)
+        if (!nested) {
+          idata[["col"]] = col[round(rescale_num(idata$by, from = range(datapoints$by), to = c(1, 100)))]
+          idata[["bg"]] = bg[round(rescale_num(idata$by, from = range(datapoints$by), to = c(1, 100)))]
+          idata = list(idata)
+        } else {
+          idata[["bg"]] = bg[round(rescale_num(idata$col, from = range(datapoints$col), to = c(1, 100)))]
+          idata[["col"]] = col[round(rescale_num(idata$col, from = range(datapoints$col), to = c(1, 100)))]
+          idata = lapply(idata, split, iby)
+          idata = do.call(function(...) Map("list", ...), idata)
+        }
+        
       } else {
         idata = lapply(idata, split, iby)
         idata = do.call(function(...) Map("list", ...), idata)
@@ -1324,9 +1356,15 @@ tinyplot.formula = function(
 
   ## placeholder for legend title
   legend_args = list(x = NULL)
+  
+  ## catch of col passed as formula (e.g., for nested colors)
+  col_fml = NULL
+  if (!is.null(col) && inherits(col, "formula")) {
+    col_fml = col
+  }
 
   ## process all formulas
-  tf = tinyformula(formula, facet)
+  tf = tinyformula(formula, facet, col_fml)
 
   ## set up model frame
   m = match.call(expand.dots = FALSE)
@@ -1374,6 +1412,16 @@ tinyplot.formula = function(
       attr(facet, "facet_nrow") = length(unique(yfacet))
     }
   }
+  
+  ## extract col (if any)
+  nested = FALSE
+  if (!is.null(col_fml)) {
+    col = tinyframe(tf$col, mf)
+    colnam = names(col)
+    col = if (length(colnam) == 1L) col[[colnam]] else interaction(col, sep = ":")
+    # col = factor(tapply(col, by, FUN = `[[`, 1))  # grab group colours
+    nested = TRUE
+  }
 
   ## nice axis and legend labels
   dens_type = (is.atomic(type) && identical(type, "density")) || (!is.atomic(type) && identical(type$name, "density"))
@@ -1391,7 +1439,9 @@ tinyplot.formula = function(
     if (is.null(ylab)) ylab = ynam
     if (is.null(xlab)) xlab = xnam
   }
-  if (!is.null(by)) {
+  if (nested) {
+    legend_args[["title"]] = if (length(colnam) == 1L) colnam else sprintf("interaction(%s)", paste(colnam, collapse = ", "))
+  } else if (!is.null(by)) {
     legend_args[["title"]] = if (length(bynam) == 1L) bynam else sprintf("interaction(%s)", paste(bynam, collapse = ", "))
   }
 
@@ -1418,6 +1468,7 @@ tinyplot.formula = function(
     lty = lty,
     lwd = lwd,
     restore.par = restore.par,
+    nested = nested,
     ...
   )
 }