Repository: arviz-devs/arviz

Issue:
Summary shows coordinates without dims as data_vars
**Describe the bug**
xarray allows to include coordinates without dimensions, which are treated incorrectly as variables by summary when converting to dataframe.

**To Reproduce**
```python
import arviz as az
import pymc3 as pm
import numpy as np

vals=np.array([[1,2,3], [4,5,6]])
with pm.Model() as m:
    sd = pm.HalfNormal("sd", sigma=10)
    mu = pm.Normal("mu", mu=0, sd=sd, shape=(2, 1))
    y = pm.Normal("y", mu=mu, sd=1, observed=vals)
    trace=pm.sample(draws=100, tune=100)

data = az.from_pymc3(trace=trace)
print(az.summary(data))

data.posterior = data.posterior.squeeze()
print(az.summary(data))
```

which prints:

```
          mean     sd  hpd_3%  hpd_97%  ...  ess_sd  ess_bulk  ess_tail  r_hat
mu[0,0]  1.968  0.562   0.838    2.935  ...   240.0     247.0     282.0   1.01
mu[1,0]  4.941  0.512   3.962    5.843  ...   629.0     620.0     350.0   1.01
sd       6.480  3.559   1.925   13.606  ...   178.0     228.0     207.0   1.01

[3 rows x 11 columns]
           mean     sd  hpd_3%  hpd_97%  ...  ess_sd  ess_bulk  ess_tail  r_hat
mu[0]     1.968  0.562   0.838    2.935  ...   240.0     247.0     282.0   1.01
mu[1]     4.941  0.512   3.962    5.843  ...   629.0     620.0     350.0   1.01
mu_dim_1  0.000  0.000   0.000    0.000  ...     0.0       0.0       0.0   0.00
sd        6.480  3.559   1.925   13.606  ...   178.0     228.0     207.0   1.01

[4 rows x 11 columns]
```

**Expected behavior**
Both datasets yield the same summary dataframe.

**Additional context**
Working on development version ArviZ.


Discussion:

Patch:
diff --git a/.travis.yml b/.travis.yml
index 5c04265..b7ef992 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,6 +2,10 @@ dist: xenial
 
 language: python
 
+branches:
+  only:
+    - master
+
 matrix:
   include:
     - name: "Python 3.6 Unit Test"
diff --git a/arviz/plots/essplot.py b/arviz/plots/essplot.py
index d099b71..ce53c60 100644
--- a/arviz/plots/essplot.py
+++ b/arviz/plots/essplot.py
@@ -1,6 +1,7 @@
 """Plot quantile or local effective sample sizes."""
 import numpy as np
 import xarray as xr
+from scipy.stats import rankdata
 
 from ..data import convert_to_dataset
 from ..stats import ess
