diff --git a/nltools/data/adjacency.py b/nltools/data/adjacency.py
index e37d22a0..984f2119 100644
--- a/nltools/data/adjacency.py
+++ b/nltools/data/adjacency.py
@@ -61,8 +61,7 @@ class Adjacency(object):
 
     '''
 
-    def __init__(self, data=None, Y=None, matrix_type=None, labels=None,
-                 **kwargs):
+    def __init__(self, data=None, Y=None, matrix_type=None, labels=[], **kwargs):
         if matrix_type is not None:
             if matrix_type.lower() not in ['distance', 'similarity', 'directed',
                                            'distance_flat', 'similarity_flat',
@@ -126,7 +125,7 @@ def __init__(self, data=None, Y=None, matrix_type=None, labels=None,
         else:
             self.Y = pd.DataFrame()
 
-        if labels is not None:
+        if labels:
             if not isinstance(labels, (list, np.ndarray)):
                 raise ValueError("Make sure labels is a list or numpy array.")
             if self.is_single_matrix:
@@ -147,7 +146,7 @@ def __init__(self, data=None, Y=None, matrix_type=None, labels=None,
                         raise ValueError("All lists of labels must be same length as shape of data.")
                     self.labels = deepcopy(labels)
         else:
-            self.labels = None
+            self.labels = []
 
     def __repr__(self):
         return ("%s.%s(shape=%s, square_shape=%s, Y=%s, is_symmetric=%s,"
@@ -162,7 +161,7 @@ def __repr__(self):
 
     def __getitem__(self, index):
         new = self.copy()
-        if isinstance(index, int):
+        if isinstance(index, (int, np.integer)):
             new.data = np.array(self.data[index, :]).squeeze()
             new.is_single_matrix = True
         else:
@@ -184,35 +183,41 @@ def __iter__(self):
 
     def __add__(self, y):
         new = deepcopy(self)
-        if isinstance(y, (int, float)):
+        if isinstance(y, (int, np.integer, float, np.floating)):
             new.data = new.data + y
-        if isinstance(y, Adjacency):
+        elif isinstance(y, Adjacency):
             if self.shape() != y.shape():
                 raise ValueError('Both Adjacency() instances need to be the '
                                  'same shape.')
             new.data = new.data + y.data
+        else:
+            raise ValueError('Can only add int, float, or Adjacency')
         return new
 
     def __sub__(self, y):
         new = deepcopy(self)
-        if isinstance(y, (int, float)):
+        if isinstance(y, (int, np.integer, float, np.floating)):
             new.data = new.data - y
-        if isinstance(y, Adjacency):
+        elif isinstance(y, Adjacency):
             if self.shape() != y.shape():
                 raise ValueError('Both Adjacency() instances need to be the '
                                  'same shape.')
             new.data = new.data - y.data
+        else:
+            raise ValueError('Can only subtract int, float, or Adjacency')
         return new
 
     def __mul__(self, y):
         new = deepcopy(self)
-        if isinstance(y, (int, float)):
+        if isinstance(y, (int, np.integer, float, np.floating)):
             new.data = new.data * y
-        if isinstance(y, Adjacency):
+        elif isinstance(y, Adjacency):
             if self.shape() != y.shape():
                 raise ValueError('Both Adjacency() instances need to be the '
                                  'same shape.')
             new.data = np.multiply(new.data, y.data)
+        else:
+            raise ValueError('Can only multiply int, float, or Adjacency')
         return new
 
     @staticmethod
@@ -330,8 +335,8 @@ def plot(self, limit=3, *args, **kwargs):
         ''' Create Heatmap of Adjacency Matrix'''
 
         if self.is_single_matrix:
-            f, a = plt.subplots(nrows=1, figsize=(7, 5))
-            if self.labels is None:
+            _, a = plt.subplots(nrows=1, figsize=(7, 5))
+            if not self.labels:
                 sns.heatmap(self.squareform(), square=True, ax=a,
                             *args, **kwargs)
             else:
@@ -341,8 +346,8 @@ def plot(self, limit=3, *args, **kwargs):
                             *args, **kwargs)
         else:
             n_subs = np.minimum(len(self), limit)
-            f, a = plt.subplots(nrows=n_subs, figsize=(7, len(self)*5))
-            if self.labels is None:
+            _, a = plt.subplots(nrows=n_subs, figsize=(7, len(self)*5))
+            if not self.labels:
                 for i in range(n_subs):
                     sns.heatmap(self[i].squareform(), square=True, ax=a[i],
                                 *args, **kwargs)
@@ -352,7 +357,7 @@ def plot(self, limit=3, *args, **kwargs):
                                 xticklabels=self.labels[i],
                                 yticklabels=self.labels[i],
                                 ax=a[i], *args, **kwargs)
-        return f
+        return
 
     def mean(self, axis=0):
         ''' Calculate mean of Adjacency
@@ -548,18 +553,18 @@ def _convert_data_similarity(data, perm_type=None, ignore_diagonal=ignore_diagon
                                     metric=metric, n_permute=n_permute,
                                     **kwargs) for x in self]
 
-    def distance(self, method='correlation', **kwargs):
+    def distance(self, metric='correlation', **kwargs):
         ''' Calculate distance between images within an Adjacency() instance.
 
         Args:
-            method: (str) type of distance metric (can use any scikit learn or
+            metric: (str) type of distance metric (can use any scikit learn or
                     sciypy metric)
 
         Returns:
             dist: (Adjacency) Outputs a 2D distance matrix.
 
         '''
-        return Adjacency(pairwise_distances(self.data, metric=method, **kwargs),
+        return Adjacency(pairwise_distances(self.data, metric=metric, **kwargs),
                          matrix_type='distance')
 
     def threshold(self, upper=None, lower=None, binarize=False):
@@ -611,7 +616,7 @@ def to_graph(self):
                 G = nx.DiGraph(self.squareform())
             else:
                 G = nx.Graph(self.squareform())
-            if self.labels is not None:
+            if self.labels:
                 labels = {x: y for x, y in zip(G.nodes, self.labels)}
                 nx.relabel_nodes(G, labels, copy=False)
             return G
@@ -687,7 +692,7 @@ def plot_label_distance(self, labels=None, ax=None):
                            palette={"Within": "lightskyblue", "Between": "red"}, ax=ax)
         f.set_ylabel('Average Distance')
         f.set_title('Average Group Distance')
-        return f
+        return
 
     def stats_label_distance(self, labels=None, n_permute=5000, n_jobs=-1):
         ''' Calculate permutation tests on within and between label distance.
@@ -745,10 +750,9 @@ def plot_silhouette(self, labels=None, ax=None, permutation_test=True,
             if len(labels) != distance.shape[0]:
                 raise ValueError('Labels must be same length as distance matrix')
 
-        (f, outAll) = plot_silhouette(distance, labels, ax=None,
+        return plot_silhouette(distance, pd.Series(labels), ax=None,
                                       permutation_test=True,
                                       n_permute=5000, **kwargs)
-        return (f, outAll)
 
     def bootstrap(self, function, n_samples=5000, save_weights=False,
                   n_jobs=-1, random_state=None, *args, **kwargs):
@@ -779,7 +783,7 @@ def bootstrap(self, function, n_samples=5000, save_weights=False,
         bootstrapped = Adjacency(bootstrapped)
         return summarize_bootstrap(bootstrapped, save_weights=save_weights)
 
-    def plot_mds(self, n_components=2, metric=True, labels_color=None,
+    def plot_mds(self, n_components=2, metric=True, labels=None, labels_color=None,
                  cmap=plt.cm.hot_r, n_jobs=-1, view=(30, 20),
                  figsize=[12, 8], ax=None, *args, **kwargs):
         ''' Plot Multidimensional Scaling
