Repository: bids-standard/pybids
Pull request #548 (2019-11-15 18:52:18)

ENH: Model spec module and associated utilities
This PR introduces a new `model_design` module containing new classes meant to serve as intermediate abstractions between the `BIDSVariableCollection` classes produces by PyBIDS and model estimation classes (e.g., those in fitlins). There are various other improvements and changes.

Summary of changes:
* New `bids.analysis.model_spec` module. The key feature here is the `GLMMSpec` class, which provides a Pythonic interface to the components of the model specified in the BIDS-StatsModels `model` block (now including variance components/random effects).
* Transformations are now managed via the new `TransformerManager`, which can be used independently of the transformations implemented in PyBIDS, and allows registration of arbitrary named transformations.
* The API in the `bids.analysis` module has changed fairly dramatically. The `DesignInfo` namedtuple is gone, along with `get_design_matrix()`. Internally, the `AnalysisNode` abstraction has been removed. This means you now get design matrices by calling `to_df()` directly on a `BIDSVariableCollection`. You can get collections from a `Step` using `get_collections()`. See the tests for examples. On the whole, this approach is slightly more verbose for the user, but much more computationally efficient (pretty much everything is now computed only when needed, whereas before it was easy to waste cycles generating structures that would never get used). I think it's also conceptually cleaner. Basically, instead of passing around lists of `ContrastInfo`, `DesignInfo`, etc., the expectation is now that the user gets the collections they need, and then loops over those explicitly and constructs design matrices, retrieve contrasts, and whatnot *inside* the loop.
* The `to_df()` method in the `BIDSVariableCollection` has undergone some changes. Complexity has been reduced somewhat, in that there are fewer options for controlling how sparse and dense variables are dealt with. Basically the idea is now that if you want to force the output to look different, you may have to explicitly call `.resample()` or `.to_dense()` on the collection before you generate the DF. Again this is a little bit more verbose, but I think it makes what's going on less opaque, and improves code clarity.
* Resampling and densification functionality have been streamlined throughout the package as described in #572. There are now collection-level `.to_dense()` and `.resample()` methods that share nearly all of the underlying logic. I'm not thrilled about having two separate functions that do almost the same thing–especially because each operation can implicitly invoke the other—but semantically, I think it makes sense to keep them separate.
* The `variables.kollekshuns` module has been renamed to `variables.collections` now that we've dropped Python 2 support and don't have to worry about reserved module names.



Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6f533b0d..c72fbd87 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -291,7 +291,7 @@ This is a major release that introduces the following features:
 	* `bids.variables.entities.py`: Classes for representing BIDS hierarchies as a graph-like structure.
 	* `bids.variables.variables.py`: Classes and functions for representing and manipulating non-imaging data read from BIDS projects (e.g., fMRI events, densely-sampled physiological measures, etc.).
 	* `bids.variables.io.py`: Tools for loading variable data from BIDS projects.
-	* `bids.variables.kollekshuns`: Containers that facilitate aggregation and manipulation of `Variable` classes.
+	* `bids.variables.collections`: Containers that facilitate aggregation and manipulation of `Variable` classes.
 * Extensions to the `BIDSLayout` class that make it easy to retrieve data/variables from the project (i.e., `Layout.get_collections`)
 * A new `auto_model` utility that generates simple BIDS-Model specifications from BIDS projects (thanks to @Shotgunosine)
 * A new `reports` module that generates methods sections from metadata in BIDS projects (thanks to @tsalo)
diff --git a/bids/analysis/__init__.py b/bids/analysis/__init__.py
index e70531e8..a9182a6f 100644
--- a/bids/analysis/__init__.py
+++ b/bids/analysis/__init__.py
@@ -1,4 +1,9 @@
 from .analysis import Analysis
 from .auto_model import auto_model
+from .transformations import TransformerManager
 
-__all__ = ['Analysis', 'auto_model']
+__all__ = [
+    'Analysis',
+    'auto_model',
+    'TransformerManager'
+]
diff --git a/bids/analysis/analysis.py b/bids/analysis/analysis.py
index db6c002d..008a2860 100644
--- a/bids/analysis/analysis.py
+++ b/bids/analysis/analysis.py
@@ -1,13 +1,18 @@
+"""BIDS-StatsModels functionality."""
+
 import json
-from io import open
-from bids.layout import BIDSLayout
-from bids.utils import matches_entities, convert_JSON
-from bids.variables import BIDSVariableCollection, merge_collections
-from . import transformations as transform
 from collections import namedtuple, OrderedDict
+from itertools import chain
+
 import numpy as np
 import pandas as pd
-from itertools import chain
+
+from bids.layout import BIDSLayout
+from bids.utils import matches_entities, convert_JSON
+from bids.variables import (BIDSVariableCollection, SparseRunVariable,
+                            merge_collections)
+from bids.analysis import transformations as tm
+from .model_spec import create_model_spec
 
 
 class Analysis(object):
@@ -42,7 +47,7 @@ class Analysis(object):
         level = index.lower()
         name_matches = list(filter(lambda x: x.name == level, self.steps))
         if not name_matches:
-            raise KeyError('There is no block with the name "%s".' % index)
+            raise KeyError('There is no step with the name "%s".' % index)
         return name_matches[0]
 
 
@@ -57,7 +62,6 @@ class Analysis(object):
         steps = self.model['steps']
         self.steps = []
         for i, step_args in enumerate(steps):
-            step_args['level'] = step_args['level'].lower()
             step = Step(self.layout, index=i, **step_args)
             self.steps.append(step)
 
@@ -69,7 +73,7 @@ class Analysis(object):
         steps : list
             Optional list of steps to set up. Each element
             must be either an int giving the index of the step in the
-            JSON config block list, or a str giving the (unique) name of
+            JSON config step list, or a str giving the (unique) name of
             the step, as specified in the JSON config. Steps that do not
             match either index or name will be skipped.
         drop_na : bool
@@ -78,60 +82,59 @@ class Analysis(object):
             from event files.
         """
 
-        # In the beginning, there was nothing
-        input_nodes = None
+        # The first Step in the sequence can't have any contrast inputs
+        input_contrasts = None
 
         # Use inputs from model, and update with kwargs
         selectors = self.model.get('input', {}).copy()
         selectors.update(kwargs)
 
-        for i, b in enumerate(self.steps):
+        for i, step in enumerate(self.steps):
 
-            # Skip any steps whose names or indexes don't match block list
-            if steps is not None and i not in steps and b.name not in steps:
+            # Skip any steps whose names or indexes don't match step list
+            if steps is not None and i not in steps and step.name not in steps:
                 continue
 
-            b.setup(input_nodes, drop_na=drop_na, **selectors)
-            input_nodes = b.output_nodes
+            step.setup(input_contrasts, drop_na=drop_na, **selectors)
+            input_contrasts = [step.get_contrasts(c)
+                               for c in step.get_collections(**selectors)]
+            input_contrasts = list(chain(*input_contrasts))
 
 
 class Step(object):
-    """Represents a single analysis block from a BIDS-Model specification.
+    """Represents a single analysis step from a BIDS-Model specification.
 
     Parameters
     ----------
     layout : :obj:`bids.layout.BIDSLayout`
         The BIDSLayout containing all project files.
     level : str
-        The BIDS keyword to use as the grouping variable; must be
-        one of ['run', 'session', 'subject', or 'dataset'].
+        The BIDS keyword to use as the grouping variable; must be one of
+        ['run', 'session', 'subject', or 'dataset'].
     index : int
-        The numerical index of the current Step within the
-        sequence of steps.
+        The numerical index of the current Step within the sequence of steps.
     name : str
-        Optional name to assign to the block. Must be specified
-        in order to enable name-based indexing in the parent Analysis.
+        Optional name to assign to the step. Must be specified in order to
+        enable name-based indexing in the parent Analysis.
     transformations : list
         List of BIDS-Model transformations to apply.
     model : dict
-        The 'model' part of the BIDS-Model block specification.
+        The 'model' part of the BIDS-StatsModels specification.
     contrasts : list
-        List of contrasts to apply to the parameter estimates
-        generated when the model is fit.
-    input_nodes : list
-        Optional list of AnalysisNodes to use as input to
-        this Step (typically, the output from the preceding Step).
+        List of contrasts to apply to the parameter estimates generated when
+        the model is fit.
+    inputs : list
+        Optional list of BIDSVariableCollections to use as input to this Step
+        (typically, the outputs from the preceding Step).
     dummy_contrasts : dict
-        Optional dictionary specifying which conditions to create
-        indicator contrasts for. Dictionary must include a
-        "type" key ('t' or 'FEMA'), and optionally a subset of "conditions".
-        This parameter is over-written by the setting
-        in setup() if the latter is passed.
+        Optional dictionary specifying which conditions to create indicator
+        contrasts for. Dictionary must include a "type" key ('t' or 'FEMA'),
+        and optionally a subset of "conditions". This parameter is over-written
+        by the setting in setup() if the latter is passed.
     """
 
     def __init__(self, layout, level, index, name=None, transformations=None,
-                 model=None, contrasts=None, input_nodes=None,
-                 dummy_contrasts=False):
+                 model=None, contrasts=None, inputs=None, dummy_contrasts=False):
         self.layout = layout
         self.level = level.lower()
         self.index = index
@@ -139,23 +142,26 @@ class Step(object):
         self.transformations = transformations or []
         self.model = model or None
         self.contrasts = contrasts or []
-        self.input_nodes = input_nodes or []
+        self.inputs = inputs or []
         self.dummy_contrasts = dummy_contrasts
-        self.output_nodes = []
+        self._collections = []
 
-    def _filter_objects(self, objects, kwargs):
-        # Keeps only objects that match target entities, and also removes those
-        # keys from the kwargs dict.
+    def _filter_collections(self, collections, kwargs):
+        # Keeps only collections that match target entities, and also removes
+        # those keys from the kwargs dict.
         valid_ents = {'task', 'subject', 'session', 'run'}
         entities = {k: kwargs.pop(k) for k in dict(kwargs) if k in valid_ents}
-        objects = [o for o in objects if o.matches_entities(entities)]
-        return (objects, kwargs)
+        collections = [c for c in collections if matches_entities(c, entities)]
+        return (collections, kwargs)
 
-    def _group_objects(self, objects):
+    def _group_objects_by_entities(self, objects):
         # Group list of objects into bins defined by all entities at current
-        # Step level or higher.
+        # Step level or higher. E.g., if the level is 'subject', the
+        # returned list will have one element per subject, where each element
+        # is a list containing all objects that belongs to that subject. Any
+        # object with a defined .entities attribute is groupable.
         if self.level == 'dataset':
-            return [objects]
+            return {'dataset': objects}
         groups = OrderedDict()
         valid_ents = ['subject', 'session', 'task', 'run']
         valid_ents = valid_ents[:(valid_ents.index(self.level) + 1)]
@@ -165,37 +171,78 @@ class Step(object):
             if key not in groups:
                 groups[key] = []
             groups[key].append(o)
-        return list(groups.values())
-
-    def _concatenate_input_nodes(self, nodes):
-        data, entities = [], []
-        for n in nodes:
-            contrasts = [c.name for c in n.contrasts]
-            row = pd.Series(np.ones(len(contrasts)), index=contrasts)
-            data.append(row)
-            entities.append(pd.Series(n.entities))
-        data = pd.concat(data, axis=1, sort=True).T
-        entities = pd.concat(entities, axis=1, sort=True).T
-        return BIDSVariableCollection.from_df(data, entities, self.level)
-
-    def setup(self, input_nodes=None, drop_na=False, **kwargs):
+        return groups
+
+    def _merge_contrast_inputs(self, inputs):
+        """ Merges a list of ContrastInfo tuples and constructs a dict mapping
+        from units of the current level to BIDSVariableCollections.
+
+        Parameters
+        ----------
+        inputs: [[ContrastInfo]]
+            List of list of ContrastInfo tuples passed from the previous Step.
+            Each element in the outer list maps to the output of a unit at the
+            previous level; each element in the inner list is a ContrastInfo
+            tuple. E.g., if contrast information is being passed from run-level
+            to subject-level, each outer element is a run.
+
+        Returns
+        -------
+        A dictionary, where the keys are the values of the entities at the
+        current level (e.g., '01', '02'...) and the values are
+        BIDSVariableCollection containing contrast information.
+
+        Notes
+        -----
+        Each output BIDSVariableCollection contains information for a single
+        unit at the present level. The variables in the collection reflect the
+        union of all contrasts found in one or more of the inputs. A value of
+        1 indicates that the contrast is present for a given row in the input;
+        0 indicates that the contrast was missing.
+        """
+
+        groups = self._group_objects_by_entities(inputs)
+
+        ent_cols = list(list(groups.values())[0][0].entities.keys())
+
+        collections = {}
+
+        for name, contrasts in groups.items():
+            # Create a DF with contrasts in rows and contrast names in columns
+            data = [{**c.entities, **{c.name: 1}} for c in contrasts]
+            data = pd.DataFrame.from_records(data)
+            # Group by all entities and sum, collapsing over rows belonging
+            # to the current unit
+            data = data.groupby(ent_cols).sum()
+            # Split the DF up into separate data and entities DFs
+            entities = data.index.to_frame(index=False)
+            data = data.reset_index(drop=True)
+            # Construct the collection
+            coll = BIDSVariableCollection.from_df(data, entities, self.level)
+            collections[name] = coll
+
+        return collections
+
+    def setup(self, inputs=None, drop_na=False, **kwargs):
         """Set up the Step and construct the design matrix.
 
         Parameters
         ----------
-        input_nodes : list
-            Optional list of Node objects produced by
-            the preceding Step in the analysis. If None, uses any inputs
-            passed in at Step initialization.
+        inputs : list
+            Optional list of BIDSVariableCollections produced as output by the
+            preceding Step in the analysis. If None, uses inputs passed at
+            initialization (if any).
         drop_na : bool
-            Boolean indicating whether or not to automatically
-            drop events that have a n/a amplitude when reading in data
-            from event files.
+            Boolean indicating whether or not to automatically drop events that
+            have a n/a amplitude when reading in data from event files.
         kwargs : dict
             Optional keyword arguments to pass onto load_variables.
         """
-        self.output_nodes = []
-        input_nodes = input_nodes or self.input_nodes or []
+        self._collections = []
+
+        # Convert input contrasts to a list of BIDSVariableCollections
+        inputs = inputs or self.inputs or []
+        input_grps = self._merge_contrast_inputs(inputs) if inputs else {}
 
         # TODO: remove the scan_length argument entirely once we switch tests
         # to use the synthetic dataset with image headers.
@@ -203,263 +250,65 @@ class Step(object):
             kwargs = kwargs.copy()
             kwargs.pop('scan_length', None)
 
+        # Now handle variables read from the BIDS dataset: read them in, filter
+        # on passed selectors, and group by unit of current level
         collections = self.layout.get_collections(self.level, drop_na=drop_na,
                                                   **kwargs)
