diff --git a/DESCRIPTION b/DESCRIPTION
index 036e7ae..6114413 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -10,20 +10,19 @@ Authors@R: c(
 Description: An S7 framework to organize and manage multiple sets of factors, 
     for instance when tracing or converting feature IDs across databases. 
     Methods for MultiFactor aim to follow factor behaviour.
+Depends: 
+    R (>= 4.1.0)
 Imports:
-    S7, 
+    S7,
     rlang,
     generics,
-    Matrix, 
-    forcats, 
+    Matrix,
+    forcats,
     igraph
-Depends: 
-    R (>= 4.1.0)
 License: GPL-3
 LazyData: false
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.3
 Suggests: 
     anansi,
     ggplot2,
@@ -36,3 +35,4 @@ Suggests:
 Config/testthat/edition: 3
 URL: https://minotau-r.github.io/MultiFactor/
 VignetteBuilder: knitr
+Config/roxygen2/version: 8.0.0
diff --git a/NAMESPACE b/NAMESPACE
index c0d7068..1624de5 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method("[","MultiFactor::LinkMap")
 S3method("[<-","MultiFactor::MultiFactor")
 S3method("[[<-","MultiFactor::MultiFactor")
+S3method("levels<-","MultiFactor::LinkMap")
 S3method("levels<-","MultiFactor::MultiFactor")
 S3method(as.igraph,"MultiFactor::LinkMap")
 S3method(as.igraph,"MultiFactor::MultiFactor")
@@ -51,8 +53,10 @@ importFrom(igraph,induced_subgraph)
 importFrom(igraph,sample_gnm)
 importFrom(igraph,shortest_paths)
 importFrom(rlang,as_label)
+importFrom(rlang,check_dots_empty0)
 importFrom(rlang,f_lhs)
 importFrom(rlang,f_rhs)
+importFrom(rlang,is_missing)
 importFrom(stats,fisher.test)
 importFrom(stats,reformulate)
 importFrom(utils,count.fields)
diff --git a/R/AllClasses.R b/R/AllClasses.R
index b809eaf..f90a375 100644
--- a/R/AllClasses.R
+++ b/R/AllClasses.R
@@ -11,6 +11,8 @@
 #'     application-specific tags.
 #' @param x `data.frame` with two named columns that can be coerced to factors.
 #'     Optionally, additional columns will be stored as metadata.
+#' @param metadata Optional `data.frame` with same number of rows as x. Contains
+#'     information about the feature link in that row.
 #' @returns a `LinkMap` object.
 #' @examples
 #' # Generate random linkage input
@@ -31,38 +33,47 @@ LinkMap <- S7::new_class(
     parent = S7::class_data.frame,
     properties = list(
         levels = S7::new_property(
-            getter = function(self) lapply(
-                S7::S7_data(self)[c(1, 2)], levels
-            )
+            getter = function(self) lapply(S7::S7_data(self), levels),
+            setter = function(self, value) {
+                x <- `class<-`(S7::S7_data(self), "data.frame")
+                S7::S7_data(self) <- .unify_levels_LinkMap(x, value)
+                return(self)
+            }
         ),
         metadata = S7::new_property(
-            getter = function(self) `class<-`(
-                S7::S7_data(self), "data.frame"
-            )[-c(1, 2)]
+            class = S7::class_data.frame,
+            getter = function(self) self@metadata
         )
     ),
-    constructor = function(x) {
-        if(S7::S7_inherits(x, LinkMap)) return(x)
+    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")
+            }
+        if(!NCOL(metadata)) {
+            metadata <- data.frame(row.names = seq_len(NROW(x)))
+        } else {
+            stopifnot(
+                "Arg 'x' must have the same number of rows as 'metadata'" =
+                    NROW(x) == NROW(metadata)
+            )
+        }
         # Factorize x
-        x[ c(1, 2)] <- lapply(x[c(1, 2)], as.factor)
-        x <- x[ !duplicated(x[, c(1, 2)]), ]
+        x[] <- lapply(x, factor)
+        i <- !duplicated(x)
+        x <- x[i, , drop = FALSE]
+        metadata <- metadata[i, , drop = FALSE]
 
-        S7::new_object(x)
-    },
+        S7::new_object(x, metadata = metadata)
+        },
     validator = function(self) {
-        if(!is.data.frame(self)) {
-            "Must be a data.frame. "
-        }
-        if(! NCOL(self) == 2L) {
-            "Must be a data.frame with exactly two columns. "
-        }
-        if(!length(colnames(self)) == 2L) {
-            "Both columns must be named. "
-        }
-        if(!all(vapply(self[c(1, 2)], is.factor, NA, USE.NAMES = FALSE))){
-            "Both columns must be factors. "
+        if( !is.data.frame(self) ) { "Must be a data.frame." }
+        if( NCOL(self) != 2L ) { "Must be a data.frame with two columns." }
+        if( length(colnames(self)) != 2L ) { "Both columns must be named." }
+        if( !all(vapply(self, is.factor, NA, USE.NAMES = FALSE)) ) {
+            "Both columns must be factors."
         }
     }
 )