@@ -787,12 +791,11 @@ def plot_mds(self, n_components=2, metric=True, labels_color=None,
             Args:
                 n_components: (int) Number of dimensions to project (can be 2 or 3)
                 metric: (bool) Perform metric or non-metric dimensional scaling; default
+                labels: (list) Can override labels stored in Adjacency Class
                 labels_color: (str) list of colors for labels, if len(1) then make all same color
                 n_jobs: (int) Number of parallel jobs
                 view: (tuple) view for 3-Dimensional plot; default (30,20)
 
-            Returns:
-                fig: returns matplotlib figure
         '''
 
         if self.matrix_type != 'distance':
@@ -801,10 +804,15 @@ def plot_mds(self, n_components=2, metric=True, labels_color=None,
             raise ValueError("MDS only works on single matrices.")
         if n_components not in [2, 3]:
             raise ValueError('Cannot plot {0}-d image'.format(n_components))
+        if labels is not None:
+            if len(labels) != self.square_shape()[0]:
+                raise ValueError("Make sure labels matches the same shape as Adjaency data")
+        else:
+            labels = self.labels
         if labels_color is not None:
-            if self.labels is None:
+            if len(labels) == 0:
                 raise ValueError("Make sure that Adjacency object has labels specified.")
-            if len(self.labels) != len(labels_color):
+            if len(labels) != len(labels_color):
                 raise ValueError("Length of labels_color must match self.labels.")
 
         # Run MDS
@@ -814,7 +822,6 @@ def plot_mds(self, n_components=2, metric=True, labels_color=None,
 
         # Create Plot
         if ax is None:  # Create axis
-            returnFig = True
             fig = plt.figure(figsize=figsize)
             if n_components == 3:
                 ax = fig.add_subplot(111, projection='3d')
@@ -830,21 +837,18 @@ def plot_mds(self, n_components=2, metric=True, labels_color=None,
 
         # Plot labels
         if labels_color is None:
-            labels_color = ['black'] * len(self.labels)
+            labels_color = ['black'] * len(labels)
         if n_components == 3:
-            for ((x, y, z), label, color) in zip(proj, self.labels, labels_color):
+            for ((x, y, z), label, color) in zip(proj, labels, labels_color):
                 ax.text(x, y, z, label, color='white', bbox=dict(facecolor=color, alpha=1, boxstyle="round,pad=0.3"))
         else:
-            for ((x, y), label, color) in zip(proj, self.labels, labels_color):
+            for ((x, y), label, color) in zip(proj, labels, labels_color):
                 ax.text(x, y, label, color='white',  # color,
                         bbox=dict(facecolor=color, alpha=1, boxstyle="round,pad=0.3"))
 
         ax.xaxis.set_visible(False)
         ax.yaxis.set_visible(False)
 
-        if returnFig:
-            return fig
-
     def distance_to_similarity(self, beta=1):
         '''Convert distance matrix to similarity matrix
 
@@ -918,14 +922,15 @@ def regress(self, X, mode='ols', **kwargs):
             stats['beta'].data, stats['t'].data, stats['p'].data = b.squeeze(), t.squeeze(), p.squeeze()
             stats['residual'] = self.copy()
             stats['residual'].data = res
+            stats['df'] = df
         else:
             raise ValueError('X must be a Design_Matrix or Adjacency Instance.')
 
         return stats
 
     def social_relations_model(self, summarize_results=True, nan_replace=True):
-        '''Estimate the social relations model from a matrix for a round-robin design
-
+        '''Estimate the social relations model from a matrix for a round-robin design.
+        
         X_{ij} = m + \alpha_i + \beta_j + g_{ij} + \episolon_{ijl}
 
         where X_{ij} is the score for person i rating person j, m is the group mean,
@@ -1133,7 +1138,7 @@ def fix_missing(data):
             if data.is_single_matrix:
                 X, coord = fix_missing(data)
             else:
-                X = []; coord = [];
+                X = []; coord = []
                 for d in data:
                     m, c = fix_missing(d)
                     X.append(m)
diff --git a/nltools/data/brain_data.py b/nltools/data/brain_data.py
index f015fd9f..7661db3e 100644
--- a/nltools/data/brain_data.py
+++ b/nltools/data/brain_data.py
@@ -38,7 +38,7 @@
 from nltools.analysis import Roc
 from nilearn.input_data import NiftiMasker
 from nilearn.plotting import plot_stat_map
-from nilearn.image import resample_img, smooth_img
+from nilearn.image import smooth_img, resample_to_img
 from nilearn.masking import intersect_masks
 from nilearn.regions import connected_regions, connected_label_regions
 from nltools.utils import (set_algorithm,
@@ -47,6 +47,7 @@
                            _bootstrap_apply_func,
                            set_decomposition_algorithm,
                            check_brain_data,
+                           check_brain_data_is_single,
                            _roi_func,
                            get_mni_from_img_resolution,
                            _df_meta_to_arr)
@@ -207,7 +208,7 @@ def __repr__(self):
 
     def __getitem__(self, index):
         new = deepcopy(self)
-        if isinstance(index, int):
+        if isinstance(index, (int, np.integer)):
             new.data = np.array(self.data[index, :]).squeeze()
         else:
             if isinstance(index, slice):
@@ -243,51 +244,59 @@ def __len__(self):
 
     def __add__(self, y):
         new = deepcopy(self)
-        if isinstance(y, (int, float)):
+        if isinstance(y, (int, np.integer, float, np.floating)):
             new.data = new.data + y
-        if isinstance(y, Brain_Data):
+        elif isinstance(y, Brain_Data):
             if self.shape() != y.shape():
                 raise ValueError("Both Brain_Data() instances need to be the "
                                  "same shape.")
             new.data = new.data + y.data
+        else:
+            raise ValueError('Can only add int, float, or Brain_Data')
         return new
 
     def __radd__(self, y):
         new = deepcopy(self)
-        if isinstance(y, (int, float)):
+        if isinstance(y, (int, np.integer, float, np.floating)):
             new.data = y + new.data
         elif isinstance(y, Brain_Data):
             if self.shape() != y.shape():
                 raise ValueError("Both Brain_Data() instances need to be the "
                                  "same shape.")
             new.data = y.data + new.data
+        else:
+            raise ValueError('Can only add int, float, or Brain_Data')
         return new
 
     def __sub__(self, y):
         new = deepcopy(self)
-        if isinstance(y, (int, float)):
+        if isinstance(y, (int, np.integer, float, np.floating)):
             new.data = new.data - y
         elif isinstance(y, Brain_Data):
             if self.shape() != y.shape():
                 raise ValueError('Both Brain_Data() instances need to be the '
                                  'same shape.')
             new.data = new.data - y.data
+        else:
+            raise ValueError('Can only add int, float, or Brain_Data')
         return new
 
     def __rsub__(self, y):
         new = deepcopy(self)
-        if isinstance(y, (int, float)):
+        if isinstance(y, (int, np.integer, float, np.floating)):
             new.data = y - new.data
         elif isinstance(y, Brain_Data):
             if self.shape() != y.shape():
                 raise ValueError('Both Brain_Data() instances need to be the '
                                  'same shape.')
             new.data = y.data - new.data
+        else:
+            raise ValueError('Can only add int, float, or Brain_Data')
         return new
 
     def __mul__(self, y):
         new = deepcopy(self)
-        if isinstance(y, (int, float)):
+        if isinstance(y, (int, np.integer, float, np.floating)):
             new.data = new.data * y
         elif isinstance(y, Brain_Data):
             if self.shape() != y.shape():
@@ -301,17 +310,21 @@ def __mul__(self, y):
                                  'images in Brain_Data instance.')
             else:
                 new.data = np.dot(new.data.T, y).T
+        else:
+            raise ValueError('Can only multiply int, float, list, or Brain_Data')
         return new
 
     def __rmul__(self, y):
         new = deepcopy(self)