-        objects = collections + input_nodes
+        collections, _ = self._filter_collections(collections, kwargs)
+        groups = self._group_objects_by_entities(collections)
 
-        objects, kwargs = self._filter_objects(objects, kwargs)
-        groups = self._group_objects(objects)
+        # Merge in the inputs
+        for key, input_ in input_grps.items():
+            if key not in groups:
+                groups[key] = []
+            groups[key].append(input_)
 
         # Set up and validate variable lists
         model = self.model or {}
         X = model.get('x', [])
 
-        for grp in groups:
-            # Split into separate lists of Collections and Nodes
-            input_nodes = [o for o in grp if isinstance(o, AnalysisNode)]
-            colls = list(set(grp) - set(input_nodes))
-
-            if input_nodes:
-                node_coll = self._concatenate_input_nodes(input_nodes)
-                colls.append(node_coll)
+        for grp, colls in groups.items():
+            coll = merge_collections(colls)
 
-            coll = merge_collections(colls) if len(colls) > 1 else colls[0]
+            colls = tm.TransformerManager().transform(coll, self.transformations)
 
-            coll = apply_transformations(coll, self.transformations)
             if X:
-                transform.Select(coll, X)
+                tm.Select(coll, X)
 
-            node = AnalysisNode(self.level, coll, self.contrasts, input_nodes,
-                                self.dummy_contrasts)
+            self._collections.append(coll)
 
-            self.output_nodes.append(node)
-
-    def get_design_matrix(self, names=None, format='long', mode='both',
-                          force=False, sampling_rate='TR', **kwargs):
-        """Get design matrix and associated information.
+    def get_collections(self, **filters):
+        """Returns BIDSVariableCollections at the current Step.
 
         Parameters
         ----------
-        names : list
-            Optional list of names of variables to include in the
-            returned design matrix. If None, all variables are included.
-        format : str
-            Whether to return the design matrix in 'long' or
-            'wide' format. Note that dense design matrices are always
-            returned in 'wide' format.
-        mode : str
-            Specifies whether to return variables in a sparse
-            representation ('sparse'), dense representation ('dense'), or
-            both ('both').
-        force : bool
-            Indicates how to handle columns not of the type
-            indicated by the mode argument. When False, variables of the
-            non-selected type will be silently ignored. When True,
-            variables will be forced to the desired representation. For
-            example, if mode='dense' and force=True, sparse variables will
-            be converted to dense variables and included in the returned
-            design matrix in the .dense attribute. The force argument is
-            ignored entirely if mode='both'.
-        sampling_rate : {'TR', 'highest'} or float
-            Sampling rate at which to
-            generate the dense design matrix. When 'TR', the repetition
-            time is used, if available, to select the sampling rate (1/TR).
-            When 'highest', all variables are resampled to the highest
-            sampling rate of any variable. The sampling rate may also be
-            specified explicitly in Hz. Has no effect on sparse design
-            matrices.
-        kwargs : dict
-            Optional keyword arguments. Includes (1) selectors used
-            to constrain which of the available nodes get returned
-            (e.g., passing subject=['01', '02'] will return design
-            information only for subjects '01' and '02'), and (2) arguments
-            passed on to each Variable's to_df() call (e.g.,
-            sampling_rate, entities, timing, etc.).
+        filters : dict
+            Optional keyword filters used to constrain which of the available
+            collections get returned (e.g., passing subject=['01', '02'] will
+            return collections for only subjects '01' and '02').
 
         Returns
         -------
-        list of DesignMatrixInfo namedtuples
-            one tuple per unit of the current
-            analysis level (e.g., if level='run', each element in the list
-            represents the design matrix for a single run).
+        list of BIDSVariableCollection instances
+            One instance per unit of the current analysis level (e.g., if
+            level='run', each element in the list represents the collection
+            for a single run).
         """
-        nodes, kwargs = self._filter_objects(self.output_nodes, kwargs)
-        return [n.get_design_matrix(names, format, mode=mode, force=force,
-                                    sampling_rate=sampling_rate, **kwargs)
-                for n in nodes]
+        return self._filter_collections(self._collections, filters)[0]
 
-    def get_contrasts(self, names=None, variables=None, **kwargs):
-        """Return contrast information for the current block.
+    def get_contrasts(self, collection, names=None, variables=None):
+        """Return contrast information at this step for the passed collection.
 
         Parameters
         ----------
+        collection : BIDSVariableCollection
+            The collection to generate/validate contrasts for.
         names : list
-            Optional list of names of contrasts to return. If
-            None (default), all contrasts are returned.
+            Optional list of names of contrasts to return. If None (default),
+            all contrasts are returned.
         variables : bool
-            Optional list of strings giving the names of
-            design matrix columns to use when generating the matrix of
-            weights.
-        kwargs : dict
-            Optional keyword arguments used to constrain which of the
-            available nodes get returned (e.g., passing subject=['01',
-            '02'] will return contrast  information only for subjects '01'
-            and '02').
-
-        Returns
-        -------
-        list
-            A list with one element per unit of the current analysis level
-            (e.g., if level='run', each element in the list representing the
-            contrast information for a single run). Each element is a list of
-            ContrastInfo namedtuples (one per contrast).
-
-        """
-        nodes, kwargs = self._filter_objects(self.output_nodes, kwargs)
-        return [n.get_contrasts(names, variables) for n in nodes]
-
-
-DesignMatrixInfo = namedtuple('DesignMatrixInfo',
-                              ('sparse', 'dense', 'entities'))
-
-
-ContrastInfo = namedtuple('ContrastInfo', ('name', 'weights', 'type',
-                                           'entities'))
-
-
-class AnalysisNode(object):
-    """A single analysis node generated within a Step.
-
-    Parameters
-    ----------
-    level : str
-        The level of the Node. Most be one of 'run', 'session',
-        'subject', or 'dataset'.
-    collection : :obj:`bids.variables.kollekshuns.BIDSVariableCollection`
-        The BIDSVariableCollection containing variables at this Node.
-    contrasts : list
-        A list of contrasts defined in the originating Step.
-    dummy_contrasts : list
-        Optional dictionary specifying which conditions to create
-        indicator contrasts for. Dictionary must include a
-        "type" key ('t' or 'FEMA'), and optionally a subset of "conditions".
-        This parameter is over-written by the setting
-        in setup() if the latter is passed.
-    """
-
-    def __init__(self, level, collection, contrasts, input_nodes=None,
-                 dummy_contrasts=None):
-        self.level = level.lower()
-        self.collection = collection
-        self._block_contrasts = contrasts
-        self.input_nodes = input_nodes
-        self.dummy_contrasts = dummy_contrasts
-        self._contrasts = None
-
-    @property
-    def entities(self):
-        return self.collection.entities
-
-    @property
-    def contrasts(self):
-        if self._contrasts is None:
-            self.get_contrasts()
-        return self._contrasts
-
-    def get_design_matrix(self, names=None, format='long', mode='both',
-                          force=False, sampling_rate='TR', **kwargs):
-        """Get design matrix and associated information.
-
-        Parameters
-        ----------
-        names : list
-            Optional list of names of variables to include in the
-            returned design matrix. If None, all variables are included.
-        format : str
-            Whether to return the design matrix in 'long' or
-            'wide' format. Note that dense design matrices are always
-            returned in 'wide' format.
-        mode : str
-            Specifies whether to return variables in a sparse
-            representation ('sparse'), dense representation ('dense'), or
-            both ('both').
-        force : bool
-            Indicates how to handle columns not of the type
-            indicated by the mode argument. When False, variables of the
-            non-selected type will be silently ignored. When True,
-            variables will be forced to the desired representation. For
-            example, if mode='dense' and force=True, sparse variables will
-            be converted to dense variables and included in the returned
-            design matrix in the .dense attribute. The force argument is
-            ignored entirely if mode='both'.
-        sampling_rate : {'TR', 'highest'} or float
-            Sampling rate at which to
-            generate the dense design matrix. When 'TR', the repetition
-            time is used, if available, to select the sampling rate (1/TR).
-            When 'highest', all variables are resampled to the highest
-            sampling rate of any variable. The sampling rate may also be
-            specified explicitly in Hz. Has no effect on sparse design
-            matrices.
-        kwargs : dict
-            Optional keyword arguments to pass onto each Variable's
-            to_df() call (e.g., sampling_rate, entities, timing, etc.).
-
-        Returns
-        -------
-        A DesignMatrixInfo namedtuple.
-        """
-        sparse_df, dense_df = None, None
-        coll = self.collection
-
-        if self.level != 'run' and mode != 'sparse':
-            mode = 'sparse'
-
-        include_sparse = include_dense = (force and mode != 'both')
-
-        if mode in ['sparse', 'both']:
-            kwargs['sparse'] = True
-            sparse_df = coll.to_df(names, format, include_dense=include_dense,
-                                   **kwargs)
-
-        if mode in ['dense', 'both']:
-            # The current implementation of pivoting to wide in
-            # BIDSVariableCollection.to_df() breaks if we don't have the
-            # temporal columns to index on, so we force their inclusion first
-            # and then drop them afterwards.
-            kwargs['timing'] = True
-            kwargs['sparse'] = False
-
-            if sampling_rate == 'TR':
-                trs = {var.run_info[0].tr for var in self.collection.variables.values()}
-                if not trs:
-                    raise ValueError("Repetition time unavailable; specify sampling_rate "
-                                     "explicitly")
-                elif len(trs) > 1:
-                    raise ValueError("Non-unique Repetition times found ({!r}); specify "
-                                     "sampling_rate explicitly")
-                sampling_rate = 1. / trs.pop()
-            elif sampling_rate == 'highest':
-                sampling_rate = None
-            dense_df = coll.to_df(names, format='wide',
-                                  include_sparse=include_sparse,
-                                  sampling_rate=sampling_rate, **kwargs)
-            if dense_df is not None:
-                dense_df = dense_df.drop(['onset', 'duration'], axis=1)
-
-        return DesignMatrixInfo(sparse_df, dense_df, self.entities)
-
-    def get_contrasts(self, names=None, variables=None):
-        """Return contrast information for the current block.
-
-        Parameters
-        ----------
-        names : list
-            Optional list of names of contrasts to return. If
-            None (default), all contrasts are returned.
-        variables : bool
-            Optional list of strings giving the names of
-            design matrix columns to use when generating the matrix of
-            weights.
+            Optional list of strings giving the names of design matrix columns
+            to use when generating the matrix of weights.
 
         Returns
         -------
@@ -469,21 +318,20 @@ class AnalysisNode(object):
         Notes
         -----
         The 'variables' argument take precedence over the natural process
-        of column selection. I.e.,
-            if a variable shows up in a contrast, but isn't named in
-            variables, it will *not* be included in the returned
+        of column selection. I.e., if a variable shows up in a contrast, but
+        isn't named in variables, it will *not* be included in the result.
         """
 
         # Verify that there are no invalid columns in the condition_lists
-        all_conds = [c['condition_list'] for c in self._block_contrasts]
+        all_conds = [c['condition_list'] for c in self.contrasts]
         all_conds = set(chain(*all_conds))
-        bad_conds = all_conds - set(self.collection.variables.keys())
+        bad_conds = all_conds - set(collection.variables.keys())
         if bad_conds:
             raise ValueError("Invalid condition names passed in one or more "
                              " contrast condition lists: %s." % bad_conds)
 
         # Construct a list of all contrasts, including dummy contrasts
-        contrasts = list(self._block_contrasts)
+        contrasts = list(self.contrasts)
 
         # Check that all contrasts have unique name
         contrast_names = [c['name'] for c in contrasts]
@@ -494,9 +342,9 @@ class AnalysisNode(object):
         if self.dummy_contrasts:
             if 'conditions' in self.dummy_contrasts:
                 conditions = [c for c in self.dummy_contrasts['conditions']
-                              if c in self.collection.variables.keys()]
+                              if c in collection.variables.keys()]
             else:
-                conditions = self.collection.variables.keys()
+                conditions = collection.variables.keys()
 
             for col_name in conditions:
                 if col_name not in contrast_names:
@@ -522,52 +370,47 @@ class AnalysisNode(object):
 
             test_type = c.get('type', ('t' if len(weights) == 1 else 'F'))
 
-            return ContrastInfo(c['name'], weights, test_type, self.entities)
+            return ContrastInfo(c['name'], weights, test_type,
+                                collection.entities)
 
-        self._contrasts = [setup_contrast(c) for c in contrasts]
+        return [setup_contrast(c) for c in contrasts]
 
-        return self._contrasts
-
-    def matches_entities(self, entities, strict=False):
-        """Determine whether current AnalysisNode matches passed entities.
+    def get_model_spec(self, collection, sampling_rate='TR'):
+        """Get a ModelSpec instance for the passed collection.
 
         Parameters
         ----------
-        entities : dict
-            Dictionary of entities to match. Keys are entity
-            names; values are single values or lists.
-        strict : bool
-            If True, _all_ entities in the current Node must
-            match in order to return True.
+        collection : BIDSVariableCollection
+            The BIDSVariableCollection to construct a model for.
+        sampling_rate : {'TR', 'highest'} or float
+            For run-level models, the sampling rate at which to generate the
+            design matrix. When 'TR', the repetition time is used, if
+            available, to select the sampling rate (1/TR). When 'highest', all
+            variables are resampled to the highest sampling rate of any
+            variable in the collection. The sampling rate may also be specified
+            explicitly in Hz. Has no effect on non-run-level collections.
+
+        Returns
+        -------
+        A bids.analysis.model_spec.ModelSpec instance.
+
+        Notes
+        -----
+        If the passed BIDSVariableCollection contains any sparse variables,
+        they will be automatically converted to dense (using the specified
+        sampling rate) before the ModelSpec is constructed. For non-run-level
+        collections, timing is irrelevant, so the design matrix is constructed
+        based on the "as-is" values found in each variable.
         """
-        return matches_entities(self, entities, strict)
+        if self.model is None:
+            raise ValueError("Cannot generate a ModelSpec instance; no "
+                             "BIDS-StatsModels model specification found "
+                             "for this step!")
 
+        if collection.level == 'run':
+            collection = collection.resample(sampling_rate, force_dense=True)
+        return create_model_spec(collection, self.model)
 
-def apply_transformations(collection, transformations, select=None):
-    """Apply all transformations to the variables in the collection.
 