@@ -115,8 +126,18 @@ MultiFactor <- S7::new_class(
     parent = S7::class_list,
     properties  = list(
         levels = S7::new_property(
-            class = S7::class_list, getter = function(self) self@levels,
-             setter = function(self, value) .set_levels_MultiFactor(self, value)
+            class = S7::class_list,
+            getter = function(self) self@levels,
+            setter = function(self, value) {
+                # If initializing, don't re-unify
+                if( !length(levels(self)) ) {
+                    self@levels <- value
+                } else {
+                    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")
@@ -127,7 +148,7 @@ MultiFactor <- S7::new_class(
             )
         )
     ),
-    constructor = function(x, levels = NULL) {
+    constructor = function(x, levels = list()) {
         # Check input
         if(is.data.frame(x)) x <- LinkMap(x)
         if(S7::S7_inherits(x, LinkMap)) x <- list(x = x)
@@ -140,7 +161,7 @@ MultiFactor <- S7::new_class(
         names(x) <- paste0("x_", seq_along(x))
         x   <- .merge_linkmaps(x)
 
-        if(is.null(levels))  levels <- .build_levels(x)
+        if( !length(levels) )  levels <- .build_levels_from_linkmap_list(x)
         x <- .unify_levels(x, levels)
 
         names(x) <- vapply(
@@ -150,7 +171,7 @@ MultiFactor <- S7::new_class(
         )
 
         S7::new_object(
-            .parent = x,
+            x,
             levels = levels
         )
     },
@@ -169,10 +190,10 @@ MultiFactor <- S7::new_class(
     if(! is.data.frame(x) ) {
         stop("Must be a data.frame. ")
     }
-    if(! NCOL(x) >= 2L ) {
+    if(! NCOL(x) == 2L ) {
         stop("Must be a data.frame with at least two key columns. ")
     }
-    if(! length(colnames(x)[seq_len(2L)]) == 2L ) {
+    if(! length(colnames(x)) == 2L ) {
         stop("Both key columns must be named. ")
     }
     return( TRUE )
@@ -208,7 +229,7 @@ MultiFactor <- S7::new_class(
     mx <- switch(mode,
                  "counts" = unlist(
                      lapply(x, function(y) {
-                         lapply(y[seq_len(2L)], function(z) length(unique(z)))
+                         lapply(y, function(z) length(unique(z)))
                      }),
                      use.names = FALSE
                  ),
@@ -231,7 +252,7 @@ MultiFactor <- S7::new_class(
 
 
 .merge_linkmaps <- function(x) {
-    all_names <- lapply(x, \(x) sort(names(x)))
+    all_names <- lapply(x, function(x) sort(names(x)))
     if(!any(duplicated(all_names))) return(x)
 
     dup_names <- unique(all_names[duplicated(all_names)])
@@ -305,7 +326,15 @@ MultiFactor <- S7::new_class(
 
 #' Given a list of linkmaps x, return a unified and sorted list of levels.
 #' @noRd
-.build_levels <- function(x) {
+.build_levels_from_df_list <- function(x) {
+    all_lvs <- unique(unlist(lapply(x, colnames), FALSE, FALSE))
+    lvs <- .gather_all_levels(x, all_lvs)
+    return(lvs)
+}
+
+#' Given a list of linkmaps x, return a unified and sorted list of levels.
+#' @noRd
+.build_levels_from_linkmap_list <- function(x) {
     all_lvs <- unique(unlist(lapply(x, colnames), FALSE, FALSE))
     lvs <- .gather_all_levels(x, all_lvs)
     return(lvs)
@@ -334,19 +363,25 @@ MultiFactor <- S7::new_class(
 #' @noRd
 .unify_levels <- function(x, levels) lapply(x, .unify_levels_LinkMap, levels)
 
+
+# TODO metadata and .data are now separate props. Use lapply for indexing?
+# FIXME
 #' @importFrom forcats lvls_expand
 #'
 .unify_levels_LinkMap <- function(x, levels) {
     x[] <- mapply(
-        forcats::lvls_expand, x[seq_len(2L)], levels[colnames(x)], SIMPLIFY = FALSE
+        forcats::lvls_expand,
+        x,
+        levels[colnames(x)],
+        SIMPLIFY = FALSE
     )
     return(x)
 }
 
 .validLinkMap <- function(x) {
     is.data.frame(x) &&
-        NCOL(x) == 2L &&
-        length(colnames(x)) == 2L &&
+        NCOL(x) >= 2L &&
+        length(colnames(x)) >= 2L &&
         all(vapply(x[seq_len(2L)], is.factor, NA, USE.NAMES = FALSE))
 }
 
diff --git a/R/LinkMap-methods.R b/R/LinkMap-methods.R
index b0f70b1..178571e 100644
--- a/R/LinkMap-methods.R
+++ b/R/LinkMap-methods.R
@@ -24,28 +24,6 @@
 #' @returns A `LinkMap`
 NULL
 
-S7::method(names, LinkMap) <- function(x) names(
-    S7::S7_data(x)[c(1, 2)]
-    )
-
-S7::method(dimnames, LinkMap) <- function(x) dimnames(
-    `class<-`(S7::S7_data(x), "data.frame")[c(1, 2)]
-    )
-
-S7::method(dim, LinkMap) <- function(x) dim(
-    `class<-`(S7::S7_data(x), "data.frame")[c(1, 2)]
-    )
-
-S7::method(print, LinkMap) <- function(x, ...) print(
-    `class<-`(S7::S7_data(x), "data.frame")[c(1, 2)], ...
-)
-
-S7::method(str, LinkMap) <- function(object, ...) str(
-    `class<-`(S7::S7_data(object), "data.frame")
-)
-
-S7::method(levels, LinkMap) <- function(x) lapply(x[c(1, 2)], levels)
-
 #' @param use.names `Boolean scalar` Should names be provided.
 #'     (Default: `TRUE`)
 #' @noRd
@@ -54,6 +32,45 @@ S7::method(nlevels, LinkMap) <- function(x, use.names = TRUE) lengths(
     levels(x), use.names
     )
 
+S7::method(levels, LinkMap) <- function(x) x@levels
+
+#' @export
+#'
+`levels<-.MultiFactor::LinkMap` <- function(x, value) {
+    x@levels <- value
+    return(x)
+}
+
+# S7::method(`levels<-`, LinkMap) <- function(x, value) {
+#     `levels<-.MultiFactor::LinkMap`(x, value)
+# }
+
+
+#' @importFrom rlang check_dots_empty0 is_missing
+#' @export
+#'
+`[.MultiFactor::LinkMap` <- function(x, i, ...) {
+    rlang::check_dots_empty0(...)
+    if(! rlang::is_missing(i))  {
+        metadata <- x@metadata
+        metadata <- `[.data.frame`(metadata, i, , drop = FALSE)
+
+        x <- `class<-`(S7::S7_data(x), "data.frame")
+        x <- `[.data.frame`(x, i, , drop = FALSE)
+
+        x <- LinkMap(x, metadata)
+        }
+    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.
 #' @name as.matrix.LinkMap
@@ -69,8 +86,8 @@ S7::method(nlevels, LinkMap) <- function(x, use.names = TRUE) lengths(
 #' @seealso [Matrix::sparseMatrix()]
 #'
 `as.matrix.MultiFactor::LinkMap` <- function(
-        x, terms = colnames(x),
-        dims = nlevels(x[terms]), dimnames = levels(x)[terms],
+        x, terms = colnames(x)[seq_len(2L)],
+        dims = nlevels(x)[terms], dimnames = levels(x)[terms],
         ...
 ) Matrix::sparseMatrix(
     i = x[[terms[1L]]], j = x[[terms[2L]]],
diff --git a/R/MultiFactor-methods.R b/R/MultiFactor-methods.R
index 44266df..17330bb 100644
--- a/R/MultiFactor-methods.R
+++ b/R/MultiFactor-methods.R
@@ -104,12 +104,10 @@ S7::method(print, MultiFactor) <- function(x, ...) {
 S7::method(levels, MultiFactor) <- function(x) {
     x@levels
 }
-
 #' @export
 `levels<-.MultiFactor::MultiFactor` <- function(x, value)
     .set_levels_MultiFactor(x, value)
 
-
 S7::method(dimnames, MultiFactor) <- function(x) {
     dimnames(x@map)
 }
@@ -155,7 +153,7 @@ S7::method(`[[`, MultiFactor) <- function(x, i) base::`[[`(S7::S7_data(x), i)
     lv_list <- lapply(x[mfs], levels)
     if(any(lms)) {
         lm_lst <- x[lms]
-        lv_list <- c(lv_list, list(.build_levels(lm_lst)))
+        lv_list <- c(lv_list, list(.build_levels_from_linkmap_list(lm_lst)))
     }
     all_lvs <- unique(unlist(lapply(lv_list, names), FALSE, FALSE))
     all_lvs <- .reduce_level_list(lv_list, all_lvs)
diff --git a/R/MultiFactor-wrangle-utils.R b/R/MultiFactor-wrangle-utils.R
index 26bd0a6..1dcc025 100644
--- a/R/MultiFactor-wrangle-utils.R
+++ b/R/MultiFactor-wrangle-utils.R
@@ -77,4 +77,15 @@ S7::method(augment, MultiFactor) <-
     return(x)
 }
 
+#' Safely get all names from a MultiFactor as a data.frame.
+#' @noRd
+#'
+.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.VALUE = c(NA_character_, NA_character_))
+        )
+    )
+}
 
diff --git a/R/as.LinkMap.R b/R/as.LinkMap.R
index 8bed800..08065ef 100644
--- a/R/as.LinkMap.R
+++ b/R/as.LinkMap.R
@@ -45,19 +45,27 @@ NULL
 
 S7::method(as.LinkMap, LinkMap) <- function(x) x
 
+
 S7::method(as.LinkMap, S7::class_data.frame) <- function(
-    x, edge.names = NULL){
+        x, edge.names = NULL
+){
+    if( NCOL(x) >= 3L ) {
+        metadata <- x[, -seq_len(2L)]
+        x <- x[, seq_len(2L)]
+    } else {
+        metadata <- data.frame(row.names = seq_len(NROW(x)))
+        }
     # Assign custom edge names
     if( !is.null(edge.names) ){
         colnames(x) <- edge.names
     }
     # Convert to LinkMap
-    x2y <- LinkMap(x)
-    return(x2y)
+    x <- LinkMap(x, metadata)
+    return(x)
 }
 
 S7::method(as.LinkMap, S7::class_list) <- function(
-    x, y = NULL, edge.names = NULL){
+    x, y = NULL, edge.names = NULL ){
     # Use list names
     if( is.null(y) ){
         y <- names(x)
diff --git a/R/graph_conversion_utils.R b/R/graph_conversion_utils.R
index 64d73bd..a3d6b91 100644
--- a/R/graph_conversion_utils.R
+++ b/R/graph_conversion_utils.R
@@ -83,14 +83,9 @@ S7::method(as.igraph, MultiFactor) <-
 #' x_df <- mf_as_graph_df(x)
 #'
 mf_as_graph_df <- function(x) {
-    x <- MultiFactor(x)
-
-    res <- as.data.frame.matrix(
-        t(vapply(x, names, c(NA_character_, NA_character_)))
-    )
-    cbind.data.frame(
-        res,
-        x@metadata
-    )
+    if(!S7::S7_inherits(x, MultiFactor)) x <- MultiFactor(x)
+    res <- as.data.frame.matrix( .all_names_in_list_mf(x) )
+    res <- cbind.data.frame(  res, x@metadata )
+    return(res)
 }
 
diff --git a/R/levels-utils.R b/R/levels-utils.R
index c160bf9..c9ec14c 100644
--- a/R/levels-utils.R
+++ b/R/levels-utils.R
@@ -6,7 +6,7 @@
 #' @noRd
 #'
 .set_levels_MultiFactor <- function(x, levels) {
-    all_lvs <- unique(unlist(lapply(x, colnames), FALSE, FALSE))
+    all_lvs <- colnames(x)
     matched_lvs <- names(levels) %in% all_lvs
 
     stopifnot( "No overlap in 'x' and 'levels'." = sum(matched_lvs) >= 1L )
@@ -17,13 +17,11 @@
     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
+        )
 
-    if(S7::S7_inherits(x, MultiFactor)) {
-        S7::S7_data(x) <- .unify_levels(S7::S7_data(x), new_levels)
-        x@levels <- new_levels
-    } else {
-        x <- MultiFactor(x, new_levels)
-    }
     return(x)
 }
 
diff --git a/R/randomMultiFactor.R b/R/randomMultiFactor.R
index 33c945a..abd1703 100644
--- a/R/randomMultiFactor.R
+++ b/R/randomMultiFactor.R
@@ -147,7 +147,7 @@ stopifnot("If provided, 'x' must be a list of two named character vectors" =
         "'sparseness' must be a proportion [0-1]. " = sparseness <= 1 &&
             sparseness > 0
     )
-    LinkMap(.randomLinkDF(
+    as.LinkMap(.randomLinkDF(
         x[[1]], x[[2]], names(x)[[1]], names(x)[[2]],
         p = 1 - sparseness
         ))
diff --git a/R/select_path.R b/R/select_path.R
index d4600d6..edbf847 100644
--- a/R/select_path.R
+++ b/R/select_path.R
@@ -208,7 +208,7 @@ termSeq <- function(terms, x) {
         )
     )
     g <- igraph::graph_from_data_frame(
-        d = t(vapply(x, names, c(NA_character_, NA_character_))),
+        d = .all_names_in_list_mf(x),
         directed = FALSE
     )
     sp <- igraph::all_shortest_paths(
diff --git a/R/subgroup_apply.R b/R/subgroup_apply.R
index d19fb89..2486cf4 100644
--- a/R/subgroup_apply.R
+++ b/R/subgroup_apply.R
@@ -18,12 +18,17 @@
 #'
 #' # Prepare data
 #'
-#' link <- anansi::kegg_link()
-#' data("FMT_data", package = "anansi")
-#' x <- FMT_KOs
+#' link <- trade_posts()
+#'
+#' # Generate small example feature table 'x'.
+#' n <- nlevels(link)[["clothing"]]
+#' x <- replicate(10, rbinom(n, rbinom(n, 100, runif(n)), runif(n)))
+#'
+#' #' # Ensure rownames correspond to the second (RHS) variable in the formula.
+#' x <- as.data.frame(x, row.names = levels(link)$clothing)
 #'
 #' # Apply arbitrary code to x based on group membership
-#' subgroup_apply(x, link, BY =  ec ~ ko, FUN = function(x) dim(x))
+#' subgroup_apply(x, link, BY = fruit ~ clothing, FUN = function(x) colSums(x))
 #'
 #' @seealso [weave()] [LinkMap()] [MultiFactor()]
 #' @export
@@ -63,7 +68,7 @@ subgroup_apply <- function( X, LINK, BY, FUN = NULL, ..., INDEX = "row.names" )
     res <- data.frame( row.index, seq_len(NROW(X)) )
 
     colnames(res) <- c(type, "row.index")
-    LinkMap(res, ...)
+    as.LinkMap(res)
 }
 
 #' Index a table
diff --git a/R/tidy.R b/R/tidy.R
index eeb505e..be9006d 100644
--- a/R/tidy.R
+++ b/R/tidy.R
@@ -14,19 +14,27 @@
 #' @importFrom Matrix which
 #' @examples
 #' # Prepare data
-#' link <- anansi::kegg_link()
-#' data("FMT_data", package = "anansi")
-#' x <- FMT_KOs
-#' subgroup_to_tbl(x, link, .path = ec ~ ko)
+#' link <- trade_posts()
+#'
+#' # Generate small example feature table 'x'.
+#' n <- nlevels(link)[["clothing"]]
+#' x <- replicate(10, rbinom(n, rbinom(n, 100, runif(n)), runif(n)))
+#'
+#' #' # Ensure rownames correspond to the second (RHS) variable in the formula.
+#' x <- as.data.frame(x, row.names = levels(link)$clothing)
+#'
+#' # Apply arbitrary code to x based on group membership
+#' subgroup_to_tbl(x, link, .path = fruit ~ clothing)
 #'
 #' @export
 #'
 subgroup_to_tbl <- function(x, link, .path, .index = "row.names") {
-   out <- .index_tbl_by(x, link, .path, .index)
+   out <- S7::S7_data( .index_tbl_by(x, link, .path, .index) )
+   class(out) <- "data.frame"
    out[[2L]] <- as.integer(out[[2L]])
    xcol <- if(.index == "row.names") row.names(x) else x[[.index]]
 
-   subgroup <- data.frame(y = out[1L], x = xcol[out[[2L]]])
+   subgroup <- data.frame(y = out[1L], x = xcol[out[[2L]]], check.rows = FALSE)
    colnames(subgroup) <- if(inherits(.path, "formula")) all.vars(.path) else .path
 
    out <- data.frame(subgroup, x[out[[2L]], ])
diff --git a/R/weave.R b/R/weave.R
index 4b89af6..18b9f04 100644
--- a/R/weave.R
+++ b/R/weave.R
@@ -66,7 +66,7 @@ S7::method(weave, MultiFactor) <- function(
   if( out.MF ) {
     res <- MultiFactor(res, lv_list)
   } else if( out.format == "LinkMap" ) {
-    res <- LinkMap(res)
+    res <- as.LinkMap(res)
   }
 
   return(res)
@@ -86,7 +86,7 @@ weave_along_path <- function(x, path, out.format = "LinkMap") {
   res <- .weave_single_path(x, path, out.format)
   # Remove duplicates
   if(out.format == "LinkMap") {
-    res <- res[ !duplicated(res[, c(1, 2)]), ]
+    res <- res[ !duplicated(res[, seq_len(2L)]), ]
   }
   res
 }
@@ -129,7 +129,7 @@ path_coverage <- function(x, path, out.format = "matrix") {
 }
 
 .weave_complex_formula_lvs <- function(x, lv_list) {
-  new_lvs <- .build_levels(x)
+  new_lvs <- .build_levels_from_linkmap_list(x)
   new_names <- names(new_lvs)
   kept <- intersect(new_names, names(lv_list))
   new_lvs <- list(lv_list[kept], new_lvs)
@@ -200,7 +200,7 @@ path_coverage <- function(x, path, out.format = "matrix") {
     )
   )
   g <- igraph::graph_from_data_frame(
-    d = t(vapply(x, names, c(NA_character_, NA_character_))),
+    d = .all_names_in_list_mf(x),
     directed = FALSE
   )
   igraph::all_shortest_paths(
diff --git a/R/zzz.R b/R/zzz.R
index 3921892..079d8e6 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -2,10 +2,9 @@
 #' @rawNamespace if (getRversion() < "4.3.0") importFrom("S7", "@")
 NULL
 
-rlang::on_load({
-    S7::methods_register()
-})
+
 
 .onLoad <- function(...) {
-    rlang::run_on_load()
+    S7::methods_register()
 }
+
diff --git a/man/LinkMap-class.Rd b/man/LinkMap-class.Rd
index 1569bb5..bc71d25 100644
--- a/man/LinkMap-class.Rd
+++ b/man/LinkMap-class.Rd
@@ -4,11 +4,14 @@
 \alias{LinkMap}
 \title{LinkMap S7 container class}
 \usage{
-LinkMap(x)
+LinkMap(x, metadata = NULL)
 }
 \arguments{
 \item{x}{\code{data.frame} with two named columns that can be coerced to factors.
 Optionally, additional columns will be stored as metadata.}
+
+\item{metadata}{Optional \code{data.frame} with same number of rows as x. Contains
+information about the feature link in that row.}
 }
 \value{
 a \code{LinkMap} object.
diff --git a/man/MultiFactor-class.Rd b/man/MultiFactor-class.Rd
index ffd2e24..9f90290 100644
--- a/man/MultiFactor-class.Rd
+++ b/man/MultiFactor-class.Rd
@@ -4,7 +4,7 @@
 \alias{MultiFactor}
 \title{MultiFactor S7 container class}
 \usage{
-MultiFactor(x, levels = NULL)
+MultiFactor(x, levels = list())
 }
 \arguments{
 \item{x}{a \code{LinkMap}, or named list of \code{LinkMap} objects.}
diff --git a/man/as.matrix.LinkMap.Rd b/man/as.matrix.LinkMap.Rd
index 712b022..046afc6 100644
--- a/man/as.matrix.LinkMap.Rd
+++ b/man/as.matrix.LinkMap.Rd
@@ -8,8 +8,8 @@ Convert a LinkMap to a sparse matrix object from the \code{Matrix} package.}
 \usage{
 \method{as.matrix}{`MultiFactor::LinkMap`}(
   x,
-  terms = colnames(x),
-  dims = nlevels(x[terms]),
+  terms = colnames(x)[seq_len(2L)],
+  dims = nlevels(x)[terms],
   dimnames = levels(x)[terms],
   ...
 )
diff --git a/man/subgroup_apply.Rd b/man/subgroup_apply.Rd
index 8c47e4b..67cd088 100644
--- a/man/subgroup_apply.Rd
+++ b/man/subgroup_apply.Rd
@@ -33,12 +33,17 @@ of row-based subsets of table \code{X}. Mimics lapply.
 
 # Prepare data
 
-link <- anansi::kegg_link()
-data("FMT_data", package = "anansi")
-x <- FMT_KOs
+link <- trade_posts()
+
+# Generate small example feature table 'x'.
+n <- nlevels(link)[["clothing"]]
+x <- replicate(10, rbinom(n, rbinom(n, 100, runif(n)), runif(n)))
+
+#' # Ensure rownames correspond to the second (RHS) variable in the formula.
+x <- as.data.frame(x, row.names = levels(link)$clothing)
 
 # Apply arbitrary code to x based on group membership
-subgroup_apply(x, link, BY =  ec ~ ko, FUN = function(x) dim(x))
+subgroup_apply(x, link, BY = fruit ~ clothing, FUN = function(x) colSums(x))
 
 }
 \seealso{
diff --git a/man/subgroup_to_tbl.Rd b/man/subgroup_to_tbl.Rd
index fe937e3..2c09656 100644
--- a/man/subgroup_to_tbl.Rd
+++ b/man/subgroup_to_tbl.Rd
@@ -27,9 +27,16 @@ Expand a table into a tidy table suitable for tidyverse-stype operations.
 }
 \examples{
 # Prepare data
-link <- anansi::kegg_link()
-data("FMT_data", package = "anansi")
-x <- FMT_KOs
-subgroup_to_tbl(x, link, .path = ec ~ ko)
+link <- trade_posts()
+
+# Generate small example feature table 'x'.
+n <- nlevels(link)[["clothing"]]
+x <- replicate(10, rbinom(n, rbinom(n, 100, runif(n)), runif(n)))
+
+#' # Ensure rownames correspond to the second (RHS) variable in the formula.
+x <- as.data.frame(x, row.names = levels(link)$clothing)
+
+# Apply arbitrary code to x based on group membership
+subgroup_to_tbl(x, link, .path = fruit ~ clothing)
 
 }
diff --git a/vignettes/MultiFactor.Rmd b/vignettes/MultiFactor.Rmd
index 0c020e9..1c04443 100644
--- a/vignettes/MultiFactor.Rmd
+++ b/vignettes/MultiFactor.Rmd
@@ -194,7 +194,7 @@ interact with its matrix representation, as well as access the component levels:
 dim(tp)
 dimnames(tp)
 levels(tp)
-lengths(levels(tp))
+nlevels(tp)
 ```
 Note that levels of the same type are automatically unified across all component
 LinkMaps with that type of level.
@@ -263,7 +263,7 @@ plot(lg)
 # Convert `LinkMap` linkage information to sparse adjacency `Matrix`
 ```{r 'Matrix-sp'}
 # Convert to an igraph object
-m <- as.matrix(fruit2clothing, dimnames = levels(fruit2clothing))
+m <- as.matrix(fruit2clothing)
 m
 ```