Repository: neurostuff/PyMARE

Issue:
Results objects break if user called fit() rather than fit_dataset()
The various results classes all currently require access to the estimator's last-fitted `Dataset`, which won't exist now following the API change to allow `fit()` to be called with numpy arrays. The solution is to either implicitly create and store a `Dataset`, which would be easier, but feels counterproductive (since the point of eliminating `Dataset` inputs from `fit()` was to cut down overhead), or to directly store references to the input arrays rather than to a `Dataset` in the estimator (in which case we have to add a bunch of code to every estimator, and lose some of the benefits of the `Dataset abstraction).

Discussion:
codecov-commenter: # [Codecov](https://codecov.io/gh/neurostuff/PyMARE/pull/104?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff) Report
> Merging [#104](https://codecov.io/gh/neurostuff/PyMARE/pull/104?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff) (7bbc570) into [master](https://codecov.io/gh/neurostuff/PyMARE/commit/541894e8b50ac9ec1b8d746d71397ced961bc1e0?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff) (541894e) will **increase** coverage by `0.13%`.
> The diff coverage is `94.11%`.

```diff
@@            Coverage Diff             @@
##           master     #104      +/-   ##
==========================================
+ Coverage   87.44%   87.57%   +0.13%     
==========================================
  Files          13       13              
  Lines         844      861      +17     
==========================================
+ Hits          738      754      +16     
- Misses        106      107       +1     
```


| [Impacted Files](https://codecov.io/gh/neurostuff/PyMARE/pull/104?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff) | Coverage Δ | |
|---|---|---|
| [pymare/estimators/estimators.py](https://codecov.io/gh/neurostuff/PyMARE/pull/104/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff#diff-cHltYXJlL2VzdGltYXRvcnMvZXN0aW1hdG9ycy5weQ==) | `83.64% <83.33%> (-0.01%)` | :arrow_down: |
| [pymare/estimators/combination.py](https://codecov.io/gh/neurostuff/PyMARE/pull/104/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff#diff-cHltYXJlL2VzdGltYXRvcnMvY29tYmluYXRpb24ucHk=) | `92.15% <100.00%> (+0.15%)` | :arrow_up: |
| [pymare/results.py](https://codecov.io/gh/neurostuff/PyMARE/pull/104/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff#diff-cHltYXJlL3Jlc3VsdHMucHk=) | `82.46% <100.00%> (+0.87%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/neurostuff/PyMARE/pull/104?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff).
> **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=neurostuff)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/neurostuff/PyMARE/pull/104?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff). Last update [541894e...7bbc570](https://codecov.io/gh/neurostuff/PyMARE/pull/104?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=neurostuff). 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=neurostuff).


Patch:
diff --git a/pymare/estimators/combination.py b/pymare/estimators/combination.py
index 8a8d6ce..f2ffabc 100644
--- a/pymare/estimators/combination.py
+++ b/pymare/estimators/combination.py
@@ -36,6 +36,9 @@ class CombinationTest(BaseEstimator):
 
     def fit(self, z, *args, **kwargs):
         """Fit the estimator to z-values."""
+        # This resets the Estimator's dataset_ attribute. fit_dataset will overwrite if called.
+        self.dataset_ = None
+
         if self.mode == "concordant":
             ose = self.__class__(mode="directed")
             p1 = ose.p_value(z, *args, **kwargs)
diff --git a/pymare/estimators/estimators.py b/pymare/estimators/estimators.py
index 7149c3f..7a2b1a9 100644
--- a/pymare/estimators/estimators.py
+++ b/pymare/estimators/estimators.py
@@ -191,7 +191,24 @@ class WeightedLeastSquares(BaseEstimator):
         self.tau2 = tau2
 
     def fit(self, y, X, v=None):
-        """Fit the estimator to data."""
+        """Fit the estimator to data.
+
+        Parameters
+        ----------
+        y : :obj:`numpy.ndarray` of shape (n, d)
+            The dependent variable(s) (y).
+        X : :obj:`numpy.ndarray` of shape (n, p)
+            The independent variable(s) (X).
+        v : :obj:`numpy.ndarray` of shape (n, d), optional
+            Sampling variances. If not provided, unit weights will be used.
+
+        Returns
+        -------
+        :obj:`~pymare.estimators.WeightedLeastSquares`
+        """
+        # This resets the Estimator's dataset_ attribute. fit_dataset will overwrite if called.
+        self.dataset_ = None
+
         if v is None:
             v = np.ones_like(y)
 
@@ -219,7 +236,24 @@ class DerSimonianLaird(BaseEstimator):
     """
 
     def fit(self, y, v, X):
-        """Fit the estimator to data."""
+        """Fit the estimator to data.
+
+        Parameters
+        ----------
+        y : :obj:`numpy.ndarray` of shape (n, d)
+            The dependent variable(s) (y).
+        v : :obj:`numpy.ndarray` of shape (n, d)
+            Sampling variances.
+        X : :obj:`numpy.ndarray` of shape (n, p)
+            The independent variable(s) (X).
+
+        Returns
+        -------
+        :obj:`~pymare.estimators.DerSimonianLaird`
+        """
+        # This resets the Estimator's dataset_ attribute. fit_dataset will overwrite if called.
+        self.dataset_ = None
+
         y = ensure_2d(y)
         v = ensure_2d(v)
 
@@ -266,7 +300,24 @@ class Hedges(BaseEstimator):
     """
 
     def fit(self, y, v, X):
-        """Fit the estimator to data."""
+        """Fit the estimator to data.
+
+        Parameters
+        ----------
+        y : :obj:`numpy.ndarray` of shape (n, d)
+            The dependent variable(s) (y).
+        v : :obj:`numpy.ndarray` of shape (n, d)
+            Sampling variances.
+        X : :obj:`numpy.ndarray` of shape (n, p)
+            The independent variable(s) (X).
+
+        Returns
+        -------
+        :obj:`~pymare.estimators.Hedges`
+        """
+        # This resets the Estimator's dataset_ attribute. fit_dataset will overwrite if called.
+        self.dataset_ = None
+
         k, p = X.shape[:2]
         _unit_v = np.ones_like(y)
         beta, inv_cov = weighted_least_squares(y, _unit_v, X, return_cov=True)
@@ -317,7 +368,24 @@ class VarianceBasedLikelihoodEstimator(BaseEstimator):
 
     @_loopable
     def fit(self, y, v, X):
-        """Fit the estimator to data."""
+        """Fit the estimator to data.
+
+        Parameters
+        ----------
+        y : :obj:`numpy.ndarray` of shape (n, d)
+            The dependent variable(s) (y).
+        v : :obj:`numpy.ndarray` of shape (n, d)
+            Sampling variances.
+        X : :obj:`numpy.ndarray` of shape (n, p)
+            The independent variable(s) (X).
+
+        Returns
+        -------
+        :obj:`~pymare.estimators.VarianceBasedLikelihoodEstimator`
+        """
+        # This resets the Estimator's dataset_ attribute. fit_dataset will overwrite if called.
+        self.dataset_ = None
+
         # use D-L estimate for initial values
         est_DL = DerSimonianLaird().fit(y, v, X).params_
         beta = est_DL["fe_params"]
@@ -393,7 +461,24 @@ class SampleSizeBasedLikelihoodEstimator(BaseEstimator):
 
     @_loopable
     def fit(self, y, n, X):
-        """Fit the estimator to data."""
+        """Fit the estimator to data.
+
+        Parameters
+        ----------
+        y : :obj:`numpy.ndarray` of shape (n, d)
+            The dependent variable(s) (y).
+        n : :obj:`numpy.ndarray` of shape (n, d)
+            Sample sizes.
+        X : :obj:`numpy.ndarray` of shape (n, p)
+            The independent variable(s) (X).
+
+        Returns
+        -------
+        :obj:`~pymare.estimators.SampleSizeBasedLikelihoodEstimator`
+        """
+        # This resets the Estimator's dataset_ attribute. fit_dataset will overwrite if called.
+        self.dataset_ = None
+
         if n.std() < np.sqrt(np.finfo(float).eps):
             raise ValueError(
                 "Sample size-based likelihood estimator cannot "
@@ -550,6 +635,9 @@ class StanMetaRegression(BaseEstimator):
         `groups` argument can be used to specify the nesting structure
         (i.e., which rows in `y`, `v`, and `X` belong to each study).
         """
+        # This resets the Estimator's dataset_ attribute. fit_dataset will overwrite if called.
+        self.dataset_ = None
+
         if y.ndim > 1 and y.shape[1] > 1:
             raise ValueError(
                 "The StanMetaRegression estimator currently does "
diff --git a/pymare/results.py b/pymare/results.py
index 5da53e3..b84b1f6 100644
--- a/pymare/results.py
+++ b/pymare/results.py
@@ -34,6 +34,12 @@ class MetaRegressionResults:
     tau2 : None or :obj:`numpy.ndarray` of shape (d,) or :obj:`float`, optional
         A 1-d array containing the estimated tau^2 value for each parallel dataset
         (or a float, for a single dataset). May be omitted by fixed-effects estimators.
+
+    Warning
+    -------
+    When an Estimator is fitted to arrays directly using the ``fit`` method, the Results object's
+    utility is limited.
+    Many methods will not work.
     """
 
     def __init__(self, estimator, dataset, fe_params, fe_cov, tau2=None):
@@ -47,7 +53,10 @@ class MetaRegressionResults:
     @lru_cache(maxsize=1)
     def fe_se(self):
         """Get fixed-effect standard error."""
-        cov = np.atleast_3d(self.fe_cov)  # 3rd dim is for parallel datasets
+        cov = self.fe_cov.copy()
+        if cov.ndim == 2:
+            cov = cov[None, :, :]
+
         return np.sqrt(np.diagonal(cov)).T
 
     @lru_cache(maxsize=16)
@@ -92,6 +101,11 @@ class MetaRegressionResults:
     def get_re_stats(self, method="QP", alpha=0.05):
         """Get random-effect statistics.
 
+        .. warning::
+
+            This method relies on the ``.dataset`` attribute, so the original Estimator must have
+            be fitted with ``fit_dataset``, not ``fit``.
+
         Parameters
         ----------
         method : {"QP"}, optional
@@ -115,6 +129,9 @@ class MetaRegressionResults:
         ----------
         .. footbibliography::
         """
+        if self.dataset is None:
+            raise ValueError("The Dataset is unavailable. This method requires a Dataset.")
+
         if method == "QP":
             n_datasets = np.atleast_2d(self.tau2).shape[1]
             if n_datasets > 10:
@@ -157,6 +174,11 @@ class MetaRegressionResults:
     def get_heterogeneity_stats(self):
         """Get heterogeneity statistics.
 
+        .. warning::
+
+            This method relies on the ``.dataset`` attribute, so the original Estimator must have
+            be fitted with ``fit_dataset``, not ``fit``.
+
         Returns
         -------
         :obj:`dict`
@@ -180,6 +202,9 @@ class MetaRegressionResults:
         ----------
         .. footbibliography::
         """
+        if self.dataset is None:
+            raise ValueError("The Dataset is unavailable. This method requires a Dataset.")
+
         v = self.estimator.get_v(self.dataset)
         q_fe = q_gen(self.dataset.y, v, self.dataset.X, 0)
         df = self.dataset.y.shape[0] - self.dataset.X.shape[1]
@@ -195,6 +220,11 @@ class MetaRegressionResults:
 
             This method only works for one-dimensional results.
 
+        .. warning::
+
+            This method relies on the ``.dataset`` attribute, so the original Estimator must have
+            be fitted with ``fit_dataset``, not ``fit``.
+
         Parameters
         ----------
         alpha : :obj:`float`, optional
@@ -217,6 +247,9 @@ class MetaRegressionResults:
                         the CI columns will be ``"ci_0.025"`` and ``"ci_0.975"``.
             =========== ==========================================================================
         """
+        if self.dataset is None:
+            raise ValueError("The Dataset is unavailable. This method requires a Dataset.")
+
         b_shape = self.fe_params.shape
         if len(b_shape) > 1 and b_shape[1] > 1:
             raise ValueError(
@@ -237,6 +270,11 @@ class MetaRegressionResults:
     def permutation_test(self, n_perm=1000):
         """Run permutation test.
 
+        .. warning::
+
+            This method relies on the ``.dataset`` attribute, so the original Estimator must have
+            be fitted with ``fit_dataset``, not ``fit``.
+
         Parameters
         ----------
         n_perm : :obj:`int`, optional
@@ -260,6 +298,9 @@ class MetaRegressionResults:
         This means that one can often set very high n_perm values (e.g., 100k) with little
         performance degradation.
         """
+        if self.dataset is None:
+            raise ValueError("The Dataset is unavailable. This method requires a Dataset.")
+
         n_obs, n_datasets = self.dataset.y.shape
         has_mods = self.dataset.X.shape[1] > 1
 
@@ -381,6 +422,11 @@ class CombinationTestResults:
     def permutation_test(self, n_perm=1000):
         """Run permutation test.
 
+        .. warning::
+
+            This method relies on the ``.dataset`` attribute, so the original Estimator must have
+            be fitted with ``fit_dataset``, not ``fit``.
+
         Parameters
         ----------
         n_perm : :obj:`int`, optional
@@ -403,6 +449,9 @@ class CombinationTestResults:
         set very high n_perm values (e.g., 100k) with little performance
         degradation.
         """
+        if self.dataset is None:
+            raise ValueError("The Dataset is unavailable. This method requires a Dataset.")
+
         n_obs, n_datasets = self.dataset.y.shape
 
         # create results arrays


Test patch:
diff --git a/pymare/tests/test_estimators.py b/pymare/tests/test_estimators.py
index c4531f2..339e163 100644
--- a/pymare/tests/test_estimators.py
+++ b/pymare/tests/test_estimators.py
@@ -18,6 +18,19 @@ def test_weighted_least_squares_estimator(dataset):
     est = WeightedLeastSquares().fit_dataset(dataset)
     results = est.summary()
     beta, tau2 = results.fe_params, results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
+    assert beta.shape == (2, 1)
+    assert isinstance(tau2, float)
+    assert fe_stats["est"].shape == (2, 1)
+    assert fe_stats["se"].shape == (2, 1)
+    assert fe_stats["ci_l"].shape == (2, 1)
+    assert fe_stats["ci_u"].shape == (2, 1)
+    assert fe_stats["z"].shape == (2, 1)
+    assert fe_stats["p"].shape == (2, 1)
+
+    # Check output values
     assert np.allclose(beta.ravel(), [-0.2725, 0.6935], atol=1e-4)
     assert tau2 == 0.0
 
@@ -35,6 +48,19 @@ def test_dersimonian_laird_estimator(dataset):
     est = DerSimonianLaird().fit_dataset(dataset)
     results = est.summary()
     beta, tau2 = results.fe_params, results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
+    assert beta.shape == (2, 1)
+    assert tau2.shape == (1,)
+    assert fe_stats["est"].shape == (2, 1)
+    assert fe_stats["se"].shape == (2, 1)
+    assert fe_stats["ci_l"].shape == (2, 1)
+    assert fe_stats["ci_u"].shape == (2, 1)
+    assert fe_stats["z"].shape == (2, 1)
+    assert fe_stats["p"].shape == (2, 1)
+
+    # Check output values
     assert np.allclose(beta.ravel(), [-0.1070, 0.7664], atol=1e-4)
     assert np.allclose(tau2, 8.3627, atol=1e-4)
 
@@ -43,9 +69,19 @@ def test_2d_DL_estimator(dataset_2d):
     """Test DerSimonianLaird estimator on 2D Dataset."""
     results = DerSimonianLaird().fit_dataset(dataset_2d).summary()
     beta, tau2 = results.fe_params, results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
     assert beta.shape == (2, 3)
     assert tau2.shape == (3,)
+    assert fe_stats["est"].shape == (2, 3)
+    assert fe_stats["se"].shape == (2, 3)
+    assert fe_stats["ci_l"].shape == (2, 3)
+    assert fe_stats["ci_u"].shape == (2, 3)
+    assert fe_stats["z"].shape == (2, 3)
+    assert fe_stats["p"].shape == (2, 3)
 
+    # Check output values
     # First and third sets are identical to previous DL test; second set is
     # randomly different.
     assert np.allclose(beta[:, 0], [-0.1070, 0.7664], atol=1e-4)
@@ -64,6 +100,19 @@ def test_hedges_estimator(dataset):
     est = Hedges().fit_dataset(dataset)
     results = est.summary()
     beta, tau2 = results.fe_params, results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
+    assert beta.shape == (2, 1)
+    assert tau2.shape == (1,)
+    assert fe_stats["est"].shape == (2, 1)
+    assert fe_stats["se"].shape == (2, 1)
+    assert fe_stats["ci_l"].shape == (2, 1)
+    assert fe_stats["ci_u"].shape == (2, 1)
+    assert fe_stats["z"].shape == (2, 1)
+    assert fe_stats["p"].shape == (2, 1)
+
+    # Check output values
     assert np.allclose(beta.ravel(), [-0.1066, 0.7704], atol=1e-4)
     assert np.allclose(tau2, 11.3881, atol=1e-4)
 
@@ -72,8 +121,17 @@ def test_2d_hedges(dataset_2d):
     """Test Hedges estimator on 2D Dataset."""
     results = Hedges().fit_dataset(dataset_2d).summary()
     beta, tau2 = results.fe_params, results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
     assert beta.shape == (2, 3)
     assert tau2.shape == (3,)
+    assert fe_stats["est"].shape == (2, 3)
+    assert fe_stats["se"].shape == (2, 3)
+    assert fe_stats["ci_l"].shape == (2, 3)
+    assert fe_stats["ci_u"].shape == (2, 3)
+    assert fe_stats["z"].shape == (2, 3)
+    assert fe_stats["p"].shape == (2, 3)
 
     # First and third sets are identical to single dim test; second set is
     # randomly different.
@@ -91,6 +149,19 @@ def test_variance_based_maximum_likelihood_estimator(dataset):
     est = VarianceBasedLikelihoodEstimator(method="ML").fit_dataset(dataset)
     results = est.summary()
     beta, tau2 = results.fe_params, results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
+    assert beta.shape == (2, 1)
+    assert tau2.shape == (1, 1)
+    assert fe_stats["est"].shape == (2, 1)
+    assert fe_stats["se"].shape == (2, 1)
+    assert fe_stats["ci_l"].shape == (2, 1)
+    assert fe_stats["ci_u"].shape == (2, 1)
+    assert fe_stats["z"].shape == (2, 1)
+    assert fe_stats["p"].shape == (2, 1)
+
+    # Check output values
     assert np.allclose(beta.ravel(), [-0.1072, 0.7653], atol=1e-4)
     assert np.allclose(tau2, 7.7649, atol=1e-4)
 
@@ -101,6 +172,19 @@ def test_variance_based_restricted_maximum_likelihood_estimator(dataset):
     est = VarianceBasedLikelihoodEstimator(method="REML").fit_dataset(dataset)
     results = est.summary()
     beta, tau2 = results.fe_params, results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
+    assert beta.shape == (2, 1)
+    assert tau2.shape == (1, 1)
+    assert fe_stats["est"].shape == (2, 1)
+    assert fe_stats["se"].shape == (2, 1)
+    assert fe_stats["ci_l"].shape == (2, 1)
+    assert fe_stats["ci_u"].shape == (2, 1)
+    assert fe_stats["z"].shape == (2, 1)
+    assert fe_stats["p"].shape == (2, 1)
+
+    # Check output values
     assert np.allclose(beta.ravel(), [-0.1066, 0.7700], atol=1e-4)
     assert np.allclose(tau2, 10.9499, atol=1e-4)
 
@@ -113,6 +197,20 @@ def test_sample_size_based_maximum_likelihood_estimator(dataset_n):
     beta = results.fe_params
     sigma2 = results.estimator.params_["sigma2"]
     tau2 = results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
+    assert beta.shape == (1, 1)
+    assert sigma2.shape == (1, 1)
+    assert tau2.shape == (1, 1)
+    assert fe_stats["est"].shape == (1, 1)
+    assert fe_stats["se"].shape == (1, 1)
+    assert fe_stats["ci_l"].shape == (1, 1)
+    assert fe_stats["ci_u"].shape == (1, 1)
+    assert fe_stats["z"].shape == (1, 1)
+    assert fe_stats["p"].shape == (1, 1)
+
+    # Check output values
     assert np.allclose(beta, [-2.0951], atol=1e-4)
     assert np.allclose(sigma2, 12.777, atol=1e-3)
     assert np.allclose(tau2, 2.8268, atol=1e-4)
@@ -126,6 +224,20 @@ def test_sample_size_based_restricted_maximum_likelihood_estimator(dataset_n):
     beta = results.fe_params
     sigma2 = results.estimator.params_["sigma2"]
     tau2 = results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
+    assert beta.shape == (1, 1)
+    assert sigma2.shape == (1, 1)
+    assert tau2.shape == (1, 1)
+    assert fe_stats["est"].shape == (1, 1)
+    assert fe_stats["se"].shape == (1, 1)
+    assert fe_stats["ci_l"].shape == (1, 1)
+    assert fe_stats["ci_u"].shape == (1, 1)
+    assert fe_stats["z"].shape == (1, 1)
+    assert fe_stats["p"].shape == (1, 1)
+
+    # Check output values
     assert np.allclose(beta, [-2.1071], atol=1e-4)
     assert np.allclose(sigma2, 13.048, atol=1e-3)
     assert np.allclose(tau2, 3.2177, atol=1e-4)
@@ -136,9 +248,19 @@ def test_2d_looping(dataset_2d):
     est = VarianceBasedLikelihoodEstimator().fit_dataset(dataset_2d)
     results = est.summary()
     beta, tau2 = results.fe_params, results.tau2
+    fe_stats = results.get_fe_stats()
+
+    # Check output shapes
     assert beta.shape == (2, 3)
     assert tau2.shape == (1, 3)
+    assert fe_stats["est"].shape == (2, 3)
+    assert fe_stats["se"].shape == (2, 3)
+    assert fe_stats["ci_l"].shape == (2, 3)
+    assert fe_stats["ci_u"].shape == (2, 3)
+    assert fe_stats["z"].shape == (2, 3)
+    assert fe_stats["p"].shape == (2, 3)
 
+    # Check output values
     # First and third sets are identical to single dim test; 2nd is different
     assert np.allclose(beta[:, 0], [-0.1072, 0.7653], atol=1e-4)
     assert np.allclose(tau2[0, 0], 7.7649, atol=1e-4)
diff --git a/pymare/tests/test_results.py b/pymare/tests/test_results.py
index 28aaacc..2baf3bd 100644
--- a/pymare/tests/test_results.py
+++ b/pymare/tests/test_results.py
@@ -33,6 +33,71 @@ def results_2d(fitted_estimator, dataset_2d):
     return est.fit_dataset(dataset_2d).summary()
 
 
+def test_meta_regression_results_from_arrays(dataset):
+    """Ensure that a MetaRegressionResults can be created from arrays.
+
+    This is a regression test for a bug that caused the MetaRegressionResults
+    to fail when Estimators were fitted to arrays instead of Datasets.
+    See https://github.com/neurostuff/PyMARE/issues/52 for more info.
+    """
+    est = DerSimonianLaird()
+    fitted_estimator = est.fit(y=dataset.y, X=dataset.X, v=dataset.v)
+    results = fitted_estimator.summary()
+    assert isinstance(results, MetaRegressionResults)
+    assert results.fe_params.shape == (2, 1)
+    assert results.fe_cov.shape == (2, 2, 1)
+    assert results.tau2.shape == (1,)
+
+    # fit overwrites dataset_ attribute with None
+    assert fitted_estimator.dataset_ is None
+    # fit_dataset overwrites it with the Dataset
+    fitted_estimator.fit_dataset(dataset)
+    assert isinstance(fitted_estimator.dataset_, Dataset)
+    # fit sets it back to None
+    fitted_estimator.fit(y=dataset.y, X=dataset.X, v=dataset.v)
+    assert fitted_estimator.dataset_ is None
+
+    # Some methods are not available if fit was used
+    results = fitted_estimator.summary()
+    with pytest.raises(ValueError):
+        results.get_re_stats()
+
+    with pytest.raises(ValueError):
+        results.get_heterogeneity_stats()
+
+    with pytest.raises(ValueError):
+        results.to_df()
+
+    with pytest.raises(ValueError):
+        results.permutation_test(1000)
+
+
+def test_combination_test_results_from_arrays(dataset):
+    """Ensure that a CombinationTestResults can be created from arrays.
+
+    This is a regression test for a bug that caused the MetaRegressionResults
+    to fail when Estimators were fitted to arrays instead of Datasets.
+    See https://github.com/neurostuff/PyMARE/issues/52 for more info.
+    """
+    fitted_estimator = StoufferCombinationTest().fit(z=dataset.y)
+    results = fitted_estimator.summary()
+    assert isinstance(results, CombinationTestResults)
+    assert results.p.shape == (1,)
+
+    # fit overwrites dataset_ attribute with None
+    assert fitted_estimator.dataset_ is None
+    # fit_dataset overwrites it with the Dataset
+    fitted_estimator.fit_dataset(dataset)
+    assert isinstance(fitted_estimator.dataset_, Dataset)
+    # fit sets it back to None
+    fitted_estimator.fit(z=dataset.y)
+    assert fitted_estimator.dataset_ is None
+
+    # Some methods are not available if fit was used
+    with pytest.raises(ValueError):
+        fitted_estimator.summary().permutation_test(1000)
+
+
 def test_meta_regression_results_init_1d(fitted_estimator):
     """Test MetaRegressionResults from 1D data."""
     est = fitted_estimator