-    Parameters
-    ----------
-    transformations : list
-        List of transformations to apply.
-    select : list
-        Optional list of names of variables to retain after all
-        transformations are applied.
-    """
-    for t in transformations:
-        kwargs = dict(t)
-        func = kwargs.pop('name')
-        cols = kwargs.pop('input', None)
-
-        if isinstance(func, str):
-            if func in ('and', 'or'):
-                func += '_'
-            if not hasattr(transform, func):
-                raise ValueError("No transformation '%s' found!" % func)
-            func = getattr(transform, func)
-            func(collection, cols, **kwargs)
-
-    if select is not None:
-        transform.Select(collection, select)
-
-    return collection
+ContrastInfo = namedtuple('ContrastInfo', ('name', 'weights', 'type',
+                                           'entities'))
diff --git a/bids/analysis/model_spec.py b/bids/analysis/model_spec.py
new file mode 100644
index 00000000..b165dc1e
--- /dev/null
+++ b/bids/analysis/model_spec.py
@@ -0,0 +1,332 @@
+from abc import ABCMeta, abstractmethod
+
+import pandas as pd
+import numpy as np
+
+from bids.variables import BIDSRunVariableCollection
+from bids.utils import convert_JSON
+
+
+def create_model_spec(collection, model):
+    kind = model.get('type', 'glm').lower()
+    SpecCls = {
+        'glm': GLMMSpec
+    }[kind]
+    return SpecCls.from_collection(collection, model)
+
+
+class ModelSpec(metaclass=ABCMeta):
+    """Base class for all ModelSpec classes."""
+    @abstractmethod
+    def from_collection(self, collection, model):
+        """Initialize from a BIDSVariableCollection instance."""
+        pass
+
+
+class GLMMSpec(ModelSpec):
+    """Generalized Linear Mixed Model specification.
+s
+    Parameters
+    ----------
+    terms : list of Term
+        A list of Term instances to include in the GLMMSpec instance.
+    X: pd.DataFrame
+        A pandas DataFrame containing the fixed effect design matrix
+        (i.e., the X matrix in the typical mixed effect formulation). Each
+        column will be internally converted to a separate Term instance.
+    Z: pd.DataFrame
+        A pandas DataFrame containing the random effect/grouping matrix
+        (i.e., the Z matrix in the typical mixed effect formulation). Columns
+        that share variance components are identified by the groups argument.
+    groups: NDArray
+        A binary 2d array with dimension k x v, where k is the number of
+        columns in Z and v is the number of distinct variance components. A
+        value of 1 indicates that the i'th of k rows is a level in the j's of
+        v variance components. If Z is passed and groups is None, it is assumed
+        that all columns in Z share the same single variance.
+    sigma: NDArray
+        A 2d array giving the covariance matrix for the variance components
+        defined in the groups argument. Has dimension v x v, where v is the
+        number of columns in groups. If None (default), no constraint is
+        imposed and the covariance is directly estimated.
+    family: str
+        The name of the family to use for the error distribution. By default,
+        gaussian.
+    link: str
+        The name of the link function to use. Default depends on family. In the
+        case of a gaussian (default family), an identity link is used.
+    priors: dict
+        Optional specification of default priors to use for new terms.
+    """
+    def __init__(self, terms=None, X=None, Z=None, groups=None, sigma=None,
+                 family=None, link=None, priors=None):
+        self.terms = {}
+        self.family = family
+        self.link = link
+        self.sigma = sigma
+
+        if priors is not None:
+            self.set_priors(priors)
+
+        if terms is not None:
+            for t in terms:
+                self.add_term(t)
+        if X is not None:
+            self.build_fixed_terms(X)
+        if Z is not None:
+            self.build_variance_components(Z, groups, sigma)
+
+    def set_priors(self, fixed=None, random=None):
+        pass
+
+    def build_fixed_terms(self, X):
+        """Build one or more fixed terms from the columns of a pandas DF.
+
+        Parameters
+        ----------
+        X : pd.DataFrame
+            A pandas DataFrame containing variables to convert to Term
+            instances. Each column is converted to a different (fixed) Term,
+            with the name taken from the column name.
+        """
+        for col in X.columns:
+            data = X.loc[:, col].values
+            cat = data.dtype.name in ('str', 'category', 'object')
+            # TODO: get default prior
+            t = Term(col, data, categorical=cat)
+            self.add_term(t)
+
+    def build_variance_components(self, Z, groups=None, sigma=None, names=None):
+        """Build one or more variance components from the columns of a binary
+        grouping matrix and variance specification.
+
+        Arguments:
+            Z (DataFrame, NDArray): A binary 2D array or pandas DataFrame. Each
+                column represents a column/predictor, each row represents an
+                observation.
+            groups (2DArray): A 2D binary array that maps the columns of Z
+                onto variance components. Has dimension n_rows(Z) x k,
+                where k is the number of distinct variance components. If None,
+                a single group over all columns of Z is assumed.
+            sigma (2DArray): A k x k 2D covariance matrix specifying the
+                covariances between variance components.
+            names (list): Optional list specifying the names of the groups. 
+        """
+        if groups is None:
+            groups = np.ones((Z.shape[1], 1))
+        n_grps = groups.shape[1]
+
+        if names is None:
+            names = getattr(groups, 'columns',
+                            ['VC{}'.format(i) for i in range(n_grps)])
+
+        # Work with array instead of DF
+        if hasattr(groups, 'values'):
+            groups = groups.values
+
+        for i in range(n_grps):
+            z_grp = Z[:, groups[:, i].astype(bool)]
+            # TODO: select default prior
+            vc = VarComp(names[i], z_grp)
+            self.add_term(vc)
+
+    def add_term(self, term):
+        """Add a new Term to the instance.
+
+        Parameters
+        ----------
+        term : Term
+            A Term instance to add to the current instance.
+        """
+        if term.name in self.terms:
+            raise ValueError("Term with name {} already exists!"
+                             .format(term.name))
+        self.terms[term.name] = term
+
+    @property
+    def X(self):
+        """Return X design matrix (i.e., fixed component of model)."""
+        if not self.fixed_terms:
+            return None
+        names, cols = zip(*[(c.name, c.values) for c in self.fixed_terms])
+        return pd.DataFrame(np.c_[cols], columns=names)
+
+    @property
+    def Z(self):
+        """Return Z design matrix (i.e., random effects/variance components).
+        """
+        if not self.variance_components:
+            return None
+        names, cols = [], []
+        for c in self.variance_components:
+            cols.append(c.values)
+            names.extend(['{}.{}'.format(c.name, i)
+                          for i in range(c.values.shape[1])])
+        return pd.DataFrame(np.concatenate(cols, axis=1), columns=names)
+
+    @property
+    def fixed_terms(self):
+        """Return a list of all available fixed effects."""
+        return [t for t in self.terms.values() if not isinstance(t, VarComp)]
+
+    @property
+    def variance_components(self):
+        """Return a list of all available variance components."""
+        return [t for t in self.terms.values() if isinstance(t, VarComp)]
+
+    @classmethod
+    def from_collection(cls, collection, model):
+        """ Initialize a GLMMSpec instance from a BIDSVariableCollection and
+        a BIDS-StatsModels JSON spec.
+
+        Parameters
+        ----------
+        collection : BIDSVariableCollection
+            A BIDSVariableCollection containing variable information.
+        model : dict
+            The "Model" section from a BIDS-StatsModel specification.
+
+        Returns
+        -------
+        A GLMMSpec instance.
+        """
+        if isinstance(collection, BIDSRunVariableCollection):
+            if not collection.all_dense():
+                raise ValueError("Input BIDSRunVariableCollection contains at "
+                                 "least one sparse variable. All variables must"
+                                 " be dense!")
+
+        kwargs = {}
+
+        # Fixed terms
+        model = convert_JSON(model)
+        names = model.get('x', [])
+        if names:
+            names = collection.match_variables(names)
+            X = collection.to_df(names).loc[:, names]
+            kwargs['X'] = X
+
+        # Variance components
+        vcs = model.get('variance_components', [])
+        Z_list = []
+        if vcs:
+            for vc in vcs:
+                # Levels can either be defined by the levels of a single
+                # categorical ("LevelsFrom") or by a set of binary variables.
+                if 'levels_from' in vc:
+                    data = collection.variables[vc['levels_from']].values
+                    Z_list.append(pd.get_dummies(data).values)
+                else:
+                    names = collection.match_variables(vc['levels'])
+                    df = collection.to_df(names).loc[:, names]
+                    Z_list.append(df.values)
+
+            Z = np.concatenate(Z_list, axis=1)
+            groups = np.zeros((Z.shape[1], len(Z_list)))
+            c = 0
+            for i, vc in enumerate(Z_list):
+                n = vc.shape[1]
+                groups[c:(c+n), i] = 1
+                c += n
+            groups = pd.DataFrame(groups, columns=[vc['name'] for vc in vcs])
+
+            kwargs['Z'] = Z
+            kwargs['groups'] = groups
+
+        error = model.get('error')
+        if error:
+            kwargs['family'] = error.get('family')
+            kwargs['link'] = error.get('link')
+
+        return GLMMSpec(**kwargs)
+
+
+class Term(object):
+    """Represents a model term.
+
+    Parameters
+    ----------
+    name : str
+        The name of the term.
+    values : iterable
+        A 1d array or other iterable containing the predictor values.
+    categorical : bool
+        Indicates whether or not the Term represents a categorical variable.
+    prior : dict
+        Optional specification of the prior distribution for the Term.
+    metadata : dict
+        Arbitrary metadata to store internally.
+    """
+    def __init__(self, name, values, categorical=False, prior=None,
+                 metadata=None):
+        self.name = name
+        self.values = values
+        self.categorical = categorical
+        self.prior = prior
+        self.metadata = metadata or {}
+
+
+class VarComp(Term):
+    """Represents a variance component/random effect.
+
+    Parameters
+    ----------
+    name : str
+        The name of the variance component.
+    values : iterable
+        A 2d binary array identifying the observations that belong to the
+        levels of the variance component. Has dimension n x k, where n is the
+        number of observed rows in the dataset and k is the number of levels
+        in the factor.
+    prior : dict
+        Optional specification of the prior distribution for the VarComp.
+    metadata : dict
+        Arbitrary metadata to store internally.
+    """
+    def __init__(self, name, values, prior=None, metadata=None):
+        super(VarComp, self).__init__(name, values, categorical=True,
+                                      prior=prior, metadata=metadata)
+        self.index_vec = self.dummies_to_vec(values)
+
+    @staticmethod
+    def dummies_to_vec(dummies):
+        """Convert dummy-coded columns to a single integer index.
+
+        Parameters
+        ----------
+        dummies : NDArray
+            2d binary array to recode as a single vector.
+
+        Notes
+        -----
+        Used for the sake of computational efficiency (i.e., to avoid lots of
+        large matrix multiplications in the backends), invert the dummy-coding
+        process and represent full-rank dummies as a vector of indices into the
+        coefficients.
+        """
+        vec = np.zeros(len(dummies), dtype=int)
+        for i in range(dummies.shape[1]):
+            vec[(dummies[:, i] == 1)] = i + 1
+        return vec
+
+
+class Prior(object):
+    '''Abstract specification of a term prior.
+
+    Parameters
+    ----------
+    name : str
+        Name of prior distribution (e.g., Normal, Bernoulli, etc.)
+    kwargs: dict
+        Optional keywords specifying the parameters of the named distribution.
+
+    Notes
+    -----
+    At present there's no controlled vocabulary of supported prior names and
+    arguments, but users implementing new Bayesian estimators are encouraged to
+    use the names used in PyMC3 (e.g., 'Normal', parameterized with mu and
+    sd arguments).
+    '''
+    def __init__(self, name, **kwargs):
+        self.name = name
+        self.kwargs = kwargs
diff --git a/bids/analysis/transformations/__init__.py b/bids/analysis/transformations/__init__.py
index 6d820bf7..7a20ad53 100644
--- a/bids/analysis/transformations/__init__.py
+++ b/bids/analysis/transformations/__init__.py
@@ -2,6 +2,7 @@ from .compute import (Sum, Product, Scale, Orthogonalize, Threshold, And, Or,
                       Not, Demean, Convolve)
 from .munge import (Split, Rename, Assign, Copy, Factor, Filter, Select,
                     Delete, DropNA, Replace, ToDense, Group, Resample)
+from .base import TransformerManager
 
 __all__ = [
     'And',
@@ -26,5 +27,6 @@ __all__ = [
     'Split',
     'Sum',
     'Threshold',
-    'ToDense'
+    'ToDense',
+    'TransformerManager'
 ]
diff --git a/bids/analysis/transformations/base.py b/bids/analysis/transformations/base.py
index 03257db7..b1e1e1ea 100644
--- a/bids/analysis/transformations/base.py
+++ b/bids/analysis/transformations/base.py
@@ -10,8 +10,9 @@ import inspect
 import numpy as np
 import pandas as pd
 
-from bids.utils import listify
+from bids.utils import listify, convert_JSON
 from bids.variables import SparseRunVariable
+from bids.analysis import transformations as pbt
 
 
 class Transformation(metaclass=ABCMeta):
@@ -372,3 +373,77 @@ class Transformation(metaclass=ABCMeta):
                 _align([c] + align_variables)
         else:
             _align(listify(variables) + align_variables)
+
+
+class TransformerManager(object):
+    """Handles registration and application of transformations to
+    BIDSVariableCollections.
+
+    Parameters
+    ----------
+    default: object
+        A module or other object containing default transformations as
+            attributes. Any named transformation not explicitly registered on
+            the TransformerManager instance is expected to be found here.
+            If None, the PyBIDS transformations module is used.
+    """
+
+    def __init__(self, default=None):
+        self.transformations = {}
+        if default is None:
+            # Default to PyBIDS transformations
+            default = pbt
+        self.default = default
+
+    def _sanitize_name(self, name):
+        """ Replace any invalid/reserved transformation names with acceptable
+        equivalents.
+
+        Parameters
+        ----------
+        name: str
+            The name of the transformation to sanitize.
+        """
+        if name in ('And', 'Or'):
+            name += '_'
+        return name
+
+    def register(self, name, func):
+        """Register a new transformation handler.
+
+        Parameters
+        ----------
+        name : str
+            The name of the transformation to handle.
+        func : callable
+            The callable to invoke when the named transformation is applied.
+        """
+        name = self._sanitize_name(name)
+        self.transformations[name] = func
+
+    def transform(self, collection, transformations):
+        """Apply all transformations to the variables in the collection.
+
+        Parameters
+        ----------
+        collection: BIDSVariableCollection
+            The BIDSVariableCollection containing variables to transform.
+        transformations : list
+            List of transformations to apply.
+        """
+        for t in transformations:
+            t = convert_JSON(t) # make sure all keys are snake case
+            kwargs = dict(t)
+            name = self._sanitize_name(kwargs.pop('name'))
+            cols = kwargs.pop('input', None)
+
+            # Check registered transformations; fall back on default module
+            func = self.transformations.get(name, None)
+            if func is None:
+                if not hasattr(self.default, name):
+                    raise ValueError("No transformation '%s' found: either "
+                                     "explicitly register a handler, or pass a"
+                                     " default module that supports it." % name)
+                func = getattr(self.default, name)
+                func(collection, cols, **kwargs)
+        return collection
diff --git a/bids/analysis/transformations/munge.py b/bids/analysis/transformations/munge.py
index cd4c1c1f..38c4c7d6 100644
--- a/bids/analysis/transformations/munge.py
+++ b/bids/analysis/transformations/munge.py
@@ -235,8 +235,8 @@ class Rename(Transformation):
     _allow_categorical = ('variables',)
 
     def _transform(self, var):
-        """Rename happens automatically in the base class, so all we need to
-        do is unset the original variable in the collection. """
+        # Rename happens automatically in the base class, so all we need to
+        # do is unset the original variable in the collection.
         self.collection.variables.pop(var.name)
         return var.values
 
diff --git a/bids/layout/layout.py b/bids/layout/layout.py
index 6a2bec56..e5f90829 100644
--- a/bids/layout/layout.py
+++ b/bids/layout/layout.py
@@ -1072,8 +1072,8 @@ class BIDSLayout(object):
             observation is handled separately, and the result is returned
             as a list.
         sampling_rate : int or str
-            If level='run', the sampling rate to pass on to the returned
-            :obj:`bids.variables.kollekshuns.BIDSRunVariableCollection`.
+            If level='run', the sampling rate to pass onto the returned
+            :obj:`bids.variables.collections.BIDSRunVariableCollection`.
         skip_empty : bool
             Whether or not to skip empty Variables (i.e., where there are no
             rows/records in a file after applying any filtering operations
@@ -1084,10 +1084,10 @@ class BIDSLayout(object):
 
         Returns
         -------
-        list of :obj:`bids.variables.kollekshuns.BIDSVariableCollection`
-            or :obj:`bids.variables.kollekshuns.BIDSVariableCollection`
+        list of :obj:`bids.variables.collections.BIDSVariableCollection`
+            or :obj:`bids.variables.collections.BIDSVariableCollection`
             A list if merge=False;
-            a single :obj:`bids.variables.kollekshuns.BIDSVariableCollection`
+            a single :obj:`bids.variables.collections.BIDSVariableCollection`
             if merge=True.
 
         """