-        if isinstance(y, (int, float)):
+        if isinstance(y, (int, np.integer, float, np.floating)):
             new.data = y * new.data
         elif isinstance(y, Brain_Data):
             if self.shape() != y.shape():
                 raise ValueError("Both Brain_Data() instances need to be the "
                                  "same shape.")
             new.data = np.multiply(y.data, new.data)
+        else:
+            raise ValueError('Can only multiply int, float, or Brain_Data')
         return new
 
     def __iter__(self):
@@ -323,16 +336,54 @@ def shape(self):
 
         return self.data.shape
 
-    def mean(self):
-        """ Get mean of each voxel across images. """
+    def mean(self, axis=0):
+        ''' Get mean of each voxel or image
+            
+            Args:
+                axis: (int) across images=0 (default), within images=1
+            
+            Returns:
+                out: (float/np.array/Brain_Data)
+
+        '''
 
         out = deepcopy(self)
-        if len(self.shape()) > 1:
-            out.data = np.mean(self.data, axis=0)
-            out.X = pd.DataFrame()
-            out.Y = pd.DataFrame()
-        else:
+        if check_brain_data_is_single(self):
             out = np.mean(self.data)
+        else:
+            if axis == 0:
+                out.data = np.mean(self.data, axis=0)
+                out.X = pd.DataFrame()
+                out.Y = pd.DataFrame()
+            elif axis == 1:
+                out = np.mean(self.data, axis=1)
+            else:
+                raise ValueError('axis must be 0 or 1')
+        return out
+
+    def median(self, axis=0):
+        ''' Get median of each voxel or image
+            
+            Args:
+                axis: (int) across images=0 (default), within images=1
+            
+            Returns:
+                out: (float/np.array/Brain_Data)
+                
+        '''
+
+        out = deepcopy(self)
+        if check_brain_data_is_single(self):
+            out = np.median(self.data)
+        else:
+            if axis == 0:
+                out.data = np.median(self.data, axis=0)
+                out.X = pd.DataFrame()
+                out.Y = pd.DataFrame()
+            elif axis == 1:
+                out = np.median(self.data, axis=1)
+            else:
+                raise ValueError('axis must be 0 or 1')
         return out
 
     def std(self):
