diff --git a/examples/02_meta-analysis/plot_run_meta-analysis.py b/examples/02_meta-analysis/plot_run_meta-analysis.py
index ecde5a6..c18af3c 100644
--- a/examples/02_meta-analysis/plot_run_meta-analysis.py
+++ b/examples/02_meta-analysis/plot_run_meta-analysis.py
@@ -14,8 +14,6 @@
 # -----------------------------------------------------------------------------
 from pprint import pprint
 
-from pymare import core, datasets, estimators
-
 ###############################################################################
 # Load the data
 # -----------------------------------------------------------------------------
@@ -24,8 +22,14 @@
 #
 # We only want to do a mean analysis, so we won't have any covariates except for
 # an intercept.
+import numpy as np
+
+from pymare import core, datasets, estimators
+
 data, meta = datasets.michael2013()
-dset = core.Dataset(data=data, y="yi", v="vi", X=None, add_intercept=True)
+data["age"] = np.random.randint(1, 100, size=data.shape[0])
+data["n_houses"] = np.random.randint(1, 10, size=data.shape[0])
+dset = core.Dataset(data=data, y="yi", v="vi", X=["age", "n_houses"], add_intercept=True)
 dset.to_df()
 
 ###############################################################################
@@ -62,6 +66,9 @@
 perm_results = results.permutation_test(n_perm=1000)
 perm_results.to_df()
 
+###############################################################################
+print(results.summary())
+
 ###############################################################################
 # References
 # -----------------------------------------------------------------------------
diff --git a/pymare/results.py b/pymare/results.py
index b84b1f6..95a7282 100644
--- a/pymare/results.py
+++ b/pymare/results.py
@@ -14,6 +14,7 @@
     az = None
 
 from pymare.stats import q_gen, q_profile
+from pymare.utils import _get_sig_code
 
 
 class MetaRegressionResults:
@@ -327,6 +328,10 @@ def permutation_test(self, n_perm=1000):
 
         exact = n_exact < n_perm
         if exact:
+            LGR.warning(
+                "The number of requested permutations is larger than the number of possible ones. "
+                f"An exact test with {n_exact} permutations will be performed."
+            )
             n_perm = n_exact
 
         # Loop over parallel datasets
@@ -379,6 +384,55 @@ def permutation_test(self, n_perm=1000):
 
         return PermutationTestResults(self, params, n_perm, exact)
 
