Repository: sappelhoff/pyprep
Pull request #88 (2021-06-15 20:52:42)

Refactor & improved documentation of NoisyChannels
<!--Thanks for contributing-->
<!--If this is your first time, please make sure to read the contributing guideline-->
<!--https://github.com/sappelhoff/pyprep/blob/master/.github/CONTRIBUTING.md-->

# PR Description

This has been sitting half-done for a while now, but I've finally got this into a good draft state! This PR is a major cleanup of the code/docs for the non-RANSAC `NoisyChannels` methods, including:

- Splitting the code for each method into small, documented chunks
- Moving he code for 50 Hz bandpass filtering to its own separate private method, so that all noisy detection methods can be run independently/repeatedly without accidental re-filtering of the data, and that we ensure the data is always filtered pre-RANSAC.
- Rewriting the docstrings for the class and all the methods to better explain what they all do, now that I've mostly wrapped my head around the logic
- Implemented a simple `madmed` function in utils to replace `robust.mad` for a mild speed boost, which also removes the entire dependency of `statsmodels`
- Refactoring `find_bad_by_correlation` to avoid division-by-zero and the related warning suppression code
- Identifying and documenting some magic numbers

Some other stuff I noticed but didn't clean up yet:
- The bad window fraction threshold is `frac_bad` for the bad-by-correlation method, but is `fraction_bad` for the bad-by-RANSAC method. Should these be made consistent?
- Technically bad-by-deviation, bad-by-correlation, and bad-by-dropout channels are all supposed to be detected and excluded prior to RANSAC, but right now the code doesn't enforce this. Should we add a chunk of code to ensure those have all been run first?
- For RANSAC, the name of the argument specifying the proportion of channels to use per sample is `sample_prop` in the docstring, but `fraction_good` in the function definition. Since `fraction_good` isn't technically accurate anymore, I'm assuming I should switch to `sample_prop` throughout?

This doesn't really move us any closer to 0.4.0, but I find code cleanup like this sort of relaxing so it was nice to pick away at between other work. Also, please let me know if any of the new code or docs is unclear!

# Merge Checklist

<!--To merge your PR we need to first take the following points into account.-->
<!--Please just leave this checklist untouched-->

- [x] the PR has been reviewed and all comments are resolved
- [x] all [CI][what-is-ci] checks pass
- [ ] (if applicable): the PR description includes the phrase `closes #<issue-number>` to [automatically close an issue][auto-close-documentation]
- [ ] (if applicable): bug fixes, new features, or [API][what-is-api] changes are documented in [whats_new.rst][whats-new-file]


[what-is-ci]: https://help.github.com/en/actions/building-and-testing-code-with-continuous-integration/about-continuous-integration
[auto-close-documentation]: https://help.github.com/en/github/managing-your-work-on-github/linking-a-pull-request-to-an-issue#linking-a-pull-request-to-an-issue-using-a-keyword
[what-is-api]: https://en.wikipedia.org/wiki/Application_programming_interface
[whats-new-file]: https://github.com/sappelhoff/pyprep/blob/master/docs/whats_new.rst