diff --git a/bids/variables/__init__.py b/bids/variables/__init__.py
index 039582b5..f6f6c217 100644
--- a/bids/variables/__init__.py
+++ b/bids/variables/__init__.py
@@ -1,6 +1,6 @@
 from .variables import (SimpleVariable, SparseRunVariable,
                         DenseRunVariable, merge_variables)
-from .kollekshuns import (BIDSVariableCollection, BIDSRunVariableCollection,
+from .collections import (BIDSVariableCollection, BIDSRunVariableCollection,
                           merge_collections)
 from .io import load_variables
 
diff --git a/bids/variables/collections.py b/bids/variables/collections.py
new file mode 100644
index 00000000..81dc0e32
--- /dev/null
+++ b/bids/variables/collections.py
@@ -0,0 +1,591 @@
+"""Classes and functions related to the management of sets of BIDSVariables."""
+from copy import copy
+import warnings
+import re
+from collections import OrderedDict
+from itertools import chain
+import fnmatch
+
+import numpy as np
+import pandas as pd
+from pandas.api.types import is_numeric_dtype
+
+from .variables import (SparseRunVariable, SimpleVariable, DenseRunVariable,
+                        merge_variables, BIDSVariable)
+from bids.utils import listify, matches_entities
+
+
+class BIDSVariableCollection(object):
+    """A container for one or more variables extracted from variable files
+    at a single level of analysis.
+
+    Parameters
+    ----------
+    variables : list
+        A list of BIDSVariables or SimpleVariables.
+
+    Notes
+    -----
+    Variables in the list must all share the same analysis level, which
+    must be one of 'session', 'subject', or 'dataset' level. For
+    run-level Variables, use the BIDSRunVariableCollection.
+    """
+
+    def __init__(self, variables):
+
+        if not variables:
+            raise ValueError("No variables were provided")
+        SOURCE_TO_LEVEL = {
+            'events': 'run',
+            'physio': 'run',
+            'stim': 'run',
+            'regressors': 'run',
+            'scans': 'session',
+            'sessions': 'subject',
+            'participants': 'dataset'
+        }
+        var_levels = set([SOURCE_TO_LEVEL[v.source] if v.source in
+                          SOURCE_TO_LEVEL else v.source for v in variables])
+
+        # TODO: relax this requirement & allow implicit merging between levels
+        if len(var_levels) > 1:
+            raise ValueError("A Collection cannot be initialized from "
+                             "variables at more than one level of analysis. "
+                             "Levels found in input variables: %s" %
+                             var_levels)
+        elif not var_levels:
+            raise ValueError(
+                "None of the provided variables matched any of the known levels, which are: %s"
+                % (', '.join(sorted(SOURCE_TO_LEVEL.values())))
+            )
+
+        self.level = list(var_levels)[0]
+        variables = self.merge_variables(variables)
+        self.variables = {v.name: v for v in variables}
+        self._index_entities()
+
+        # Container for variable groups (see BIDS-StatsModels spec)--maps from
+        # group names to lists of variables.
+        self.groups = {}
+
+    @staticmethod
+    def merge_variables(variables, **kwargs):
+        """Concatenates Variables along row axis.
+
+        Parameters
+        ----------
+        variables : list
+            List of Variables to merge. Variables can have
+            different names (and all Variables that share a name will be
+            concatenated together).
+
+        Returns
+        -------
+        list
+            A list of Variables.
+        """
+        var_dict = OrderedDict()
+        for v in variables:
+            if v.name not in var_dict:
+                var_dict[v.name] = []
+            var_dict[v.name].append(v)
+        return [merge_variables(vars_, **kwargs)
+                for vars_ in list(var_dict.values())]
+
+    def to_df(self, variables=None, format='wide', fillna=np.nan,
+              entities=True, timing=True):
+        """Merge variables into a single pandas DataFrame.
+
+        Parameters
+        ----------
+        variables : list of str or BIDSVariable
+            Optional list of variables or variable names to retain. If strings
+            are passed, each one gives the name of a variable in the current
+            collection. If BIDSVariables are passed, they will be used as-is.
+            If None, all variables are returned. Strings and BIDSVariables
+            cannot be mixed in the list.
+        format : {'wide', 'long'}
+            Whether to return a DataFrame in 'wide' or 'long' format. In 'wide'
+            format, each row is defined by a unique entity combination, and
+            each variable is in a separate column. In 'long' format, each row
+            is a unique combination of entities and variable names, and a
+            single 'amplitude' column provides the value.
+        fillna : value
+            Replace missing values with the specified value.
+        entities : bool
+            Whether or not to include a column for each entity.
+        timing : bool
+            Whether or not to include onset and duration columns.
+
+        Returns
+        -------
+        :obj:`pandas.DataFrame`
+            A pandas DataFrame.
+        """
+
+        if variables is None:
+            variables = list(self.variables.keys())
+
+        # Can receive already-selected Variables from sub-classes
+        if not isinstance(variables[0], BIDSVariable):
+            variables = [v for v in self.variables.values()
+                         if v.name in variables]
+
+        # Convert all variables to separate DFs.
+        # Note: bad things can happen if we pass the conditions, entities, and
+        # timing flags through to the individual variables and then do
+        # concat/reshaping operations. So instead, we set them all to True
+        # temporarily, do what we need to, then drop them later if needed.
+        dfs = [v.to_df(True, True, timing=True) for v in variables]
+
+        # Always concatenate along row axis (for format='wide', we'll pivot).
+        df = pd.concat(dfs, axis=0, sort=True)
+
+        ind_cols = list(set(df.columns) - {'condition', 'amplitude'})
+        ent_cols = list(set(ind_cols) - {'onset', 'duration'})
+
+        if format == 'long':
+            df = df.reset_index(drop=True).fillna(fillna)
+        else:
+            df['amplitude'] = df['amplitude'].fillna('n/a')
+            df = df.pivot_table(index=ind_cols, columns='condition',
+                                values='amplitude', aggfunc='first')
+            df = df.reset_index().replace('n/a', fillna)
+            df.columns.name = None
+
+        # Drop any columns we don't want
+        if not timing:
+            df.drop(['onset', 'duration'], axis=1, inplace=True)
+        if not entities:
+            df.drop(ent_cols, axis=1, inplace=True)
+
+        return df
+
+    @classmethod
+    def from_df(cls, data, entities=None, source='contrast'):
+        """Create a Collection from a pandas DataFrame.
+
+        Parameters
+        ----------
+        df : :obj:`pandas.DataFrame`
+            The DataFrame to convert to a Collection. Each
+            column will be converted to a SimpleVariable.
+        entities : :obj:`pandas.DataFrame`
+            An optional second DataFrame containing
+            entity information.
+        source : str
+            The value to set as the source for all Variables.
+
+        Returns
+        -------
+        BIDSVariableCollection
+        """
+        variables = []
+        for col in data.columns:
+            _data = pd.DataFrame(data[col].values, columns=['amplitude'])
+            if entities is not None:
+                _data = pd.concat([_data, entities], axis=1, sort=True)
+            variables.append(SimpleVariable(name=col, data=_data, source=source))
+        return BIDSVariableCollection(variables)
+
+    def clone(self):
+        """Returns a shallow copy of the current instance, except that all
+        variables are deep-cloned.
+        """
+        clone = copy(self)
+        clone.variables = {k: v.clone() for (k, v) in self.variables.items()}
+        return clone
+
+
+    def _index_entities(self):
+        """Sets current instance's entities based on the existing index.
+
+        Notes
+        -----
+        Only entity key/value pairs common to all rows in all contained
+        Variables are returned. E.g., if a Collection contains Variables
+        extracted from runs 1, 2 and 3 from subject '01', the returned dict
+        will be {'subject': '01'}; the runs will be excluded as they vary
+        across the Collection contents.
+        """
+        all_ents = pd.DataFrame.from_records(
+            [v.entities for v in self.variables.values()])
+        constant = all_ents.apply(lambda x: x.nunique() == 1)
+        if constant.empty:
+            self.entities = {}
+        else:
+            keep = all_ents.columns[constant]
+            ents = {k: all_ents[k].dropna().iloc[0] for k in keep}
+            self.entities = {k: v for k, v in ents.items() if pd.notnull(v)}
+
+    def __getitem__(self, var):
+        if var in self.variables:
+            return self.variables[var]
+        keys = list(self.variables.keys())
+        raise ValueError("No variable named '{}' found in this collection. "
+                         "Available names are {}.".format(var, keys))
+
+    def __setitem__(self, var, obj):
+        # Ensure name matches collection key, but raise warning if needed.
+        if obj.name != var:
+            warnings.warn("The provided key to use in the collection ('%s') "
+                          "does not match the passed Column object's existing "
+                          "name ('%s'). The Column name will be set to match "
+                          "the provided key." % (var, obj.name))
+            obj.name = var
+        self.variables[var] = obj
+
+    def match_variables(self, pattern, return_type='name', match_type='unix'):
+        """Return columns whose names match the provided pattern.
+
+        Parameters
+        ----------
+        pattern : str, list
+            One or more regex patterns to match all variable names against.
+        return_type : {'name', 'variable'}
+            What to return. Must be one of:
+            'name': Returns a list of names of matching variables.
+            'variable': Returns a list of Variable objects whose names
+            match.
+        match_type : str
+            Matching approach to use. Either 'regex' (full-blown regular
+                expression matching) or 'unix' (unix-style pattern matching
+                via the fnmatch module).
+
+        Returns
+        -------
+        A list of all matching variables or variable names
+        """
+        pattern = listify(pattern)
+        results = []
+        for patt in pattern:
+            if match_type.lower().startswith('re'):
+                patt = re.compile(patt)
+                vars_ = [v for v in self.variables.keys() if patt.search(v)]
+            else:
+                vars_ = fnmatch.filter(list(self.variables.keys()), patt)
+            if return_type.startswith('var'):
+                vars_ = [self.variables[v] for v in vars_]
+            results.extend(vars_)
+        return results
+
+
+class BIDSRunVariableCollection(BIDSVariableCollection):
+    """A container for one or more RunVariables--i.e., Variables that have a
+    temporal dimension.
+
+    Parameters
+    ----------
+    variables : list
+        A list of SparseRunVariable and/or DenseRunVariable.
+    sampling_rate : float
+        Sampling rate (in Hz) to use when working with
+        dense representations of variables. If None, defaults to 10.
+
+    Notes
+    -----
+    Variables in the list must all be at the 'run' level. For other
+    levels (session, subject, or dataset), use the
+    BIDSVariableCollection.
+    """
+
+    def __init__(self, variables, sampling_rate=None):
+        # Don't put the default value in signature because None is passed from
+        # several places and we don't want multiple conflicting defaults.
+        if sampling_rate:
+            if isinstance(sampling_rate, str):
+                raise ValueError("Sampling rate must be numeric.")
+        self.sampling_rate = sampling_rate or 10
+        super(BIDSRunVariableCollection, self).__init__(variables)
+
+    def get_dense_variables(self, variables=None):
+        """Returns a list of all stored DenseRunVariables."""
+        if variables is None:
+            variables = set(self.variables.keys())
+        return [v for v in self.variables.values()
+                if isinstance(v, DenseRunVariable) and v.name in variables]
+
+    def get_sparse_variables(self, variables=None):
+        """Returns a list of all stored SparseRunVariables."""
+        if variables is None:
+            variables = set(self.variables.keys())
+        return [v for v in self.variables.values()
+                if isinstance(v, SparseRunVariable) and v.name in variables]
+
+    def all_dense(self):
+        return len(self.get_dense_variables()) == len(self.variables)
+
+    def all_sparse(self):
+        return len(self.get_sparse_variables()) == len(self.variables)
+
+    def _get_sampling_rate(self, sampling_rate):
+        """Parse sampling rate argument and return appropriate value."""
+
+        if sampling_rate is None:
+            return self.sampling_rate
+
+        if isinstance(sampling_rate, (float, int)):
+            return sampling_rate
+
+        if sampling_rate == 'TR':
+            trs = {var.run_info[0].tr for var in self.variables.values()}
+            if not trs:
+                raise ValueError("Repetition time unavailable; specify "
+                                    "sampling_rate in Hz explicitly or set to"
+                                    " 'highest'.")
+            elif len(trs) > 1:
+                raise ValueError("Non-unique Repetition times found "
+                                    "({!r}); specify sampling_rate explicitly"
+                                    .format(trs))
+            return 1. / trs.pop()
+
+        if sampling_rate.lower() == 'highest':
+            dense_vars = self.get_dense_variables()
+            # If no dense variables are available, fall back on instance SR
+            if not dense_vars:
+                return self.sampling_rate
+            var_srs = [v.sampling_rate for v in dense_vars]
+            if len(var_srs) == 1:
+                return var_srs[0]
+            return max(*var_srs)
+
+        raise ValueError("Invalid sampling_rate value '{}' provided. Must be "
+                         "a float, None, 'TR', or 'highest'."
+                         .format(sampling_rate))
+
+    def _densify_and_resample(self, sampling_rate=None, variables=None,
+                              resample_dense=False, force_dense=False,
+                              in_place=False, kind='linear'):
+
+        sampling_rate = self._get_sampling_rate(sampling_rate)
+
+        _dense, _sparse = [], []
+
+        # Filter variables and sort by class
+        for name, var in self.variables.items():
+            if variables is not None and name not in variables:
+                continue
+            if isinstance(var, DenseRunVariable):
+                _dense.append(var)
+            else:
+                _sparse.append(var)
+
+        _variables = {}
+
+        if force_dense:
+            for v in _sparse:
+                if is_numeric_dtype(v.values):
+                    _variables[v.name] = v.to_dense(sampling_rate)
+
+        if resample_dense:
+            for v in _dense:
+                _variables[v.name] = v.resample(sampling_rate, kind=kind)
+
+        coll = self if in_place else self.clone()
+
+        if in_place:
+            coll.variables.update(_variables)
+        else:
+            coll.variables = _variables
+
+        coll.sampling_rate = sampling_rate
+        return coll
+
+    def to_dense(self, sampling_rate=None, variables=None, in_place=False,
+                 kind='linear'):
+        """ Convert all contained SparseRunVariables to DenseRunVariables.
+
+        Parameters
+        ----------
+        sampling_rate : None, {'TR', 'highest'}, float
+            Sampling rate to use when densifying sparse variables. If None,
+            uses the currently stored instance value. If 'TR', the repetition
+            time is used, if available, to select the sampling rate (1/TR).
+            If 'highest', all variables are resampled to the highest sampling
+            rate of any of the existing dense variables. The sampling rate may
+            also be specified explicitly in Hz as a float.
+        variables : list
+            Optional list of names of Variables to resample. If None, all
+            variables are resampled.
+        in_place : bool
+            When True, all variables are overwritten in-place.
+            When False, returns resampled versions of all variables.
+        kind : str
+            Argument to pass to scipy's interp1d; indicates the kind of
+            interpolation approach to use. See interp1d docs for valid values.
+
+        Returns
+        -------
+        A BIDSVariableCollection (if in_place is False).
+
+        Notes
+        -----
+        Categorical variables are ignored.
+        """
+        return self._densify_and_resample(sampling_rate, variables,
+                                          resample_dense=False,
+                                          in_place=in_place, kind=kind,
+                                          force_dense=True)
+
+    def resample(self, sampling_rate=None, variables=None, force_dense=False,
+                 in_place=False, kind='linear'):
+        """Resample all dense variables (and optionally, sparse ones) to the
+        specified sampling rate.
+
+        Parameters
+        ----------
+        sampling_rate : int or float
+            Target sampling rate (in Hz). If None, uses the instance value.
+        variables : list
+            Optional list of names of Variables to resample. If None, all
+            variables are resampled.
+        force_dense : bool
+            if True, all sparse variables will be forced to dense.
+        in_place : bool
+            When True, all variables are overwritten in-place.
+            When False, returns resampled versions of all variables.
+        kind : str
+            Argument to pass to scipy's interp1d; indicates the kind of
+            interpolation approach to use. See interp1d docs for valid values.
+
+        Returns
+        -------
+        A BIDSVariableCollection (if in_place is False).
+        """
+        return self._densify_and_resample(sampling_rate, variables,
+                                   force_dense=force_dense, in_place=in_place,
+                                   kind=kind, resample_dense=True)
+
+    def to_df(self, variables=None, format='wide', fillna=np.nan,
+              sampling_rate='highest', include_sparse=True, include_dense=True,
+              entities=True, timing=True):
+        """Merge columns into a single pandas DataFrame.
+
+        Parameters
+        ----------
+        variables : list
+            Optional list of variable names to retain;
+            if None, all variables are written out.
+        format : str
+            Whether to return a DataFrame in 'wide' or 'long'
+            format. In 'wide' format, each row is defined by a unique
+            onset/duration, and each variable is in a separate column. In
+            'long' format, each row is a unique combination of onset,
+            duration, and variable name, and a single 'amplitude' column
+            provides the value.
+        fillna : value
+            Replace missing values with the specified value.
+        sampling_rate : float
+            Specifies the sampling rate to use for all variables in the event
+            that resampling needs to be performed (i.e., if some variables are
+            sparse, or if dense variables have different sampling rates). Must
+            be one of 'TR', 'highest', None, or a float (specifying the rate in
+            Hz). If None, uses the instance sampling rate (10 Hz by default).
+        include_sparse : bool
+            Whether or not to include sparse variables in the output.
+        include_dense : bool
+            Whether or not to include dense variables in the output.
+        entities : bool
+            Whether or not to include a column for each entity.
+        timing : bool
+            Whether or not to include onset and duration columns.
+
+        Returns
+        -------
+        :obj:`pandas.DataFrame`
+            A pandas DataFrame.
+
+        Notes
+        -----
+        The precise format of the resulting DataFrame depends on the variables
+        contained in the current instance. If all variables are sparse, the
+        output will also be sparse--i.e., the events in the DataFrame may have
+        non-uniform timing. If at least one dense variable is present, and the
+        user has not explicitly excluded dense variables (by setting
+        include_dense=False), all selected variables will be implicitly
+        converted to dense using the specified `sampling_rate` (if provided).
+        To avoid unexpected behavior, we recommend converting mixed collections
+        to all-dense form explicitly via the `to_dense()` or `resample()`
+        methods before calling `to_df()`.
+        """
+
+        if not include_sparse and not include_dense:
+            raise ValueError("You can't exclude both dense and sparse "
+                             "variables! That leaves nothing!")
+
+        _vars = []
+
+        if include_sparse:
+            _vars += self.get_sparse_variables(variables)
+
+        if include_dense:
+            _vars += self.get_dense_variables(variables)
+
+        if not _vars:
+            raise ValueError("No variables were selected for output.")
+
+        # If all variables are sparse/simple, we can pass them as-is. Otherwise
+        # we first force all variables to dense via .resample().
+        if all(isinstance(v, SimpleVariable) for v in _vars):
+            variables = _vars
+        else:
+            sampling_rate = sampling_rate or self.sampling_rate
+            var_names = [v.name for v in _vars]
+            collection = self.resample(sampling_rate, variables=var_names,
+                                       force_dense=include_sparse)
+            variables = list(collection.variables.values())
+
+        return super().to_df(variables, format, fillna, entities=entities,
+                             timing=timing)
+
+
+def merge_collections(collections, sampling_rate='highest', output_level=None):
+    """Merge two or more collections at the same level of analysis.
+
+    Parameters
+    ----------
+    collections : list
+        List of Collections to merge.
+    sampling_rate : int or str
+        Sampling rate to use if it becomes necessary
+        to resample DenseRunVariables. Either an integer or 'highest' (see
+        merge_variables docstring for further explanation).
+    output_level : str, optional
+        Assign a new level (e.g., 'run', 'subject', etc.) to the merged
+        collection. If None, the current level is retained.
+
+    Returns
+    -------
+    BIDSVariableCollection or BIDSRunVariableCollection
+        Result type depends on the type of the input collections.
+    """
+    collections = listify(collections)
+    if len(collections) == 1:
+        return collections[0]
+
+    levels = set([c.level for c in collections])
+    if len(levels) > 1:
+        raise ValueError("At the moment, it's only possible to merge "
+                         "Collections at the same level of analysis. You "
+                         "passed collections at levels: %s." % levels)
+
+    variables = list(chain(*[c.variables.values() for c in collections]))
+    cls = collections[0].__class__
+
+    variables = cls.merge_variables(variables, sampling_rate=sampling_rate)
+
+    if isinstance(collections[0], BIDSRunVariableCollection):
+        # 'auto' was renamed to 'highest' circa 0.10, but check for both
+        if sampling_rate in {'auto', 'highest'}:
+            rates = [var.sampling_rate for var in variables
+                     if isinstance(var, DenseRunVariable)]
+
+            sampling_rate = rates[0] if rates else None
+
+        return cls(variables, sampling_rate)
+
+    # For non-run collections, we may need to set a different output level
+    coll = cls(variables)
+    if output_level is not None:
+        coll.level = output_level
+    return coll
diff --git a/bids/variables/entities.py b/bids/variables/entities.py
index f4b27928..b3d77499 100644
--- a/bids/variables/entities.py
+++ b/bids/variables/entities.py
@@ -2,9 +2,12 @@
 
 from itertools import chain
 from collections import namedtuple
-from . import kollekshuns as clc
+
 import pandas as pd
 
+from . import collections as clc
+from bids.utils import matches_entities
+
 
 class Node(object):
     """Base class for objects that represent a single object in the BIDS
@@ -124,7 +127,7 @@ class NodeIndex(object):
 
         for n in nodes:
             var_set = list(n.variables.values())
-            var_set = [v for v in var_set if v.matches_entities(entities)]
+            var_set = [v for v in var_set if matches_entities(v, entities)]
             if names is not None:
                 var_set = [v for v in var_set if v.name in names]
             # Additional filtering on Variables past run level, because their
diff --git a/bids/variables/kollekshuns.py b/bids/variables/kollekshuns.py
deleted file mode 100644
index 190ab8c5..00000000
--- a/bids/variables/kollekshuns.py
+++ /dev/null
@@ -1,455 +0,0 @@
-"""Classes and functions related to the management of sets of BIDSVariables.
-
-Why 'kollekshuns'? Because 'collections' would conflict with the standard lib
-module of the same name on Python 2. We could go with something sensible but
-verbose like 'variable_collections', but that would make far too much sense.
-"""
-
-from copy import copy
-import warnings
-import re
-from collections import OrderedDict
-from itertools import chain
-import fnmatch
-
-import numpy as np
-import pandas as pd
-from pandas.api.types import is_numeric_dtype
-
-from .variables import (SparseRunVariable, SimpleVariable, DenseRunVariable,
-                        merge_variables, BIDSVariable)
-from bids.utils import listify, matches_entities
-
-
-class BIDSVariableCollection(object):
-    """A container for one or more variables extracted from variable files
-    at a single level of analysis.
-
-    Parameters
-    ----------
-    variables : list
-        A list of BIDSVariables or SimpleVariables.
-
-    Notes
-    -----
-    Variables in the list must all share the same analysis level, which
-    must be one of 'session', 'subject', or 'dataset' level. For
-    run-level Variables, use the BIDSRunVariableCollection.
-    """
-
-    def __init__(self, variables):
-
-        if not variables:
-            raise ValueError("No variables were provided")
-        SOURCE_TO_LEVEL = {
-            'events': 'run',
-            'physio': 'run',
-            'stim': 'run',
-            'regressors': 'run',
-            'scans': 'session',
-            'sessions': 'subject',
-            'participants': 'dataset'
-        }
-        var_levels = set([SOURCE_TO_LEVEL[v.source] if v.source in
-                          SOURCE_TO_LEVEL else v.source for v in variables])
-
-        # TODO: relax this requirement & allow implicit merging between levels
-        if len(var_levels) > 1:
-            raise ValueError("A Collection cannot be initialized from "
-                             "variables at more than one level of analysis. "
-                             "Levels found in input variables: %s" %
-                             var_levels)
-        elif not var_levels:
-            raise ValueError(
-                "None of the provided variables matched any of the known levels, which are: %s"
-                % (', '.join(sorted(SOURCE_TO_LEVEL.values())))
-            )
-
-        self.level = list(var_levels)[0]
-        variables = self.merge_variables(variables)
-        self.variables = {v.name: v for v in variables}
-        self._index_entities()
-
-        # Container for variable groups (see BIDS-StatsModels spec)--maps from
-        # group names to lists of variables.
-        self.groups = {}
-
-    @staticmethod
-    def merge_variables(variables, **kwargs):
-        """Concatenates Variables along row axis.
-
-        Parameters
-        ----------
-        variables : list
-            List of Variables to merge. Variables can have
-            different names (and all Variables that share a name will be
-            concatenated together).
-
-        Returns
-        -------
-        list
-            A list of Variables.
-        """
-        var_dict = OrderedDict()
-        for v in variables:
-            if v.name not in var_dict:
-                var_dict[v.name] = []
-            var_dict[v.name].append(v)
-        return [merge_variables(vars_, **kwargs)
-                for vars_ in list(var_dict.values())]
-
-    def to_df(self, variables=None, format='wide', fillna=np.nan, **kwargs):
-        """Merge variables into a single pandas DataFrame.
-
-        Parameters
-        ----------
-        variables : list
-            Optional list of column names to retain; if None,
-            all variables are returned.
-        format : {'wide', 'long'}
-            Whether to return a DataFrame in 'wide' or 'long'
-            format. In 'wide' format, each row is defined by a unique
-            entity combination, and each variable is in a separate column.
-            In 'long' format, each row is a unique combination of entities
-            and variable names, and a single 'amplitude' column provides
-            the value.
-        fillna : value
-            Replace missing values with the specified value.
-        kwargs : dict
-            Optional keyword arguments to pass onto each Variable's
-            to_df() call (e.g., condition, entities, and timing).
-
-        Returns
-        -------
-        :obj:`pandas.DataFrame`
-            A pandas DataFrame.
-        """
-
-        if variables is None:
-            variables = list(self.variables.keys())
-
-        # Can receive already-selected Variables from sub-classes
-        if not isinstance(variables[0], BIDSVariable):
-            variables = [v for v in self.variables.values()
-                         if v.name in variables]
-
-        dfs = [v.to_df(**kwargs) for v in variables]
-        df = pd.concat(dfs, axis=0, sort=True)
-
-        if format == 'long':
-            return df.reset_index(drop=True).fillna(fillna)
-
-        ind_cols = list(set(df.columns) - {'condition', 'amplitude'})
-
-        df['amplitude'] = df['amplitude'].fillna('n/a')
-        df = df.pivot_table(index=ind_cols, columns='condition',
-                            values='amplitude', aggfunc='first')
-        df = df.reset_index().replace('n/a', fillna)
-        df.columns.name = None
-        return df
-
-    @classmethod
-    def from_df(cls, data, entities=None, source='contrast'):
-        """Create a Collection from a pandas DataFrame.
-
-        Parameters
-        ----------
-        data : :obj:`pandas.DataFrame`
-            The DataFrame to convert to a Collection. Each
-            column will be converted to a SimpleVariable.
-        entities : :obj:`pandas.DataFrame`
-            An optional second DataFrame containing
-            entity information.
-        source : str
-            The value to set as the source for all Variables.
-
-        Returns
-        -------
-        BIDSVariableCollection
-        """
-        variables = []
-        for col in data.columns:
-            _data = pd.DataFrame(data[col].values, columns=['amplitude'])
-            if entities is not None:
-                _data = pd.concat([_data, entities], axis=1, sort=True)
-            variables.append(SimpleVariable(name=col, data=_data, source=source))
-        return BIDSVariableCollection(variables)
-
-    def clone(self):
-        """Returns a shallow copy of the current instance, except that all
-        variables are deep-cloned.
-        """
-        clone = copy(self)
-        clone.variables = {k: v.clone() for (k, v) in self.variables.items()}
-        return clone
-
-    def matches_entities(self, entities, strict=False):
-        """Checks whether current Collection's entities match the input. """
-        return matches_entities(self, entities, strict)
-
-    def _index_entities(self):
-        """Sets current instance's entities based on the existing index.
-
-        Notes
-        -----
-        Only entity key/value pairs common to all rows in all contained
-        Variables are returned. E.g., if a Collection contains Variables
-        extracted from runs 1, 2 and 3 from subject '01', the returned dict
-        will be {'subject': '01'}; the runs will be excluded as they vary
-        across the Collection contents.
-        """
-        all_ents = pd.DataFrame.from_records(
-            [v.entities for v in self.variables.values()])
-        constant = all_ents.apply(lambda x: x.nunique() == 1)
-        if constant.empty:
-            self.entities = {}
-        else:
-            keep = all_ents.columns[constant]
-            ents = {k: all_ents[k].dropna().iloc[0] for k in keep}
-            self.entities = {k: v for k, v in ents.items() if pd.notnull(v)}
-
-    def __getitem__(self, var):
-        return self.variables[var]
-
-    def __setitem__(self, var, obj):
-        # Ensure name matches collection key, but raise warning if needed.
-        if obj.name != var:
-            warnings.warn("The provided key to use in the collection ('%s') "
-                          "does not match the passed Column object's existing "
-                          "name ('%s'). The Column name will be set to match "
-                          "the provided key." % (var, obj.name))
-            obj.name = var
-        self.variables[var] = obj
-
-    def match_variables(self, pattern, return_type='name', match_type='unix'):
-        """Return columns whose names match the provided pattern.
-
-        Parameters
-        ----------
-        pattern : str
-            A regex pattern to match all variable names against.
-        return_type : {'name', 'variable'}
-            What to return. Must be one of:
-            'name': Returns a list of names of matching variables.
-            'variable': Returns a list of Variable objects whose names
-            match.
-        match_type : str
-            Matching approach to use. Either 'regex' (full-blown regular
-                expression matching) or 'unix' (unix-style pattern matching
-                via the fnmatch module).
-        """
-        if match_type.lower().startswith('re'):
-            pattern = re.compile(pattern)
-            vars_ = [v for v in self.variables.keys() if pattern.search(v)]
-        else:
-            vars_ = fnmatch.filter(list(self.variables.keys()), pattern)
-
-        return vars_ if return_type == 'name' \
-            else [self.variables[v] for v in vars_]
-
-
-class BIDSRunVariableCollection(BIDSVariableCollection):
-    """A container for one or more RunVariables--i.e., Variables that have a
-    temporal dimension.
-
-    Parameters
-    ----------
-    variables : list
-        A list of SparseRunVariable and/or DenseRunVariable.
-    sampling_rate : float
-        Sampling rate (in Hz) to use when working with
-        dense representations of variables. If None, defaults to 10.
-
-    Notes
-    -----
-    Variables in the list must all be at the 'run' level. For other
-    levels (session, subject, or dataset), use the
-    BIDSVariableCollection.
-    """
-
-    def __init__(self, variables, sampling_rate=None):
-        # Don't put the default value in signature because None is passed from
-        # several places and we don't want multiple conflicting defaults.
-        if sampling_rate:
-            if isinstance(sampling_rate, str):
-                raise ValueError("Sampling rate must be numeric.")
-        self.sampling_rate = sampling_rate or 10
-        super(BIDSRunVariableCollection, self).__init__(variables)
-
-    def _none_dense(self):
-        return all([isinstance(v, SimpleVariable)
-                    for v in self.variables.values()])
-
-    def _all_dense(self):
-        return all([isinstance(v, DenseRunVariable)
-                    for v in self.variables.values()])
-
-    def resample(self, sampling_rate=None, variables=None, force_dense=False,
-                 in_place=False, kind='linear'):
-        """Resample all dense variables (and optionally, sparse ones) to the
-        specified sampling rate.
-
-        Parameters
-        ----------
-        sampling_rate : int or float
-            Target sampling rate (in Hz). If None,
-            uses the instance sampling rate.
-        variables : list
-            Optional list of Variables to resample. If None,
-            all variables are resampled.
-        force_dense : bool
-            if True, all sparse variables will be forced to
-            dense.
-        in_place : bool
-            When True, all variables are overwritten in-place.
-            When False, returns resampled versions of all variables.
-        kind : str
-            Argument to pass to scipy's interp1d; indicates the
-            kind of interpolation approach to use. See interp1d docs for
-            valid values.
-        """
-
-        # Store old sampling rate-based variables
-        sampling_rate = sampling_rate or self.sampling_rate
-
-        _variables = {}
-
-        for name, var in self.variables.items():
-            if variables is not None and name not in variables:
-                continue
-            if isinstance(var, SparseRunVariable):
-                if force_dense and is_numeric_dtype(var.values):
-                    var = var.to_dense()
-                    _variables[name] = var
-                else:
-                    continue
-
-            _var = var.resample(sampling_rate,
-                                inplace=in_place,
-                                kind=kind)
-            if not in_place:  # None if in_place; no update needed
-                _variables[name] = _var
-
-        if in_place:  # Replace densified variables
-            for k, v in _variables.items():
-                self.variables[k] = v
-            self.sampling_rate = sampling_rate
-        else:
-            return _variables
-
-    def to_df(self, variables=None, format='wide', sparse=True,
-              sampling_rate=None, include_sparse=True, include_dense=True,
-              **kwargs):
-        """Merge columns into a single pandas DataFrame.
-
-        Parameters
-        ----------
-        variables : list
-            Optional list of variable names to retain;
-            if None, all variables are written out.
-        format : str
-            Whether to return a DataFrame in 'wide' or 'long'
-            format. In 'wide' format, each row is defined by a unique
-            onset/duration, and each variable is in a separate column. In
-            'long' format, each row is a unique combination of onset,
-            duration, and variable name, and a single 'amplitude' column
-            provides the value.
-        sparse : bool
-            If True, variables will be kept in a sparse
-            format provided they are all internally represented as such.
-            If False, a dense matrix (i.e., uniform sampling rate for all
-            events) will be exported. Will be ignored if at least one
-            variable is dense.
-        sampling_rate : float
-            If a dense matrix is written out, the
-            sampling rate (in Hz) to use for downsampling. Defaults to the
-            value currently set in the instance.
-        kwargs : dict
-            Optional keyword arguments to pass onto each Variable's
-            to_df() call (e.g., condition, entities, and timing).
-        include_sparse : bool
-            Whether or not to include sparse Variables.
-        include_dense : bool
-            Whether or not to include dense Variables.
-
-        Returns
-        -------
-        :obj:`pandas.DataFrame`
-            A pandas DataFrame.
-        """
-
-        if not include_sparse and not include_dense:
-            raise ValueError("You can't exclude both dense and sparse "
-                             "variables! That leaves nothing!")
-
-        if variables is None:
-            variables = list(self.variables.keys())
-
-        if not include_sparse:
-            variables = [v for v in variables if
-                         isinstance(self.variables[v], DenseRunVariable)]
-
-        if not include_dense:
-            variables = [v for v in variables if not
-                         isinstance(self.variables[v], DenseRunVariable)]
-
-        if not variables:
-            return None
-
-        _vars = [self.variables[v] for v in variables]
-        if sparse and all(isinstance(v, SimpleVariable) for v in _vars):
-            variables = _vars
-
-        else:
-            sampling_rate = sampling_rate or self.sampling_rate
-
-            # Make sure all variables have the same sampling rate
-            variables = list(self.resample(sampling_rate, variables,
-                                           force_dense=True,
-                                           in_place=False).values())
-
-        return super(BIDSRunVariableCollection, self).to_df(variables, format,
-                                                            **kwargs)
-
-
-def merge_collections(collections, force_dense=False, sampling_rate='auto'):
-    """Merge two or more collections at the same level of analysis.
-
-    Parameters
-    ----------
-    collections : list
-        List of Collections to merge.
-    sampling_rate : int or str
-        Sampling rate to use if it becomes necessary
-        to resample DenseRunVariables. Either an integer or 'auto' (see
-        merge_variables docstring for further explanation).
-
-    Returns
-    -------
-    BIDSVariableCollection or BIDSRunVariableCollection
-        Result type depends on the type of the input collections.
-    """
-    if len(listify(collections)) == 1:
-        return collections
-
-    levels = set([c.level for c in collections])
-    if len(levels) > 1:
-        raise ValueError("At the moment, it's only possible to merge "
-                         "Collections at the same level of analysis. You "
-                         "passed collections at levels: %s." % levels)
-
-    variables = list(chain(*[c.variables.values() for c in collections]))
-    cls = collections[0].__class__
-
-    variables = cls.merge_variables(variables, sampling_rate=sampling_rate)
-
-    if isinstance(collections[0], BIDSRunVariableCollection):
-        if sampling_rate == 'auto':
-            rates = [var.sampling_rate for var in variables
-                     if isinstance(var, DenseRunVariable)]
-
-            sampling_rate = rates[0] if rates else None
-
-        return cls(variables, sampling_rate)
-
-    return cls(variables)
diff --git a/bids/variables/variables.py b/bids/variables/variables.py
index 8bf3dba1..b4afc9f2 100644
--- a/bids/variables/variables.py
+++ b/bids/variables/variables.py
@@ -1,14 +1,16 @@
 """ Classes for representing BIDS variables. """
 
-import numpy as np
-import pandas as pd
 import math
 import warnings
 from copy import deepcopy
 from abc import abstractmethod, ABCMeta
-from bids.utils import listify
 from itertools import chain
 from functools import reduce
+
+import numpy as np
+import pandas as pd
+
+from bids.utils import listify
 from bids.utils import matches_entities
 
 class BIDSVariable(metaclass=ABCMeta):
@@ -23,7 +25,7 @@ class BIDSVariable(metaclass=ABCMeta):
         self.name = name
         self.values = values
         self.source = source
-        self._index_entities()
+        self.entities = self._extract_entities()
 
     def clone(self, data=None, **kwargs):
         """Clone (deep copy) the current column, optionally replacing its
@@ -221,12 +223,8 @@ class BIDSVariable(metaclass=ABCMeta):
 
         return data.reset_index(drop=True)
 
-    def matches_entities(self, entities, strict=False):
-        """Checks whether current Variable's entities match the input. """
-        return matches_entities(self, entities, strict)
-
-    def _index_entities(self):
-        """Returns a dict of entities for the current Variable.
+    def _extract_entities(self):
+        """Returns a dict of all non-varying entities for the current Variable.
 
         Notes
         -----
@@ -238,10 +236,10 @@ class BIDSVariable(metaclass=ABCMeta):
         """
         constant = self.index.apply(lambda x: x.nunique() == 1)
         if constant.empty:
-            self.entities = {}
+            return {}
         else:
             keep = self.index.columns[constant]
-            self.entities = {k: self.index[k].dropna().iloc[0] for k in keep}
+            return {k: self.index[k].dropna().iloc[0] for k in keep}
 
 
 class SimpleVariable(BIDSVariable):
@@ -373,44 +371,65 @@ class SparseRunVariable(SimpleVariable):
         Parameters
         ----------
         sampling_rate : float or None
-            Sampling rate (in Hz) to use when constructing the DenseRunVariable.
+            Sampling rate (in Hz) to use when constructing the DenseRunVariable
 
         Returns
         -------
         DenseRunVariable
         """
-        if sampling_rate is None:
-            # Cast onsets and durations to milliseconds
-            onsets = np.round(self.onset * 1000).astype(int)
-            durations = np.round(self.duration * 1000).astype(int)
-            sampling_rate = 1000. / reduce(math.gcd, [*onsets, *durations])
-
-        duration = int(math.ceil(sampling_rate * self.get_duration()))
+        # Cast onsets and durations to milliseconds
+        onsets = np.round(self.onset * 1000).astype(int)
+        durations = np.round(self.duration * 1000).astype(int)
+        gcd = np.gcd.reduce(np.r_[onsets, durations])
+        bin_sr = 1000. / gcd
+
+        # never use a computed SR smaller than the requested one, because
+        # when events are widely-spaced and timing is very regular, this can
+        # result in a nasty loss of precision in the resampling step.
+        if sampling_rate is not None:
+            bin_sr = max(bin_sr, sampling_rate)
+
+        duration = int(math.ceil(bin_sr * self.get_duration()))
         ts = np.zeros(duration, dtype=self.values.dtype)
 
-        onsets = np.round(self.onset * sampling_rate).astype(int)
-        durations = np.round(self.duration * sampling_rate).astype(int)
+        onsets = np.round(self.onset * bin_sr).astype(int)
+        durations = np.round(self.duration * bin_sr).astype(int)
 
         run_i, start, last_ind = 0, 0, 0
         for i, val in enumerate(self.values.values):
             if onsets[i] < last_ind:
-                start += self.run_info[run_i].duration * sampling_rate
+                start += self.run_info[run_i].duration * bin_sr
                 run_i += 1
             _onset = int(start + onsets[i])
             _offset = int(_onset + durations[i])
             if _onset >= duration:
-                warnings.warn("The onset time of a variable seems to exceed the runs"
-                              "duration, hence runs are incremented by one internally.")
+                warnings.warn("The onset time of a variable seems to exceed "
+                              "the runs duration, hence runs are incremented "
+                              "by one internally.")
             ts[_onset:_offset] = val
             last_ind = onsets[i]
 
         run_info = list(self.run_info)
-        return DenseRunVariable(
+        dense_var = DenseRunVariable(
             name=self.name,
             values=ts,
             run_info=run_info,
             source=self.source,
-            sampling_rate=sampling_rate)
+            sampling_rate=bin_sr)
+
+        if sampling_rate is not None and bin_sr != sampling_rate:
+            dense_var.resample(sampling_rate, inplace=True)
+
+        return dense_var
+
+    def _extract_entities(self):
+        # Get all entities common to all runs. The super method already does
+        # this for entities that show up in filenames, so we just add the
+        # ones that show up in the RunInfo tuples, as those include metadata.
+        ent_items = [run.entities.items() for run in self.run_info]
+        entities = reduce(lambda x, y: x & y, ent_items, ent_items[0])
+        base_ents = super()._extract_entities()
+        return dict(entities, **base_ents)
 
     @classmethod
     def _merge(cls, variables, name, **kwargs):
@@ -513,7 +532,6 @@ class DenseRunVariable(BIDSVariable):
         if sampling_rate == self.sampling_rate:
             return
 
-        old_sr = self.sampling_rate
         n = len(self.index)
 
         self.index = self._build_entity_index(self.run_info, sampling_rate)
@@ -579,8 +597,8 @@ class DenseRunVariable(BIDSVariable):
                     msg = ("Cannot merge DenseRunVariables (%s) with different"
                            " sampling rates (%s). Either specify an integer "
                            "sampling rate to use for all variables, or set "
-                           "sampling_rate='auto' to use the highest sampling "
-                           "rate found." % (name, rates))
+                           "sampling_rate='highest' to use the highest sampling"
+                           " rate found." % (name, rates))
                     raise ValueError(msg)
 
         variables = [v.resample(sampling_rate) for v in variables]
