Coerce tau2 to an array in permutation_test (#94)

* Coerce tau2 to an array in permutation_test.

* Update plot_run_meta-analysis.py

* Update plot_run_meta-analysis.py

* Update plot_run_meta-analysis.py

examples/02_meta-analysis/plot_run_meta-analysis.py

@@ -31,20 +31,37 @@
 ###############################################################################
 # Now we fit a model
 # -----------------------------------------------------------------------------
+# You must first initialize the estimator, after which you can use
+# :meth:`~pymare.estimators.estimators.BaseEstimator.fit` to fit the model to
+# numpy arrays, or
+# :meth:`~pymare.estimators.estimators.BaseEstimator.fit_dataset` to fit it to
+# a :class:`~pymare.core.Dataset`.
+#
+# The :meth:`~pymare.estimators.estimators.BaseEstimator.summary` function
+# will return a :class:`~pymare.results.MetaRegressionResults` object,
+# which contains the results of the analysis.
 est = estimators.WeightedLeastSquares().fit_dataset(dset)
 results = est.summary()
 results.to_df()
 
 ###############################################################################
 # We can also extract some useful information from the results object
 # -----------------------------------------------------------------------------
-# Here we'll take a look at the heterogeneity statistics.
+# The :meth:`~pymare.results.MetaRegressionResults.get_heterogeneity_stats`
+# method will calculate heterogeneity statistics.
 pprint(results.get_heterogeneity_stats())
 
 ###############################################################################
-# We can also estimate the confidence interval for :math:`\tau^2`.
+# The :meth:`~pymare.results.MetaRegressionResults.get_re_stats` method will
+# estimate the confidence interval for :math:`\tau^2`.
 pprint(results.get_re_stats())
 
+###############################################################################
+# The :meth:`~pymare.results.MetaRegressionResults.permutation_test` method
+# will run a permutation test to estimate more accurate p-values.
+perm_results = results.permutation_test(n_perm=1000)
+perm_results.to_df()
+
 ###############################################################################
 # References
 # -----------------------------------------------------------------------------

pymare/results.py

@@ -116,18 +116,18 @@ def get_re_stats(self, method="QP", alpha=0.05):
         .. footbibliography::
         """
         if method == "QP":
-            n_iters = np.atleast_2d(self.tau2).shape[1]
-            if n_iters > 10:
+            n_datasets = np.atleast_2d(self.tau2).shape[1]
+            if n_datasets > 10:
                 warn(
                     "Method 'QP' is not parallelized; it may take a while to "
-                    "compute CIs for {} parallel tau^2 values.".format(n_iters)
+                    f"compute CIs for {n_datasets} parallel tau^2 values."
                 )
 
             # Make sure we have an estimate of v if it wasn't observed
             v = self.estimator.get_v(self.dataset)
 
             cis = []
-            for i in range(n_iters):
+            for i in range(n_datasets):
                 args = {
                     "y": self.dataset.y[:, i],
                     "v": v[:, i],
@@ -266,6 +266,11 @@ def permutation_test(self, n_perm=1000):
         fe_stats = self.get_fe_stats()
         re_stats = self.get_re_stats()
 
+        # Ensure that tau2 is an array
+        tau2 = re_stats["tau^2"]
+        if not isinstance(tau2, (list, tuple, np.ndarray)):
+            tau2 = np.full(n_datasets, tau2)
+
         # create results arrays
         fe_p = np.zeros_like(self.fe_params)
         rfx = self.tau2 is not None
@@ -323,7 +328,7 @@ def permutation_test(self, n_perm=1000):
                 fe_obs = fe_obs[:, None]
             fe_p[:, i] = (np.abs(fe_obs) < np.abs(params["fe_params"])).mean(1)
             if rfx:
-                abs_obs = np.abs(re_stats["tau^2"][i])
+                abs_obs = np.abs(tau2[i])
                 tau_p[i] = (abs_obs < np.abs(params["tau2"])).mean()
 
         # p-values can't be smaller than 1/n_perm