@@ -445,14 +496,14 @@ def plot(self, limit=5, anatomical=None, view='axial', threshold_upper=None, thr
                 anatomical = get_mni_from_img_resolution(self, img_type='plot')
 
             if self.data.ndim == 1:
-                f, a = plt.subplots(nrows=1, figsize=(15, 2))
+                _, a = plt.subplots(nrows=1, figsize=(15, 2))
                 plot_stat_map(self.to_nifti(), anatomical,
                               cut_coords=range(-40, 50, 10), display_mode='z',
                               black_bg=True, colorbar=True, draw_cross=False,
                               axes=a, **kwargs)
             else:
                 n_subs = np.minimum(self.data.shape[0], limit)
-                f, a = plt.subplots(nrows=n_subs, figsize=(15, len(self) * 2))
+                _, a = plt.subplots(nrows=n_subs, figsize=(15, len(self) * 2))
                 for i in range(n_subs):
                     plot_stat_map(self[i].to_nifti(), anatomical,
                                   cut_coords=range(-40, 50, 10),
@@ -462,7 +513,7 @@ def plot(self, limit=5, anatomical=None, view='axial', threshold_upper=None, thr
                                   draw_cross=False,
                                   axes=a[i],
                                   **kwargs)
-            return f
+            return
         elif view in ['glass', 'mni', 'full']:
             if self.data.ndim == 1:
                 return plot_brain(self, how=view, thr_upper=threshold_upper, thr_lower=threshold_lower, **kwargs)
@@ -846,11 +897,11 @@ def flatten_array(data):
             raise ValueError('Method must be one of: correlation, dot_product, cosine')
         return flatten_array(pexp)
 
-    def distance(self, method='euclidean', **kwargs):
+    def distance(self, metric='euclidean', **kwargs):
         """ Calculate distance between images within a Brain_Data() instance.
 
             Args:
-                method: (str) type of distance metric (can use any scikit learn or
+                metric: (str) type of distance metric (can use any scikit learn or
                         sciypy metric)
 
             Returns:
@@ -858,7 +909,7 @@ def distance(self, method='euclidean', **kwargs):
 
         """
 
-        return Adjacency(pairwise_distances(self.data, metric=method, **kwargs),
+        return Adjacency(pairwise_distances(self.data, metric=metric, **kwargs),
                          matrix_type='Distance')
 
     def multivariate_similarity(self, images, method='ols'):
@@ -1183,67 +1234,114 @@ def predict_multi(self, algorithm=None, cv_dict=None, method='searchlight', rois
             out = Brain_Data(out, mask=self.mask)
         return out
 
-    def apply_mask(self, mask):
+    def apply_mask(self, mask, resample_mask_to_brain=False):
         """ Mask Brain_Data instance
 
+        Note target data will be resampled into the same space as the mask. If you would like the mask
+        resampled into the Brain_Data space, then set resample_mask_to_brain=True.
+
         Args:
-            mask: (Brain_Data or nifti object) mask to apply to Brain_Data object
+            mask: (Brain_Data or nifti object) mask to apply to Brain_Data object.
+            resample_mask_to_brain: (bool) Will resample mask to brain space before applying mask (default=False).
 
         Returns:
             masked: (Brain_Data) masked Brain_Data object
 
         """
 
-        if isinstance(mask, Brain_Data):
-            mask = mask.to_nifti()  # convert to nibabel
+        masked = deepcopy(self)
+        mask = check_brain_data(mask)
+        if not check_brain_data_is_single(mask):
+            raise ValueError('Mask must be a single image')
+
+        if check_brain_data_is_single(self):
+            n_vox = len(self)
+        else:
+            n_vox = self.shape()[1]
 
-        if not isinstance(mask, nib.Nifti1Image):
-            if isinstance(mask, six.string_types):
-                if os.path.isfile(mask):
-                    mask = nib.load(mask)
-                if not ((self.mask.get_affine() == mask.get_affine()).all()) & (self.mask.shape[0:3] == mask.shape[0:3]):
-                    mask = resample_img(mask, target_affine=self.mask.get_affine(), target_shape=self.mask.shape)
-            else:
-                raise ValueError("Mask is not a nibabel instance, Brain_Data "
-                                 "instance, or a valid file name.")
+        if resample_mask_to_brain: 
+            mask = resample_to_img(mask.to_nifti(), masked.to_nifti())
+            mask = check_brain_data(mask, masked.mask)
 
-        masked = deepcopy(self)
-        nifti_masker = NiftiMasker(mask_img=mask)
-        masked.data = nifti_masker.fit_transform(self.to_nifti())
-        masked.nifti_masker = nifti_masker
+        nifti_masker = NiftiMasker(mask_img=mask.to_nifti()).fit()
+
+        if n_vox == len(mask):
+            if check_brain_data_is_single(masked):
+                masked.data = masked.data[mask.data.astype(bool)]
+            else:
+                masked.data = masked.data[:, mask.data.astype(bool)]
+            masked.nifti_masker = nifti_masker
+        else:
+            masked.data = nifti_masker.fit_transform(masked.to_nifti())
+            masked.nifti_masker = nifti_masker
         if (len(masked.shape()) > 1) & (masked.shape()[0] == 1):
             masked.data = masked.data.flatten()
         return masked
 
-    def extract_roi(self, mask, method='mean'):
+    def extract_roi(self, mask, metric='mean', n_components=None):
         """ Extract activity from mask
 
         Args:
-            mask: (nifiti) nibabel mask can be binary or numbered for
+            mask: (nifti) nibabel mask can be binary or numbered for
                   different rois
-            method: type of extraction method (default=mean)
+            metric: type of extraction method ['mean', 'median', 'pca'], (default=mean)
                     NOTE: Only mean currently works!
+            n_components: if metric='pca', number of components to return (takes any input into sklearn.Decomposition.PCA)
 
         Returns:
             out: mean within each ROI across images
 
         """
+
+        metrics = ['mean','median','pca']
+
         mask = check_brain_data(mask)
         ma = mask.copy()
 
-        if method != 'mean':
-            raise ValueError('Only mean is currently implemented.')
+        if metric not in metrics:
+            raise NotImplementedError
 
         if len(np.unique(ma.data)) == 2:
-            out = np.mean(self.data[:, np.where(ma.data)].squeeze(), axis=1)
+            masked = self.apply_mask(ma)
+            if check_brain_data_is_single(masked):
+                if metric == 'mean':
+                    out = masked.mean()
+                elif metric == 'median':
+                    out = masked.median()
+                else:
+                    raise ValueError('Not possible to run PCA on a single image')
+            else:
+                if metric == 'mean':
+                    out = masked.mean(axis=1)
+                elif metric == 'median':
+                    out = masked.median(axis=1)
+                else:
+                    output = masked.decompose(algorithm='pca', n_components=n_components, axis='images')
+                    out = output['weights'].T
         elif len(np.unique(ma.data)) > 2:
             # make sure each ROI id is an integer
             ma.data = np.round(ma.data).astype(int)
             all_mask = expand_mask(ma)
-            out = []
-            for i in range(all_mask.shape()[0]):
-                out.append(np.mean(self.data[:, np.where(all_mask[i].data)].squeeze(), axis=1))
-            out = np.array(out)
+            if check_brain_data_is_single(self):
+                if metric == 'mean':
+                    out = np.array([self.apply_mask(m).mean() for m in all_mask])
+                elif metric == 'median':
+                    out = np.array([self.apply_mask(m).median() for m in all_mask])
+                else:
+                    raise ValueError('Not possible to run PCA on a single image')
+            else:
+                if metric == 'mean':
+                    out = np.array([self.apply_mask(m).mean(axis=1) for m in all_mask])
+                elif metric == 'median':
+                    out = np.array([self.apply_mask(m).median(axis=1) for m in all_mask])
+                else:
+                    out = []
+                    for m in all_mask:
+                        masked = self.apply_mask(m)
+                        output = masked.decompose(algorithm='pca', n_components=n_components, axis='images')
+                        out.append(output['weights'].T)
+        else:
+            raise ValueError('Mask must be binary or integers')
         return out
 
     def icc(self, icc_type='icc2'):
@@ -1654,7 +1752,7 @@ def decompose(self, algorithm='pca', axis='voxels', n_components=None,
 
         Args:
             algorithm: (str) Algorithm to perform decomposition
-                        types=['pca','ica','nnmf','fa']
+                        types=['pca','ica','nnmf','fa','dictionary','kernelpca']
             axis: dimension to decompose ['voxels','images']
             n_components: (int) number of components. If None then retain
                         as many as possible.
diff --git a/nltools/mask.py b/nltools/mask.py
index 17380aa7..8399ae7d 100644
--- a/nltools/mask.py
+++ b/nltools/mask.py
@@ -178,9 +178,9 @@ def roi_to_brain(data, mask_x):
     must correspond to ROI numbers.
 
     This is useful for populating a parcellation scheme by a vector of Values
-    
+
     Args:
-        data: Pandas series or dataframe of ROI by observation
+        data: Pandas series, dataframe, list, np.array of ROI by observation
         mask_x: an expanded binary mask
     Returns:
         out: (Brain_Data) Brain_Data instance where each ROI is now populated
@@ -197,6 +197,32 @@ def series_to_brain(data, mask_x):
             raise ValueError('Data must have the same number of rows as mask has ROIs.')
         return Brain_Data([mask_x[x]*data[x] for x in data.keys()]).sum()
 
+    if not isinstance(data, (pd.Series, pd.DataFrame)):
+        if isinstance(data, list):
+            if len(data) != len(mask_x):
+                raise ValueError('Data must have the same number of rows as mask has ROIs.')
+            else:
+                data = pd.Series(data)
+        elif isinstance(data, np.ndarray):
+            if len(data.shape) == 1:
+                if len(data) != len(mask_x):
+                    raise ValueError('Data must have the same number of rows as mask has ROIs.')
+                else:
+                    data = pd.Series(data)
+            elif len(data.shape) == 2:
+                data = pd.DataFrame(data)
+                if data.shape[0] != len(mask_x):
+                    if data.shape[1] == len(mask_x):
+                        data = data.T
+                    else:
+                        raise ValueError('Data must have the same number of rows as rois in mask')
+            else:
+                raise NotImplementedError
+
+        else:
+            raise NotImplementedError
+
+
     if len(mask_x) != data.shape[0]:
         raise ValueError('Data must have the same number of rows as mask has ROIs.')
 
@@ -204,3 +230,5 @@ def series_to_brain(data, mask_x):
         return series_to_brain(data, mask_x)
     elif isinstance(data, pd.DataFrame):
         return Brain_Data([series_to_brain(data[x], mask_x) for x in data.keys()])
+    else:
+        raise ValueError("Data must be a pandas series or data frame.")
diff --git a/nltools/plotting.py b/nltools/plotting.py
index 623cc68f..d221a45e 100644
--- a/nltools/plotting.py
+++ b/nltools/plotting.py
@@ -26,10 +26,14 @@
 import seaborn as sns
 import matplotlib.pyplot as plt
 import numpy as np
+from numpy.fft import fft, fftfreq
 from nltools.stats import two_sample_permutation, one_sample_permutation
 from nilearn.plotting import plot_glass_brain, plot_stat_map, view_img, view_img_on_surf
 from nltools.prefs import MNI_Template, resolve_mni_path
 from nltools.utils import attempt_to_import
+from ipywidgets import BoundedFloatText, BoundedIntText
+from ipywidgets import interact
+import sklearn
 import warnings
 import os
 
@@ -350,7 +354,7 @@ def dist_from_hyperplane_plot(stats_output):
     """
 
     if "dist_from_hyperplane_xval" in stats_output.columns:
-        fig = sns.factorplot(
+        sns.factorplot(
             "subject_id",
             "dist_from_hyperplane_xval",
             hue="Y",
@@ -358,7 +362,7 @@ def dist_from_hyperplane_plot(stats_output):
             kind="point",
         )
     else:
-        fig = sns.factorplot(
+        sns.factorplot(
             "subject_id",
             "dist_from_hyperplane_all",
             hue="Y",
@@ -368,7 +372,7 @@ def dist_from_hyperplane_plot(stats_output):
     plt.xlabel("Subject", fontsize=16)
     plt.ylabel("Distance from Hyperplane", fontsize=16)
     plt.title("Classification", fontsize=18)
-    return fig
+    return 
 
 
 def scatterplot(stats_output):
@@ -383,13 +387,13 @@ def scatterplot(stats_output):
     """
 
     if "yfit_xval" in stats_output.columns:
-        fig = sns.lmplot("Y", "yfit_xval", data=stats_output)
+        sns.lmplot("Y", "yfit_xval", data=stats_output)
     else:
-        fig = sns.lmplot("Y", "yfit_all", data=stats_output)
+        sns.lmplot("Y", "yfit_all", data=stats_output)
     plt.xlabel("Y", fontsize=16)
     plt.ylabel("Predicted Value", fontsize=16)
     plt.title("Prediction", fontsize=18)
-    return fig
+    return 
 
 
 def probability_plot(stats_output):
@@ -403,13 +407,13 @@ def probability_plot(stats_output):
 
     """
     if "Probability_xval" in stats_output.columns:
-        fig = sns.lmplot("Y", "Probability_xval", data=stats_output, logistic=True)
+        sns.lmplot("Y", "Probability_xval", data=stats_output, logistic=True)
     else:
-        fig = sns.lmplot("Y", "Probability_all", data=stats_output, logistic=True)
+        sns.lmplot("Y", "Probability_all", data=stats_output, logistic=True)
     plt.xlabel("Y", fontsize=16)
     plt.ylabel("Predicted Probability", fontsize=16)
     plt.title("Prediction", fontsize=18)
-    return fig
+    return 
 
     # # and plot the result
     # plt.figure(1, figsize=(4, 3))
@@ -435,13 +439,13 @@ def roc_plot(fpr, tpr):
 
     """
 
-    fig = plt.figure()
+    plt.figure()
     plt.plot(fpr, tpr, color="red", linewidth=3)
     # fig = sns.tsplot(tpr,fpr,color='red',linewidth=3)
     plt.xlabel("(1 - Specificity)", fontsize=16)
     plt.ylabel("Sensitivity", fontsize=16)
     plt.title("ROC Plot", fontsize=18)
-    return fig
+    return
 
 
 def plot_stacked_adjacency(adjacency1, adjacency2, normalize=True, **kwargs):
@@ -604,9 +608,9 @@ def plot_between_label_distance(
             ]["Distance"].mean()
 
     if ax is None:
-        f, ax = plt.subplots(1)
+        _, ax = plt.subplots(1)
     else:
-        f = plt.figure()
+        plt.figure()
 
     if permutation_test:
         mn_dist_out = pd.DataFrame(
@@ -640,10 +644,10 @@ def plot_between_label_distance(
             ax=ax,
             cbar=False,
         )
-        return (f, out, within_dist_out, mn_dist_out, p_dist_out)
+        return (out, within_dist_out, mn_dist_out, p_dist_out)
     else:
-        f = sns.heatmap(within_dist_out, ax=ax, square=True, **kwargs)
-        return (f, out, within_dist_out)
+        sns.heatmap(within_dist_out, ax=ax, square=True, **kwargs)
+        return (out, within_dist_out)
 
 
 def plot_silhouette(
@@ -703,9 +707,9 @@ def plot_silhouette(
     # Plot
     with sns.axes_style("white"):
         if ax is None:
-            f, ax = plt.subplots(1, figsize=figsize)
+            _, ax = plt.subplots(1, figsize=figsize)
         else:
-            f = plt.plot(figsize=figsize)
+            plt.plot(figsize=figsize)
     x_lower = 10
     labelX = []
     for labelInd in range(n_clusters):
@@ -749,7 +753,56 @@ def plot_silhouette(
             else:
                 temp["p"] = 999
             outAll = outAll.append(temp)
-        return (f, outAll)
+        return outAll
     else:
-        return f
+        return
+
+def component_viewer(output, tr=2.0):
+    ''' This a function to interactively view the results of a decomposition analysis
+
+    Args:
+        output: (dict) output dictionary from running Brain_data.decompose()
+        tr: (float) repetition time of data
+    '''
+
+    def component_inspector(component, threshold):
+        '''This a function to be used with ipywidgets to interactively view a decomposition analysis
+
+            Make sure you have tr and output assigned to variables.
+
+            Example:
+
+                from ipywidgets import BoundedFloatText, BoundedIntText
+                from ipywidgets import interact
+
+                tr = 2.4
+                output = data_filtered_smoothed.decompose(algorithm='ica', n_components=30, axis='images', whiten=True)
+
+                interact(component_inspector, component=BoundedIntText(description='Component', value=0, min=0, max=len(output['components'])-1),
+                      threshold=BoundedFloatText(description='Threshold', value=2.0, min=0, max=4, step=.1))
+
+        '''
+        _, ax = plt.subplots(nrows=3, figsize=(12,8))
+        thresholded = (output['components'][component] - output['components'][component].mean())*(1/output['components'][component].std())
+        thresholded.data[np.abs(thresholded.data) <= threshold] = 0
+        plot_stat_map(thresholded.to_nifti(), cut_coords=range(-40, 70, 10),
+                      display_mode='z', black_bg=True, colorbar=True, annotate=False,
+                      draw_cross=False, axes=ax[0])
+        if isinstance(output['decomposition_object'], (sklearn.decomposition.PCA)):
+            var_exp = output['decomposition_object'].explained_variance_ratio_[component]
+            ax[0].set_title(f"Component: {component}/{len(output['components'])}, Variance Explained: {var_exp:2.2}", fontsize=18)
+        else:
+            ax[0].set_title(f"Component: {component}/{len(output['components'])}", fontsize=18)
+
+        ax[1].plot(output['weights'][:, component], linewidth=2, color='red')
+        ax[1].set_ylabel('Intensity (AU)', fontsize=18)
+        ax[1].set_title(f'Timecourse (TR={tr})', fontsize=16)
+        y = fft(output['weights'][:, component])
+        f = fftfreq(len(y), d=tr)
+        ax[2].plot(f[f > 0], np.abs(y)[f > 0]**2)
+        ax[2].set_ylabel('Power', fontsize=18)
+        ax[2].set_xlabel('Frequency (Hz)', fontsize=16)
+
+    interact(component_inspector, component=BoundedIntText(description='Component', value=0, min=0, max=len(output['components'])-1),
+              threshold=BoundedFloatText(description='Threshold', value=2.0, min=0, max=4, step=.1))
 
diff --git a/nltools/stats.py b/nltools/stats.py
index ffa07dfe..f248e322 100644
--- a/nltools/stats.py
+++ b/nltools/stats.py
@@ -989,6 +989,7 @@ def regress(X, Y, mode='ols', stats='full', **kwargs):
     if mode == 'ols' or mode == 'robust':
 
         b = np.dot(np.linalg.pinv(X), Y)
+        
         # Return betas and stop other computations if that's all that's requested
         if stats == 'betas':
             return b.squeeze()
@@ -1006,7 +1007,10 @@ def regress(X, Y, mode='ols', stats='full', **kwargs):
             axis_func = [_robust_estimator, 0, res, X, robust_estimator, nlags]
             stderr = np.apply_along_axis(*axis_func)
 
-        t = b / stderr
+        # Then only compute t-stats at voxels where the standard error is at least .000001
+        t = np.zeros_like(b)
+        t[stderr > 1.e-6] = b[stderr > 1.e-6] / stderr[stderr > 1.e-6]
+
         # Return betas and ts and stop other computations if that's all that's requested
         if stats == 'tstats':
             return b.squeeze(), t.squeeze()
@@ -1076,7 +1080,7 @@ def regress_permutation(X, Y, n_permute=5000, tail=2, random_state=None, verbose
     # We could optionally Save (X.T * X)^-1 * X.T so we dont have to invert each permutation, but this would require not relying on regress() and because the second-level design mat is probably on the small side we might not actually save that much time
     # inv = np.linalg.pinv(X)
 
-    for i in range(n_permute):
+    for _ in range(n_permute):
         _, _t = regress(func(X.values), Y, stats='tstats', **kwargs)
         if tail == 2:
             p += np.abs(_t) >= np.abs(t)
@@ -1251,7 +1255,7 @@ def align(data, method='deterministic_srm', n_features=None, axis=0,
 
 def procrustes(data1, data2):
     '''Procrustes analysis, a similarity test for two data sets.
-
+    
     Each input matrix is a set of points or vectors (the rows of the matrix).
     The dimension of the space is the number of columns of each matrix. Given
     two identically sized matrices, procrustes standardizes both such that:
diff --git a/nltools/tests/test_adjacency.py b/nltools/tests/test_adjacency.py
index 1247d1f7..4ecbbb2d 100644
--- a/nltools/tests/test_adjacency.py
+++ b/nltools/tests/test_adjacency.py
@@ -14,22 +14,18 @@ def test_type_single(sim_adjacency_single):
     assert dat_single2.matrix_type == 'similarity'
     assert sim_adjacency_single.issymmetric
 
-
 def test_type_directed(sim_adjacency_directed):
     assert not sim_adjacency_directed.issymmetric
 
-
 def test_length(sim_adjacency_multiple):
     assert len(sim_adjacency_multiple) == sim_adjacency_multiple.data.shape[0]
     assert len(sim_adjacency_multiple[0]) == 1
 
-
 def test_indexing(sim_adjacency_multiple):
     assert len(sim_adjacency_multiple[0]) == 1
     assert len(sim_adjacency_multiple[0:4]) == 4
     assert len(sim_adjacency_multiple[0, 2, 3]) == 3
 
-
 def test_arithmetic(sim_adjacency_directed):
     assert(sim_adjacency_directed+5).data[0] == sim_adjacency_directed.data[0]+5
     assert(sim_adjacency_directed-.5).data[0] == sim_adjacency_directed.data[0]-.5
@@ -39,16 +35,13 @@ def test_arithmetic(sim_adjacency_directed):
     assert np.all(np.isclose((sim_adjacency_directed*2 - sim_adjacency_directed).data,
                              sim_adjacency_directed.data))
 
-
 def test_copy(sim_adjacency_multiple):
     assert np.all(sim_adjacency_multiple.data == sim_adjacency_multiple.copy().data)
 
-
 def test_squareform(sim_adjacency_multiple):
     assert len(sim_adjacency_multiple.squareform()) == len(sim_adjacency_multiple)
     assert sim_adjacency_multiple[0].squareform().shape == sim_adjacency_multiple[0].square_shape()
 
-
 def test_write_multiple(sim_adjacency_multiple, tmpdir):
     sim_adjacency_multiple.write(os.path.join(str(tmpdir.join('Test.csv'))),
                                  method='long')
@@ -187,15 +180,12 @@ def test_bootstrap(sim_adjacency_multiple):
 
 
 def test_plot(sim_adjacency_multiple):
-    f = sim_adjacency_multiple[0].plot()
-    assert isinstance(f, plt.Figure)
-    f = sim_adjacency_multiple.plot()
-    assert isinstance(f, plt.Figure)
+    sim_adjacency_multiple[0].plot()
+    sim_adjacency_multiple.plot()
 
 
 def test_plot_mds(sim_adjacency_single):
-    f = sim_adjacency_single.plot_mds()
-    assert isinstance(f, plt.Figure)
+    sim_adjacency_single.plot_mds()
 
 
 def test_similarity_conversion(sim_adjacency_single):
@@ -231,7 +221,6 @@ def test_regression():
     assert np.allclose(np.array([out['Group1'], out['Group2']]), np.array([1, 0]), rtol=1e-01)  # np.allclose(np.sum(stats['beta']-np.array([1,2,3])),0)
 
 
-
 def test_social_relations_model():
     data = Adjacency(np.array([[np.nan, 8, 5, 10],
                     [7, np.nan, 7, 6],
diff --git a/nltools/tests/test_brain_data.py b/nltools/tests/test_brain_data.py
index 302225e6..f42193d2 100644
--- a/nltools/tests/test_brain_data.py
+++ b/nltools/tests/test_brain_data.py
@@ -1,4 +1,5 @@
 import os
+import pytest
 import numpy as np
 import nibabel as nb
 import pandas as pd
@@ -7,7 +8,7 @@
                           Adjacency,
                           Groupby)
 from nltools.stats import threshold, align
-from nltools.mask import create_sphere
+from nltools.mask import create_sphere, roi_to_brain
 from nltools.utils import get_resource_path
 from nltools.mask import expand_mask
 # from nltools.prefs import MNI_Template
@@ -69,20 +70,29 @@ def test_load(tmpdir):
 def test_shape(sim_brain_data):
     assert sim_brain_data.shape() == shape_2d
 
-
 def test_mean(sim_brain_data):
     assert sim_brain_data.mean().shape()[0] == shape_2d[1]
-
+    assert sim_brain_data.mean().shape()[0] == shape_2d[1]
+    assert len(sim_brain_data.mean(axis=1)) == shape_2d[0]
+    with pytest.raises(ValueError):
+        sim_brain_data.mean(axis='1')
+    assert isinstance(sim_brain_data[0].mean(), (float, np.floating))
+
+def test_median(sim_brain_data):
+    assert sim_brain_data.median().shape()[0] == shape_2d[1]
+    assert sim_brain_data.median().shape()[0] == shape_2d[1]
+    assert len(sim_brain_data.median(axis=1)) == shape_2d[0]
+    with pytest.raises(ValueError):
+        sim_brain_data.median(axis='1')
+    assert isinstance(sim_brain_data[0].median(), (float, np.floating))
 
 def test_std(sim_brain_data):
     assert sim_brain_data.std().shape()[0] == shape_2d[1]
 
-
 def test_sum(sim_brain_data):
     s = sim_brain_data.sum()
     assert s.shape() == sim_brain_data[1].shape()
 
-
 def test_add(sim_brain_data):
     new = sim_brain_data + sim_brain_data
     assert new.shape() == shape_2d
@@ -118,25 +128,23 @@ def test_indexing(sim_brain_data):
     assert d.shape[0:3] == shape_3d
     assert Brain_Data(d)
 
-
 def test_concatenate(sim_brain_data):
     out = Brain_Data([x for x in sim_brain_data])
     assert isinstance(out, Brain_Data)
     assert len(out) == len(sim_brain_data)
 
-
 def test_append(sim_brain_data):
     assert sim_brain_data.append(sim_brain_data).shape()[0] == shape_2d[0]*2
 
-
 def test_ttest(sim_brain_data):
     out = sim_brain_data.ttest()
     assert out['t'].shape()[0] == shape_2d[1]
-    distance = sim_brain_data.distance(method='correlation')
+
+def test_distance(sim_brain_data):
+    distance = sim_brain_data.distance(metric='correlation')
     assert isinstance(distance, Adjacency)
     assert distance.square_shape()[0] == shape_2d[0]
 
-
 def test_regress(sim_brain_data):
     sim_brain_data.X = pd.DataFrame({'Intercept': np.ones(len(sim_brain_data.Y)),
                                      'X1': np.array(sim_brain_data.Y).flatten()}, index=None)
@@ -163,7 +171,6 @@ def test_regress(sim_brain_data):
     tt = threshold(out['t'][i], out['p'][i], .05)
     assert isinstance(tt, Brain_Data)
 
-
 def test_randomise(sim_brain_data):
     sim_brain_data.X = pd.DataFrame({'Intercept': np.ones(len(sim_brain_data.Y))})
 
@@ -192,28 +199,49 @@ def test_apply_mask(sim_brain_data):
     assert isinstance(s1, nb.Nifti1Image)
     masked_dat = sim_brain_data.apply_mask(s1)
     assert masked_dat.shape()[1] == np.sum(s1.get_data() != 0)
+    masked_dat = sim_brain_data.apply_mask(s1, resample_mask_to_brain=True)
+    assert masked_dat.shape()[1] == np.sum(s1.get_data() != 0)
 
 
 def test_extract_roi(sim_brain_data):
     mask = create_sphere([12, 10, -8], radius=10)
-    assert len(sim_brain_data.extract_roi(mask)) == shape_2d[0]
-
+    assert len(sim_brain_data.extract_roi(mask, metric='mean')) == shape_2d[0]
+    assert len(sim_brain_data.extract_roi(mask, metric='median')) == shape_2d[0]
+    n_components = 2
+    assert sim_brain_data.extract_roi(mask, metric='pca', n_components=n_components).shape == (n_components, shape_2d[0])
+    with pytest.raises(NotImplementedError):
+        sim_brain_data.extract_roi(mask, metric='p')
+
+    assert isinstance(sim_brain_data[0].extract_roi(mask, metric='mean'), (float, np.floating))
+    assert isinstance(sim_brain_data[0].extract_roi(mask, metric='median'), (float, np.floating))
+    with pytest.raises(ValueError):
+        sim_brain_data[0].extract_roi(mask, metric='pca')
+    with pytest.raises(NotImplementedError):
+        sim_brain_data[0].extract_roi(mask, metric='p')
+
+    s1 = create_sphere([15, 10, -8], radius=10)
+    s2 = create_sphere([-15, 10, -8], radius=10)
+    s3 = create_sphere([0, -15, -8], radius=10)
+    masks = Brain_Data([s1, s2, s3])
+    mask = roi_to_brain([1,2,3], masks)
+    assert len(sim_brain_data[0].extract_roi(mask, metric='mean')) == len(masks)
+    assert len(sim_brain_data[0].extract_roi(mask, metric='median')) == len(masks)
+    assert sim_brain_data.extract_roi(mask, metric='mean').shape == (len(masks), shape_2d[0])
+    assert sim_brain_data.extract_roi(mask, metric='median').shape == (len(masks), shape_2d[0])
+    assert len(sim_brain_data.extract_roi(mask, metric='pca', n_components=n_components)) == len(masks)
 
 def test_r_to_z(sim_brain_data):
     z = sim_brain_data.r_to_z()
     assert z.shape() == sim_brain_data.shape()
 
-
 def test_copy(sim_brain_data):
     d_copy = sim_brain_data.copy()
     assert d_copy.shape() == sim_brain_data.shape()
 
-
 def test_detrend(sim_brain_data):
     detrend = sim_brain_data.detrend()
     assert detrend.shape() == sim_brain_data.shape()
 
-
 def test_standardize(sim_brain_data):
     s = sim_brain_data.standardize()
     assert s.shape() == sim_brain_data.shape()
@@ -237,7 +265,6 @@ def test_groupby_aggregate(sim_brain_data):
     assert isinstance(mn, Brain_Data)
     assert len(mn.shape()) == 1
 
-
 def test_threshold():
     s1 = create_sphere([12, 10, -8], radius=10)
     s2 = create_sphere([22, -2, -22], radius=10)
diff --git a/nltools/tests/test_mask.py b/nltools/tests/test_mask.py
index 652ea80e..524f79b2 100644
--- a/nltools/tests/test_mask.py
+++ b/nltools/tests/test_mask.py
@@ -1,5 +1,7 @@
-from nltools.mask import create_sphere
+from nltools.mask import create_sphere, roi_to_brain
+from nltools.data import Brain_Data
 import numpy as np
+import pandas as pd
 
 
 def test_create_sphere():
@@ -11,3 +13,32 @@ def test_create_sphere():
     assert np.sum(a.get_data()) >= 553  # 571
     a = create_sphere(radius=10, coordinates=[[0, 0, 0], [15, 0, 25]])
     assert np.sum(a.get_data()) >= 1013  # 1051
+
+
+def test_roi_to_brain():
+    s1 = create_sphere([15, 10, -8], radius=10)
+    s2 = create_sphere([-15, 10, -8], radius=10)
+    s3 = create_sphere([0, -15, -8], radius=10)
+    masks = Brain_Data([s1,s2,s3])
+
+    d = [1,2,3]
+    m = roi_to_brain(d, masks)
+    assert np.all([np.any(m.data==x) for x in d])
+
+    d = pd.Series([1.1, 2.1, 3.1])
+    m = roi_to_brain(d, masks)
+    assert np.all([np.any(m.data==x) for x in d])
+
+    d = np.array([1, 2, 3])
+    m = roi_to_brain(d, masks)
+    assert np.all([np.any(m.data==x) for x in d])
+
+    d = pd.DataFrame([np.ones(10)*x for x in [1, 2, 3]])
+    m = roi_to_brain(d, masks)
+    assert len(m) == d.shape[1]
+    assert np.all([np.any(m[0].data==x) for x in d[0]])
+
+    d = np.array([np.ones(10)*x for x in [1, 2, 3]])
+    m = roi_to_brain(d, masks)
+    assert len(m) == d.shape[1]
+    assert np.all([np.any(m[0].data==x) for x in d[0]])
\ No newline at end of file
diff --git a/nltools/tests/test_stats.py b/nltools/tests/test_stats.py
index 63d8cd11..785851a8 100644
--- a/nltools/tests/test_stats.py
+++ b/nltools/tests/test_stats.py
@@ -140,7 +140,7 @@ def test_align():
     assert len(data) == len(out['transformation_matrix'])
     assert data[0].shape == out['common_model'].shape
     transformed = np.dot(data[0].T, out['transformation_matrix'][0])
-    np.testing.assert_almost_equal(0, np.sum(out['transformed'][0]-transformed.T))
+    np.testing.assert_almost_equal(np.sum(out['transformed'][0]-transformed.T), 0, decimal=5)
     assert len(out['isc']) == out['transformed'][0].shape[0]
 
     out = align(data, method='probabilistic_srm')
@@ -148,7 +148,7 @@ def test_align():
     assert len(data) == len(out['transformation_matrix'])
     assert data[0].shape == out['common_model'].shape
     transformed = np.dot(data[0].T, out['transformation_matrix'][0])
-    np.testing.assert_almost_equal(0, np.sum(out['transformed'][0]-transformed.T))
+    np.testing.assert_almost_equal(np.sum(out['transformed'][0]-transformed.T), 0, decimal=5)
     assert len(out['isc']) == out['transformed'][0].shape[0]
 
     out2 = align(data, method='procrustes')
@@ -158,7 +158,7 @@ def test_align():
     assert len(data) == len(out2['disparity'])
     centered = data[0].T-np.mean(data[0].T, 0)
     transformed = (np.dot(centered/np.linalg.norm(centered), out2['transformation_matrix'][0])*out2['scale'][0])
-    np.testing.assert_almost_equal(0, np.sum(out2['transformed'][0]-transformed.T))
+    np.testing.assert_almost_equal(np.sum(out2['transformed'][0]-transformed.T), 0, decimal=5)
     assert out['transformed'][0].shape == out2['transformed'][0].shape
     assert out['transformation_matrix'][0].shape == out2['transformation_matrix'][0].shape
     assert len(out['isc']) == out['transformed'][0].shape[0]
@@ -170,7 +170,7 @@ def test_align():
     assert len(data) == len(out['transformation_matrix'])
     assert data[0].shape() == out['common_model'].shape()
     transformed = np.dot(d1.data, out['transformation_matrix'][0])
-    np.testing.assert_almost_equal(0, np.sum(out['transformed'][0].data-transformed))
+    np.testing.assert_almost_equal(np.sum(out['transformed'][0].data-transformed), 0, decimal=5)
     assert len(out['isc']) == out['transformed'][0].shape()[1]
 
     out = align(data, method='probabilistic_srm')
@@ -178,7 +178,7 @@ def test_align():
     assert len(data) == len(out['transformation_matrix'])
     assert data[0].shape() == out['common_model'].shape()
     transformed = np.dot(d1.data, out['transformation_matrix'][0])
-    np.testing.assert_almost_equal(0, np.sum(out['transformed'][0].data-transformed))
+    np.testing.assert_almost_equal(np.sum(out['transformed'][0].data-transformed), 0, decimal=5)
     assert len(out['isc']) == out['transformed'][0].shape()[1]
 
     out2 = align(data, method='procrustes')
@@ -188,7 +188,7 @@ def test_align():
     assert len(data) == len(out2['disparity'])
     centered = data[0].data-np.mean(data[0].data, 0)
     transformed = (np.dot(centered/np.linalg.norm(centered), out2['transformation_matrix'][0])*out2['scale'][0])
-    np.testing.assert_almost_equal(0, np.sum(out2['transformed'][0].data-transformed))
+    np.testing.assert_almost_equal(np.sum(out2['transformed'][0].data-transformed), 0, decimal=5)
     assert out['transformed'][0].shape() == out2['transformed'][0].shape()
     assert out['transformation_matrix'][0].shape == out2['transformation_matrix'][0].shape
     assert len(out['isc']) == out['transformed'][0].shape()[1]
@@ -208,7 +208,7 @@ def test_align():
     assert len(data) == len(out['transformation_matrix'])
     assert data[0].shape == out['common_model'].shape
     transformed = np.dot(data[0], out['transformation_matrix'][0])
-    np.testing.assert_almost_equal(0, np.sum(out['transformed'][0]-transformed))
+    np.testing.assert_almost_equal(np.sum(out['transformed'][0]-transformed), 0, decimal=5)
     assert len(out['isc']) == out['transformed'][0].shape[1]
 
     out = align(data, method='probabilistic_srm', axis=1)
@@ -216,7 +216,7 @@ def test_align():
     assert len(data) == len(out['transformation_matrix'])
     assert data[0].shape == out['common_model'].shape
     transformed = np.dot(data[0], out['transformation_matrix'][0])
-    np.testing.assert_almost_equal(0, np.sum(out['transformed'][0]-transformed))
+    np.testing.assert_almost_equal(np.sum(out['transformed'][0]-transformed), 0, decimal=5)
     assert len(out['isc']) == out['transformed'][0].shape[1]
 
     out2 = align(data, method='procrustes', axis=1)
@@ -226,7 +226,7 @@ def test_align():
     assert len(data) == len(out2['disparity'])
     centered = data[0]-np.mean(data[0], 0)
     transformed = (np.dot(centered/np.linalg.norm(centered), out2['transformation_matrix'][0])*out2['scale'][0])
-    np.testing.assert_almost_equal(0, np.sum(out2['transformed'][0]-transformed))
+    np.testing.assert_almost_equal(np.sum(out2['transformed'][0]-transformed), 0, decimal=5)
     assert out['transformed'][0].shape == out2['transformed'][0].shape
     assert out['transformation_matrix'][0].shape == out2['transformation_matrix'][0].shape
     assert len(out['isc']) == out['transformed'][0].shape[1]
@@ -238,7 +238,7 @@ def test_align():
     assert len(data) == len(out['transformation_matrix'])
     assert data[0].shape() == out['common_model'].shape()
     transformed = np.dot(d1.data.T, out['transformation_matrix'][0])
-    np.testing.assert_almost_equal(0, np.sum(out['transformed'][0].data-transformed.T))
+    np.testing.assert_almost_equal(np.sum(out['transformed'][0].data-transformed.T), 0, decimal=5)
     assert len(out['isc']) == out['transformed'][0].shape()[0]
 
     out = align(data, method='probabilistic_srm', axis=1)
@@ -246,7 +246,7 @@ def test_align():
     assert len(data) == len(out['transformation_matrix'])
     assert data[0].shape() == out['common_model'].shape()
     transformed = np.dot(d1.data.T, out['transformation_matrix'][0])
-    np.testing.assert_almost_equal(0, np.sum(out['transformed'][0].data-transformed.T))
+    np.testing.assert_almost_equal(np.sum(out['transformed'][0].data-transformed.T), 0, decimal=5)
     assert len(out['isc']) == out['transformed'][0].shape()[0]
 
     out2 = align(data, method='procrustes', axis=1)
@@ -256,7 +256,7 @@ def test_align():
     assert len(data) == len(out2['disparity'])
     centered = data[0].data.T-np.mean(data[0].data.T, 0)
     transformed = (np.dot(centered/np.linalg.norm(centered), out2['transformation_matrix'][0])*out2['scale'][0])
-    np.testing.assert_almost_equal(0, np.sum(out2['transformed'][0].data-transformed.T))
+    np.testing.assert_almost_equal(np.sum(out2['transformed'][0].data-transformed.T), 0, decimal=5)
     assert out['transformed'][0].shape() == out2['transformed'][0].shape()
     assert out['transformation_matrix'][0].shape == out2['transformation_matrix'][0].shape
     assert len(out['isc']) == out['transformed'][0].shape()[0]
diff --git a/nltools/tests/test_utils.py b/nltools/tests/test_utils.py
new file mode 100644
index 00000000..6ece6c53
--- /dev/null
+++ b/nltools/tests/test_utils.py
@@ -0,0 +1,17 @@
+
+from nltools.utils import check_brain_data, check_brain_data_is_single
+from nltools.mask import create_sphere
+from nltools.data import Brain_Data
+import numpy as np
+
+def test_check_brain_data(sim_brain_data):
+    mask = Brain_Data(create_sphere([15, 10, -8], radius=10))
+    a = check_brain_data(sim_brain_data)
+    assert isinstance(a, Brain_Data)
+    b = check_brain_data(sim_brain_data, mask=mask)
+    assert isinstance(b, Brain_Data)
+    assert b.shape()[1] == np.sum(mask.data==1)
+
+def test_check_brain_data_is_single(sim_brain_data):
+    assert not check_brain_data_is_single(sim_brain_data)
+    assert check_brain_data_is_single(sim_brain_data[0])
\ No newline at end of file
diff --git a/nltools/utils.py b/nltools/utils.py
index 780dcecb..899d15d6 100644
--- a/nltools/utils.py
+++ b/nltools/utils.py
@@ -179,7 +179,8 @@ def set_decomposition_algorithm(algorithm, n_components=None, *args, **kwargs):
     Args:
         algorithm: The decomposition algorithm to use. Either a string or an
                     (uninitialized) scikit-learn decomposition object.
-                    If string must be one of 'pca','nnmf', ica','fa'
+                    If string must be one of 'pca','nnmf', ica','fa',
+                    'dictionary', 'kernelpca'.
         kwargs: Additional keyword arguments to pass onto the scikit-learn
                 clustering object.
 
@@ -201,8 +202,9 @@ def load_class(import_string):
         'pca': 'sklearn.decomposition.PCA',
         'ica': 'sklearn.decomposition.FastICA',
         'nnmf': 'sklearn.decomposition.NMF',
-        'fa': 'sklearn.decomposition.FactorAnalysis'
-        }
+        'fa': 'sklearn.decomposition.FactorAnalysis',
+        'dictionary': 'sklearn.decomposition.DictionaryLearning',
+        'kernelpca': 'sklearn.decomposition.KernelPCA'}
 
     if algorithm in algs.keys():
         alg = load_class(algs[algorithm])
@@ -287,17 +289,35 @@ def check_square_numpy_matrix(data):
     return data
 
 
-def check_brain_data(data):
+def check_brain_data(data, mask=None):
     '''Check if data is a Brain_Data Instance.'''
     from nltools.data import Brain_Data
 
     if not isinstance(data, Brain_Data):
         if isinstance(data, nib.Nifti1Image):
-            data = Brain_Data(data)
+            data = Brain_Data(data, mask=mask)
         else:
             raise ValueError("Make sure data is a Brain_Data instance.")
+    else:
+        if mask is not None:
+            data = data.apply_mask(mask)
     return data
 
+def check_brain_data_is_single(data):
+    '''Logical test if Brain_Data instance is a single image
+    
+    Args:
+        data: brain data
+    
+    Returns:
+        (bool)
+    
+    '''
+    data = check_brain_data(data)
+    if len(data.shape()) > 1:
+        return False
+    else:
+        return True
 
 def _roi_func(brain, roi, algorithm, cv_dict, **kwargs):
     '''Brain_Data.predict_multi() helper function'''
@@ -322,7 +342,7 @@ def generate_jitter(n_trials, mean_time=5, min_time=2, max_time=12, atol=.2):
 
     Returns:
         data: (np.array) jitter for each trial
-
+ 
     '''
 
     def generate_data(n_trials, scale=5, min_time=2, max_time=12):
diff --git a/nltools/version.py b/nltools/version.py
index d9e8e4ca..4152a263 100644
--- a/nltools/version.py
+++ b/nltools/version.py
@@ -1,4 +1,4 @@
 """Specifies current version of nltools to be used by setup.py and __init__.py
 """
 
-__version__ = '0.3.16'
+__version__ = '0.3.17'