@@ -609,9 +627,10 @@ def merge_variables(variables, name=None, **kwargs):
         Optional keyword arguments to pass onto the class-specific merge() call.
         Possible args:
             - sampling_rate (int, str): The sampling rate to use if resampling
-              of DenseRunVariables is necessary for harmonization. If 'auto',
-              the highest sampling rate found will be used. This argument is
-              only used when passing DenseRunVariables in the variables list.
+              of DenseRunVariables is necessary for harmonization. If
+              'highest', the highest sampling rate found will be used. This
+              argument is only used when passing DenseRunVariables in the
+              variables list.
 
     Returns
     -------
diff --git a/doc/api.rst b/doc/api.rst
index a6c63861..f9452704 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -36,7 +36,6 @@ API Reference
    bids.analysis.auto_model
    bids.analysis.analysis.Step
    bids.analysis.analysis.AnalysisNode
-   bids.analysis.analysis.apply_transformations
 
 .. currentmodule:: bids
 
@@ -73,7 +72,7 @@ API Reference
    bids.variables.merge_variables
    bids.variables.io
    bids.variables.entities
-   bids.variables.kollekshuns
+   bids.variables.collections
    bids.variables.variables
 
 .. currentmodule:: bids


Test patch:
diff --git a/bids/analysis/tests/test_analysis.py b/bids/analysis/tests/test_analysis.py
index 3c990403..0b893a0c 100644
--- a/bids/analysis/tests/test_analysis.py
+++ b/bids/analysis/tests/test_analysis.py
@@ -1,10 +1,13 @@
 from os.path import join
