diff --git a/NAMESPACE b/NAMESPACE
index 1624de5..345b48d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -10,7 +10,9 @@ S3method(as.igraph,"MultiFactor::MultiFactor")
 S3method(as.matrix,"MultiFactor::LinkMap")
 S3method(augment,"MultiFactor::MultiFactor")
 S3method(c,"MultiFactor::MultiFactor")
+S3method(stack,"MultiFactor::MultiFactor")
 S3method(subset,"MultiFactor::MultiFactor")
+S3method(unique,"MultiFactor::LinkMap")
 export(LinkMap)
 export(MultiFactor)
 export(as.LinkMap)
@@ -33,6 +35,7 @@ importFrom(Matrix,rowSums)
 importFrom(Matrix,sparseMatrix)
 importFrom(Matrix,t)
 importFrom(Matrix,which)
+importFrom(S7,"prop<-")
 importFrom(S7,S7_dispatch)
 importFrom(S7,new_class)
 importFrom(S7,new_property)
@@ -54,6 +57,7 @@ importFrom(igraph,sample_gnm)
 importFrom(igraph,shortest_paths)
 importFrom(rlang,as_label)
 importFrom(rlang,check_dots_empty0)
+importFrom(rlang,dots_list)
 importFrom(rlang,f_lhs)
 importFrom(rlang,f_rhs)
 importFrom(rlang,is_missing)
@@ -62,4 +66,5 @@ importFrom(stats,reformulate)
 importFrom(utils,count.fields)
 importFrom(utils,data)
 importFrom(utils,download.file)
+importFrom(utils,stack)
 importMethodsFrom(Matrix,"%&%")