+    def summary(self):
+        """Print a summary of the results.
+
+        Returns
+        -------
+        :obj:`str`
+            A summary of the results, modeled after ``metafor``'s results.
+        """
+        n_obs, n_datasets = self.dataset.y.shape
+        df = n_obs - self.dataset.X.shape[1]
+
+        if n_datasets > 1:
+            raise ValueError(
+                "More than one set of results found! "
+                "A summary table cannot be displayed for multidimensional results at the moment."
+            )
+
+        re_stats = self.get_re_stats()
+        heterogeneity_stats = self.get_heterogeneity_stats()
+
+        # Extract info
+        tau2 = re_stats["tau^2"]
+        if isinstance(tau2, (list, tuple, np.ndarray)):
+            tau2 = tau2[0]
+
+        tau2_se = "n/a"
+        tau2_estimator = self.estimator.__class__.__name__
+
+        model_results = self.to_df()
+        model_results[""] = model_results["p-value"].apply(_get_sig_code)
+
+        pattern = f"""Random-Effects Model (k = {n_obs}; tau^2 estimator: {tau2_estimator})
+
+tau^2 (estimated amount of total heterogeneity): {tau2:.04f} (SE = {tau2_se})
+tau (square root of estimated tau^2 value):      {np.sqrt(tau2):.04f}
+I^2 (total heterogeneity / total variability):   {heterogeneity_stats['I^2'][0]:.2f}%
+H^2 (total variability / sampling variability):  {heterogeneity_stats['H'][0] ** 2:.2f}
+
+Test for Heterogeneity:
+Q(df = {df}) = {heterogeneity_stats['Q'][0]:.04f}, p-val = {heterogeneity_stats['p(Q)'][0]:.04f}
+
+Model Results:
+
+{model_results.to_string(index=False, float_format="%.4f", na_rep="n/a")}
+
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"""
+        return pattern
+
 
 class CombinationTestResults:
     """Container for results generated by p-value combination methods.
@@ -460,6 +514,11 @@ def permutation_test(self, n_perm=1000):
         # Calculate # of permutations and determine whether to use exact test
         n_exact = 2**n_obs
         if n_exact < n_perm:
+            LGR.warning(
+                "The number of requested permutations is larger than the number of possible ones. "
+                f"An exact test with {n_exact} permutations will be performed."
+            )
+
             perms = np.array(list(itertools.product([-1, 1], repeat=n_obs))).T
             exact = True
             n_perm = n_exact
@@ -497,6 +556,44 @@ def permutation_test(self, n_perm=1000):
 
         return PermutationTestResults(self, {"fe_p": p_p}, n_perm, exact)
 
+    def summary(self):
+        """Print a summary of the results.
+
+        Returns
+        -------
+        :obj:`str`
+            A summary of the results, modeled after ``metafor``'s results.
+        """
+        n_obs, n_datasets = self.dataset.y.shape
+        df = n_obs - self.dataset.X.shape[1]
+
+        if n_datasets > 1:
+            raise ValueError(
+                "More than one set of results found! "
+                "A summary table cannot be displayed for multidimensional results at the moment."
+            )
+
+        model_results = self.to_df()
+        model_results[""] = model_results["p-value"].apply(_get_sig_code)
+
+        pattern = f"""Random-Effects Model (k = {n_obs}; tau^2 estimator: {tau2_estimator})
+
+tau^2 (estimated amount of total heterogeneity): {tau2:.04f} (SE = {tau2_se})
+tau (square root of estimated tau^2 value):      {np.sqrt(tau2):.04f}
+I^2 (total heterogeneity / total variability):   {heterogeneity_stats['I^2'][0]:.2f}%
+H^2 (total variability / sampling variability):  {heterogeneity_stats['H'][0] ** 2:.2f}
+
+Test for Heterogeneity:
+Q(df = {df}) = {heterogeneity_stats['Q'][0]:.04f}, p-val = {heterogeneity_stats['p(Q)'][0]:.04f}
+
+Model Results:
+
+{model_results.to_string(index=False, float_format="%.4f", na_rep="n/a")}
+
+---
+Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"""
+        return pattern
+
 
 class PermutationTestResults:
     """Lightweight container to hold and display permutation test results."""
diff --git a/pymare/tests/test_results.py b/pymare/tests/test_results.py
index 2baf3bd..b85fee1 100644
--- a/pymare/tests/test_results.py
+++ b/pymare/tests/test_results.py
@@ -108,6 +108,7 @@ def test_meta_regression_results_init_1d(fitted_estimator):
     assert results.fe_params.shape == (2, 1)
     assert results.fe_cov.shape == (2, 2, 1)
     assert results.tau2.shape == (1,)
+    assert isinstance(results.summary(), str)
 
 
 def test_meta_regression_results_init_2d(results_2d):
@@ -116,6 +117,8 @@ def test_meta_regression_results_init_2d(results_2d):
     assert results_2d.fe_params.shape == (2, 3)
     assert results_2d.fe_cov.shape == (2, 2, 3)
     assert results_2d.tau2.shape == (1, 3)
+    with pytest.raises(ValueError):
+        results_2d.summary()
 
 
 def test_mrr_fe_se(results, results_2d):
diff --git a/pymare/utils.py b/pymare/utils.py
index 1e11cb7..af58366 100644
--- a/pymare/utils.py
+++ b/pymare/utils.py
@@ -19,3 +19,16 @@ def _listify(obj):
     This provides a simple way to accept flexible arguments.
     """
     return obj if isinstance(obj, (list, tuple, type(None), np.ndarray)) else [obj]
+
+
+def _get_sig_code(p_value):
+    if p_value < 0.001:
+        return "***"
+    elif p_value < 0.01:
+        return "**"
+    elif p_value < 0.05:
+        return "*"
+    elif p_value < 0.1:
+        return "."
+    else:
+        return " "