+
+import numpy as np
+import pytest
+
 from bids.analysis import Analysis
-from bids.analysis.analysis import ContrastInfo, DesignMatrixInfo
+from bids.analysis.analysis import ContrastInfo
 from bids.layout import BIDSLayout
 from bids.tests import get_test_data_path
-import numpy as np
-import pytest
+from bids.variables import BIDSVariableCollection
 
 
 @pytest.fixture
@@ -17,92 +20,42 @@ def analysis():
     return analysis
 
 
-def test_design_matrix_info(analysis):
-    result = analysis['run'].get_design_matrix(subject=['01', '02', '03'])
-    for dmi in result:
-        assert isinstance(dmi, DesignMatrixInfo)
-        assert dmi._fields == ('sparse', 'dense', 'entities')
-        assert hasattr(dmi.sparse, 'shape')
-        assert dmi.dense is None
-        assert isinstance(dmi.entities, dict)
-
-
-def test_get_design_matrix_arguments(analysis):
-    kwargs = dict(run=1, subject='01', sparse=True)
-    result = analysis['run'].get_design_matrix(**kwargs)
-    result = result[0]
-    assert result.sparse.shape == (172, 9)
-    assert result.dense is None
-
-    kwargs = dict(run=1, subject='01', mode='dense', force=False)
-    result = analysis['run'].get_design_matrix(**kwargs)[0]
-    assert result.sparse is None
-    assert result.dense is None
-
-    kwargs = dict(run=1, subject='01', mode='dense', force=True,
-                  sampling_rate='highest')
-    result = analysis['run'].get_design_matrix(**kwargs)[0]
-    assert result.sparse is None
-    assert result.dense.shape == (4800, 10)
-
-    kwargs = dict(run=1, subject='01', mode='dense', force=True,
-                  sampling_rate='TR')
-    result = analysis['run'].get_design_matrix(**kwargs)[0]
-    assert result.sparse is None
-    assert result.dense.shape == (240, 10)
-
-    kwargs = dict(run=1, subject='01', mode='dense', force=True,
-                  sampling_rate=0.5)
-    result = analysis['run'].get_design_matrix(**kwargs)[0]
-    assert result.sparse is None
-    assert result.dense.shape == (240, 10)
-
-    # format='long' should be ignored for dense output
-    kwargs = dict(run=1, subject='01', mode='dense', force=True,
-                  format='long', entities=False)
-    result = analysis['run'].get_design_matrix(**kwargs)[0]
-    assert result.sparse is None
-    assert result.dense.shape == (240, 1)
-
-    kwargs = dict(run=1, subject='01', mode='sparse', format='wide',
-                  entities=False)
-    result = analysis['run'].get_design_matrix(**kwargs)[0]
-    assert result.dense is None
-    assert result.sparse.shape == (86, 4)
-
-
 def test_first_level_sparse_design_matrix(analysis):