diff --git a/R/AllClasses.R b/R/AllClasses.R
index f90a375..34bee20 100644
--- a/R/AllClasses.R
+++ b/R/AllClasses.R
@@ -48,10 +48,13 @@ LinkMap <- S7::new_class(
     constructor = function(x, metadata = NULL) {
         # Check input
         stopifnot(.check_input_df(x))
+
+
         if(S7::S7_inherits(x, LinkMap)) {
             if( !NCOL(metadata) ) { metadata <- x@metadata }
             x <- `class<-`(S7::S7_data(x), "data.frame")
-            }
+        }
+        x <- `row.names<-.data.frame`(x, NULL)
         if(!NCOL(metadata)) {
             metadata <- data.frame(row.names = seq_len(NROW(x)))
         } else {
@@ -136,8 +139,7 @@ MultiFactor <- S7::new_class(
                     self <- .set_levels_MultiFactor(self, value)
                 }
                 return(self)
-            },
-            default = quote(as.list(colnames(self)))
+            }
             ),
         map = S7::new_property(
             getter = function(self) .mapMultiFactor(self, mode = "counts")
@@ -164,11 +166,7 @@ MultiFactor <- S7::new_class(
         if( !length(levels) )  levels <- .build_levels_from_linkmap_list(x)
         x <- .unify_levels(x, levels)
 
-        names(x) <- vapply(
-            x,
-            function(x) paste(names(x), collapse = "2"),
-            FUN.VALUE = "", USE.NAMES = FALSE
-        )
+        names(x) <- vapply(x, .linkmap2name, FUN.VALUE = "", USE.NAMES = FALSE)
 
         S7::new_object(
             x,
@@ -289,7 +287,7 @@ MultiFactor <- S7::new_class(
 
         # Filter feature ids in each df to only universally shared ones.
         x[ii] <- lapply(x[ii], function(df) {
-            return(df[df[[j]] %in% keep, ])
+            return(df[df[[j]] %in% keep ])
         })
     }
     return(x)
@@ -361,7 +359,25 @@ MultiFactor <- S7::new_class(
 #' Given a list of linkmaps x and named list of chars levels, unify all levels
 #' across x.
 #' @noRd
-.unify_levels <- function(x, levels) lapply(x, .unify_levels_LinkMap, levels)
+.unify_levels <- function(x, levels) {
+    is_s7 <- all(vapply(x, S7::S7_inherits, LinkMap, FUN.VALUE = FALSE))
+    if(is_s7) {
+        res <- lapply(x, .unify_levels_LinkMap, levels)
+    } else {
+        res <- lapply(x, .unify_levels_data.frame, levels)
+    }
+    return(res)
+}
+
+#' @importFrom forcats lvls_expand
+#'
+.unify_levels_data.frame <- function(x, levels) {
+    x <- mapply(
+        forcats::lvls_expand, x,
+        levels, SIMPLIFY = FALSE
+    )
+    return(x)
+}
 
 
 # TODO metadata and .data are now separate props. Use lapply for indexing?
@@ -369,19 +385,21 @@ MultiFactor <- S7::new_class(
 #' @importFrom forcats lvls_expand
 #'
 .unify_levels_LinkMap <- function(x, levels) {
+
     x[] <- mapply(
         forcats::lvls_expand,
         x,
         levels[colnames(x)],
         SIMPLIFY = FALSE
     )
+    #S7::S7_data(x) <- `class<-`(old, "data.frame")
     return(x)
 }
 
 .validLinkMap <- function(x) {
     is.data.frame(x) &&
-        NCOL(x) >= 2L &&
-        length(colnames(x)) >= 2L &&
-        all(vapply(x[seq_len(2L)], is.factor, NA, USE.NAMES = FALSE))
+        NCOL(x) == 2L &&
+        length(colnames(x)) == 2L &&
+        all(vapply(x, is.factor, NA, USE.NAMES = FALSE))
 }
 
diff --git a/R/LinkMap-methods.R b/R/LinkMap-methods.R
index 178571e..0c152f1 100644
--- a/R/LinkMap-methods.R
+++ b/R/LinkMap-methods.R
@@ -41,10 +41,13 @@ S7::method(levels, LinkMap) <- function(x) x@levels
     return(x)
 }
 
-# S7::method(`levels<-`, LinkMap) <- function(x, value) {
-#     `levels<-.MultiFactor::LinkMap`(x, value)
-# }
-
+#' @export
+#'
+`unique.MultiFactor::LinkMap` <- function(x, incomparables = FALSE, ...) {
+  if (!isFALSE(incomparables))
+    .NotYetUsed("incomparables != FALSE")
+  x[! duplicated(x) ]
+}
 
 #' @importFrom rlang check_dots_empty0 is_missing
 #' @export
@@ -63,13 +66,6 @@ S7::method(levels, LinkMap) <- function(x) x@levels
     return(x)
 }
 
-S7::method(`[`, MultiFactor) <- function(x, i) {
-
-    MultiFactor(base::`[`(S7::S7_data(x), i))
-}
-
-S7::method(`[[`, MultiFactor) <- function(x, i) base::`[[`(S7::S7_data(x), i)
-
 
 #' @title Convert a LinkMap to a sparse matrix.
 #' Convert a LinkMap to a sparse matrix object from the `Matrix` package.
@@ -86,7 +82,7 @@ S7::method(`[[`, MultiFactor) <- function(x, i) base::`[[`(S7::S7_data(x), i)
 #' @seealso [Matrix::sparseMatrix()]
 #'
 `as.matrix.MultiFactor::LinkMap` <- function(
-        x, terms = colnames(x)[seq_len(2L)],
+        x, terms = colnames(x),
         dims = nlevels(x)[terms], dimnames = levels(x)[terms],
         ...
 ) Matrix::sparseMatrix(
@@ -94,3 +90,19 @@ S7::method(`[[`, MultiFactor) <- function(x, i) base::`[[`(S7::S7_data(x), i)
     dims = dims, dimnames = dimnames,
     ...
 )
+
+
+
+##### LinkMap utils
+
+#' @noRd
+.formula2name <- function(.path) {
+   var_list <- .path_parse_formula(.path)
+   var_vctr <- vapply(var_list, paste, collapse = ".", FUN.VALUE = character(1L))
+   var_name <- paste(var_vctr, collapse = "2")
+   return(var_name)
+}
+
+.colnames2name <- function(x) paste(x, collapse = "2")
+
+.linkmap2name <- function(x) paste(names(x), collapse = "2")
diff --git a/R/MultiFactor-methods.R b/R/MultiFactor-methods.R
index 17330bb..398b2a5 100644
--- a/R/MultiFactor-methods.R
+++ b/R/MultiFactor-methods.R
@@ -105,8 +105,9 @@ S7::method(levels, MultiFactor) <- function(x) {
     x@levels
 }
 #' @export
-`levels<-.MultiFactor::MultiFactor` <- function(x, value)
+`levels<-.MultiFactor::MultiFactor` <- function(x, value) {
     .set_levels_MultiFactor(x, value)
+}
 
 S7::method(dimnames, MultiFactor) <- function(x) {
     dimnames(x@map)
@@ -126,6 +127,7 @@ S7::method(`[[`, MultiFactor) <- function(x, i) base::`[[`(S7::S7_data(x), i)
 
 })
 
+
 #' @export
 #'
 `[<-.MultiFactor::MultiFactor` <- function(x, i, value) {
diff --git a/R/MultiFactor-wrangle-utils.R b/R/augment-subset.R
similarity index 77%
rename from R/MultiFactor-wrangle-utils.R
rename to R/augment-subset.R
index 1dcc025..3ab4e9d 100644
--- a/R/MultiFactor-wrangle-utils.R
+++ b/R/augment-subset.R
@@ -1,6 +1,6 @@
 #' Tools to modify MultiFactors
-#' @name MultiFactor-wrangle-methods
-#' @rdname MultiFactor-wrangle-methods
+#' @name augment-subset
+#' @rdname augment-subset
 #' @description
 #' Generates a new `MultiFactor` object by cross-referencing the elements of a
 #'     given `MultiFactor`.
@@ -14,12 +14,14 @@
 #' # Generate a random MultiFactor
 #' x <- randomMultiFactor()
 #'
-#' # Use augment to build upon the same MultiFactor.
+#' # Use augment to chain together operations like weave and stack, in order.
 #' augment(x,
 #'     weave(x, a ~ c),
+#'     stack(x, a + b ~ c + d),
 #'     weave(x, d ~ f)
 #' )
-#' # Setting a LinkMap to NULL deletes it from the MultiFactor
+#'
+#' # Setting a LinkMap to NULL by name deletes it from the MultiFactor
 #' augment(x, a2b = NULL )
 #'
 NULL
@@ -33,7 +35,7 @@ NULL
     if(drop.unmatched) x <- .trimMultiFactor(x)
     if(is.null(subset)) return(x)
     if(by_path){
-        subset <- unlist(.path_terms(subset))
+        subset <- unlist(.path_parse(subset))
         stopifnot("Argument `subset` must be length 2 if by_path` is TRUE" =
                       length(subset) == 2L)
         subset <- termSeq(subset, x)
@@ -58,18 +60,21 @@ S7::method(augment, MultiFactor) <-
     function(x, ...) `augment.MultiFactor::MultiFactor`(x, ...)
 
 #' @export
-#' @rdname MultiFactor-wrangle-methods
+#' @rdname augment-subset
 #' @name augment.MultiFactor
+#' @importFrom rlang dots_list
 #'
 `augment.MultiFactor::MultiFactor` <- function(x, ...) {
     old_names <- rownames(x)
-    dots <- list(...)
+    dots <- rlang::dots_list(...)
     dot_names <- names(dots)
-
     for ( i in seq_along(dot_names) ) {
         res <- eval(dots[[i]])
-        if( dot_names[i] %in% old_names ) {
-            x[[dot_names[i]]] <- res
+        res_name <- dot_names[i]
+        if(res_name == "") { res_name <- .linkmap2name(res) }
+
+        if( res_name %in% old_names ) {
+            x[[res_name]] <- res
         } else {
             x <- c(x, res)
         }
@@ -83,7 +88,7 @@ S7::method(augment, MultiFactor) <-
 .all_names_in_list_mf <- function(x) {
     as.data.frame.matrix(
         t(vapply(X = x,
-                 FUN = function(i) return(names(i)[seq_len(2L)]),
+                 FUN = function(i) return(names(i)),
                  FUN.VALUE = c(NA_character_, NA_character_))
         )
     )
diff --git a/R/levels-utils.R b/R/levels-utils.R
index c9ec14c..2e9dc86 100644
--- a/R/levels-utils.R
+++ b/R/levels-utils.R
@@ -2,6 +2,9 @@
 #' @param x MultiFactor or appropriately formatted list
 #' @param levels A `named list of character vectors`, to be used as
 #'     replacement levels.
+#' @param merge A boolean. Whether to merge or overwrite (default) overlapping
+#'     levels
+#' @importFrom S7 prop<-
 #' @returns a MultiFactor with updated levels.
 #' @noRd
 #'
@@ -9,18 +12,11 @@
     all_lvs <- colnames(x)
     matched_lvs <- names(levels) %in% all_lvs
 
-    stopifnot( "No overlap in 'x' and 'levels'." = sum(matched_lvs) >= 1L )
-
     new_levels <- levels[matched_lvs]
-
-    unchanged <- ! all_lvs %in% names(levels)
-    if( sum(unchanged) >= 1L ) {
-        new_levels <- c( new_levels, .gather_all_levels(x, all_lvs[unchanged]) )
-    }
     S7::S7_data(x) <- .unify_levels(
         `class<-`(S7::S7_data(x), "data.frame"),
         new_levels
-        )
+    )
 
     return(x)
 }
diff --git a/R/path-utils.R b/R/path-utils.R
new file mode 100644
index 0000000..421b1ff
--- /dev/null
+++ b/R/path-utils.R
@@ -0,0 +1,108 @@
+.check_path <- function(.path) {
+    # Initialize output with defaults
+    res <- c(info = "minimal", vars = "ordinary", class = "character")
+    #Some grace for select_path output
+    if(is.list(.path) && length(.path) == 1L) {.path <- .path[[1L]] }
+
+    # Defenses
+    stopifnot("'.path' requires at least two variables." = length(.path) >= 2L )
+    stopifnot(
+        "'.path' must be a formula or a (list of) character vector(s)." =
+            inherits(.path, c("character", "formula", "list"))
+    )
+    if(inherits(.path, "formula")) {
+        # Formula case
+        res["class"] <- "formula"
+        # TODO
+        # FIXME
+        if(.path_function_is_detailed(.path)) {
+            res["info"] <- "detailed"
+        }
+    } else {
+        # Character case
+        if(is.list(.path)) {
+            stopifnot(
+                "'.path' list elements must all be character vectors." =
+                    all(vapply(.path, is.character, FUN.VALUE = FALSE))
+            )
+            if( any(lengths(.path) >= 2L) ) { res["vars"] <- "complex" }
+
+        }
+        if(length(.path) >= 3L ) res["info"] <- "detailed"
+
+    }
+    return(res)
+}
+
+.path_ordinary_to_full <- function(x, .path) {
+    if(inherits(.path, "formula")) {
+    all_terms <- unlist(strsplit(rlang::as_label(.path), " ~ ", fixed = TRUE))
+    } else {
+        all_terms <- .path
+    }
+    terms <- all_terms[c(1L, length(all_terms))]
+    include <- all_terms[-c(1, length(all_terms))]
+    if(!length(include)) { include <- NULL}
+    full_path <- .select_path( x, terms, include)
+    return(full_path)
+}
+
+.build_path <- function(x, .path, pc) {
+    if(pc["class"] == "factor") {
+
+    }
+    if(pc["vars"] == "complex") .path_prep_complex(.path)
+}
+
+#' Standardize terms
+#' @returns a length 2 character vector of y, x.
+#' @noRd
+#'
+.path_parse <- function(.path) {
+
+    if(inherits(.path, "formula")) return(.path_parse_formula(.path))
+
+    if(inherits(.path, "character")) {
+        if(length(.path == 2L)) return(.path) else stop(
+            "Length of '.path' is must be exactly 2 using character input. ",
+            "Use formula syntax for more control."
+        )
+    }
+}
+
+#' @returns BOOL
+#' @noRd
+#' @importFrom rlang as_label
+#'
+.path_function_is_detailed <- function(.path) {
+    # Returns TRUE if more than one "~" seen.
+    length(unlist(strsplit(rlang::as_label(.path), "~", fixed = TRUE))) > 2L
+}
+
+#' @importFrom stats reformulate
+#' @noRd
+#' @returns a list of step-wise formulae.
+.path_prep_complex <- function(.path) {
+    all_terms <- unlist(strsplit(rlang::as_label(.path), "~", fixed = TRUE))
+    lapply(
+        seq_len(length(all_terms) -1L),
+        function(i) stats::reformulate(all_terms[i+1L], all_terms[i])
+        )
+    }
+
+
+#' Standardize terms
+#' @importFrom rlang f_lhs f_rhs
+#' @returns a list of length 2 containing character vectors of y, x.
+#' @noRd
+#'
+.path_parse_formula <- function(.path) {
+    stopifnot( "'.path' must be a formula." = inherits(.path, "formula"))
+    y_vars <- rlang::f_lhs(.path)
+    x_vars <- rlang::f_rhs(.path)
+    if( is.null(y_vars) || is.null(x_vars) ) {
+        stop("Neither formula side can be empty.")
+    }
+    lapply(list(y_vars, x_vars), all.vars)
+}
+
diff --git a/R/select_path.R b/R/select_path.R
index edbf847..c7ecda6 100644
--- a/R/select_path.R
+++ b/R/select_path.R
@@ -27,7 +27,7 @@
 select_path <- function(
         x, .path, include = NULL, exclude = NULL, exact = NULL, as.edges = FALSE
 ) {
-    paths <- .select_path(x, .path_terms(.path), include, exclude, exact)
+    paths <- .select_path(x, .path_parse(.path), include, exclude, exact)
 
     if( as.edges ) paths <- lapply(paths, .V_path_as_E_path)
     return(paths)
@@ -46,6 +46,7 @@ select_path <- function(
         paths <- apply( term.grid, 1L, .apply_shortest_ps, g )
         paths <- lapply( unlist(paths, FALSE, FALSE), as_ids )
     }
+    return(paths)
 }
 
 .subset_paths <- function(g, include, exclude, exact) {
@@ -174,23 +175,12 @@ select_path <- function(
 }
 
 
-.weave_mult <- function(x, terms, out.format) {
-    res <- apply(
-    expand.grid(terms), 1L, .weave_ordinary_terms_df,
-    x = x, out.format = out.format, simplify = FALSE
-    )
-    if(out.format == "LinkMap") {
-        cn <- vapply(lapply(terms, unique), paste, collapse = ".", "")
-        res <- do.call( rbind.data.frame, lapply(res, `colnames<-`, cn) )
-        res <- LinkMap(res)
-    }
-    return(res)
-    }
+
 
 
 
 .weave_ordinary_terms_df <- function(df, x, out.format) .weave_ordinary_terms(
-    x, c(df[seq_len(2L)]), out.format
+    x, c(df), out.format
 )
 
 
diff --git a/R/stack.R b/R/stack.R
new file mode 100644
index 0000000..1c0a88a
--- /dev/null
+++ b/R/stack.R
@@ -0,0 +1,70 @@
+#' Combine levels across several LinkMaps in a MultiFactor
+#' @name stack
+#' @rdname stack.MultiFactor
+#' @description
+#' Generates a new `LinkMap` object by cross-referencing the elements of a
+#'     given `MultiFactor`. Elements can be merged by including several names,
+#'     separated by the plus (`+`) sign. See examples.
+#' @param x a `MultiFactor`
+#' @param .path a `formula` of length 2 with with levels to be merged separated
+#'     by the plus (`+`) sign. Optionally, a list with two character vectors,
+#'     signifying the variables to be combined at the left and right hand side,
+#'     respectively.
+#' @param ... Additional arguments (unused.)
+#' @param out.format `Character scalar`. One of `'LinkMap'`, `'matrix'`.
+#' @returns a `LinkMap` or `sparse Matrix`.
+#' @examples
+#' x <- randomMultiFactor()
+#' # Merge variables with "+" operator, new names get concatenated with ".":
+#' stack(x, b ~ c + d)
+#'
+NULL
+
+#' @export
+#'
+S7::method(stack, MultiFactor) <- function(
+        x, .path, out.format = c("LinkMap", "matrix"), ...
+        ) `stack.MultiFactor::MultiFactor`(x, .path, out.format, ...)
+
+#' @importFrom utils stack
+#' @export
+#' @rdname stack.MultiFactor
+#'
+`stack.MultiFactor::MultiFactor` <- function(
+        x, .path, out.format = c("LinkMap", "matrix"), ...
+) {
+    out.format <- match.arg(out.format, c("LinkMap", "matrix"))
+    .p_check <- .check_path(.path)
+    stopifnot(
+        "stack does not support '.path' with multiple tildes " =
+            .p_check["info"] == "minimal"
+        )
+    if( .p_check["class"] == "formula" ) {
+        terms <- .path_parse_formula(.path)
+    }  else {
+        terms <- .path
+    }
+    stopifnot(
+        "At least one side of '.path' must include several variables." =
+            any( lengths(terms) != 1L)
+        )
+    res <- .stack_terms(x, terms, out.format = "LinkMap")
+    if(out.format == "matrix") {
+        res <- `as.matrix.MultiFactor::LinkMap`(res)
+    }
+    return(res)
+}
+
+.stack_terms <- function(x, terms, out.format) {
+    res <- apply(
+        expand.grid(terms), 1L, .weave_ordinary_terms_df,
+        x = x, out.format = out.format, simplify = FALSE
+    )
+    if(out.format == "LinkMap") {
+        cn <- vapply(lapply(terms, unique), paste, collapse = ".", "")
+        res <- do.call( rbind.data.frame, lapply(res, `colnames<-`, cn) )
+        res <- LinkMap(res)
+    }
+    return(res)
+}
+
diff --git a/R/subgroup_apply.R b/R/subgroup_apply.R
index 2486cf4..9a67743 100644
--- a/R/subgroup_apply.R
+++ b/R/subgroup_apply.R
@@ -91,7 +91,7 @@ subgroup_apply <- function( X, LINK, BY, FUN = NULL, ..., INDEX = "row.names" )
     # Ensure link is a MultiFactor
     link <- MultiFactor(link)
 
-    terms <- .path_terms(.path)
+    terms <- .path_parse(.path)
     X <- .index_tbl(X, type = terms[[2L]], .i = .i)
     terms[[2L]] <- colnames(X)[[2L]]
     X <- MultiFactor(X)
diff --git a/R/weave-coverage.R b/R/weave-coverage.R
index 03eec1a..a20dde0 100644
--- a/R/weave-coverage.R
+++ b/R/weave-coverage.R
@@ -22,7 +22,7 @@
 #' )
 #'
 #' # Now enrich test input
-#' weave(x, a ~ b ~ c) |>
+#' weave(x, a ~ b ~ c)
 #'     #test_set_enrichment()
 #' @noRd
 test_set_enrichment <- function(
diff --git a/R/weave-formula-utils.R b/R/weave-formula-utils.R
deleted file mode 100644
index 99ee4f3..0000000
--- a/R/weave-formula-utils.R
+++ /dev/null
@@ -1,89 +0,0 @@
-.check_by_terms <- function(.path) {
-    res <- c("single", "character")
-    #Some grace for select_path output
-    if(is.list(.path)) {
-        stopifnot(
-            "'.path' must be a formula or character vector." = length(.path) == 1L
-        )
-        .path <- .path[[1L]]
-    }
-    stopifnot(
-        "'.path' must be a character vector or a formula." =
-            inherits(.path, c("character", "formula"))
-    )
-    if(is.character(.path)) {
-        if(length(.path) < 2L) stop(
-            "Length of '.path' is must be at least 2 using character input. "
-        )
-        if(length(.path) >= 3L ) res[1L] <- "full"
-    }
-    if(inherits(.path, "formula")) {
-        res[2L] <- "formula"
-        if(.path_is_complex(.path)) {
-            res[1L] <- "complex"
-        }
-    }
-    return(res)
-}
-
-
-#' Standardize terms
-#' @returns a length 2 character vector of y, x.
-#' @noRd
-#'
-.path_terms <- function(.path) {
-
-
-    if(inherits(.path, "formula")) return(.weave_parse_formula(.path))
-
-    if(inherits(.path, "character")) {
-        if(length(.path == 2L)) return(.path) else stop(
-            "Length of '.path' is must be exactly 2 using character input. ",
-            "Use formula syntax for more control."
-        )
-    }
-}
-
-#' @returns BOOL
-#' @noRd
-#' @importFrom rlang as_label
-#'
-.path_is_complex <- function(.path) {
-    # Only support complex through formula
-    if(is.character(.path)) return(FALSE)
-
-    stopifnot(
-        "'.path' must be a character vector or a formula." =
-            inherits(.path, c("character", "formula"))
-    )
-    # Returns TRUE if more than one "~" seen.
-    length(unlist(strsplit(rlang::as_label(.path), "~", fixed = TRUE))) > 2L
-}
-
-#' @importFrom stats reformulate
-#' @noRd
-#' @returns a list of step-wise formulae.
-.path_prep_complex_call <- function(.path) {
-    all_terms <- unlist(strsplit(rlang::as_label(.path), "~", fixed = TRUE))
-    lapply(
-        seq_len(length(all_terms) -1L),
-        function(i) stats::reformulate(all_terms[i+1L], all_terms[i])
-        )
-    }
-
-
-#' Standardize terms
-#' @importFrom rlang f_lhs f_rhs
-#' @returns a list of length 2 containing character vectors of y, x.
-#' @noRd
-#'
-.weave_parse_formula <- function(.path) {
-    stopifnot( "'.path' must be a formula." = inherits(.path, "formula"))
-    y_vars <- rlang::f_lhs(.path)
-    x_vars <- rlang::f_rhs(.path)
-    if( is.null(y_vars) || is.null(x_vars) ) {
-        stop("Neither formula side can be empty.")
-    }
-    lapply(list(y_vars, x_vars), all.vars)
-}
-
diff --git a/R/weave.R b/R/weave.R
index 18b9f04..c57e046 100644
--- a/R/weave.R
+++ b/R/weave.R
@@ -21,57 +21,33 @@
 #' weave(x, b ~ c)
 #' weave(x, b ~ a, out.format = "matrix")
 #'
-#' # Merge variables with "+" operator, new names get concatenated with ".":
-#' weave(x, b ~ c + d)
-#'
-#' # Control intermediate variable types "~", returning a MultiFactor:
-#' weave(x, a ~ c ~ e )
-#'
-#' # Combine merging and intermediate stops:
-#' weave(x, a ~ b + c ~ d + e ~ f )
-#'
 NULL
 
 #' @export
 #'
 S7::method(weave, MultiFactor) <- function(
-    x, .path, out.format = c("LinkMap", "matrix", "MultiFactor"),
+    x, .path, out.format = c("LinkMap", "matrix"),
     include = NULL, exclude = NULL, exact = NULL
 ) {
-  out.format <- match.arg(out.format, c("LinkMap", "matrix", "MultiFactor"))
+  out.format <- match.arg(out.format, c("LinkMap", "matrix"))
   lv_list <- levels(x)
-  out.MF <- out.format == "MultiFactor"
-  if(out.MF) {
-    out.format <- "LinkMap"
-  }
-  .path_type <- .check_by_terms(.path)
 
-  if(.path_type[1] == "complex") {
-    out.MF <- TRUE
-    res <- .weave_complex_formula(x, .path, "LinkMap")
-    lv_list <- .weave_complex_formula_lvs(res, lv_list)
-  }
-  if(.path_type[1] == "full") {
-    res <- .weave_single_path(x, .path, out.format)
-  }
-  if( .path_type[1L] == "single" ) {
-    terms <- if( .path_type[2L] == "formula" ) .weave_parse_formula(.path) else .path
+  .p_check <- .check_path(.path)
 
-    if ( all( lengths(terms) == 1L) ) {
-      res <- .weave_ordinary_terms(x, unlist(terms), out.format)
-    } else {
-      res <- .weave_mult(x, terms, out.format)
-    }
+  if(.p_check["vars"] == "complex") {
+    stop("weave() '.path' cannot contain '+'. Use stack() to prepare input.")
   }
-  if( out.MF ) {
-    res <- MultiFactor(res, lv_list)
-  } else if( out.format == "LinkMap" ) {
-    res <- as.LinkMap(res)
+    full_path <- .path_ordinary_to_full(x, .path)[[1L]]
+    res <- .weave_full_path(x, full_path, out.format)
+
+  if( out.format == "LinkMap" ) {
+    res <- LinkMap(res)
   }
 
   return(res)
 }
 
+
 weave_along_path <- function(x, path, out.format = "LinkMap") {
   # tolerate single path result in list
   if(length(path) == 1L) path <- path[[1L]]
@@ -83,7 +59,7 @@ weave_along_path <- function(x, path, out.format = "LinkMap") {
     "All entries in 'path' must be found in colnames(x)." =
       all( path %in% colnames(x) )
   )
-  res <- .weave_single_path(x, path, out.format)
+  res <- .weave_full_path(x, path, out.format)
   # Remove duplicates
   if(out.format == "LinkMap") {
     res <- res[ !duplicated(res[, seq_len(2L)]), ]
@@ -124,7 +100,7 @@ path_coverage <- function(x, path, out.format = "matrix") {
 }
 
 .weave_complex_formula <- function(x, .path, out.format) {
-  .path_list <- .path_prep_complex_call(.path)
+  .path_list <- .path_prep_complex(.path)
   res <- lapply( .path_list, weave, x = x, out.format = out.format )
 }
 
@@ -149,7 +125,7 @@ path_coverage <- function(x, path, out.format = "matrix") {
   lv_list <- levels(x)
   # Determine required ids in order, only keep relevant elements of link.
   all_terms <- .select_path(x, terms, include = NULL, exclude = NULL, exact = NULL)
-  res <- lapply(all_terms, .weave_single_path, x = x, out.format = "LinkMap")
+  res <- lapply(all_terms, .weave_full_path, x = x, out.format = "LinkMap")
   res <- do.call(
     rbind.data.frame,
     c(res, make.row.names = FALSE, stringsAsFactors = TRUE)
@@ -164,7 +140,7 @@ path_coverage <- function(x, path, out.format = "matrix") {
 }
 
 
-.weave_single_path <- function(x, all_terms, out.format) {
+.weave_full_path <- function(x, all_terms, out.format) {
   x <- subsetByPath(x, all_terms)
   terms <- all_terms[c(1L, length(all_terms))]
   # Construct dictionary
diff --git a/man/as.matrix.LinkMap.Rd b/man/as.matrix.LinkMap.Rd
index 046afc6..bce7a43 100644
--- a/man/as.matrix.LinkMap.Rd
+++ b/man/as.matrix.LinkMap.Rd
@@ -8,7 +8,7 @@ Convert a LinkMap to a sparse matrix object from the \code{Matrix} package.}
 \usage{
 \method{as.matrix}{`MultiFactor::LinkMap`}(
   x,
-  terms = colnames(x)[seq_len(2L)],
+  terms = colnames(x),
   dims = nlevels(x)[terms],
   dimnames = levels(x)[terms],
   ...
diff --git a/man/MultiFactor-wrangle-methods.Rd b/man/augment-subset.Rd
similarity index 72%
rename from man/MultiFactor-wrangle-methods.Rd
rename to man/augment-subset.Rd
index 9b0d55c..ca5b7d0 100644
--- a/man/MultiFactor-wrangle-methods.Rd
+++ b/man/augment-subset.Rd
@@ -1,7 +1,7 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/MultiFactor-wrangle-utils.R
-\name{MultiFactor-wrangle-methods}
-\alias{MultiFactor-wrangle-methods}
+% Please edit documentation in R/augment-subset.R
+\name{augment-subset}
+\alias{augment-subset}
 \alias{augment.MultiFactor}
 \alias{augment.MultiFactor::MultiFactor}
 \title{Tools to modify MultiFactors}
@@ -27,12 +27,14 @@ require(generics)
 # Generate a random MultiFactor
 x <- randomMultiFactor()
 
-# Use augment to build upon the same MultiFactor.
+# Use augment to chain together operations like weave and stack, in order.
 augment(x,
     weave(x, a ~ c),
+    stack(x, a + b ~ c + d),
     weave(x, d ~ f)
 )
-# Setting a LinkMap to NULL deletes it from the MultiFactor
+
+# Setting a LinkMap to NULL by name deletes it from the MultiFactor
 augment(x, a2b = NULL )
 
 }
diff --git a/man/stack.MultiFactor.Rd b/man/stack.MultiFactor.Rd
new file mode 100644
index 0000000..a007d88
--- /dev/null
+++ b/man/stack.MultiFactor.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/stack.R
+\name{stack}
+\alias{stack}
+\alias{stack.MultiFactor::MultiFactor}
+\title{Combine levels across several LinkMaps in a MultiFactor}
+\usage{
+\method{stack}{`MultiFactor::MultiFactor`}(x, .path, out.format = c("LinkMap", "matrix"), ...)
+}
+\arguments{
+\item{x}{a \code{MultiFactor}}
+
+\item{.path}{a \code{formula} of length 2 with with levels to be merged separated
+by the plus (\code{+}) sign. Optionally, a list with two character vectors,
+signifying the variables to be combined at the left and right hand side,
+respectively.}
+
+\item{out.format}{\verb{Character scalar}. One of \code{'LinkMap'}, \code{'matrix'}.}
+
+\item{...}{Additional arguments (unused.)}
+}
+\value{
+a \code{LinkMap} or \verb{sparse Matrix}.
+}
+\description{
+Generates a new \code{LinkMap} object by cross-referencing the elements of a
+given \code{MultiFactor}. Elements can be merged by including several names,
+separated by the plus (\code{+}) sign. See examples.
+}
+\examples{
+x <- randomMultiFactor()
+# Merge variables with "+" operator, new names get concatenated with ".":
+stack(x, b ~ c + d)
+
+}
diff --git a/man/weave-methods.Rd b/man/weave-methods.Rd
index f848712..9b42f68 100644
--- a/man/weave-methods.Rd
+++ b/man/weave-methods.Rd
@@ -31,13 +31,4 @@ x <- randomMultiFactor()
 weave(x, b ~ c)
 weave(x, b ~ a, out.format = "matrix")
 
-# Merge variables with "+" operator, new names get concatenated with ".":
-weave(x, b ~ c + d)
-
-# Control intermediate variable types "~", returning a MultiFactor:
-weave(x, a ~ c ~ e )
-
-# Combine merging and intermediate stops:
-weave(x, a ~ b + c ~ d + e ~ f )
-
 }
diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml
index a0afbe2..f8442f2 100644
--- a/pkgdown/_pkgdown.yml
+++ b/pkgdown/_pkgdown.yml
@@ -19,6 +19,7 @@ reference:
 - contents:
   - weave
   - weave-methods
+  - stack
   - select_path
   - subgroup_apply
   - subgroup_to_tbl
@@ -26,7 +27,7 @@ reference:
 - contents:
   - MultiFactor
   - MultiFactor-methods
-  - MultiFactor-wrangle-methods
+  - augment-subset
   - LinkMap
   - LinkMap-methods
   - nlevels
diff --git a/tests/testthat/test-weave-formula-utils.R b/tests/testthat/test-weave-formula-utils.R
index df0a8ed..93c21d3 100644
--- a/tests/testthat/test-weave-formula-utils.R
+++ b/tests/testthat/test-weave-formula-utils.R
@@ -1,6 +1,6 @@
 test_that("complex formula is parsing works", {
     f <- a ~ b + c  ~ d + e + f ~ g
-    res <- lapply(.path_prep_complex_call(f), deparse1)
+    res <- lapply(.path_prep_complex(f), deparse1)
     expect_identical(
         res, list("a ~ b + c", "b + c ~ d + e + f", "d + e + f ~ g")
     )