@@ -296,7 +297,7 @@ def plot_ess(
 
             values = data[var_name].sel(**selection).values.flatten()
             mask = idata.sample_stats[rug_kind].values.flatten()
-            values = np.argsort(values)[mask]
+            values = rankdata(values)[mask]
             rug_space = np.max(x) * rug_kwargs.pop("space")
             rug_x, rug_y = values / (len(mask) - 1), np.zeros_like(values) - rug_space
             ax_.plot(rug_x, rug_y, **rug_kwargs)
diff --git a/arviz/plots/forestplot.py b/arviz/plots/forestplot.py
index 85a1ba3..610f7b2 100644
--- a/arviz/plots/forestplot.py
+++ b/arviz/plots/forestplot.py
@@ -208,7 +208,10 @@ def plot_forest(
         idx += 1
 
     for ax in axes:
-        ax.grid(False)
+        if kind == "ridgeplot":
+            ax.grid(False)
+        else:
+            ax.grid(False, axis="y")
         # Remove ticklines on y-axes
         ax.tick_params(axis="y", left=False, right=False)
 
diff --git a/arviz/plots/mcseplot.py b/arviz/plots/mcseplot.py
index e45cde1..b93fca6 100644
--- a/arviz/plots/mcseplot.py
+++ b/arviz/plots/mcseplot.py
@@ -1,6 +1,7 @@
 """Plot quantile MC standard error."""
 import numpy as np
 import xarray as xr
+from scipy.stats import rankdata
 
 from ..data import convert_to_dataset
 from ..stats import mcse
@@ -194,7 +195,7 @@ def plot_mcse(
             rug_kwargs.setdefault("markersize", rug_kwargs.pop("ms", 2 * _markersize))
 
             mask = idata.sample_stats[rug_kind].values.flatten()
-            values = np.argsort(values)[mask]
+            values = rankdata(values)[mask]
             y_min, y_max = ax_.get_ylim()
             y_min = y_min if errorbar else 0
             rug_space = (y_max - y_min) * rug_kwargs.pop("space")
diff --git a/arviz/plots/ppcplot.py b/arviz/plots/ppcplot.py
index fe5a483..e0b1abf 100644
--- a/arviz/plots/ppcplot.py
+++ b/arviz/plots/ppcplot.py
@@ -36,6 +36,7 @@ def plot_ppc(
     animated=False,
     animation_kwargs=None,
     legend=True,
+    ax=None,
 ):
     """
     Plot for posterior predictive checks.
@@ -99,6 +100,8 @@ def plot_ppc(
         Keywords passed to `animation.FuncAnimation`.
     legend : bool
         Add legend to figure. By default True.
+    ax : axes
+        Matplotlib axes. Defaults to None.
 
     Returns
     -------
@@ -247,9 +250,22 @@ def plot_ppc(
         figsize, textsize, rows, cols
     )
 
-    fig, axes = _create_axes_grid(length_plotters, rows, cols, figsize=figsize)
+    if ax is None:
+        fig, axes = _create_axes_grid(length_plotters, rows, cols, figsize=figsize)
+    else:
+        axes = np.ravel(ax)
+        if len(axes) != length_plotters:
+            raise ValueError(
+                "Found {} variables to plot but {} axes instances. They must be equal.".format(
+                    length_plotters, len(axes)
+                )
+            )
+        if animated:
+            fig = axes[0].get_figure()
+            if not all([ax.get_figure() is fig for ax in axes]):
+                raise ValueError("All axes must be on the same figure for animation to work")
 
-    for i, ax in enumerate(axes):
+    for i, ax_i in enumerate(axes):
         var_name, selection, obs_vals = obs_plotters[i]
         pp_var_name, _, pp_vals = pp_plotters[i]
         dtype = posterior_predictive[pp_var_name].dtype.kind
@@ -263,7 +279,7 @@ def plot_ppc(
             plot_kwargs = {"color": "C5", "alpha": alpha, "linewidth": 0.5 * linewidth}
             if dtype == "i":
                 plot_kwargs["drawstyle"] = "steps-pre"
-            ax.plot([], color="C5", label="Posterior predictive {}".format(pp_var_name))
+            ax_i.plot([], color="C5", label="Posterior predictive {}".format(pp_var_name))
 
             if dtype == "f":
                 plot_kde(
@@ -271,14 +287,14 @@ def plot_ppc(
                     label="Observed {}".format(var_name),
                     plot_kwargs={"color": "k", "linewidth": linewidth, "zorder": 3},
                     fill_kwargs={"alpha": 0},
-                    ax=ax,
+                    ax=ax_i,
                     legend=legend,
                 )
             else:
                 bins = get_bins(obs_vals)
                 hist, bin_edges = np.histogram(obs_vals, bins=bins, density=True)
                 hist = np.concatenate((hist[:1], hist))
-                ax.plot(
+                ax_i.plot(
                     bin_edges,
                     hist,
                     label="Observed {}".format(var_name),
@@ -290,7 +306,7 @@ def plot_ppc(
 
             if animated:
                 animate, init = _set_animation(
-                    pp_sampled_vals, ax, dtype=dtype, kind=kind, plot_kwargs=plot_kwargs
+                    pp_sampled_vals, ax_i, dtype=dtype, kind=kind, plot_kwargs=plot_kwargs
                 )
 
             else:
@@ -308,7 +324,7 @@ def plot_ppc(
                         hist = np.concatenate((hist[:1], hist))
                         pp_densities.extend([bin_edges, hist])
 
-                ax.plot(*pp_densities, **plot_kwargs)
+                ax_i.plot(*pp_densities, **plot_kwargs)
 
             if mean:
                 if dtype == "f":
@@ -321,7 +337,7 @@ def plot_ppc(
                             "zorder": 2,
                         },
                         label="Posterior predictive mean {}".format(pp_var_name),
-                        ax=ax,
+                        ax=ax_i,
                         legend=legend,
                     )
                 else:
@@ -329,7 +345,7 @@ def plot_ppc(
                     bins = get_bins(vals)
                     hist, bin_edges = np.histogram(vals, bins=bins, density=True)
                     hist = np.concatenate((hist[:1], hist))
-                    ax.plot(
+                    ax_i.plot(
                         bin_edges,
                         hist,
                         color="C0",
@@ -339,12 +355,12 @@ def plot_ppc(
                         linestyle="--",
                         drawstyle=plot_kwargs["drawstyle"],
                     )
-            ax.tick_params(labelsize=xt_labelsize)
-            ax.set_yticks([])
+            ax_i.tick_params(labelsize=xt_labelsize)
+            ax_i.set_yticks([])
 
         elif kind == "cumulative":
             drawstyle = "default" if dtype == "f" else "steps-pre"
-            ax.plot(
+            ax_i.plot(
                 *_empirical_cdf(obs_vals),
                 color="k",
                 linewidth=linewidth,
@@ -355,7 +371,7 @@ def plot_ppc(
             if animated:
                 animate, init = _set_animation(
                     pp_sampled_vals,
-                    ax,
+                    ax_i,
                     kind=kind,
                     alpha=alpha,
                     drawstyle=drawstyle,
@@ -370,12 +386,12 @@ def plot_ppc(
                     pp_x, pp_density = _empirical_cdf(vals)
                     pp_densities.extend([pp_x, pp_density])
 
-                ax.plot(
+                ax_i.plot(
                     *pp_densities, alpha=alpha, color="C5", drawstyle=drawstyle, linewidth=linewidth
                 )
-            ax.plot([], color="C5", label="Posterior predictive {}".format(pp_var_name))
+            ax_i.plot([], color="C5", label="Posterior predictive {}".format(pp_var_name))
             if mean:
-                ax.plot(
+                ax_i.plot(
                     *_empirical_cdf(pp_vals.flatten()),
                     color="C0",
                     linestyle="--",
@@ -383,7 +399,7 @@ def plot_ppc(
                     drawstyle=drawstyle,
                     label="Posterior predictive mean {}".format(pp_var_name)
                 )
-            ax.set_yticks([0, 0.5, 1])
+            ax_i.set_yticks([0, 0.5, 1])
 
         elif kind == "scatter":
             if mean:
@@ -397,7 +413,7 @@ def plot_ppc(
                             "zorder": 3,
                         },
                         label="Posterior predictive mean {}".format(pp_var_name),
-                        ax=ax,
+                        ax=ax_i,
                         legend=legend,
                     )
                 else:
@@ -405,7 +421,7 @@ def plot_ppc(
                     bins = get_bins(vals)
                     hist, bin_edges = np.histogram(vals, bins=bins, density=True)
                     hist = np.concatenate((hist[:1], hist))
-                    ax.plot(
+                    ax_i.plot(
                         bin_edges,
                         hist,
                         color="C0",
@@ -416,7 +432,7 @@ def plot_ppc(
                         drawstyle="steps-pre",
                     )
 
-            _, limit = ax.get_ylim()
+            _, limit = ax_i.get_ylim()
             limit *= 1.05
             y_rows = np.linspace(0, limit, num_pp_samples + 1)
             jitter_scale = y_rows[1] - y_rows[0]
@@ -426,7 +442,7 @@ def plot_ppc(
             obs_yvals = np.zeros_like(obs_vals, dtype=np.float64)
             if jitter:
                 obs_yvals += np.random.uniform(low=scale_low, high=scale_high, size=len(obs_vals))
-            ax.plot(
+            ax_i.plot(
                 obs_vals,
                 obs_yvals,
                 "o",
@@ -440,7 +456,7 @@ def plot_ppc(
             if animated:
                 animate, init = _set_animation(
                     pp_sampled_vals,
-                    ax,
+                    ax_i,
                     kind=kind,
                     height=y_rows.mean() * 0.5,
                     markersize=markersize,
@@ -452,25 +468,25 @@ def plot_ppc(
                     yvals = np.full_like(vals, y, dtype=np.float64)
                     if jitter:
                         yvals += np.random.uniform(low=scale_low, high=scale_high, size=len(vals))
-                    ax.plot(
+                    ax_i.plot(
                         vals, yvals, "o", zorder=2, color="C5", markersize=markersize, alpha=alpha
                     )
 
-            ax.plot([], "C5o", label="Posterior predictive {}".format(pp_var_name))
+            ax_i.plot([], "C5o", label="Posterior predictive {}".format(pp_var_name))
 
-            ax.set_yticks([])
+            ax_i.set_yticks([])
 
         if var_name != pp_var_name:
             xlabel = "{} / {}".format(var_name, pp_var_name)
         else:
             xlabel = var_name
-        ax.set_xlabel(make_label(xlabel, selection), fontsize=ax_labelsize)
+        ax_i.set_xlabel(make_label(xlabel, selection), fontsize=ax_labelsize)
 
         if legend:
             if i == 0:
-                ax.legend(fontsize=xt_labelsize * 0.75)
+                ax_i.legend(fontsize=xt_labelsize * 0.75)
             else:
-                ax.legend([])
+                ax_i.legend([])
 
     if animated:
         ani = animation.FuncAnimation(
diff --git a/arviz/stats/__init__.py b/arviz/stats/__init__.py
index cea761f..27b69a9 100644
--- a/arviz/stats/__init__.py
+++ b/arviz/stats/__init__.py
@@ -16,7 +16,6 @@ __all__ = [
     "r2_score",
     "summary",
     "waic",
-    "effective_sample_size",
     "ELPDData",
     "ess",
     "rhat",
diff --git a/arviz/stats/diagnostics.py b/arviz/stats/diagnostics.py
index 70b08d2..054e588 100644
--- a/arviz/stats/diagnostics.py
+++ b/arviz/stats/diagnostics.py
@@ -18,7 +18,7 @@ from .stats_utils import (
 from ..data import convert_to_dataset
 from ..utils import _var_names, conditional_jit, conditional_vect, Numba, _numba_var
 
-__all__ = ["bfmi", "effective_sample_size", "ess", "rhat", "mcse", "geweke"]
+__all__ = ["bfmi", "ess", "rhat", "mcse", "geweke"]
 
 
 def bfmi(data):
@@ -62,72 +62,6 @@ def bfmi(data):
     return _bfmi(dataset.energy)
 
 
-def effective_sample_size(data, *, var_names=None, method="bulk", relative=False, prob=None):
-    r"""Calculate estimate of the effective sample size.
-
-    Function deprecated. Use `arviz.ess`.
-
-    Parameters
-    ----------
-    data : obj
-        Any object that can be converted to an az.InferenceData object.
-        Refer to documentation of az.convert_to_dataset for details.
-        For ndarray: shape = (chain, draw).
-        For n-dimensional ndarray transform first to dataset with az.convert_to_dataset.
-    var_names : list
-        Names of variables to include in the effective_sample_size_mean report
-    method : str
-        Select ess method. Valid methods are:
-
-        - "bulk"
-        - "tail"     # prob, optional
-        - "quantile" # prob
-        - "mean" (old ess)
-        - "sd"
-        - "median"
-        - "mad" (mean absolute deviance)
-        - "z_scale"
-        - "folded"
-        - "identity"
-
-    relative : bool
-        Return relative ess
-        `ress = ess / N`
-    prob : float, optional
-        probability value for "tail" and "quantile" ess functions.
-
-    Returns
-    -------
-    xarray.Dataset
-        Return the effective sample size for mean, :math:`\hat{N}_{eff}`
-
-    Notes
-    -----
-    The basic ess diagnostic is computed by:
-
-    .. math:: \hat{N}_{eff} = \frac{MN}{\hat{\tau}}
-    .. math:: \hat{\tau} = -1 + 2 \sum_{t'=0}^K \hat{P}_t'
-
-    where :math:`\hat{\rho}_t` is the estimated _autocorrelation at lag t, and T
-    is the first odd positive integer for which the sum
-    :math:`\hat{\rho}_{T+1} + \hat{\rho}_{T+1}` is negative.
-
-    The current implementation is similar to Stan, which uses Geyer's initial monotone sequence
-    criterion (Geyer, 1992; Geyer, 2011).
-
-    References
-    ----------
-    * Vehtari et al. (2019) see https://arxiv.org/abs/1903.08008
-    * https://mc-stan.org/docs/2_18/reference-manual/effective-sample-size-section.html
-      Section 15.4.2
-    * Gelman et al. BDA (2014) Formula 11.8
-    """
-    warnings.warn(
-        "Function `arviz.effective_sample_size` is deprecated. Use `arviz.ess`", DeprecationWarning
-    )
-    return ess(data, var_names=var_names, method=method, relative=relative, prob=prob)
-
-
 def ess(data, *, var_names=None, method="bulk", relative=False, prob=None):
     r"""Calculate estimate of the effective sample size.
 
diff --git a/arviz/stats/stats.py b/arviz/stats/stats.py
index 95d23b6..ef54c9e 100644
--- a/arviz/stats/stats.py
+++ b/arviz/stats/stats.py
@@ -398,6 +398,9 @@ def loo(data, pointwise=False, reff=None, scale="deviance"):
         - `log` : 1 * log-score (after Vehtari et al. (2017))
         - `negative_log` : -1 * (log-score)
 
+        A higher log-score (or a lower deviance) indicates a model with better predictive
+        accuracy.
+
     Returns
     -------
     pandas.Series with the following rows:
@@ -976,7 +979,9 @@ def summary(
         )
     metrics.extend(extra_metrics)
     metric_names.extend(extra_metric_names)
-    joined = xr.concat(metrics, dim="metric").assign_coords(metric=metric_names)
+    joined = (
+        xr.concat(metrics, dim="metric").assign_coords(metric=metric_names).reset_coords(drop=True)
+    )
 
     if fmt.lower() == "wide":
         dfs = []
@@ -1045,6 +1050,9 @@ def waic(data, pointwise=False, scale="deviance"):
         - `log` : 1 * log-score
         - `negative_log` : -1 * (log-score)
 
+        A higher log-score (or a lower deviance) indicates a model with better predictive
+        accuracy.
+
     Returns
     -------
     Series with the following rows:
@@ -1112,10 +1120,11 @@ def waic(data, pointwise=False, scale="deviance"):
     warn_mg = False
     if np.any(vars_lpd > 0.4):
         warnings.warn(
-            """For one or more samples the posterior variance of the log predictive
-        densities exceeds 0.4. This could be indication of WAIC starting to fail see
-        http://arxiv.org/abs/1507.04544 for details
-        """
+            (
+                "For one or more samples the posterior variance of the log predictive "
+                "densities exceeds 0.4. This could be indication of WAIC starting to fail. \n"
+                "See http://arxiv.org/abs/1507.04544 for details"
+            )
         )
         warn_mg = True
 
diff --git a/doc/api.rst b/doc/api.rst
index 46e6ab4..92a965c 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -64,7 +64,6 @@ Diagnostics
     bfmi
     geweke
     ess
-    effective_sample_size
     rhat
     mcse
 


Test patch:
diff --git a/arviz/tests/helpers.py b/arviz/tests/helpers.py
index 812b421..7b0be01 100644
--- a/arviz/tests/helpers.py
+++ b/arviz/tests/helpers.py
@@ -9,6 +9,9 @@ import logging
 import pytest
 import numpy as np
 
+from ..data import from_dict
+
+
 _log = logging.getLogger(__name__)
 
 
@@ -34,6 +37,126 @@ def chains():
     return 2
 
 
+def create_model(seed=10):
+    """Create model with fake data."""
+    np.random.seed(seed)
+    nchains = 4
+    ndraws = 500
+    data = {
+        "J": 8,
+        "y": np.array([28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0]),
+        "sigma": np.array([15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0]),
+    }
+    posterior = {
+        "mu": np.random.randn(nchains, ndraws),
+        "tau": abs(np.random.randn(nchains, ndraws)),
+        "eta": np.random.randn(nchains, ndraws, data["J"]),
+        "theta": np.random.randn(nchains, ndraws, data["J"]),
+    }
+    posterior_predictive = {"y": np.random.randn(nchains, ndraws, len(data["y"]))}
+    sample_stats = {
+        "energy": np.random.randn(nchains, ndraws),
+        "diverging": np.random.randn(nchains, ndraws) > 0.90,
+        "max_depth": np.random.randn(nchains, ndraws) > 0.90,
+        "log_likelihood": np.random.randn(nchains, ndraws, data["J"]),
+    }
+    prior = {
+        "mu": np.random.randn(nchains, ndraws) / 2,
+        "tau": abs(np.random.randn(nchains, ndraws)) / 2,
+        "eta": np.random.randn(nchains, ndraws, data["J"]) / 2,
+        "theta": np.random.randn(nchains, ndraws, data["J"]) / 2,
+    }
+    prior_predictive = {"y": np.random.randn(nchains, ndraws, len(data["y"])) / 2}
+    sample_stats_prior = {
+        "energy": np.random.randn(nchains, ndraws),
+        "diverging": (np.random.randn(nchains, ndraws) > 0.95).astype(int),
+    }
+    model = from_dict(
+        posterior=posterior,
+        posterior_predictive=posterior_predictive,
+        sample_stats=sample_stats,
+        prior=prior,
+        prior_predictive=prior_predictive,
+        sample_stats_prior=sample_stats_prior,
+        observed_data={"y": data["y"]},
+        dims={"y": ["obs_dim"], "log_likelihood": ["obs_dim"]},
+        coords={"obs_dim": range(data["J"])},
+    )
+    return model
+
+
+def create_multidimensional_model(seed=10):
+    """Create model with fake data."""
+    np.random.seed(seed)
+    nchains = 4
+    ndraws = 500
+    ndim1 = 5
+    ndim2 = 7
+    data = {
+        "y": np.random.normal(size=(ndim1, ndim2)),
+        "sigma": np.random.normal(size=(ndim1, ndim2)),
+    }
+    posterior = {
+        "mu": np.random.randn(nchains, ndraws),
+        "tau": abs(np.random.randn(nchains, ndraws)),
+        "eta": np.random.randn(nchains, ndraws, ndim1, ndim2),
+        "theta": np.random.randn(nchains, ndraws, ndim1, ndim2),
+    }
+    posterior_predictive = {"y": np.random.randn(nchains, ndraws, ndim1, ndim2)}
+    sample_stats = {
+        "energy": np.random.randn(nchains, ndraws),
+        "diverging": np.random.randn(nchains, ndraws) > 0.90,
+        "log_likelihood": np.random.randn(nchains, ndraws, ndim1, ndim2),
+    }
+    prior = {
+        "mu": np.random.randn(nchains, ndraws) / 2,
+        "tau": abs(np.random.randn(nchains, ndraws)) / 2,
+        "eta": np.random.randn(nchains, ndraws, ndim1, ndim2) / 2,
+        "theta": np.random.randn(nchains, ndraws, ndim1, ndim2) / 2,
+    }
+    prior_predictive = {"y": np.random.randn(nchains, ndraws, ndim1, ndim2) / 2}
+    sample_stats_prior = {
+        "energy": np.random.randn(nchains, ndraws),
+        "diverging": (np.random.randn(nchains, ndraws) > 0.95).astype(int),
+    }
+    model = from_dict(
+        posterior=posterior,
+        posterior_predictive=posterior_predictive,
+        sample_stats=sample_stats,
+        prior=prior,
+        prior_predictive=prior_predictive,
+        sample_stats_prior=sample_stats_prior,
+        observed_data={"y": data["y"]},
+        dims={"y": ["dim1", "dim2"], "log_likelihood": ["dim1", "dim2"]},
+        coords={"dim1": range(ndim1), "dim2": range(ndim2)},
+    )
+    return model
+
+
+@pytest.fixture(scope="module")
+def models():
+    """Fixture containing 2 mock inference data instances for testing."""
+    # blank line to keep black and pydocstyle happy
+
+    class Models:
+        model_1 = create_model(seed=10)
+        model_2 = create_model(seed=11)
+
+    return Models()
+
+
+@pytest.fixture(scope="module")
+def multidim_models():
+    """Fixture containing 2 mock inference data instances with multidimensional data for testing."""
+    # blank line to keep black and pydocstyle happy
+
+    class Models:
+        model_1 = create_multidimensional_model(seed=10)
+        model_2 = create_multidimensional_model(seed=11)
+
+    return Models()
+
+
 def check_multiple_attrs(test_dict, parent):
     """Perform multiple hasattr checks on InferenceData objects.
 
diff --git a/arviz/tests/test_diagnostics.py b/arviz/tests/test_diagnostics.py
index 3dc59e6..7851b2e 100644
--- a/arviz/tests/test_diagnostics.py
+++ b/arviz/tests/test_diagnostics.py
@@ -9,7 +9,7 @@ import pytest
 
 from ..data import load_arviz_data, from_cmdstan
 from ..plots.plot_utils import xarray_var_iter
-from ..stats import bfmi, rhat, ess, mcse, geweke, effective_sample_size
+from ..stats import bfmi, rhat, ess, mcse, geweke
 from ..stats.diagnostics import (
     ks_summary,
     _ess,
@@ -273,11 +273,6 @@ class TestDiagnostics:
         assert ess_hat > n_low
         assert ess_hat < n_high
 
-    def test_old_effective_sample_size(self):
-        with pytest.deprecated_call():
-            old_ess = effective_sample_size(np.random.randn(4, 100), method="bulk")
-        assert old_ess
-
     @pytest.mark.parametrize(
         "method",
         (
diff --git a/arviz/tests/test_plots.py b/arviz/tests/test_plots.py
index d0148e6..bb3c2a5 100644
--- a/arviz/tests/test_plots.py
+++ b/arviz/tests/test_plots.py
@@ -9,7 +9,13 @@ import pytest
 
 from ..data import from_dict, load_arviz_data
 from ..stats import compare, loo, waic
-from .helpers import eight_schools_params  # pylint: disable=unused-import
+from .helpers import (  # pylint: disable=unused-import
+    eight_schools_params,
+    models,
+    create_model,
+    multidim_models,
+    create_multidimensional_model,
+)
 from ..plots import (
     plot_density,
     plot_trace,
@@ -39,120 +45,6 @@ np.random.seed(0)
 os.environ["ARVIZ_LOAD"] = "EAGER"
 
 
-def create_model(seed=10):
-    """Create model with fake data."""
-    np.random.seed(seed)
-    nchains = 4
-    ndraws = 500
-    data = {
-        "J": 8,
-        "y": np.array([28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0]),
-        "sigma": np.array([15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0]),
-    }
-    posterior = {
-        "mu": np.random.randn(nchains, ndraws),
-        "tau": abs(np.random.randn(nchains, ndraws)),
-        "eta": np.random.randn(nchains, ndraws, data["J"]),
-        "theta": np.random.randn(nchains, ndraws, data["J"]),
-    }
-    posterior_predictive = {"y": np.random.randn(nchains, ndraws, len(data["y"]))}
-    sample_stats = {
-        "energy": np.random.randn(nchains, ndraws),
-        "diverging": np.random.randn(nchains, ndraws) > 0.90,
-        "max_depth": np.random.randn(nchains, ndraws) > 0.90,
-        "log_likelihood": np.random.randn(nchains, ndraws, data["J"]),
-    }
-    prior = {
-        "mu": np.random.randn(nchains, ndraws) / 2,
-        "tau": abs(np.random.randn(nchains, ndraws)) / 2,
-        "eta": np.random.randn(nchains, ndraws, data["J"]) / 2,
-        "theta": np.random.randn(nchains, ndraws, data["J"]) / 2,
-    }
-    prior_predictive = {"y": np.random.randn(nchains, ndraws, len(data["y"])) / 2}
-    sample_stats_prior = {
-        "energy": np.random.randn(nchains, ndraws),
-        "diverging": (np.random.randn(nchains, ndraws) > 0.95).astype(int),
-    }
-    model = from_dict(
-        posterior=posterior,
-        posterior_predictive=posterior_predictive,
-        sample_stats=sample_stats,
-        prior=prior,
-        prior_predictive=prior_predictive,
-        sample_stats_prior=sample_stats_prior,
-        observed_data={"y": data["y"]},
-        dims={"y": ["obs_dim"], "log_likelihood": ["obs_dim"]},
-        coords={"obs_dim": range(data["J"])},
-    )
-    return model
-
-
-def create_multidimensional_model(seed=10):
-    """Create model with fake data."""
-    np.random.seed(seed)
-    nchains = 4
-    ndraws = 500
-    ndim1 = 5
-    ndim2 = 7
-    data = {
-        "y": np.random.normal(size=(ndim1, ndim2)),
-        "sigma": np.random.normal(size=(ndim1, ndim2)),
-    }
-    posterior = {
-        "mu": np.random.randn(nchains, ndraws),
-        "tau": abs(np.random.randn(nchains, ndraws)),
-        "eta": np.random.randn(nchains, ndraws, ndim1, ndim2),
-        "theta": np.random.randn(nchains, ndraws, ndim1, ndim2),
-    }
-    posterior_predictive = {"y": np.random.randn(nchains, ndraws, ndim1, ndim2)}
-    sample_stats = {
-        "energy": np.random.randn(nchains, ndraws),
-        "diverging": np.random.randn(nchains, ndraws) > 0.90,
-        "log_likelihood": np.random.randn(nchains, ndraws, ndim1, ndim2),
-    }
-    prior = {
-        "mu": np.random.randn(nchains, ndraws) / 2,
-        "tau": abs(np.random.randn(nchains, ndraws)) / 2,
-        "eta": np.random.randn(nchains, ndraws, ndim1, ndim2) / 2,
-        "theta": np.random.randn(nchains, ndraws, ndim1, ndim2) / 2,
-    }
-    prior_predictive = {"y": np.random.randn(nchains, ndraws, ndim1, ndim2) / 2}
-    sample_stats_prior = {
-        "energy": np.random.randn(nchains, ndraws),
-        "diverging": (np.random.randn(nchains, ndraws) > 0.95).astype(int),
-    }
-    model = from_dict(
-        posterior=posterior,
-        posterior_predictive=posterior_predictive,
-        sample_stats=sample_stats,
-        prior=prior,
-        prior_predictive=prior_predictive,
-        sample_stats_prior=sample_stats_prior,
-        observed_data={"y": data["y"]},
-        dims={"y": ["dim1", "dim2"], "log_likelihood": ["dim1", "dim2"]},
-        coords={"dim1": range(ndim1), "dim2": range(ndim2)},
-    )
-    return model
-
-
-@pytest.fixture(scope="module")
-def models():
-    class Models:
-        model_1 = create_model(seed=10)
-        model_2 = create_model(seed=11)
-
-    return Models()
-
-
-@pytest.fixture(scope="module")
-def multidim_models():
-    class Models:
-        model_1 = create_multidimensional_model(seed=10)
-        model_2 = create_multidimensional_model(seed=11)
-
-    return Models()
-
-
 @pytest.fixture(scope="function", autouse=True)
 def clean_plots(request, save_figs):
     """Close plots after each test, optionally save if --save is specified during test invocation"""
@@ -649,6 +541,25 @@ def test_plot_ppc_bad(models, kind):
         plot_ppc(models.model_1, num_pp_samples="bad_val")
 
 
+@pytest.mark.parametrize("kind", ["density", "cumulative", "scatter"])
+def test_plot_ppc_ax(models, kind, fig_ax):
+    """Test ax argument of plot_ppc."""
+    _, ax = fig_ax
+    axes = plot_ppc(models.model_1, kind=kind, ax=ax)
+    assert axes[0] is ax
+
+
+def test_plot_ppc_bad_ax(models, fig_ax):
+    _, ax = fig_ax
+    _, ax2 = plt.subplots(1, 2)
+    with pytest.raises(ValueError, match="same figure"):
+        plot_ppc(
+            models.model_1, ax=[ax, *ax2], flatten=[], coords={"obs_dim": [1, 2, 3]}, animated=True
+        )
+    with pytest.raises(ValueError, match="2 axes"):
+        plot_ppc(models.model_1, ax=ax2)
+
+
 @pytest.mark.parametrize("var_names", (None, "mu", ["mu", "tau"]))
 def test_plot_violin(models, var_names):
     axes = plot_violin(models.model_1, var_names=var_names)
diff --git a/arviz/tests/test_stats.py b/arviz/tests/test_stats.py
index 397c2e4..1f6c9f2 100644
--- a/arviz/tests/test_stats.py
+++ b/arviz/tests/test_stats.py
@@ -1,4 +1,4 @@
-# pylint: disable=redefined-outer-name
+# pylint: disable=redefined-outer-name, no-member
 import os
 from copy import deepcopy
 import numpy as np
@@ -23,7 +23,7 @@ from ..stats import (
 )
 from ..stats.stats import _gpinv
 from ..utils import Numba
-from .helpers import check_multiple_attrs
+from .helpers import check_multiple_attrs, multidim_models  # pylint: disable=unused-import
 
 os.environ["ARVIZ_LOAD"] = "EAGER"
 
@@ -69,8 +69,12 @@ def test_r2_score_multivariate():
 
 
 @pytest.mark.parametrize("method", ["stacking", "BB-pseudo-BMA", "pseudo-BMA"])
-def test_compare_same(centered_eight, method):
-    data_dict = {"first": centered_eight, "second": centered_eight}
+@pytest.mark.parametrize("multidim", [True, False])
+def test_compare_same(centered_eight, multidim_models, method, multidim):
+    if multidim:
+        data_dict = {"first": multidim_models.model_1, "second": multidim_models.model_1}
+    else:
+        data_dict = {"first": centered_eight, "second": centered_eight}
 
     weight = compare(data_dict, method=method)["weight"]
     assert_allclose(weight[0], weight[1])
@@ -95,6 +99,17 @@ def test_compare_different(centered_eight, non_centered_eight, ic, method, scale
     assert_allclose(np.sum(weight), 1.0)
 
 
+@pytest.mark.parametrize("ic", ["waic", "loo"])
+@pytest.mark.parametrize("method", ["stacking", "BB-pseudo-BMA", "pseudo-BMA"])
+def test_compare_different_multidim(multidim_models, ic, method):
+    model_dict = {"model_1": multidim_models.model_1, "model_2": multidim_models.model_2}
+    weight = compare(model_dict, ic=ic, method=method)["weight"]
+
+    # this should hold because the same seed is always used
+    assert weight["model_1"] >= weight["model_2"]
+    assert_allclose(np.sum(weight), 1.0)
+
+
 def test_compare_different_size(centered_eight, non_centered_eight):
     centered_eight = deepcopy(centered_eight)
     centered_eight.posterior = centered_eight.posterior.drop("Choate", "school")
@@ -110,10 +125,9 @@ def test_compare_different_size(centered_eight, non_centered_eight):
 
 
 @pytest.mark.parametrize("var_names_expected", ((None, 10), ("mu", 1), (["mu", "tau"], 2)))
-def test_summary_var_names(var_names_expected):
+def test_summary_var_names(centered_eight, var_names_expected):
     var_names, expected = var_names_expected
-    centered = load_arviz_data("centered_eight")
-    summary_df = summary(centered, var_names=var_names)
+    summary_df = summary(centered_eight, var_names=var_names)
     assert len(summary_df.index) == expected
 
 
@@ -202,10 +216,17 @@ def test_summary_bad_unpack_order(centered_eight, order):
 
 
 @pytest.mark.parametrize("scale", ["deviance", "log", "negative_log"])
-def test_waic(centered_eight, scale):
+@pytest.mark.parametrize("multidim", (True, False))
+def test_waic(centered_eight, multidim_models, scale, multidim):
     """Test widely available information criterion calculation"""
-    assert waic(centered_eight, scale=scale) is not None
-    assert waic(centered_eight, pointwise=True, scale=scale) is not None
+    if multidim:
+        assert waic(multidim_models.model_1, scale=scale) is not None
+        waic_pointwise = waic(multidim_models.model_1, pointwise=True, scale=scale)
+    else:
+        assert waic(centered_eight, scale=scale) is not None
+        waic_pointwise = waic(centered_eight, pointwise=True, scale=scale)
+    assert waic_pointwise is not None
+    assert "waic_i" in waic_pointwise
 
 
 def test_waic_bad(centered_eight):
@@ -246,8 +267,20 @@ def test_waic_print(centered_eight, scale):
     assert waic_data == waic_pointwise
 
 
-def test_loo(centered_eight):
-    assert loo(centered_eight) is not None
+@pytest.mark.parametrize("scale", ["deviance", "log", "negative_log"])
+@pytest.mark.parametrize("multidim", (True, False))
+def test_loo(centered_eight, multidim_models, scale, multidim):
+    """Test approximate leave one out criterion calculation"""
+    if multidim:
+        assert loo(multidim_models.model_1, scale=scale) is not None
+        loo_pointwise = loo(multidim_models.model_1, pointwise=True, scale=scale)
+    else:
+        assert loo(centered_eight, scale=scale) is not None
+        loo_pointwise = loo(centered_eight, pointwise=True, scale=scale)
+    assert loo_pointwise is not None
+    assert "loo_i" in loo_pointwise
+    assert "pareto_k" in loo_pointwise
+    assert "loo_scale" in loo_pointwise
 
 
 def test_loo_one_chain(centered_eight):
@@ -257,16 +290,6 @@ def test_loo_one_chain(centered_eight):
     assert loo(centered_eight) is not None
 
 
-@pytest.mark.parametrize("scale", ["deviance", "log", "negative_log"])
-def test_loo_pointwise(centered_eight, scale):
-    """Test pointwise loo with different scales."""
-    loo_results = loo(centered_eight, scale=scale, pointwise=True)
-    assert loo_results is not None
-    assert hasattr(loo_results, "loo_scale")
-    assert hasattr(loo_results, "pareto_k")
-    assert hasattr(loo_results, "loo_i")
-
-
 def test_loo_bad(centered_eight):
     with pytest.raises(TypeError):
         loo(np.random.randn(2, 10))
@@ -287,12 +310,16 @@ def test_loo_warning(centered_eight):
     centered_eight = deepcopy(centered_eight)
     # make one of the khats infinity
     centered_eight.sample_stats["log_likelihood"][:, :, 1] = 10
-    with pytest.warns(UserWarning):
+    with pytest.warns(UserWarning) as record:
         assert loo(centered_eight, pointwise=True) is not None
+    assert len(record) == 1
+    assert "Estimated shape parameter" in str(record[0].message)
     # make all of the khats infinity
-    centered_eight.sample_stats["log_likelihood"][:, :, :] = 0
-    with pytest.warns(UserWarning):
+    centered_eight.sample_stats["log_likelihood"][:, :, :] = 1
+    with pytest.warns(UserWarning) as record:
         assert loo(centered_eight, pointwise=True) is not None
+    assert len(record) == 1
+    assert "Estimated shape parameter" in str(record[0].message)
 
 
 @pytest.mark.parametrize("scale", ["deviance", "log", "negative_log"])
@@ -305,10 +332,9 @@ def test_loo_print(centered_eight, scale):
     assert loo_data == loo_pointwise[: len(loo_data)]
 
 
-def test_psislw():
-    data = load_arviz_data("centered_eight")
-    pareto_k = loo(data, pointwise=True, reff=0.7)["pareto_k"]
-    log_likelihood = data.sample_stats.log_likelihood  # pylint: disable=no-member
+def test_psislw(centered_eight):
+    pareto_k = loo(centered_eight, pointwise=True, reff=0.7)["pareto_k"]
+    log_likelihood = centered_eight.sample_stats.log_likelihood  # pylint: disable=no-member
     log_likelihood = log_likelihood.stack(samples=("chain", "draw"))
     assert_allclose(pareto_k, psislw(-log_likelihood, 0.7)[1])
 
@@ -388,6 +414,46 @@ def test_loo_pit(centered_eight, args):
     assert np.all((loo_pit_data >= 0) & (loo_pit_data <= 1))
 
 
+@pytest.mark.parametrize(
+    "args",
+    [
+        {"y": "y"},
+        {"y": "y", "y_hat": "y"},
+        {"y": "arr", "y_hat": "y"},
+        {"y": "y", "y_hat": "arr"},
+        {"y": "arr", "y_hat": "arr"},
+        {"y": "y", "y_hat": "y", "log_weights": "arr"},
+        {"y": "arr", "y_hat": "y", "log_weights": "arr"},
+        {"y": "y", "y_hat": "arr", "log_weights": "arr"},
+        {"idata": False},
+    ],
+)
+def test_loo_pit_multidim(multidim_models, args):
+    y = args.get("y", None)
+    y_hat = args.get("y_hat", None)
+    log_weights = args.get("log_weights", None)
+    idata = multidim_models.model_1
+    y_arr = idata.observed_data.y
+    y_hat_arr = idata.posterior_predictive.y.stack(samples=("chain", "draw"))
+    log_like = idata.sample_stats.log_likelihood.stack(samples=("chain", "draw"))
+    n_samples = len(log_like.samples)
+    ess_p = ess(idata.posterior, method="mean")
+    reff = np.hstack([ess_p[v].values.flatten() for v in ess_p.data_vars]).mean() / n_samples
+    log_weights_arr = psislw(-log_like, reff=reff)[0]
+
+    if args.get("idata", True):
+        if y == "arr":
+            y = y_arr
+        if y_hat == "arr":
+            y_hat = y_hat_arr
+        if log_weights == "arr":
+            log_weights = log_weights_arr
+        loo_pit_data = loo_pit(idata=idata, y=y, y_hat=y_hat, log_weights=log_weights)
+    else:
+        loo_pit_data = loo_pit(idata=None, y=y_arr, y_hat=y_hat_arr, log_weights=log_weights_arr)
+    assert np.all((loo_pit_data >= 0) & (loo_pit_data <= 1))
+
+
 @pytest.mark.parametrize("input_type", ["idataarray", "idatanone_ystr", "yarr_yhatnone"])
 def test_loo_pit_bad_input(centered_eight, input_type):
     """Test incompatible input combinations."""