-    result = analysis['run'].get_design_matrix(subject=['01'])
-    assert len(result) == 3
-    df = result[0].sparse
+    collections = analysis['run'].get_collections(subject=['01'])
+    assert len(collections) == 3
+    df = collections[0].to_df(format='long')
     assert df.shape == (172, 9)
     assert df['condition'].nunique() == 2
-    assert set(result[0][0].columns) == {'amplitude', 'onset', 'duration',
+    assert set(df.columns) == {'amplitude', 'onset', 'duration',
                                          'condition', 'subject', 'run',
                                          'task', 'datatype', 'suffix'}
 
 
 def test_post_first_level_sparse_design_matrix(analysis):
 
-    result = analysis['session'].get_design_matrix(entities=False)
-    assert len(result) == 2
-    assert len(result[0]) == 3
-    assert result[0].sparse.shape == (9, 2)
-    assert result[0].entities == {
+    collections = analysis['session'].get_collections()
+    assert len(collections) == 2
+    result = collections[0].to_df(format='long')
+    assert result.shape == (9, 11)
+    result = collections[0].to_df(format='long', entities=False)
+    assert result.shape == (9, 2)
+    entities = {
         'subject': '01',
         'task': 'mixedgamblestask',
         'datatype': 'func',
         'suffix': 'bold'}
+    assert not set(entities.keys()) - set(collections[0].entities.keys())
+    assert not set(entities.values()) - set(collections[0].entities.values())
 
     # Participant level and also check integer-based indexing
-    result = analysis['participant'].get_design_matrix()
-    assert len(result) == 2
+    collections = analysis['participant'].get_collections()
+    assert len(collections) == 2
     assert analysis[2].name == 'participant'
 
     # Dataset level
-    result = analysis['group'].get_design_matrix()
-    assert len(result) == 1
-    data = result[0].sparse
+    collections = analysis['group'].get_collections()
+    assert len(collections) == 1
+    data = collections[0].to_df(format='long')
     assert len(data) == 10
     assert data['subject'].nunique() == 2
 
@@ -112,11 +65,18 @@ def test_post_first_level_sparse_design_matrix(analysis):
 
     # Calling an invalid level name should raise an exception
     with pytest.raises(KeyError):
-        result = analysis['nonexistent_name'].get_design_matrix()
+        result = analysis['nonexistent_name'].to_df()
+
+
+def test_step_get_collections(analysis):
+    collections = analysis['run'].get_collections(subject='01')
+    assert len(collections) == 3
+    assert isinstance(collections[0], BIDSVariableCollection)
 
 
 def test_contrast_info(analysis):
-    contrast_lists = analysis['run'].get_contrasts(subject='01')
+    colls = analysis['run'].get_collections(subject='01')
+    contrast_lists = [analysis['run'].get_contrasts(c) for c in colls]
     assert len(contrast_lists) == 3
     for cl in contrast_lists:
         assert len(cl) == 3
@@ -133,8 +93,9 @@ def test_contrast_info(analysis):
 
 def test_contrast_info_with_specified_variables(analysis):
     varlist = ['RT', 'dummy']
-    contrast_lists = analysis['run'].get_contrasts(subject='01',
-                                                   variables=varlist)
+    colls = analysis['run'].get_collections(subject='01')
+    contrast_lists = [analysis['run'].get_contrasts(c, variables=varlist)
+                      for c in colls]
     assert len(contrast_lists) == 3
     for cl in contrast_lists:
         assert len(cl) == 3
@@ -149,8 +110,9 @@ def test_contrast_info_with_specified_variables(analysis):
 
 
 def test_contrast_info_F_contrast(analysis):
-    contrast_lists = analysis['run'].get_contrasts(subject='01',
-                                                   names=["crummy-F"])
+    colls = analysis['run'].get_collections(subject='01')
+    contrast_lists = [analysis['run'].get_contrasts(c, names=["crummy-F"])
+                      for c in colls]
     assert len(contrast_lists) == 3
     for cl in contrast_lists:
         assert len(cl) == 1
@@ -164,22 +126,37 @@ def test_contrast_info_F_contrast(analysis):
 
 
 def test_dummy_contrasts(analysis):
-    names = [c.name for c in analysis['run'].get_contrasts(subject='01')[0]]
-    session = analysis['session'].get_contrasts(subject='01')[0]
+    collection = analysis['run'].get_collections(subject='01')[0]
+    names = [c.name for c in analysis['run'].get_contrasts(collection)]
+
+    collections = analysis['session'].get_collections()
+    collection = analysis['session'].get_collections(subject='01')[0]
+    session = analysis['session'].get_contrasts(collection)
     for cl in session:
         assert cl.type == 'FEMA'
         assert cl.name in names
 
-    participant = analysis['participant'].get_contrasts(subject='01')[0]
+    collection = analysis['participant'].get_collections(subject='01')[0]
+    participant = analysis['participant'].get_contrasts(collection)
     assert len(participant) == 3
     for cl in participant:
         assert cl.type == 'FEMA'
         assert cl.name in names
 
-    group = analysis['group'].get_contrasts()[0]
+    collection = analysis['group'].get_collections()[0]
+    group = analysis['group'].get_contrasts(collection)
     group_names = []
     for cl in group:
         assert cl.type == 't'
         group_names.append(cl.name)
 
     assert set(names) < set(group_names)
+
+
+def test_get_model_spec(analysis):
+    collection = analysis['run'].get_collections(subject='01', run=1)[0]
+    model_spec = analysis['run'].get_model_spec(collection, 'TR')
+    assert model_spec.__class__.__name__ == 'GLMMSpec'
+    assert model_spec.X.shape == (240, 1)
+    assert model_spec.Z is None
+    assert {'RT'} == set(model_spec.terms.keys())
diff --git a/bids/analysis/tests/test_model_spec.py b/bids/analysis/tests/test_model_spec.py
new file mode 100644
index 00000000..e9301549
--- /dev/null
+++ b/bids/analysis/tests/test_model_spec.py
@@ -0,0 +1,67 @@
+from os.path import join
+import json
+
+import pandas as pd
+import numpy as np
+import pytest
+
+from bids.variables import BIDSVariableCollection
+from bids.analysis.model_spec import (GLMMSpec, Term, VarComp)
+from bids.analysis import TransformerManager
+from bids.layout import BIDSLayout
+from bids.tests import get_test_data_path
+
+
+@pytest.fixture(scope='module')
+def collection():
+    """Create a collection at dataset level."""
+    df = pd.DataFrame({
+        'subject': np.repeat(np.arange(10), 4),
+        'run': np.tile(np.arange(4), 10),
+        'age': np.repeat(np.random.normal(40, 10, 10).astype(int), 4)
+    })
+    entities = df[['subject', 'run']]
+    return BIDSVariableCollection.from_df(df, entities=entities)
+
+
+def test_bids_variable_collection_to_model_design(collection):
+    layout_path = join(get_test_data_path(), 'ds005')
+    json_file = join(layout_path, 'models', 'ds-005_type-mfx_model.json')
+    spec = json.load(open(json_file, 'r'))
+    tm = TransformerManager()
+    collection = tm.transform(collection, spec['Steps'][1]['Transformations'])
+    md = GLMMSpec.from_collection(collection, spec['Steps'][1]['Model'])
+
+    assert len(md.terms) == 2
+    assert md.terms['age'].values.shape == (40, )
+    assert md.terms['SubjectSlopes'].values.shape == (40, 10)
+    index_vec = md.terms['SubjectSlopes'].index_vec
+    assert index_vec.shape == (40,)
+    assert np.array_equal(np.sort(np.unique(index_vec)), np.arange(10) + 1)
+    assert md.Z.shape == (40, 10)
+    assert md.Z.columns[0] == 'SubjectSlopes.0'
+
+
+def test_glmmspec_empty_init():
+    md = GLMMSpec()
+    assert md.terms == {}
+    assert md.X is None
+    assert md.Z is None
+    assert md.family is None
+    assert md.family is None
+    assert md.sigma is None
+
+
+def test_fixed_term_init():
+    t = Term('dummy', np.random.normal(size=20), categorical=True)
+    assert t.name == 'dummy'
+    assert t.values.shape == (20,)
+    assert t.categorical
+
+
+def test_var_comp_init():
+    Z = np.repeat(np.eye(20), 5, axis=0)
+    t = VarComp('sigma', Z)
+    assert t.name == 'sigma'
+    assert t.values.shape == (100, 20)
+    assert t.categorical
diff --git a/bids/analysis/tests/test_transformations.py b/bids/analysis/tests/test_transformations.py
index 97899e3a..cc3126dd 100644
--- a/bids/analysis/tests/test_transformations.py
+++ b/bids/analysis/tests/test_transformations.py
@@ -2,7 +2,7 @@
 from bids.analysis import transformations as transform
 from bids.variables import SparseRunVariable
 from bids.variables.entities import RunInfo
-from bids.variables.kollekshuns import BIDSRunVariableCollection
+from bids.variables.collections import BIDSRunVariableCollection
 from bids.layout import BIDSLayout
 import math
 import pytest
diff --git a/bids/tests/data/ds005/models/ds-005_type-mfx_model.json b/bids/tests/data/ds005/models/ds-005_type-mfx_model.json
new file mode 100644
index 00000000..c9358505
--- /dev/null
+++ b/bids/tests/data/ds005/models/ds-005_type-mfx_model.json
@@ -0,0 +1,42 @@
+{
+  "Name": "mfx_model",
+  "Description": "simple test MFX model",
+  "Steps": [
+    {
+      "Name": "run",
+      "Level": "Run",
+      "Model": {
+        "X": ["RT", "trial_type"]
+      },
+      "Transformations": [
+        {
+          "Name": "Scale",
+          "Input": "RT"
+        }
+      ],
+      "DummyContrasts": {
+        "Conditions": ["trial_type", "RT"],
+        "Type": "t"
+      }
+    },
+    {
+      "Name": "group",
+      "Level": "Dataset",
+      "Transformations": [
+        {
+          "Name": "Factor",
+          "Input": "subject"
+        }
+      ],
+      "Model": {
+        "X": ["age"],
+        "VarianceComponents": [
+          {
+            "Name": "SubjectSlopes",
+            "Levels": "subject*"
+          }
+        ]
+      }
+    }
+  ]
+}
diff --git a/bids/variables/tests/test_collections.py b/bids/variables/tests/test_collections.py
index d5ee9509..69ff67ef 100644
--- a/bids/variables/tests/test_collections.py
+++ b/bids/variables/tests/test_collections.py
@@ -1,17 +1,20 @@
-from bids.layout import BIDSLayout
-import pytest
 from os.path import join, dirname, abspath
+
+import pytest
+
+from bids.layout import BIDSLayout
 from bids.tests import get_test_data_path
 from bids.variables import (DenseRunVariable, SparseRunVariable,
                             merge_collections)
-
+from bids.variables.entities import RunInfo
 
 @pytest.fixture(scope="module")
 def run_coll():
     path = join(get_test_data_path(), 'ds005')
     layout = BIDSLayout(path)
+    # Limit to a few subjects to reduce test running time
     return layout.get_collections('run', types=['events'], merge=True,
-                                  scan_length=480)
+                                  scan_length=480, subject=['01', '02', '04'])
 
 
 @pytest.fixture(scope="module")
@@ -19,7 +22,7 @@ def run_coll_list():
     path = join(get_test_data_path(), 'ds005')
     layout = BIDSLayout(path)
     return layout.get_collections('run', types=['events'], merge=False,
-                                  scan_length=480)
+                                  scan_length=480, subject=['01', '02', '04'])
 
 
 def test_run_variable_collection_init(run_coll):