codecov-commenter: # [Codecov](https://codecov.io/gh/sappelhoff/pyprep/pull/88?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) Report
> Merging [#88](https://codecov.io/gh/sappelhoff/pyprep/pull/88?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) (3101557) into [master](https://codecov.io/gh/sappelhoff/pyprep/commit/fb4958fa0e6e0b04d7f2bf039725cb5d954d84b0?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) (fb4958f) will **decrease** coverage by `0.00%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/sappelhoff/pyprep/pull/88/graphs/tree.svg?width=650&height=150&src=pr&token=WniFfAHs8z&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff)](https://codecov.io/gh/sappelhoff/pyprep/pull/88?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff)

```diff
@@            Coverage Diff             @@
##           master      #88      +/-   ##
==========================================
- Coverage   98.39%   98.39%   -0.01%     
==========================================
  Files           7        7              
  Lines         686      685       -1     
==========================================
- Hits          675      674       -1     
  Misses         11       11              
```


| [Impacted Files](https://codecov.io/gh/sappelhoff/pyprep/pull/88?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) | Coverage Δ | |
|---|---|---|
| [pyprep/find\_noisy\_channels.py](https://codecov.io/gh/sappelhoff/pyprep/pull/88/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff#diff-cHlwcmVwL2ZpbmRfbm9pc3lfY2hhbm5lbHMucHk=) | `100.00% <100.00%> (ø)` | |
| [pyprep/utils.py](https://codecov.io/gh/sappelhoff/pyprep/pull/88/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff#diff-cHlwcmVwL3V0aWxzLnB5) | `99.29% <100.00%> (+0.02%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/sappelhoff/pyprep/pull/88?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/sappelhoff/pyprep/pull/88?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff). Last update [fb4958f...3101557](https://codecov.io/gh/sappelhoff/pyprep/pull/88?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff).

yjmantilla: 
These are only my opinions:

>The bad window fraction threshold is frac_bad for the bad-by-correlation method, but is fraction_bad for the bad-by-RANSAC method. Should these be made consistent?

I say it would be nice since it is both explicit and consistent.

>Technically bad-by-deviation, bad-by-correlation, and bad-by-dropout channels are all supposed to be detected and excluded prior to RANSAC, but right now the code doesn't enforce this. Should we add a chunk of code to ensure those have all been run first?

I think RANSAC as a module should be standalone like it is in autoreject. If anything, include the clarification on the docstring. This reminds me of the ICA analysis in mne: low frequency drifts should be removed but that is said as a warning, not as an absolute requirement inside the code.

![image](https://user-images.githubusercontent.com/36543115/122128209-f66fb980-cdf9-11eb-9b77-80586da025a2.png)

>For RANSAC, the name of the argument specifying the proportion of channels to use per sample is sample_prop in the docstring, but fraction_good in the function definition. Since fraction_good isn't technically accurate anymore, I'm assuming I should switch to sample_prop throughout?

As before, that would be nice indeed. Like being consistent should be something we look for.

Patch:
diff --git a/pyprep/find_noisy_channels.py b/pyprep/find_noisy_channels.py
index 6b289b4..4f2dc9d 100644
--- a/pyprep/find_noisy_channels.py
+++ b/pyprep/find_noisy_channels.py
@@ -3,18 +3,42 @@ import mne
 import numpy as np
 from mne.utils import check_random_state
 from scipy import signal
-from statsmodels import robust
 
 from pyprep.ransac import find_bad_by_ransac
 from pyprep.removeTrend import removeTrend
-from pyprep.utils import _mat_iqr, _mat_quantile, filter_design
+from pyprep.utils import _filter_design, _mad, _mat_iqr, _mat_quantile
 
 
 class NoisyChannels:
-    """Implements the functionality of the `findNoisyChannels` function.
+    """Detect bad channels in an EEG recording using a range of methods.
 
-    It is a part of the PREP (preprocessing pipeline) for EEG data recorded using
-    10-20 montage style described in [1]_.
+    This class provides a number of methods for detecting bad channels across a
+    full-session EEG recording. Specifically, this class implements all of the
+    noisy channel detection methods used in the PREP pipeline, as described in [1]_.
+    The detection methods in this class can be run independently, or can be run
+    all at once using the :meth:`~.find_all_bads` method.
+
+    At present, only EEG channels are supported and any non-EEG channels in the
+    provided data will be ignored.
+
+    Parameters
+    ----------
+    raw : mne.io.Raw
+        An MNE Raw object to check for bad EEG channels.
+    do_detrend : bool, optional
+        Whether or not low-frequency (<1.0 Hz) trends should be removed from the
+        EEG signal prior to bad channel detection. This should always be set to
+        ``True`` unless the signal has already had low-frequency trends removed.
+        Defaults to ``True``.
+    random_state : {int, None, np.random.RandomState}, optional
+        The seed to use for random number generation within RANSAC. This can be
+        ``None``, an integer, or a :class:`~numpy.random.RandomState` object.
+        If ``None``, a random seed will be obtained from the operating system.
+        Defaults to ``None``.
+    matlab_strict : bool, optional
+        Whether or not PyPREP should strictly follow MATLAB PREP's internal
+        math, ignoring any improvements made in PyPREP over the original code
+        (see :ref:`matlab-diffs` for more details). Defaults to ``False``.
 
     References
     ----------
@@ -25,26 +49,6 @@ class NoisyChannels:
     """
 
     def __init__(self, raw, do_detrend=True, random_state=None, matlab_strict=False):
-        """Initialize the class.
-
-        Parameters
-        ----------
-        raw : mne.io.Raw
-            The MNE raw object.
-        do_detrend : bool, optional
-            Whether or not to remove a trend from the data upon initializing the
-            `NoisyChannels` object. Defaults to ``True``.
-        random_state : {int, None, np.random.RandomState}, optional
-            The random seed at which to initialize the class. If random_state
-            is an int, it will be used as a seed for RandomState.
-            If ``None``, the seed will be obtained from the operating system
-            (see RandomState for details). Default is ``None``.
-        matlab_strict : bool, optional
-            Whether or not PyPREP should strictly follow MATLAB PREP's internal
-            math, ignoring any improvements made in PyPREP over the original code
-            (see :ref:`matlab-diffs` for more details). Defaults to ``False``.
-
-        """
         # Make sure that we got an MNE object
         assert isinstance(raw, mne.io.BaseRaw)
 
@@ -92,65 +96,73 @@ class NoisyChannels:
         # Make a subset of the data containing only usable EEG channels
         self.usable_idx = np.isin(ch_names, bads_by_nan_flat, invert=True)
         self.EEGData = self.raw_mne.get_data(picks=ch_names[self.usable_idx])
-        self.EEGData_beforeFilt = self.EEGData
+        self.EEGFiltered = None
 
         # Get usable EEG channel names & channel counts
         self.ch_names_new = np.asarray(ch_names[self.usable_idx])
         self.n_chans_new = len(self.ch_names_new)
 
+    def _get_filtered_data(self):
+        """Apply a [1 Hz - 50 Hz] bandpass filter to the EEG signal.
+
+        Only applied if the sample rate is above 100 Hz to avoid violating the
+        Nyquist theorem.
+
+        """
+        if self.sample_rate <= 100:
+            return self.EEGData.copy()
+
+        bandpass_filter = _filter_design(
+            N_order=100,
+            amp=np.array([1, 1, 0, 0]),
+            freq=np.array([0, 90 / self.sample_rate, 100 / self.sample_rate, 1]),
+        )
+        EEG_filt = np.zeros_like(self.EEGData)
+        for i in range(EEG_filt.shape[0]):
+            EEG_filt[i, :] = signal.filtfilt(bandpass_filter, 1, self.EEGData[i, :])
+
+        return EEG_filt
+
     def get_bads(self, verbose=False):
-        """Get a list of all bad channels.
+        """Get a list of all channels currently flagged as bad.
 
-        This function makes a list of all the bad channels and prints them if verbose
-        is True.
+        Note that this method does not perform any bad channel detection itself,
+        and only reports channels already detected as bad by other methods.
 
         Parameters
         ----------
         verbose : bool
-            If verbose, print a summary of bad channels.
+            If ``True``, a summary of the channels currently flagged as by bad per
+            category is printed. Defaults to ``False``.
+
+        Returns
+        -------
+        bads : list
+            THe names of all bad channels detected by any method so far.
 
         """
-        bads = (
-            self.bad_by_nan
-            + self.bad_by_flat
-            + self.bad_by_deviation
-            + self.bad_by_hf_noise
-            + self.bad_by_SNR
-            + self.bad_by_correlation
-            + self.bad_by_dropout
-            + self.bad_by_ransac
-        )
-        bads = list(set(bads))
+        bads = {
+            "n/a": self.bad_by_nan,
+            "flat": self.bad_by_flat,
+            "deviation": self.bad_by_deviation,
+            "hf noise": self.bad_by_hf_noise,
+            "correl": self.bad_by_correlation,
+            "SNR": self.bad_by_SNR,
+            "dropout": self.bad_by_dropout,
+            "RANSAC": self.bad_by_ransac,
+        }
+
+        all_bads = set()
+        for bad_chs in bads.values():
+            all_bads.update(bad_chs)
 
         if verbose:
-            print("Found {} uniquely bad channels.".format(len(bads)))
-            print("\n{} by n/a: {}".format(len(self.bad_by_nan), self.bad_by_nan))
-            print("\n{} by flat: {}".format(len(self.bad_by_flat), self.bad_by_flat))
-            print(
-                "\n{} by deviation: {}".format(
-                    len(self.bad_by_deviation), self.bad_by_deviation
-                )
-            )
-            print(
-                "\n{} by hf noise: {}".format(
-                    len(self.bad_by_hf_noise), self.bad_by_hf_noise
-                )
-            )
-            print(
-                "\n{} by correl: {}".format(
-                    len(self.bad_by_correlation), self.bad_by_correlation
-                )
-            )
-            print("\n{} by SNR {}".format(len(self.bad_by_SNR), self.bad_by_SNR))
-            print(
-                "\n{} by dropout: {}".format(
-                    len(self.bad_by_dropout), self.bad_by_dropout
-                )
-            )
-            print(
-                "\n{} by ransac: {}".format(len(self.bad_by_ransac), self.bad_by_ransac)
-            )
-        return bads
+            out = f"Found {len(all_bads)} uniquely bad channels:\n"
+            for bad_type, bad_chs in bads.items():
+                out += f"\n{len(bad_chs)} by {bad_type}: {bad_chs}\n"
+            print(out)
+
+        return list(all_bads)
 
     def find_all_bads(self, ransac=True, channel_wise=False, max_chunk_size=None):
         """Call all the functions to detect bad channels.
@@ -183,8 +195,10 @@ class NoisyChannels:
             effect. Defaults to ``None``.
 
         """
-        self.find_bad_by_nan_flat()
+        # NOTE: Bad-by-NaN/flat is already run during init, no need to re-run here
         self.find_bad_by_deviation()
+        self.find_bad_by_hfnoise()
+        self.find_bad_by_correlation()
         self.find_bad_by_SNR()
         if ransac:
             self.find_bad_by_ransac(
@@ -192,7 +206,16 @@ class NoisyChannels:
             )
 
     def find_bad_by_nan_flat(self):
-        """Detect channels that appear flat or have NaN values."""
+        """Detect channels than contain NaN values or have near-flat signals.
+
+        A channel is considered flat if its standard deviation or its median
+        absolute deviation from the median (MAD) are below 1e-9 microvolts.
+
+        This method is run automatically when a ``NoisyChannels`` object is
+        initialized, preventing flat or NaN-containing channels from interfering
+        with the detection of other types of bad channels.
+
+        """
         # Get all EEG channels from original copy of data
         EEGData = self.raw_mne.get_data()
 
@@ -202,7 +225,7 @@ class NoisyChannels:
 
         # Detect channels with flat or extremely weak signals
         FLAT_THRESHOLD = 1e-15  # corresponds to 10e-10 µV in MATLAB PREP
-        flat_by_mad = robust.mad(EEGData, axis=1, c=1) < FLAT_THRESHOLD
+        flat_by_mad = _mad(EEGData, axis=1) < FLAT_THRESHOLD
         flat_by_stdev = np.std(EEGData, axis=1) < FLAT_THRESHOLD
         flat_channel_mask = flat_by_mad | flat_by_stdev
         flat_channels = self.ch_names_original[flat_channel_mask]
@@ -212,234 +235,272 @@ class NoisyChannels:
         self.bad_by_flat = flat_channels.tolist()
 
     def find_bad_by_deviation(self, deviation_threshold=5.0):
-        """Robust z-score of the robust standard deviation for each channel is calculated.
+        """Detect channels with abnormally high or low overall amplitudes.
 
-        Channels having a z-score greater than 5 are detected as bad.
+        A channel is considered "bad-by-deviation" if its amplitude deviates
+        considerably from the median channel amplitude, as calculated using a
+        robust Z-scoring method and the given deviation threshold.
+
+        Amplitude Z-scores are calculated using the formula
+        ``(channel_amplitude - median_amplitude) / amplitude_sd``, where
+        channel amplitudes are calculated using a robust outlier-resistant estimate
+        of the signals' standard deviations (IQR scaled to units of SD), and the
+        amplitude SD is the IQR-based SD of those amplitudes.
 
         Parameters
         ----------
-        deviation_threshold : float
-            z-score threshold above which channels will be labelled bad.
+        deviation_threshold : float, optional
+            The minimum absolute z-score of a channel for it to be considered
+            bad-by-deviation. Defaults to ``5.0``.
 
         """
-        deviation_channel_mask = [False] * (self.n_chans_original)
-        robust_channel_deviation = np.zeros(self.n_chans_original, dtype=np.float64)
-        channel_deviation = np.zeros(self.n_chans_new)
-        for i in range(0, self.n_chans_new):
-            channel_deviation[i] = 0.7413 * _mat_iqr(self.EEGData[i, :])
-        channel_deviationSD = 0.7413 * _mat_iqr(channel_deviation)
-        channel_deviationMedian = np.nanmedian(channel_deviation)
-        robust_channel_deviation[self.usable_idx] = np.divide(
-            np.subtract(channel_deviation, channel_deviationMedian), channel_deviationSD
-        )
-        for i in range(0, self.n_chans_original):
-            deviation_channel_mask[i] = abs(
-                robust_channel_deviation[i]
-            ) > deviation_threshold or np.isnan(robust_channel_deviation[i])
+        IQR_TO_SD = 0.7413  # Scales units of IQR to units of SD, assuming normality
+        # Reference: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/IQR.html
+
+        # Get channel amplitudes and the median / robust SD of those amplitudes
+        chan_amplitudes = _mat_iqr(self.EEGData, axis=1) * IQR_TO_SD
+        amp_sd = _mat_iqr(chan_amplitudes) * IQR_TO_SD
+        amp_median = np.nanmedian(chan_amplitudes)
+
+        # Calculate robust Z-scores for the channel amplitudes
+        amplitude_zscore = np.zeros(self.n_chans_original)
+        amplitude_zscore[self.usable_idx] = (chan_amplitudes - amp_median) / amp_sd
+
+        # Flag channels with amplitudes that deviate excessively from the median
+        abnormal_amplitude = np.abs(amplitude_zscore) > deviation_threshold
+        deviation_channel_mask = np.isnan(amplitude_zscore) | abnormal_amplitude
+
+        # Update names of bad channels by excessive deviation & save additional info
         deviation_channels = self.ch_names_original[deviation_channel_mask]
         self.bad_by_deviation = deviation_channels.tolist()
         self._extra_info["bad_by_deviation"].update(
             {
-                "median_channel_deviation": channel_deviationMedian,
-                "channel_deviation_sd": channel_deviationSD,
-                "robust_channel_deviations": robust_channel_deviation,
+                "median_channel_amplitude": amp_median,
+                "channel_amplitude_sd": amp_sd,
+                "robust_channel_deviations": amplitude_zscore,
             }
         )
 
     def find_bad_by_hfnoise(self, HF_zscore_threshold=5.0):
-        """Determine noise of channel through high frequency ratio.
+        """Detect channels with abnormally high amounts of high-frequency noise.
 
-        Noisiness of the channel is determined by finding the ratio of the median
-        absolute deviation of high frequency to low frequency components.
+        The noisiness of a channel is defined as the amplitude of its
+        high-frequency (>50 Hz) components divided by its overall amplitude.
+        A channel is considered "bad-by-high-frequency-noise" if its noisiness
+        is considerably higher than the median channel noisiness, as determined
+        by a robust Z-scoring method and the given Z-score threshold.
 
-        Low pass 50 Hz filter is used to separate the frequency components.
-        A robust z-score is then calculated relative to all the channels.
+        Due to the Nyquist theorem, this method will only attempt bad channel
+        detection if the sample rate of the given signal is above 100 Hz.
 
         Parameters
         ----------
-        HF_zscore_threshold : float
-            z-score threshold above which channels would be labelled as bad.
+        HF_zscore_threshold : float, optional
+            The minimum noisiness z-score of a channel for it to be considered
+            bad-by-high-frequency-noise. Defaults to ``5.0``.
 
         """
-        data_tmp = np.transpose(self.EEGData)
-        dimension = np.shape(data_tmp)
-        zscore_HFNoise = np.zeros(self.n_chans_original, dtype=np.float64)
+        MAD_TO_SD = 1.4826  # Scales units of MAD to units of SD, assuming normality
+        # Reference: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/mad.html
+
+        if self.EEGFiltered is None:
+            self.EEGFiltered = self._get_filtered_data()
+
+        # Set default values for noise parameters
+        noise_median, noise_sd = (0, 1)
+        noise_zscore = np.zeros(self.n_chans_original)
+
+        # If sample rate is high enough, calculate ratio of > 50 Hz amplitude to
+        # < 50 Hz amplitude for each channel and get robust z-scores of values
         if self.sample_rate > 100:
-            EEG_filt = np.zeros((dimension[0], dimension[1]))
-            bandpass_filter = filter_design(
-                N_order=100,
-                amp=np.array([1, 1, 0, 0]),
-                freq=np.array([0, 90 / self.sample_rate, 100 / self.sample_rate, 1]),
-            )
-            for i in range(0, dimension[1]):
-                EEG_filt[:, i] = signal.filtfilt(bandpass_filter, 1, data_tmp[:, i])
             noisiness = np.divide(
-                robust.mad(np.subtract(data_tmp, EEG_filt), c=1),
-                robust.mad(EEG_filt, c=1),
-            )
-            noisiness_median = np.nanmedian(noisiness)
-            noiseSD = (
-                np.median(np.absolute(np.subtract(noisiness, np.median(noisiness))))
-                * 1.4826
-            )
-            zscore_HFNoise[self.usable_idx] = np.divide(
-                np.subtract(noisiness, noisiness_median), noiseSD
+                _mad(self.EEGData - self.EEGFiltered, axis=1),
+                _mad(self.EEGFiltered, axis=1),
             )
-            hf_mask = np.isnan(zscore_HFNoise) | (zscore_HFNoise > HF_zscore_threshold)
-            HFNoise_channels = self.ch_names_original[hf_mask]
-        else:
-            EEG_filt = data_tmp
-            noisiness_median = 0
-            noiseSD = 1
-            HFNoise_channels = np.asarray([])
-        self.EEGData_beforeFilt = data_tmp
-        self.EEGData = np.transpose(EEG_filt)
-        self.bad_by_hf_noise = HFNoise_channels.tolist()
+            noise_median = np.nanmedian(noisiness)
+            noise_sd = np.median(np.abs(noisiness - noise_median)) * MAD_TO_SD
+            noise_zscore[self.usable_idx] = (noisiness - noise_median) / noise_sd
+
+        # Flag channels with much more high-frequency noise than the median channel
+        hf_mask = np.isnan(noise_zscore) | (noise_zscore > HF_zscore_threshold)
+        hf_noise_channels = self.ch_names_original[hf_mask]
+
+        # Update names of high-frequency noise channels & save additional info
+        self.bad_by_hf_noise = hf_noise_channels.tolist()
         self._extra_info["bad_by_hf_noise"].update(
             {
-                "median_channel_noisiness": noisiness_median,
-                "channel_noisiness_sd": noiseSD,
-                "hf_noise_zscores": zscore_HFNoise,
+                "median_channel_noisiness": noise_median,
+                "channel_noisiness_sd": noise_sd,
+                "hf_noise_zscores": noise_zscore,
             }
         )
 
     def find_bad_by_correlation(
         self, correlation_secs=1.0, correlation_threshold=0.4, frac_bad=0.01
     ):
-        """Find correlation between the low frequency components of the EEG below 50 Hz.
-
-        Correlation is done using a sliding non-overlapping time window.
-        The maximum absolute correlation is as the 98th percentile of the absolute
-        values of the correlations with the other channels
-        If the maximum correlation is less than 0.4 then the channel is designated as
-        bad by correlation.
+        """Detect channels that sometimes don't correlate with any other channels.
+
+        Channel correlations are calculated by splitting the recording into
+        non-overlapping windows of time (default: 1 second), getting the absolute
+        correlations of each usable channel with every other usable channel for
+        each window, and then finding the highest correlation each channel has
+        with another channel for each window (by taking the 98th percentile of
+        the absolute correlations).
+
+        A correlation window is considered "bad" for a channel if its maximum
+        correlation with another channel is below the provided correlation
+        threshold (default: ``0.4``). A channel is considered "bad-by-correlation"
+        if its fraction of bad correlation windows is above the bad fraction
+        threshold (default: ``0.01``).
+
+        This method also detects channels with intermittent dropouts (i.e.,
+        regions of flat signal). A channel is considered "bad-by-dropout" if
+        its fraction of correlation windows with a completely flat signal is
+        above the bad fraction threshold (default: ``0.01``).
 
         Parameters
         ----------
-        correlation_secs : float
-            length of the correlation time window (default: 1 secs).
-        correlation_threshold : float
-            correlation threshold below which channel is marked bad.
-        frac_bad : float
-            percentage of data windows in which the correlation threshold was
-            not surpassed and if a channel gets a value of greater than 1%, it
-            is designated bad. Notice that if `correlation_secs` is high, and
-            thus the number of windows is low, the default `frac_bad` may
-            be too extreme. For example if only 60 windows are available
-            for a dataset, only a single bad window would be needed to
-            classify as bad.
+        correlation_secs : float, optional
+            The length (in seconds) of each correlation window. Defaults to ``1.0``.
+        correlation_threshold : float, optional
+            The lowest maximum inter-channel correlation for a channel to be
+            considered "bad" within a given window. Defaults to ``0.4``.
+        frac_bad : float, optional
+            The minimum proportion of bad windows for a channel to be considered
+            "bad-by-correlation" or "bad-by-dropout". Defaults to ``0.01`` (1% of
+            all windows).
 
         """
-        self.find_bad_by_hfnoise()  # since filtering is performed there
-        usable_idx = self.usable_idx
-        correlation_frames = correlation_secs * self.sample_rate
-        correlation_window = np.arange(correlation_frames)
-        correlation_offsets = np.arange(
-            1, (self.n_samples - correlation_frames), correlation_frames
-        )
-        w_correlation = len(correlation_offsets)
-        channel_correlation = np.ones((w_correlation, self.n_chans_original))
-        noiselevels = np.zeros((w_correlation, self.n_chans_original))
-        channel_deviations = np.zeros((w_correlation, self.n_chans_original))
-        drop = np.zeros((w_correlation, self.n_chans_original))
-        len_correlation_window = len(correlation_window)
-        EEGData = np.transpose(self.EEGData)
-        EEG_new_win = np.reshape(
-            np.transpose(EEGData[0 : len_correlation_window * w_correlation, :]),
-            (self.n_chans_new, len_correlation_window, w_correlation),
-            order="F",
-        )
-        data_win = np.reshape(
-            np.transpose(
-                self.EEGData_beforeFilt[0 : len_correlation_window * w_correlation, :]
-            ),
-            (self.n_chans_new, len_correlation_window, w_correlation),
-            order="F",
-        )
-        for k in range(0, w_correlation):
-            eeg_portion = np.transpose(np.squeeze(EEG_new_win[:, :, k]))
-            data_portion = np.transpose(np.squeeze(data_win[:, :, k]))
-            with np.errstate(invalid="ignore"):  # suppress divide-by-zero warnings
-                window_correlation = np.corrcoef(np.transpose(eeg_portion))
-            abs_corr = np.abs(
-                np.subtract(window_correlation, np.diag(np.diag(window_correlation)))
-            )
-            channel_correlation[k, usable_idx] = _mat_quantile(abs_corr, 0.98, axis=0)
-            with np.errstate(invalid="ignore"):  # suppress divide-by-zero warnings
-                noiselevels[k, usable_idx] = np.divide(
-                    robust.mad(np.subtract(data_portion, eeg_portion), c=1),
-                    robust.mad(eeg_portion, c=1),
-                )
-            channel_deviations[k, usable_idx] = 0.7413 * _mat_iqr(data_portion, axis=0)
-        for i in range(0, w_correlation):
-            for j in range(0, self.n_chans_original):
-                drop[i, j] = np.int(
-                    np.isnan(channel_correlation[i, j]) or np.isnan(noiselevels[i, j])
-                )
-                if drop[i, j] == 1:
-                    channel_correlation[i, j] = 0
-                    noiselevels[i, j] = 0
+        IQR_TO_SD = 0.7413  # Scales units of IQR to units of SD, assuming normality
+        # Reference: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/IQR.html
+
+        if self.EEGFiltered is None:
+            self.EEGFiltered = self._get_filtered_data()
+
+        # Determine the number and size (in frames) of correlation windows
+        win_size = int(correlation_secs * self.sample_rate)
+        win_offsets = np.arange(1, (self.n_samples - win_size), win_size)
+        win_count = len(win_offsets)
+
+        # Initialize per-window arrays for each type of noise info calculated below
+        max_correlations = np.ones((win_count, self.n_chans_original))
+        dropout = np.zeros((win_count, self.n_chans_original), dtype=bool)
+        noiselevels = np.zeros((win_count, self.n_chans_original))
+        channel_amplitudes = np.zeros((win_count, self.n_chans_original))
+
+        for w in range(win_count):
+            # Get both filtered and unfiltered data for the current window
+            start, end = (w * win_size, (w + 1) * win_size)
+            eeg_filtered = self.EEGFiltered[:, start:end]
+            eeg_raw = self.EEGData[:, start:end]
+
+            # Get channel amplitude info for the window
+            usable = self.usable_idx.copy()
+            channel_amplitudes[w, usable] = _mat_iqr(eeg_raw, axis=1) * IQR_TO_SD
+
+            # Check for any channel dropouts (flat signal) within the window
+            eeg_amplitude = _mad(eeg_filtered, axis=1)
+            dropout[w, usable] = eeg_amplitude == 0
+
+            # Exclude any dropout chans from further calculations (avoids div-by-zero)
+            usable[usable] = eeg_amplitude > 0
+            eeg_raw = eeg_raw[eeg_amplitude > 0, :]
+            eeg_filtered = eeg_filtered[eeg_amplitude > 0, :]
+            eeg_amplitude = eeg_amplitude[eeg_amplitude > 0]
+
+            # Get high-frequency noise ratios for the window
+            high_freq_amplitude = _mad(eeg_raw - eeg_filtered, axis=1)
+            noiselevels[w, usable] = high_freq_amplitude / eeg_amplitude
+
+            # Get inter-channel correlations for the window
+            win_correlations = np.corrcoef(eeg_filtered)
+            abs_corr = np.abs(win_correlations - np.diag(np.diag(win_correlations)))
+            max_correlations[w, usable] = _mat_quantile(abs_corr, 0.98, axis=0)
+            max_correlations[w, dropout[w, :]] = 0  # Set dropout correlations to 0
 
         # Flag channels with above-threshold fractions of bad correlation windows
-        maximum_correlations = np.transpose(channel_correlation)
-        thresholded_correlations = maximum_correlations < correlation_threshold
-        thresholded_correlations = thresholded_correlations.astype(int)
-        fraction_BadCorrelationWindows = np.mean(thresholded_correlations, axis=1)
-        bad_correlation_mask = fraction_BadCorrelationWindows > frac_bad
+        thresholded_correlations = max_correlations < correlation_threshold
+        fraction_bad_corr_windows = np.mean(thresholded_correlations, axis=0)
+        bad_correlation_mask = fraction_bad_corr_windows > frac_bad
         bad_correlation_channels = self.ch_names_original[bad_correlation_mask]
-        self.bad_by_correlation = bad_correlation_channels.tolist()
 
         # Flag channels with above-threshold fractions of drop-out windows
-        drop_out = np.transpose(drop)
-        fraction_BadDropOutWindows = np.mean(drop_out, axis=1)
-        dropout_mask = fraction_BadDropOutWindows > frac_bad
+        fraction_dropout_windows = np.mean(dropout, axis=0)
+        dropout_mask = fraction_dropout_windows > frac_bad
         dropout_channels = self.ch_names_original[dropout_mask]
-        self.bad_by_dropout = dropout_channels.tolist()
 
+        # Update names of low-correlation/dropout channels & save additional info
+        self.bad_by_correlation = bad_correlation_channels.tolist()
+        self.bad_by_dropout = dropout_channels.tolist()
         self._extra_info["bad_by_correlation"] = {
-            "max_correlations": maximum_correlations,
-            "median_max_correlations": np.median(maximum_correlations, axis=1),
-            "bad_window_fractions": fraction_BadCorrelationWindows,
+            "max_correlations": np.transpose(max_correlations),
+            "median_max_correlations": np.median(max_correlations, axis=0),
+            "bad_window_fractions": fraction_bad_corr_windows,
         }
         self._extra_info["bad_by_dropout"] = {
-            "dropouts": drop_out,
-            "bad_window_fractions": fraction_BadDropOutWindows,
+            "dropouts": np.transpose(dropout.astype(np.int8)),
+            "bad_window_fractions": fraction_dropout_windows,
         }
-        self._extra_info["bad_by_deviation"]["channel_deviations"] = channel_deviations
+        self._extra_info["bad_by_deviation"]["channel_amplitudes"] = channel_amplitudes
         self._extra_info["bad_by_hf_noise"]["noise_levels"] = noiselevels
 
     def find_bad_by_SNR(self):
-        """Determine the channels that fail both by correlation and HF noise."""
-        self.find_bad_by_correlation()
-        set_hf = set(self.bad_by_hf_noise)
-        set_correlation = set(self.bad_by_correlation)
-        self.bad_by_SNR = list(set_correlation.intersection(set_hf))
-        return None
+        """Detect channels that have a low signal-to-noise ratio.
+
+        Channels are considered "bad-by-SNR" if they are bad by both high-frequency
+        noise and bad by low correlation.
+
+        """
+        # Get names of bad-by-HF-noise and bad-by-correlation channels
+        if not len(self._extra_info["bad_by_hf_noise"]) > 1:
+            self.find_bad_by_hfnoise()
+        if not len(self._extra_info["bad_by_correlation"]):
+            self.find_bad_by_correlation()
+        bad_by_hf = set(self.bad_by_hf_noise)
+        bad_by_corr = set(self.bad_by_correlation)
+
+        # Flag channels bad by both HF noise and low correlation as bad by low SNR
+        self.bad_by_SNR = list(bad_by_corr.intersection(bad_by_hf))
 
     def find_bad_by_ransac(
         self,
         n_samples=50,
-        fraction_good=0.25,
+        sample_prop=0.25,
         corr_thresh=0.75,
-        fraction_bad=0.4,
+        frac_bad=0.4,
         corr_window_secs=5.0,
         channel_wise=False,
         max_chunk_size=None,
     ):
-        """Detect channels that are not predicted well by other channels.
-
-        This method is a wrapper for the :func:`ransac.find_bad_by_ransac`
+        """Detect channels that are predicted poorly by other channels.
+
+        This method uses a random sample consensus approach (RANSAC, see [1]_,
+        and a short discussion in [2]_) to try and predict what the signal should
+        be for each channel based on the signals and spatial locations of other
+        currently-good channels. RANSAC correlations are calculated by splitting
+        the recording into non-overlapping windows of time (default: 5 seconds)
+        and correlating each channel's RANSAC-predicted signal with its actual
+        signal within each window.
+
+        A RANSAC window is considered "bad" for a channel if its predicted signal
+        vs. actual signal correlation falls below the given correlation threshold
+        (default: ``0.75``). A channel is considered "bad-by-RANSAC" if its fraction
+        of bad RANSAC windows is above the given threshold (default: ``0.4``).
+
+        Due to its random sampling component, the channels identified as
+        "bad-by-RANSAC" may vary slightly between calls of this method.
+        Additionally, bad channels may vary between different montages given that
+        RANSAC's signal predictions are based on the spatial coordinates of each
+        electrode.
+
+        This method is a wrapper for the :func:`~ransac.find_bad_by_ransac`
         function.
 
-        Here, a ransac approach (see [1]_, and a short discussion in [2]_) is
-        adopted to predict a "clean EEG" dataset. After identifying clean EEG
-        channels through the other methods, the clean EEG dataset is
-        constructed by repeatedly sampling a small subset of clean EEG channels
-        and interpolation the complete data. The median of all those
-        repetitions forms the clean EEG dataset. In a second step, the original
-        and the ransac predicted data are correlated and channels, which do not
-        correlate well with themselves across the two datasets are considered
-        `bad_by_ransac`.
+        .. warning:: For optimal performance, RANSAC requires that channels bad by
+                     deviation, correlation, and/or dropout have already been
+                     flagged. Otherwise RANSAC will attempt to use those channels
+                     when making signal predictions, decreasing accuracy and thus
+                     increasing the likelihood of false positives.
 
         Parameters
         ----------
@@ -456,7 +517,7 @@ class NoisyChannels:
             The minimum predicted vs. actual signal correlation for a channel to
             be considered good within a given RANSAC window. Defaults
             to ``0.75``.
-        fraction_bad : float, optional
+        frac_bad : float, optional
             The minimum fraction of bad (i.e., below-threshold) RANSAC windows
             for a channel to be considered bad-by-RANSAC. Defaults to ``0.4``.
         corr_window_secs : float, optional
@@ -480,7 +541,7 @@ class NoisyChannels:
 
         References
         ----------
-        .. [1] Fischler, M.A., Bolles, R.C. (1981). Random rample consensus: A
+        .. [1] Fischler, M.A., Bolles, R.C. (1981). Random sample consensus: A
             Paradigm for Model Fitting with Applications to Image Analysis and
             Automated Cartography. Communications of the ACM, 24, 381-395
         .. [2] Jas, M., Engemann, D.A., Bekhti, Y., Raimondo, F., Gramfort, A.
@@ -488,19 +549,22 @@ class NoisyChannels:
             Data. NeuroImage, 159, 417-429
 
         """
+        if self.EEGFiltered is None:
+            self.EEGFiltered = self._get_filtered_data()
+
         exclude_from_ransac = (
             self.bad_by_correlation + self.bad_by_deviation + self.bad_by_dropout
         )
         self.bad_by_ransac, ch_correlations_usable = find_bad_by_ransac(
-            self.EEGData,
+            self.EEGFiltered,
             self.sample_rate,
             self.ch_names_new,
             self.raw_mne._get_channel_positions()[self.usable_idx, :],
             exclude_from_ransac,
             n_samples,
-            fraction_good,
+            sample_prop,
             corr_thresh,
-            fraction_bad,
+            frac_bad,
             corr_window_secs,
             channel_wise,
             max_chunk_size,
diff --git a/pyprep/ransac.py b/pyprep/ransac.py
index 615d258..1f53d3c 100644
--- a/pyprep/ransac.py
+++ b/pyprep/ransac.py
@@ -8,9 +8,9 @@ from pyprep.utils import (
     _correlate_arrays,
     _get_random_subset,
     _mat_round,
-    print_progress,
-    split_list,
-    verify_free_ram,
+    _print_progress,
+    _split_list,
+    _verify_free_ram,
 )
 
 
@@ -21,9 +21,9 @@ def find_bad_by_ransac(
     chn_pos,
     exclude,
     n_samples=50,
-    fraction_good=0.25,
+    sample_prop=0.25,
     corr_thresh=0.75,
-    fraction_bad=0.4,
+    frac_bad=0.4,
     corr_window_secs=5.0,
     channel_wise=False,
     max_chunk_size=None,
@@ -69,7 +69,7 @@ def find_bad_by_ransac(
     corr_thresh : float, optional
         The minimum predicted vs. actual signal correlation for a channel to
         be considered good within a given RANSAC window. Defaults to ``0.75``.
-    fraction_bad : float, optional
+    frac_bad : float, optional
         The minimum fraction of bad (i.e., below-threshold) RANSAC windows for a
         channel to be considered bad-by-RANSAC. Defaults to ``0.4``.
     corr_window_secs : float, optional
@@ -110,7 +110,7 @@ def find_bad_by_ransac(
 
     References
     ----------
-    .. [1] Fischler, M.A., Bolles, R.C. (1981). Random rample consensus: A
+    .. [1] Fischler, M.A., Bolles, R.C. (1981). Random sample consensus: A
         Paradigm for Model Fitting with Applications to Image Analysis and
         Automated Cartography. Communications of the ACM, 24, 381-395
     .. [2] Jas, M., Engemann, D.A., Bekhti, Y., Raimondo, F., Gramfort, A.
@@ -133,13 +133,13 @@ def find_bad_by_ransac(
     # Check if we have enough remaining channels
     # after exclusion of bad channels
     n_chans = data.shape[0]
-    n_pred_chns = int(np.around(fraction_good * n_chans))
+    n_pred_chns = int(np.around(sample_prop * n_chans))
 
     if n_pred_chns <= 3:
-        sample_pct = int(fraction_good * 100)
+        sample_pct = int(sample_prop * 100)
         e = "Too few channels in the original data to reliably perform RANSAC "
         e += "(minimum {0} for a sample size of {1}%)."
-        raise IOError(e.format(int(np.floor(4.0 / fraction_good)), sample_pct))
+        raise IOError(e.format(int(np.floor(4.0 / sample_prop)), sample_pct))
     elif n_chans_good < (n_pred_chns + 1):
         e = "Too many noisy channels in the data to reliably perform RANSAC "
         e += "(only {0} good channels remaining, need at least {1})."
@@ -148,10 +148,10 @@ def find_bad_by_ransac(
     # Before running, make sure we have enough memory when using the
     # smallest possible chunk size
     if channel_wise:
-        verify_free_ram(data, n_samples, 1)
+        _verify_free_ram(data, n_samples, 1)
     else:
         window_size = int(sample_rate * corr_window_secs)
-        verify_free_ram(data[:, :window_size], n_samples, n_chans_good)
+        _verify_free_ram(data[:, :window_size], n_samples, n_chans_good)
 
     # Generate random channel picks for each RANSAC sample
     random_ch_picks = []
@@ -207,7 +207,7 @@ def find_bad_by_ransac(
     # If not using window-wise RANSAC, do channel-wise RANSAC
     while mem_error and channel_wise:
         try:
-            channel_chunks = split_list(job, chunk_size)
+            channel_chunks = _split_list(job, chunk_size)
             total_chunks = len(channel_chunks)
             current = 1
             for chunk in channel_chunks:
@@ -247,7 +247,7 @@ def find_bad_by_ransac(
     frac_bad_corr_windows = np.mean(thresholded_correlations, axis=0)
 
     # find the corresponding channel names and return
-    bad_ransac_channels_idx = np.argwhere(frac_bad_corr_windows > fraction_bad)
+    bad_ransac_channels_idx = np.argwhere(frac_bad_corr_windows > frac_bad)
     bad_ransac_channels_name = complete_chn_labs[bad_ransac_channels_idx.astype(int)]
     bad_by_ransac = [i[0] for i in bad_ransac_channels_name]
     print("\nRANSAC done!")
@@ -332,7 +332,7 @@ def _ransac_by_window(data, interpolation_mats, win_size, win_count, matlab_stri
     for window in range(win_count):
 
         # Print RANSAC progress in 10% increments
-        print_progress(window + 1, win_count, every=0.1)
+        _print_progress(window + 1, win_count, every=0.1)
 
         # Get the current window of EEG data
         start = window * win_size
@@ -509,7 +509,7 @@ def _predict_median_signals_channelwise(
     n_timepts = data.shape[1]
 
     # Before running, make sure we have enough memory
-    verify_free_ram(data, ransac_samples, chunk_size)
+    _verify_free_ram(data, ransac_samples, chunk_size)
 
     # Memory seems to be fine ...
     # Make the predictions
diff --git a/pyprep/utils.py b/pyprep/utils.py
index a84698e..d6cbdca 100644
--- a/pyprep/utils.py
+++ b/pyprep/utils.py
@@ -2,7 +2,6 @@
 import math
 from cmath import sqrt
 
-import mne
 import numpy as np
 import scipy.interpolate
 from psutil import virtual_memory
@@ -17,10 +16,6 @@ def _set_diff(list1, list2):
     return list(set(list1) - set(list2))
 
 
-def _intersect(list1, list2):
-    return list(set(list1).intersection(set(list2)))
-
-
 def _mat_round(x):
     """Round a number to the nearest whole number.
 
@@ -301,7 +296,7 @@ def _correlate_arrays(a, b, matlab_strict=False):
     Notes
     -----
     In MATLAB PREP, RANSAC channel predictions are correlated with actual data
-    using a non-standard method: essentialy, it uses the standard Pearson
+    using a non-standard method: essentially, it uses the standard Pearson
     correlation formula but without subtracting the channel means from each channel
     before calculating sums of squares, i.e.,::
 
@@ -325,7 +320,41 @@ def _correlate_arrays(a, b, matlab_strict=False):
         return np.diag(np.corrcoef(a, b)[:n_chan, n_chan:])
 
 
-def filter_design(N_order, amp, freq):
+def _mad(x, axis=None):
+    """Calculate median absolute deviations from the median (MAD) for an array.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        A 1-D or 2-D numeric array to summarize.
+    axis : {int, tuple of int, None}, optional
+        Axis along which MADs should be calculated. If ``None``, the MAD will
+        be calculated for the full input array. Defaults to ``None``.
+
+    Returns
+    -------
+    mad : scalar or np.ndarray
+        If no axis is specified, returns the MAD for the full input array as a
+        single numeric value. Otherwise, returns an ``np.ndarray`` containing
+        the MAD for each index along the specified axis.
+
+    """
+    # Ensure array is either 1D or 2D
+    x = np.asarray(x)
+    if x.ndim > 2:
+        e = "Only 1D and 2D arrays are supported (input has {0} dimensions)"
+        raise ValueError(e.format(x.ndim))
+
+    # Calculate the median absolute deviation from the median
+    med = np.median(x, axis=axis)
+    if axis == 1:
+        med = med.reshape(-1, 1)  # Transposes array to allow subtraction below
+    mad = np.median(np.abs(x - med), axis=axis)
+
+    return mad
+
+
+def _filter_design(N_order, amp, freq):
     """Create FIR low-pass filter for EEG data using frequency sampling method.
 
     Parameters
@@ -373,7 +402,7 @@ def filter_design(N_order, amp, freq):
     return kernel
 
 
-def split_list(mylist, chunk_size):
+def _split_list(mylist, chunk_size):
     """Split list in chunks.
 
     Parameters
@@ -395,7 +424,7 @@ def split_list(mylist, chunk_size):
     ]
 
 
-def print_progress(current, end, start=None, stepsize=1, every=0.1):
+def _print_progress(current, end, start=None, stepsize=1, every=0.1):
     """Print the current progress in a loop.
 
     Parameters
@@ -433,70 +462,7 @@ def print_progress(current, end, start=None, stepsize=1, every=0.1):
             print("{0}%...".format(int(pct)), end=" ", flush=True)
 
 
-def make_random_mne_object(
-    ch_names,
-    ch_types,
-    times,
-    sfreq,
-    n_freq_comps=5,
-    freq_range=[10, 60],
-    scale=1e-6,
-    RNG=np.random.RandomState(1337),
-):
-    """Make a random MNE object to use for testing.
-
-    Parameters
-    ----------
-    ch_names : list
-        names of channels
-    ch_types : list
-        types of channels
-    times : np.ndarray, shape (1,)
-        Time vector to use.
-    sfreq : float
-        Sampling frequency associated with the time vector.
-    n_freq_comps : int
-        Number of signal components summed to make a signal.
-    freq_range : list, len==2
-        Signals will contain freqs from this range.
-    scale : float
-        Scaling factor applied to the signal in volts. For example 1e-6 to
-        get microvolts.
-    RNG : np.random.RandomState
-        Random state seed.
-
-    Returns
-    -------
-    raw : mne.io.Raw
-        The mne object for performing the tests.
-    n_freq_comps : int
-        Number of random frequency components to introduce.
-    freq_range : tuple
-        The low (`freq_range[0]`) and high (`freq_range[1]`) endpoints of
-        a frequency range from which random draws will be made to
-        introduce frequency components in the test data.
-    """
-    n_chans = len(ch_names)
-    signal_len = times.shape[0]
-    # Make a random signal
-    signal = np.zeros((n_chans, signal_len))
-    low = freq_range[0]
-    high = freq_range[1]
-    for chan in range(n_chans):
-        # Each channel signal is a sum of random freq sine waves
-        for freq_i in range(n_freq_comps):
-            freq = RNG.randint(low, high, signal_len)
-            signal[chan, :] += np.sin(2 * np.pi * times * freq)
-
-    signal *= scale  # scale
-
-    # Make mne object
-    info = mne.create_info(ch_names=ch_names, sfreq=sfreq, ch_types=ch_types)
-    raw = mne.io.RawArray(signal, info)
-    return raw, n_freq_comps, freq_range
-
-
-def verify_free_ram(data, n_samples, n_channels, max_prop=0.95):
+def _verify_free_ram(data, n_samples, n_channels, max_prop=0.95):
     """Check if enough memory is free to run RANSAC with the given parameters.
 
     Parameters
diff --git a/setup.cfg b/setup.cfg
index c109064..74d5e86 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -31,7 +31,6 @@ python_requires = >= 3.7
 install_requires =
   numpy >= 1.14.1
   scipy >= 1.0.0
-  statsmodels >= 0.8.0
   mne >= 0.20.0
   psutil >= 5.4.3
 packages = find:


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 626c0e6..6ef31f9 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,5 +1,6 @@
 """Configure tests."""
 import mne
+import numpy as np
 import pytest
 from mne.datasets import eegbci
 
@@ -66,3 +67,66 @@ def raw_clean(montage):
     raw.set_montage(montage)
 
     return raw
+
+
+def make_random_mne_object(
+    ch_names,
+    ch_types,
+    times,
+    sfreq,
+    n_freq_comps=5,
+    freq_range=[10, 60],
+    scale=1e-6,
+    RNG=np.random.RandomState(1337),
+):
+    """Make a random MNE object to use for testing.
+
+    Parameters
+    ----------
+    ch_names : list
+        names of channels
+    ch_types : list
+        types of channels
+    times : np.ndarray, shape (1,)
+        Time vector to use.
+    sfreq : float
+        Sampling frequency associated with the time vector.
+    n_freq_comps : int
+        Number of signal components summed to make a signal.
+    freq_range : list, len==2
+        Signals will contain freqs from this range.
+    scale : float
+        Scaling factor applied to the signal in volts. For example 1e-6 to
+        get microvolts.
+    RNG : np.random.RandomState
+        Random state seed.
+
+    Returns
+    -------
+    raw : mne.io.Raw
+        The mne object for performing the tests.
+    n_freq_comps : int
+        Number of random frequency components to introduce.
+    freq_range : tuple
+        The low (`freq_range[0]`) and high (`freq_range[1]`) endpoints of
+        a frequency range from which random draws will be made to
+        introduce frequency components in the test data.
+    """
+    n_chans = len(ch_names)
+    signal_len = times.shape[0]
+    # Make a random signal
+    signal = np.zeros((n_chans, signal_len))
+    low = freq_range[0]
+    high = freq_range[1]
+    for chan in range(n_chans):
+        # Each channel signal is a sum of random freq sine waves
+        for freq_i in range(n_freq_comps):
+            freq = RNG.randint(low, high, signal_len)
+            signal[chan, :] += np.sin(2 * np.pi * times * freq)
+
+    signal *= scale  # scale
+
+    # Make mne object
+    info = mne.create_info(ch_names=ch_names, sfreq=sfreq, ch_types=ch_types)
+    raw = mne.io.RawArray(signal, info)
+    return raw, n_freq_comps, freq_range
diff --git a/tests/test_matprep_compare.py b/tests/test_matprep_compare.py
index 02a31c2..9b8b722 100644
--- a/tests/test_matprep_compare.py
+++ b/tests/test_matprep_compare.py
@@ -171,26 +171,26 @@ class TestCompareNoisyChannels(object):
         # Gather PyPREP deviation info and MATLAB equivalents
         mat = matprep_noisy
         matprep_info = {
-            "median_channel_deviation": mat["channelDeviationMedian"],
-            "channel_deviation_sd": mat["channelDeviationSD"],
+            "median_channel_amplitude": mat["channelDeviationMedian"],
+            "channel_amplitude_sd": mat["channelDeviationSD"],
             "robust_channel_deviations": mat["robustChannelDeviation"],
-            "channel_deviations": mat["channelDeviations"],
+            "channel_amplitudes": mat["channelDeviations"],
         }
         pyprep_info = pyprep_noisy._extra_info["bad_by_deviation"]
 
-        # Compare overall medians and SDs for channel deviations
-        median_dev = "median_channel_deviation"
-        dev_sd = "channel_deviation_sd"
-        assert np.isclose(pyprep_info[median_dev] * 1e6, matprep_info[median_dev])
-        assert np.isclose(pyprep_info[dev_sd] * 1e6, matprep_info[dev_sd])
+        # Compare overall medians and SDs for channel amplitudes
+        median_amp = "median_channel_amplitude"
+        amp_sd = "channel_amplitude_sd"
+        assert np.isclose(pyprep_info[median_amp] * 1e6, matprep_info[median_amp])
+        assert np.isclose(pyprep_info[amp_sd] * 1e6, matprep_info[amp_sd])
 
-        # Compare robust deviations across channels
+        # Compare robust amplitude deviations across channels
         dev_by_chan = "robust_channel_deviations"
         assert np.allclose(pyprep_info[dev_by_chan], matprep_info[dev_by_chan])
 
-        # Compare windows of channel deviations across recording
-        chan_devs = "channel_deviations"
-        assert np.allclose(pyprep_info[chan_devs] * 1e6, matprep_info[chan_devs].T)
+        # Compare windows of channel amplitudes across recording
+        chan_amps = "channel_amplitudes"
+        assert np.allclose(pyprep_info[chan_amps] * 1e6, matprep_info[chan_amps].T)
 
         # Compare names of bad-by-deviation channels
         assert pyprep_noisy.bad_by_deviation == matprep_noisy["bads"]["by_deviation"]
diff --git a/tests/test_prep_pipeline.py b/tests/test_prep_pipeline.py
index f2b730b..48ebb54 100644
--- a/tests/test_prep_pipeline.py
+++ b/tests/test_prep_pipeline.py
@@ -6,7 +6,8 @@ import pytest
 import scipy.io as sio
 
 from pyprep.prep_pipeline import PrepPipeline
-from pyprep.utils import make_random_mne_object
+
+from .conftest import make_random_mne_object
 
 
 @pytest.mark.usefixtures("raw", "montage")
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 191e27a..0913200 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -7,10 +7,11 @@ from pyprep.utils import (
     _correlate_arrays,
     _eeglab_create_highpass,
     _get_random_subset,
+    _mad,
     _mat_iqr,
     _mat_quantile,
     _mat_round,
-    print_progress,
+    _print_progress,
 )
 
 
@@ -138,39 +139,56 @@ def test_eeglab_create_highpass():
     NOTE: EEGLAB values were obtained using breakpoints in ``pop_eegfiltnew``,
     since filter creation and data filtering are both done in the same function.
     Values here are first 4 values of the array ``b`` which contains the FIR
-    filter coefficents used by the function.
+    filter coefficients used by the function.
 
     """
-    # Compare initial FIR filter coefficents with EEGLAB
+    # Compare initial FIR filter coefficients with EEGLAB
     expected_vals = [5.3691e-5, 5.4165e-5, 5.4651e-5, 5.5149e-5]
     actual_vals = _eeglab_create_highpass(cutoff=1.0, srate=256)[:4]
     assert all(np.isclose(expected_vals, actual_vals, atol=0.001))
 
-    # Compare middle FIR filter coefficent with EEGLAB
+    # Compare middle FIR filter coefficient with EEGLAB
     vals = _eeglab_create_highpass(cutoff=1.0, srate=256)
     expected_val = 0.9961
     actual_val = vals[len(vals) // 2]
     assert np.isclose(expected_val, actual_val, atol=0.001)
 
 
+def test_mad():
+    """Test the median absolute deviation from the median (MAD) function."""
+    # Generate test data
+    tst = np.array([[1, 2, 3, 4, 8], [80, 10, 20, 30, 40], [100, 200, 800, 300, 400]])
+    expected = np.asarray([1, 10, 100])
+
+    # Compare output to expected results
+    assert all(np.equal(_mad(tst, axis=1), expected))
+    assert all(np.equal(_mad(tst.T, axis=0), expected))
+    assert _mad(tst) == 28  # Matches robust.mad from statsmodels
+
+    # Test exception with > 2-D arrays
+    tst = np.random.rand(3, 3, 3)
+    with pytest.raises(ValueError):
+        _mad(tst, axis=0)
+
+
 def test_print_progress(capsys):
     """Test the function for printing progress updates within a loop."""
     # Test printing start value
-    print_progress(1, 20)
+    _print_progress(1, 20)
     captured = capsys.readouterr()
     assert captured.out == "Progress: "
 
     # Test printing end values
     iterations = 27
     for i in range(iterations):
-        print_progress(i + 1, iterations, every=0.2)
+        _print_progress(i + 1, iterations, every=0.2)
     captured = capsys.readouterr()
     assert captured.out == "Progress: 20%... 40%... 60%... 80%... 100%\n"
 
     # Test printing of updates at right times
     iterations = 176
     for i in range(iterations):
-        print_progress(i + 1, iterations)
+        _print_progress(i + 1, iterations)
         if (i + 1) == 17:
             captured = capsys.readouterr()
             assert captured.out == "Progress: "
@@ -183,7 +201,7 @@ def test_print_progress(capsys):
     iterations = 25
     start = 5
     for i in range(start, iterations + 1):
-        print_progress(i, iterations, start=start)
+        _print_progress(i, iterations, start=start)
         if i == 6:
             captured = capsys.readouterr()
             assert captured.out == "Progress: "