@@ -27,21 +30,60 @@ def test_run_variable_collection_init(run_coll):
     assert run_coll.sampling_rate == 10
 
 
+def test_run_variable_collection_sparse_variable_accessors(run_coll):
+    coll = run_coll.clone()
+    assert coll.get_sparse_variables()
+    assert coll.all_sparse()
+    coll.variables['RT'] = coll.variables['RT'].to_dense(1)
+    assert not coll.all_sparse()
+    assert len(coll.get_sparse_variables()) + 1 == len(coll.variables)
+
+
+def test_run_variable_collection_dense_variable_accessors(run_coll):
+    coll = run_coll.clone()
+    coll.variables['RT'] = coll.variables['RT'].to_dense(1)
+    assert not coll.all_dense()
+    assert len(coll.get_dense_variables()) == 1
+    for k, v in coll.variables.items():
+        if k == 'RT':
+            continue
+        coll.variables[k] = v.to_dense(1)
+    assert coll.all_dense()
+
+
+def test_run_variable_collection_get_sampling_rate(run_coll):
+    coll = run_coll.clone()
+    assert coll._get_sampling_rate(None) == 10
+    assert coll._get_sampling_rate('TR') == 0.5
+    coll.variables['RT'].run_info[0] = RunInfo({}, 200, 10, None)
+    with pytest.raises(ValueError) as exc:
+        coll._get_sampling_rate('TR')
+        assert exc.value.message.startswith('Non-unique')
+    assert coll._get_sampling_rate('highest') == 10
+    coll.variables['RT1'] = coll.variables['RT'].to_dense(5.)
+    coll.variables['RT2'] = coll.variables['RT'].to_dense(12.)
+    assert coll._get_sampling_rate('highest') == 12.
+    assert coll._get_sampling_rate(20) == 20
+    with pytest.raises(ValueError) as exc:
+        coll._get_sampling_rate('BLARGH')
+        assert exc.value.message.startswith('Invalid')
+
+
 def test_resample_run_variable_collection(run_coll):
     run_coll = run_coll.clone()
     resampled = run_coll.resample()
-    assert not resampled  # Empty because all variables are sparse
+    assert not resampled.variables  # Empty because all variables are sparse
 
-    resampled = run_coll.resample(force_dense=True)
+    resampled = run_coll.resample(force_dense=True).variables
     assert len(resampled) == 7
     assert all([isinstance(v, DenseRunVariable) for v in resampled.values()])
     assert len(set([v.sampling_rate for v in resampled.values()])) == 1
-    targ_len = 480 * 16 * 3 * 10
+    targ_len = 480 * 3 * 3 * 10
     assert all([len(v.values) == targ_len for v in resampled.values()])
 
     sr = 20
-    resampled = run_coll.resample(sr, force_dense=True)
-    targ_len = 480 * 16 * 3 * sr
+    resampled = run_coll.resample(sr, force_dense=True).variables
+    targ_len = 480 * 3 * 3 * sr
     assert all([len(v.values) == targ_len for v in resampled.values()])
 
     run_coll.resample(sr, force_dense=True, in_place=True)
@@ -53,34 +95,139 @@ def test_resample_run_variable_collection(run_coll):
     assert all([isinstance(v, DenseRunVariable) for v in vars_])
 
 
-def test_run_variable_collection_to_df(run_coll):
+def test_run_variable_collection_to_df_all_sparse_vars(run_coll):
     run_coll = run_coll.clone()
 
-    # All variables sparse, wide format
+    timing_cols = {'onset', 'duration'}
+    entity_cols = {'subject', 'run', 'task',  'suffix', 'datatype'}
+    cond_names = {'PTval', 'RT', 'gain', 'loss', 'parametric gain', 'respcat',
+                  'respnum', 'trial_type'}
+    condition = {'condition'}
+    ampl = {'amplitude'}
+
+    # Fails because all vars are sparse
+    with pytest.raises(ValueError):
+        df = run_coll.to_df(include_sparse=False)
+
+    # Fails because no such variables exist
+    with pytest.raises(ValueError):
+        df = run_coll.to_df(variables=['rubadubdub', '999'])
+
+    # Wide format
     df = run_coll.to_df()
-    assert df.shape == (4096, 15)
-    wide_cols = {'onset', 'duration', 'subject', 'run', 'task',
-                 'PTval', 'RT', 'gain', 'loss', 'parametric gain', 'respcat',
-                 'respnum', 'trial_type', 'suffix', 'datatype'}
-    assert set(df.columns) == wide_cols
+    events_per_sub = 256
+    assert df.shape == (events_per_sub * 3, 15) 
+    assert set(df.columns) == timing_cols.union(entity_cols, cond_names)
 
-    # All variables sparse, wide format
+    # Wide format, selecting variables by name
+    df = run_coll.to_df(format='wide', variables=['RT', 'PTval'])
+    assert df.shape == (events_per_sub * 3, 9)
+    assert set(df.columns) == timing_cols.union(entity_cols, {'RT', 'PTval'})
+
+    # Long format
     df = run_coll.to_df(format='long')
-    assert df.shape == (32768, 9)
-    long_cols = {'amplitude', 'duration', 'onset', 'condition', 'run',
-                 'task', 'subject', 'suffix', 'datatype'}
-    assert set(df.columns) == long_cols
+    assert df.shape == (events_per_sub * 3 * 8, 9)
+    assert set(df.columns) == timing_cols.union(entity_cols, condition, ampl)
+
+    # Long format, selecting variables by name
+    df = run_coll.to_df(format='long', variables=['RT', 'PTval'])
+    assert df.shape == (events_per_sub * 3 * 2, 9)
+    assert set(df.columns) == timing_cols.union(entity_cols, condition, ampl)
+
+    # Wide format without entity columnns
+    df = run_coll.to_df(format='wide', entities=False)
+    assert df.shape == (events_per_sub * 3, 10)
+    assert set(df.columns) == timing_cols | cond_names
+
+    # Wide format without timing columns
+    df = run_coll.to_df(format='wide', timing=False)
+    assert df.shape == (events_per_sub * 3, 13)
+    assert set(df.columns) == entity_cols | cond_names
+
+
+def test_run_variable_collection_to_df_all_dense_vars(run_coll):
+
+    timing_cols = {'onset', 'duration'}
+    entity_cols = {'subject', 'run', 'task',  'suffix', 'datatype'}
+    cond_names = {'PTval', 'RT', 'gain', 'loss', 'parametric gain', 'respcat',
+                  'respnum', 'trial_type'}
+    md_names = {'TaskName', 'RepetitionTime', 'extension', 'SliceTiming'}
+    condition = {'condition'}
+    ampl = {'amplitude'}
+
+    # First few tests assume uniform sampling rate for all columns
+    unif_coll = run_coll.to_dense(sampling_rate=10)
+
+    # Fails because all vars are dense
+    with pytest.raises(ValueError):
+        unif_coll.to_df(include_dense=False)
+
+    # Fails because no such variables exist
+    with pytest.raises(ValueError):
+        df = unif_coll.to_df(variables=['rubadubdub', '999'])
+
+    # Wide format
+    df = unif_coll.to_df()
+    rows_per_var = 3 * 3 * 480 * 10  # subjects x runs x time x sampling rate
+    assert df.shape == (rows_per_var, 18)
+    cols = (timing_cols | entity_cols | cond_names | md_names) - {'trial_type'}
+    assert set(df.columns) == cols
+
+    # Wide format with variable selection
+    df = unif_coll.to_df(variables=['RT', 'gain'])
+    assert df.shape == (rows_per_var, 13)
+    cols = (timing_cols | entity_cols | {'RT', 'gain'} | md_names) - {'trial_type'}
+    assert set(df.columns) == cols
+
+    # Long format
+    df = unif_coll.to_df(format='long')
+    assert df.shape == (rows_per_var * 7, 13)
+    cols = timing_cols | entity_cols | condition | ampl | md_names
+    assert set(df.columns) == cols
+
+    # Long format with variable selection
+    df = unif_coll.to_df(format='long', variables=['RT', 'gain'])
+    assert df.shape == (rows_per_var * 2, 13)
+    cols = timing_cols | entity_cols | condition | ampl | md_names
+    assert set(df.columns) == cols
+
+    # Test resampling to TR
+    df = unif_coll.to_df(sampling_rate='TR')
+    n_rows = int(480 * 3 * 3 / 2)
+    assert df.shape == (n_rows, 18)
+    cols = (timing_cols | entity_cols | cond_names | md_names) - {'trial_type'}
+    assert set(df.columns) == cols
+
+    # Test resampling to highest when variable sampling rates diverge
+    unif_coll.resample(variables='RT', sampling_rate=12, in_place=True)
+    df = unif_coll.to_df(sampling_rate='highest')
+    n_rows = int(rows_per_var * 12 / 10)
+    assert df.shape == (n_rows, 18)
+
+
+def test_run_variable_collection_to_df_mixed_vars(run_coll):
+    coll = run_coll.clone()
+    coll.to_dense(10, variables=['RT', 'loss','gain'], in_place=True)
+
+    # Only sparse variables
+    df = coll.to_df(include_dense=False)
+    rows = 3 * 256
+    assert df.shape == (rows, 12)
+    assert 'RT' not in df.columns
+    assert 'respcat' in df.columns
 
-    # All variables dense, wide format
-    df = run_coll.to_df(sparse=False)
-    assert df.shape == (230400, 18)
-    extra_cols = {'TaskName', 'RepetitionTime', 'extension', 'SliceTiming'}
-    assert set(df.columns) == (wide_cols | extra_cols) - {'trial_type'}
+    # Only dense variables
+    df = coll.to_df(include_sparse=False)
+    rows = 3 * 3 * 480 * 10
+    assert df.shape == (rows, 14)
+    assert 'RT' in df.columns
+    assert 'respcat' not in df.columns
 
-    # All variables dense, wide format
-    df = run_coll.to_df(sparse=False, format='long')
-    assert df.shape == (1612800, 13)
-    assert set(df.columns) == (long_cols | extra_cols)
+    # Everything
+    df = coll.to_df(sampling_rate=5)
+    rows = 3 * 3 * 480 * 5
+    assert df.shape == (rows, 18)
+    assert not {'RT', 'respcat'} - set(df.columns)
 
 
 def test_merge_collections(run_coll, run_coll_list):
@@ -94,16 +241,16 @@ def test_merge_collections(run_coll, run_coll_list):
 def test_get_collection_entities(run_coll_list):
     coll = run_coll_list[0]
     ents = coll.entities
-    assert {'run', 'task', 'subject', 'suffix', 'datatype'} == set(ents.keys())
+    assert not ({'run', 'task', 'subject', 'suffix', 'datatype'} - set(ents.keys()))
 
     merged = merge_collections(run_coll_list[:3])
     ents = merged.entities
-    assert {'task', 'subject', 'suffix', 'datatype'} == set(ents.keys())
+    assert not ({'task', 'subject', 'suffix', 'datatype'} - set(ents.keys()))
     assert ents['subject'] == '01'
 
     merged = merge_collections(run_coll_list[3:6])
     ents = merged.entities
-    assert {'task', 'subject', 'suffix', 'datatype'} == set(ents.keys())
+    assert not ({'task', 'subject', 'suffix', 'datatype'} - set(ents.keys()))
     assert ents['subject'] == '02'
 
 
diff --git a/bids/variables/tests/test_variables.py b/bids/variables/tests/test_variables.py
index 5f42a766..54083989 100644
--- a/bids/variables/tests/test_variables.py
+++ b/bids/variables/tests/test_variables.py
@@ -66,14 +66,14 @@ def test_sparse_run_variable_to_dense(layout1):
     index = load_variables(layout1, types='events', scan_length=480)
     runs = index.get_nodes('run', {'subject': ['01', '02']})
 
-    for i, run in enumerate(runs):
+    for _, run in enumerate(runs):
         var = run.variables['RT']
         dense = var.to_dense(20)
 
-        # Check that all unique values are identical
-        sparse_vals = set(np.unique(var.values.values)) | {0}
+        # All sparse values must occur in the dense variable
         dense_vals = set(np.unique(dense.values.values))
-        assert sparse_vals == dense_vals
+        sparse_vals = set(np.unique(var.values.values)) | {0}
+        assert not (sparse_vals - dense_vals)
 
         assert len(dense.values) > len(var.values)
         assert isinstance(dense, DenseRunVariable)
@@ -135,9 +135,12 @@ def test_densify_merged_variables(layout1):
     for i in range(len(runs)):
         onset = i * n_rows
         offset = onset + n_rows
-        run_vals = vars_[i].to_dense(SR).values
-        dense_vals = dense.values.iloc[onset:offset].reset_index(drop=True)
-        assert dense_vals.equals(run_vals)
+        # resampled values from split vs. merged will be trivially different
+        # at boundaries due to interpolation, so leave a buffer
+        run_vals = vars_[i].to_dense(SR).values.values.ravel()
+        dense_vals = dense.values.values.ravel()[onset:offset]
+        assert run_vals.shape == dense_vals.shape
+        assert np.array_equal(run_vals, dense_vals)
 
 
 def test_merge_simple_variables(layout2):