Repository: arviz-devs/arviz

Issue:
plot_compare needs a legend 
There is a lot going on in the `plot_compare` output. There should be a legend displayed by default, as most users will not be intimately familiar with Statistical Rethinking. There should also be a better description in the function docstring describing what gets plotted.

Discussion:
yannmclatchie: This looks like a good start, I particularly like the interpretable title and the future idea of a level-based guide through the workflow.

I wonder if both the x-axis label and legend are necessary if they repeat the fact that they show ELPD? The ELPD difference is still somewhat confusing to me, since it is shown on the same scale and axis as the raw ELPDs. If a user is looking at the credible interval of the ELPD difference, it is likely that they are looking to see if zero is included in that interval, which is not immediately obvious from this plot. 

What may be clearer is if we produce multiple separate plots showing:
1. The ELPD distributions of the ranked models with the same dashed line for the "best" model (the current default plot but without the differences) with no legend and the currently suggested x-axis label, and;
2. A similar plot of just the difference in ELPD distributions with another dashed line at zero (something of an equivalent to [`plot.vsel` in `projpred`](https://mc-stan.org/projpred/reference/plot.vsel.html) with `deltas=TRUE`) which is only shown alongside the raw ELPD plot if the user passes `plot_deltas = True` for instance, again without legend and with an appropriate x-axis label.

Let me know what you think :pray:
aloctavodia: @yannmclatchie thanks for the suggestions, I incorporated those related to the docstring.
Regarding the plot I would like to keep this plot as close as possible to the one used in Statistical Rethinking. But maybe we can add a second plot like the one you are suggesting. Even when the scale is on the "raw ELPD" values, the ELPD difference are computed with respect to the top-ranked model and thus the vertical dashed line can be used as the reference point.



cfonnesbeck: The "new default" in the first pane is a big improvement, so thanks! I'm not as enthusiastic about the second one--dual y-axes don't work all that well IMO. 

An even simpler version would be to remove the difference from the plot altogether and add it as a text annotation to the relevant point. 

... OR allow users to plot the differences *instead* of the raw ELPD scores, if that is what they are interested in. 
aloctavodia: Thanks for the feedback. I will work on refining the "new default". We already have a keyword for removing the differences. I agree we can add a rcParam.
aloctavodia: @OriolAbril can you provide help about the errors building the docs?
OriolAbril: Should be fixed now. The "See Also" section has a special syntax. elements must be valid sphinx references and one can optionally add a description to the link after a ` : `. If the description is too long to fit in one line it needs to be on the line below _and indented_. Ref: https://numpydoc.readthedocs.io/en/latest/format.html#see-also

Patch:
diff --git a/arviz/plots/backends/bokeh/compareplot.py b/arviz/plots/backends/bokeh/compareplot.py
index 45e09b5..5fa9604 100644
--- a/arviz/plots/backends/bokeh/compareplot.py
+++ b/arviz/plots/backends/bokeh/compareplot.py
@@ -1,5 +1,7 @@
 """Bokeh Compareplot."""
 from bokeh.models import Span
+from bokeh.models.annotations import Title, Legend
+
 
 from ...plot_utils import _scale_fig_size
 from .. import show_layout
@@ -9,6 +11,8 @@ from . import backend_kwarg_defaults, create_axes_grid
 def plot_compare(
     ax,
     comp_df,
+    legend,
+    title,
     figsize,
     plot_ic_diff,
     plot_standard_error,
@@ -43,8 +47,10 @@ def plot_compare(
 
     yticks_pos = list(yticks_pos)
 
+    labels = []
+
     if plot_ic_diff:
-        ax.yaxis.ticker = yticks_pos
+        ax.yaxis.ticker = yticks_pos[::2]
         ax.yaxis.major_label_overrides = {
             dtype(key): value
             for key, value in zip(yticks_pos, yticks_labels)
@@ -63,7 +69,7 @@ def plot_compare(
             err_ys.append((y, y))
 
         # plot them
-        ax.triangle(
+        dif_tri = ax.triangle(
             comp_df[information_criterion].iloc[1:],
             yticks_pos[1::2],
             line_color=plot_kwargs.get("color_dse", "grey"),
@@ -71,7 +77,9 @@ def plot_compare(
             line_width=2,
             size=6,
         )
-        ax.multi_line(err_xs, err_ys, line_color=plot_kwargs.get("color_dse", "grey"))
+        dif_line = ax.multi_line(err_xs, err_ys, line_color=plot_kwargs.get("color_dse", "grey"))
+
+        labels.append(("ELPD difference", [dif_tri, dif_line]))
 
     else:
         ax.yaxis.ticker = yticks_pos[::2]
@@ -79,7 +87,7 @@ def plot_compare(
             key: value for key, value in zip(yticks_pos[::2], yticks_labels)
         }
 
-    ax.circle(
+    elpd_circ = ax.circle(
         comp_df[information_criterion],
         yticks_pos[::2],
         line_color=plot_kwargs.get("color_ic", "black"),
@@ -87,6 +95,7 @@ def plot_compare(
         line_width=2,
         size=6,
     )
+    elpd_label = [elpd_circ]
 
     if plot_standard_error:
         # create the coordinates for the errorbars
@@ -98,18 +107,22 @@ def plot_compare(
             err_ys.append((y, y))
 
         # plot them
-        ax.multi_line(err_xs, err_ys, line_color=plot_kwargs.get("color_ic", "black"))
+        elpd_line = ax.multi_line(err_xs, err_ys, line_color=plot_kwargs.get("color_ic", "black"))
+        elpd_label.append(elpd_line)
+
+    labels.append(("ELPD", elpd_label))
+
+    scale = comp_df["scale"][0]
 
     if insample_dev:
-        scale = comp_df[f"{information_criterion}_scale"][0]
-        p_ic = comp_df[f"p_{information_criterion}"]
+        p_ic = comp_df[f"p_{information_criterion.split('_')[1]}"]
         if scale == "log":
             correction = p_ic
         elif scale == "negative_log":
             correction = -p_ic
         elif scale == "deviance":
             correction = -(2 * p_ic)
-        ax.circle(
+        insample_circ = ax.circle(
             comp_df[information_criterion] + correction,
             yticks_pos[::2],
             line_color=plot_kwargs.get("color_insample_dev", "black"),
@@ -117,6 +130,7 @@ def plot_compare(
             line_width=2,
             size=6,
         )
+        labels.append(("In-sample ELPD", [insample_circ]))
 
     vline = Span(
         location=comp_df[information_criterion].iloc[0],
@@ -128,12 +142,21 @@ def plot_compare(
 
     ax.renderers.append(vline)
 
-    scale_col = information_criterion + "_scale"
-    if scale_col in comp_df:
-        scale = comp_df[scale_col].iloc[0].capitalize()
-    else:
-        scale = "Deviance"
-    ax.xaxis.axis_label = scale
+    if legend:
+        legend = Legend(items=labels, orientation="vertical", location="top_right")
+        ax.add_layout(legend, "above")
+        ax.legend.click_policy = "hide"
+
+    if title:
+        _title = Title()
+        _title.text = f"Model comparison\n{'higher' if scale == 'log' else 'lower'} is better"
+        ax.title = _title
+
+    if scale == "negative_log":
+        scale = "-log"
+
+    ax.xaxis.axis_label = f"{information_criterion} ({scale})"
+    ax.yaxis.axis_label = "ranked models"
     ax.y_range._property_values["start"] = -1 + step  # pylint: disable=protected-access
     ax.y_range._property_values["end"] = 0 - step  # pylint: disable=protected-access
 
diff --git a/arviz/plots/backends/matplotlib/compareplot.py b/arviz/plots/backends/matplotlib/compareplot.py
index 7c92d36..ee7f6ac 100644
--- a/arviz/plots/backends/matplotlib/compareplot.py
+++ b/arviz/plots/backends/matplotlib/compareplot.py
@@ -8,6 +8,8 @@ from . import backend_kwarg_defaults, backend_show, create_axes_grid
 def plot_compare(
     ax,
     comp_df,
+    legend,
+    title,
     figsize,
     plot_ic_diff,
     plot_standard_error,
@@ -41,29 +43,15 @@ def plot_compare(
     if ax is None:
         _, ax = create_axes_grid(1, backend_kwargs=backend_kwargs)
 
-    if plot_ic_diff:
-        ax.set_yticks(yticks_pos)
-        ax.errorbar(
-            x=comp_df[information_criterion].iloc[1:],
-            y=yticks_pos[1::2],
-            xerr=comp_df.dse[1:],
-            color=plot_kwargs.get("color_dse", "grey"),
-            fmt=plot_kwargs.get("marker_dse", "^"),
-            mew=linewidth,
-            elinewidth=linewidth,
-        )
-
-    else:
-        ax.set_yticks(yticks_pos[::2])
-
     if plot_standard_error:
         ax.errorbar(
             x=comp_df[information_criterion],
             y=yticks_pos[::2],
             xerr=comp_df.se,
+            label="ELPD",
             color=plot_kwargs.get("color_ic", "k"),
             fmt=plot_kwargs.get("marker_ic", "o"),
-            mfc="None",
+            mfc=plot_kwargs.get("marker_fc", "white"),
             mew=linewidth,
             lw=linewidth,
         )
@@ -71,16 +59,35 @@ def plot_compare(
         ax.plot(
             comp_df[information_criterion],
             yticks_pos[::2],
+            label="ELPD",
             color=plot_kwargs.get("color_ic", "k"),
             marker=plot_kwargs.get("marker_ic", "o"),
-            mfc="None",
+            mfc=plot_kwargs.get("marker_fc", "white"),
             mew=linewidth,
             lw=0,
+            zorder=3,
         )
 
+    if plot_ic_diff:
+        ax.set_yticks(yticks_pos)
+        ax.errorbar(
+            x=comp_df[information_criterion].iloc[1:],
+            y=yticks_pos[1::2],
+            xerr=comp_df.dse[1:],
+            label="ELPD difference",
+            color=plot_kwargs.get("color_dse", "grey"),
+            fmt=plot_kwargs.get("marker_dse", "^"),
+            mew=linewidth,
+            elinewidth=linewidth,
+        )
+
+    else:
+        ax.set_yticks(yticks_pos[::2])
+
+    scale = comp_df["scale"][0]
+
     if insample_dev:
-        scale = comp_df[f"{information_criterion}_scale"][0]
-        p_ic = comp_df[f"p_{information_criterion}"]
+        p_ic = comp_df[f"p_{information_criterion.split('_')[1]}"]
         if scale == "log":
             correction = p_ic
         elif scale == "negative_log":
@@ -90,6 +97,7 @@ def plot_compare(
         ax.plot(
             comp_df[information_criterion] + correction,
             yticks_pos[::2],
+            label="In-sample ELPD",
             color=plot_kwargs.get("color_insample_dev", "k"),
             marker=plot_kwargs.get("marker_insample_dev", "o"),
             mew=linewidth,
@@ -102,13 +110,20 @@ def plot_compare(
         color=plot_kwargs.get("color_ls_min_ic", "grey"),
         lw=linewidth,
     )
+    if legend:
+        ax.legend(bbox_to_anchor=(1.01, 1), loc="upper left", ncol=1, fontsize=ax_labelsize)
 
-    scale_col = information_criterion + "_scale"
-    if scale_col in comp_df:
-        scale = comp_df[scale_col].iloc[0].capitalize()
-    else:
-        scale = "Deviance"
-    ax.set_xlabel(scale, fontsize=ax_labelsize)
+    if title:
+        ax.set_title(
+            f"Model comparison\n{'higher' if scale == 'log' else 'lower'} is better",
+            fontsize=ax_labelsize,
+        )
+
+    if scale == "negative_log":
+        scale = "-log"
+
+    ax.set_xlabel(f"{information_criterion} ({scale})", fontsize=ax_labelsize)
+    ax.set_ylabel("ranked models", fontsize=ax_labelsize)
     ax.set_yticklabels(yticks_labels)
     ax.set_ylim(-1 + step, 0 - step)
     ax.tick_params(labelsize=xt_labelsize)
diff --git a/arviz/plots/compareplot.py b/arviz/plots/compareplot.py
index e26705b..70b2836 100644
--- a/arviz/plots/compareplot.py
+++ b/arviz/plots/compareplot.py
@@ -8,10 +8,12 @@ from .plot_utils import get_plotting_function
 
 def plot_compare(
     comp_df,
-    insample_dev=True,
+    insample_dev=False,
     plot_standard_error=True,
     plot_ic_diff=True,
     order_by_rank=True,
+    legend=True,
+    title=True,
     figsize=None,
     textsize=None,
     labeller=None,
@@ -24,13 +26,14 @@ def plot_compare(
     """
     Summary plot for model comparison.
 
-    This plot is in the style of the one used in the book Statistical Rethinking (Chapter 6)
-    by Richard McElreath.
+    Models are compared based on their expected log pointwise predictive density (ELPD),
+    the ELPD is estimated either by Pareto smoothed importance sampling leave-one-out
+    cross-validation (LOO) or using the widely applicable information criterion (WAIC).
+    We recommend LOO in line with the work presented by Vehtari et al. (2016) available
+    here: https://arxiv.org/abs/1507.04544.
 
-    Notes
-    -----
-    Defaults to comparing Leave-one-out (psis-loo) if present in comp_df column,
-    otherwise compares Widely Applicable Information Criterion (WAIC)
+    This plot is in the style of the one used in the book Statistical Rethinking by
+    Richard McElreath.Chapter 6 in the first edition or 7 in the second.
 
 
     Parameters
@@ -38,17 +41,22 @@ def plot_compare(
     comp_df : pd.DataFrame
         Result of the :func:`arviz.compare` method
     insample_dev : bool, optional
-        Plot in-sample deviance, that is the value of the information criteria without the
-        penalization given by the effective number of parameters (pIC). Defaults to True
+        Plot in-sample ELPD, that is the value of the information criteria without the
+        penalization given by the effective number of parameters (p_loo or p_waic).
+        Defaults to False
     plot_standard_error : bool, optional
-        Plot the standard error of the information criteria estimate. Defaults to True
+        Plot the standard error of the ELPD. Defaults to True
     plot_ic_diff : bool, optional
-        Plot standard error of the difference in information criteria between each model
+        Plot standard error of the difference in ELPD between each model
         and the top-ranked model. Defaults to True
     order_by_rank : bool
         If True (default) ensure the best model is used as reference.
+    legend : bool
+        Add legend to figure. By default True.
     figsize : tuple, optional
         If None, size is (6, num of models) inches
+    title : bool:
+        Show a tittle with a description of how to interpret the plot. Defaults to True.
     textsize: float
         Text size scaling factor for labels, titles and lines. If None it will be autoscaled based
         on ``figsize``.
@@ -94,12 +102,12 @@ def plot_compare(
         >>>                  'Non-centered 8 schools': az.load_arviz_data('non_centered_eight')})
         >>> az.plot_compare(model_compare)
 
-    Plot standard error and information criteria difference only
+    Include the in-sample ELDP
 
     .. plot::
         :context: close-figs
 
-        >>> az.plot_compare(model_compare, insample_dev=False)
+        >>> az.plot_compare(model_compare, insample_dev=True)
 
     """
     if plot_kwargs is None:
@@ -119,7 +127,7 @@ def plot_compare(
     else:
         yticks_labels = labels
 
-    _information_criterion = ["loo", "waic"]
+    _information_criterion = ["elpd_loo", "elpd_waic"]
     column_index = [c.lower() for c in comp_df.columns]
     for information_criterion in _information_criterion:
         if information_criterion in column_index:
@@ -136,6 +144,8 @@ def plot_compare(
     compareplot_kwargs = dict(
         ax=ax,
         comp_df=comp_df,
+        legend=legend,
+        title=title,
         figsize=figsize,
         plot_ic_diff=plot_ic_diff,
         plot_standard_error=plot_standard_error,
diff --git a/arviz/stats/stats.py b/arviz/stats/stats.py
index 1b23e67..8995d01 100644
--- a/arviz/stats/stats.py
+++ b/arviz/stats/stats.py
@@ -61,19 +61,19 @@ def compare(
     scale: Optional[ScaleKeyword] = None,
     var_name: Optional[str] = None,
 ):
-    r"""Compare models based on PSIS-LOO `loo` or WAIC `waic` cross-validation.
+    r"""Compare models based on  their expected log pointwise predictive density (ELPD).
 
-    LOO is leave-one-out (PSIS-LOO `loo`) cross-validation and
-    WAIC is the widely applicable information criterion.
-    Read more theory here - in a paper by some of the leading authorities
-    on model selection dx.doi.org/10.1111/1467-9868.00353
+    The ELPD is estimated either by Pareto smoothed importance sampling leave-one-out
+    cross-validation (LOO) or using the widely applicable information criterion (WAIC).
+    We recommend loo. Read more theory here - in a paper by some of the
+    leading authorities on model comparison dx.doi.org/10.1111/1467-9868.00353
 
     Parameters
     ----------
     compare_dict: dict of {str: InferenceData or ELPDData}
         A dictionary of model names and :class:`arviz.InferenceData` or ``ELPDData``.
     ic: str, optional
-        Information Criterion (PSIS-LOO `loo` or WAIC `waic`) used to compare models. Defaults to
+        Method to estimate the ELPD, available options are "loo" or "waic". Defaults to
         ``rcParams["stats.information_criterion"]``.
     method: str, optional
         Method used to estimate the weights for each model. Available options are:
@@ -112,29 +112,29 @@ def compare(
 
     Returns
     -------
-    A DataFrame, ordered from best to worst model (measured by information criteria).
+    A DataFrame, ordered from best to worst model (measured by the ELPD).
     The index reflects the key with which the models are passed to this function. The columns are:
     rank: The rank-order of the models. 0 is the best.
-    IC: Information Criteria (PSIS-LOO `loo` or WAIC `waic`).
-        Higher IC indicates higher out-of-sample predictive fit ("better" model). Default LOO.
-        If `scale` is `deviance` or `negative_log` smaller IC indicates
+    elpd: ELPD estimated either using (PSIS-LOO-CV `elpd_loo` or WAIC `elpd_waic`).
+        Higher ELPD indicates higher out-of-sample predictive fit ("better" model).
+        If `scale` is `deviance` or `negative_log` smaller values indicates
         higher out-of-sample predictive fit ("better" model).
     pIC: Estimated effective number of parameters.
-    dIC: Relative difference between each IC (PSIS-LOO `loo` or WAIC `waic`)
-          and the lowest IC (PSIS-LOO `loo` or WAIC `waic`).
-          The top-ranked model is always 0.
+    elpd_diff: The difference in ELPD between two models.
+        If more than two models are compared, the difference is computed relative to the
+        top-ranked model, that always has a elpd_diff of 0.
     weight: Relative weight for each model.
         This can be loosely interpreted as the probability of each model (among the compared model)
         given the data. By default the uncertainty in the weights estimation is considered using
         Bayesian bootstrap.
-    SE: Standard error of the IC estimate.
+    SE: Standard error of the ELPD estimate.
         If method = BB-pseudo-BMA these values are estimated using Bayesian bootstrap.
-    dSE: Standard error of the difference in IC between each model and the top-ranked model.
+    dSE: Standard error of the difference in ELPD between each model and the top-ranked model.
         It's always 0 for the top-ranked model.
-    warning: A value of 1 indicates that the computation of the IC may not be reliable.
+    warning: A value of 1 indicates that the computation of the ELPD may not be reliable.
         This could be indication of WAIC/LOO starting to fail see
         http://arxiv.org/abs/1507.04544 for details.
-    scale: Scale used for the IC.
+    scale: Scale used for the ELPD.
 
     Examples
     --------
@@ -148,7 +148,7 @@ def compare(
            ...: compare_dict = {"non centered": data1, "centered": data2}
            ...: az.compare(compare_dict)
 
-    Compare the models using LOO-CV, returning the IC in log scale and calculating the
+    Compare the models using PSIS-LOO-CV, returning the ELPD in log scale and calculating the
     weights using the stacking method.
 
     .. ipython::
@@ -157,8 +157,10 @@ def compare(
 
     See Also
     --------
-    loo : Compute the Pareto Smoothed importance sampling Leave One Out cross-validation.
-    waic : Compute the widely applicable information criterion.
+    loo :
+        Compute the ELPD using the Pareto smoothed importance sampling Leave-one-out
+        cross-validation method.
+    waic : Compute the ELPD using the widely applicable information criterion.
     plot_compare : Summary plot for model comparison.
 
     References
@@ -171,42 +173,40 @@ def compare(
     try:
         (ics_dict, scale, ic) = _calculate_ics(compare_dict, scale=scale, ic=ic, var_name=var_name)
     except Exception as e:
-        raise e.__class__("Encountered error in ic computation of compare.") from e
+        raise e.__class__("Encountered error in ELPD computation of compare.") from e
     names = list(ics_dict.keys())
     if ic == "loo":
         df_comp = pd.DataFrame(
             index=names,
             columns=[
                 "rank",
-                "loo",
+                "elpd_loo",
                 "p_loo",
-                "d_loo",
+                "elpd_diff",
                 "weight",
                 "se",
                 "dse",
                 "warning",
-                "loo_scale",
+                "scale",
             ],
             dtype=np.float_,
         )
-        scale_col = "loo_scale"
     elif ic == "waic":
         df_comp = pd.DataFrame(
             index=names,
             columns=[
                 "rank",
-                "waic",
+                "elpd_waic",
                 "p_waic",
-                "d_waic",
+                "elpd_diff",
                 "weight",
                 "se",
                 "dse",
                 "warning",
-                "waic_scale",
+                "scale",
             ],
             dtype=np.float_,
         )
-        scale_col = "waic_scale"
     else:
         raise NotImplementedError(f"The information criterion {ic} is not supported.")
 
@@ -224,12 +224,11 @@ def compare(
     if method.lower() not in ["stacking", "bb-pseudo-bma", "pseudo-bma"]:
         raise ValueError(f"The method {method}, to compute weights, is not supported.")
 
-    ic_se = f"{ic}_se"
     p_ic = f"p_{ic}"
     ic_i = f"{ic}_i"
 
     ics = pd.DataFrame.from_dict(ics_dict, orient="index")
-    ics.sort_values(by=ic, inplace=True, ascending=ascending)
+    ics.sort_values(by=f"elpd_{ic}", inplace=True, ascending=ascending)
     ics[ic_i] = ics[ic_i].apply(lambda x: x.values.flatten())
 
     if method.lower() == "stacking":
@@ -267,7 +266,7 @@ def compare(
         )
 
         weights = w_fuller(weights["x"])
-        ses = ics[ic_se]
+        ses = ics["se"]
 
     elif method.lower() == "bb-pseudo-bma":
         rows, cols, ic_i_val = _ic_matrix(ics, ic_i)
@@ -286,10 +285,10 @@ def compare(
         ses = pd.Series(z_bs.std(axis=0), index=names)  # pylint: disable=no-member
 
     elif method.lower() == "pseudo-bma":
-        min_ic = ics.iloc[0][ic]
-        z_rv = np.exp((ics[ic] - min_ic) / scale_value)
+        min_ic = ics.iloc[0][f"elpd_{ic}"]
+        z_rv = np.exp((ics[f"elpd_{ic}"] - min_ic) / scale_value)
         weights = z_rv / np.sum(z_rv)
-        ses = ics[ic_se]
+        ses = ics["se"]
 
     if np.any(weights):
         min_ic_i_val = ics[ic_i].iloc[0]
@@ -305,19 +304,19 @@ def compare(
             weight = weights[idx]
             df_comp.at[val] = (
                 idx,
-                res[ic],
+                res[f"elpd_{ic}"],
                 res[p_ic],
                 d_ic,
                 weight,
                 std_err,
                 d_std_err,
                 res["warning"],
-                res[scale_col],
+                res["scale"],
             )
 
     df_comp["rank"] = df_comp["rank"].astype(int)
     df_comp["warning"] = df_comp["warning"].astype(bool)
-    return df_comp.sort_values(by=ic, ascending=ascending)
+    return df_comp.sort_values(by=f"elpd_{ic}", ascending=ascending)
 
 
 def _ic_matrix(ics, ic_i):
@@ -343,7 +342,7 @@ def _calculate_ics(
     ic: Optional[ICKeyword] = None,
     var_name: Optional[str] = None,
 ):
-    """Calculate loo and waic information criteria only if necessary.
+    """Calculate LOO or WAIC only if necessary.
 
     It always calls the ic function with ``pointwise=True``.
 
@@ -382,22 +381,22 @@ def _calculate_ics(
     precomputed_scale = None
     if precomputed_elpds:
         _, arbitrary_elpd = precomputed_elpds.popitem()
-        precomputed_ic = arbitrary_elpd.index[0]
-        precomputed_scale = arbitrary_elpd[f"{precomputed_ic}_scale"]
+        precomputed_ic = arbitrary_elpd.index[0].split("_")[1]
+        precomputed_scale = arbitrary_elpd["scale"]
         raise_non_pointwise = False
         if not f"{precomputed_ic}_i" in arbitrary_elpd:
             raise_non_pointwise = True
         if precomputed_elpds:
             if not all(
-                elpd_data.index[0] == precomputed_ic for elpd_data in precomputed_elpds.values()
+                elpd_data.index[0].split("_")[1] == precomputed_ic
+                for elpd_data in precomputed_elpds.values()
             ):
                 raise ValueError(
                     "All information criteria to be compared must be the same "
                     "but found both loo and waic."
                 )
             if not all(
-                elpd_data[f"{precomputed_ic}_scale"] == precomputed_scale
-                for elpd_data in precomputed_elpds.values()
+                elpd_data["scale"] == precomputed_scale for elpd_data in precomputed_elpds.values()
             ):
                 raise ValueError("All information criteria to be compared must use the same scale")
             if (
@@ -734,15 +733,15 @@ def loo(data, pointwise=None, var_name=None, reff=None, scale=None):
     Returns
     -------
     ELPDData object (inherits from :class:`pandas.Series`) with the following row/attributes:
-    loo: approximated expected log pointwise predictive density (elpd)
-    loo_se: standard error of loo
+    elpd: approximated expected log pointwise predictive density (elpd)
+    se: standard error of the elpd
     p_loo: effective number of parameters
     shape_warn: bool
         True if the estimated shape parameter of
         Pareto distribution is greater than 0.7 for one or more samples
     loo_i: array of pointwise predictive accuracy, only if pointwise True
     pareto_k: array of Pareto shape values, only if pointwise True
-    loo_scale: scale of the loo results
+    scale: scale of the elpd
 
         The returned object has a custom print method that overrides pd.Series method.
 
@@ -856,22 +855,22 @@ def loo(data, pointwise=None, var_name=None, reff=None, scale=None):
                 scale,
             ],
             index=[
-                "loo",
-                "loo_se",
+                "elpd_loo",
+                "se",
                 "p_loo",
                 "n_samples",
                 "n_data_points",
                 "warning",
                 "loo_i",
                 "pareto_k",
-                "loo_scale",
+                "scale",
             ],
         )
 
     else:
         return ELPDData(
             data=[loo_lppd, loo_lppd_se, p_loo, n_samples, n_data_points, warn_mg, scale],
-            index=["loo", "loo_se", "p_loo", "n_samples", "n_data_points", "warning", "loo_scale"],
+            index=["elpd_loo", "se", "p_loo", "n_samples", "n_data_points", "warning", "scale"],
         )
 
 
@@ -1591,14 +1590,14 @@ def waic(data, pointwise=None, var_name=None, scale=None, dask_kwargs=None):
     Returns
     -------
     ELPDData object (inherits from :class:`pandas.Series`) with the following row/attributes:
-    waic: approximated expected log pointwise predictive density (elpd)
-    waic_se: standard error of waic
+    elpd_waic: approximated expected log pointwise predictive density (elpd)
+    se: standard error of the elpd
     p_waic: effective number parameters
     var_warn: bool
         True if posterior variance of the log predictive densities exceeds 0.4
     waic_i: :class:`~xarray.DataArray` with the pointwise predictive accuracy,
             only if pointwise=True
-    waic_scale: scale of the reported waic results
+    scale: scale of the elpd
 
         The returned object has a custom print method that overrides pd.Series method.
 
@@ -1691,14 +1690,14 @@ def waic(data, pointwise=None, var_name=None, scale=None, dask_kwargs=None):
                 scale,
             ],
             index=[
-                "waic",
-                "waic_se",
+                "elpd_waic",
+                "se",
                 "p_waic",
                 "n_samples",
                 "n_data_points",
                 "warning",
                 "waic_i",
-                "waic_scale",
+                "scale",
             ],
         )
     else:
@@ -1706,12 +1705,12 @@ def waic(data, pointwise=None, var_name=None, scale=None, dask_kwargs=None):
             data=[waic_sum, waic_se, p_waic, n_samples, n_data_points, warn_mg, scale],
             index=[
                 "waic",
-                "waic_se",
+                "se",
                 "p_waic",
                 "n_samples",
                 "n_data_points",
                 "warning",
-                "waic_scale",
+                "scale",
             ],
         )
 
diff --git a/arviz/stats/stats_utils.py b/arviz/stats/stats_utils.py
index 6ea004e..4a7d3a6 100644
--- a/arviz/stats/stats_utils.py
+++ b/arviz/stats/stats_utils.py
@@ -462,12 +462,12 @@ class ELPDData(pd.Series):  # pylint: disable=too-many-ancestors
 
     def __str__(self):
         """Print elpd data in a user friendly way."""
-        kind = self.index[0]
+        kind = self.index[0].split("_")[1]
 
         if kind not in ("loo", "waic"):
             raise ValueError("Invalid ELPDData object")
 
-        scale_str = SCALE_DICT[self[f"{kind}_scale"]]
+        scale_str = SCALE_DICT[self["scale"]]
         padding = len(scale_str) + len(kind) + 1
         base = BASE_FMT.format(padding, padding - 2)
         base = base.format(


Test patch:
diff --git a/arviz/tests/base_tests/test_plots_bokeh.py b/arviz/tests/base_tests/test_plots_bokeh.py
index 39c9869..c7de725 100644
--- a/arviz/tests/base_tests/test_plots_bokeh.py
+++ b/arviz/tests/base_tests/test_plots_bokeh.py
@@ -339,12 +339,12 @@ def test_plot_compare_no_ic(models):
     model_compare = compare({"Model 1": models.model_1, "Model 2": models.model_2})
 
     # Drop column needed for plotting
-    model_compare = model_compare.drop("loo", axis=1)
+    model_compare = model_compare.drop("elpd_loo", axis=1)
     with pytest.raises(ValueError) as err:
         plot_compare(model_compare, backend="bokeh", show=False)
 
     assert "comp_df must contain one of the following" in str(err.value)
-    assert "['loo', 'waic']" in str(err.value)
+    assert "['elpd_loo', 'elpd_waic']" in str(err.value)
 
 
 def test_plot_ecdf_basic():
diff --git a/arviz/tests/base_tests/test_plots_matplotlib.py b/arviz/tests/base_tests/test_plots_matplotlib.py
index a23428f..8752c55 100644
--- a/arviz/tests/base_tests/test_plots_matplotlib.py
+++ b/arviz/tests/base_tests/test_plots_matplotlib.py
@@ -1127,7 +1127,7 @@ def test_plot_posterior_skipna_combinedims():
 
 
 @pytest.mark.parametrize(
-    "kwargs", [{"insample_dev": False}, {"plot_standard_error": False}, {"plot_ic_diff": False}]
+    "kwargs", [{"insample_dev": True}, {"plot_standard_error": False}, {"plot_ic_diff": False}]
 )
 def test_plot_compare(models, kwargs):
     model_compare = compare({"Model 1": models.model_1, "Model 2": models.model_2})
@@ -1141,12 +1141,12 @@ def test_plot_compare_no_ic(models):
     model_compare = compare({"Model 1": models.model_1, "Model 2": models.model_2})
 
     # Drop column needed for plotting
-    model_compare = model_compare.drop("loo", axis=1)
+    model_compare = model_compare.drop("elpd_loo", axis=1)
     with pytest.raises(ValueError) as err:
         plot_compare(model_compare)
 
     assert "comp_df must contain one of the following" in str(err.value)
-    assert "['loo', 'waic']" in str(err.value)
+    assert "['elpd_loo', 'elpd_waic']" in str(err.value)
 
 
 @pytest.mark.parametrize(
diff --git a/arviz/tests/base_tests/test_rcparams.py b/arviz/tests/base_tests/test_rcparams.py
index 71eb6b4..62ad421 100644
--- a/arviz/tests/base_tests/test_rcparams.py
+++ b/arviz/tests/base_tests/test_rcparams.py
@@ -280,10 +280,10 @@ def test_data_load():
 def test_stats_information_criterion(models):
     rcParams["stats.information_criterion"] = "waic"
     df_comp = compare({"model1": models.model_1, "model2": models.model_2})
-    assert "waic" in df_comp.columns
+    assert "elpd_waic" in df_comp.columns
     rcParams["stats.information_criterion"] = "loo"
     df_comp = compare({"model1": models.model_1, "model2": models.model_2})
-    assert "loo" in df_comp.columns
+    assert "elpd_loo" in df_comp.columns
 
 
 def test_http_type_request(monkeypatch):
diff --git a/arviz/tests/base_tests/test_stats.py b/arviz/tests/base_tests/test_stats.py
index 0990cc2..7c9bf77 100644
--- a/arviz/tests/base_tests/test_stats.py
+++ b/arviz/tests/base_tests/test_stats.py
@@ -233,7 +233,7 @@ def test_compare_multiple_obs(multivariable_log_likelihood, centered_eight, non_
     }
     with pytest.raises(TypeError, match="several log likelihood arrays"):
         get_log_likelihood(compare_dict["problematic"])
-    with pytest.raises(TypeError, match="error in ic computation"):
+    with pytest.raises(TypeError, match="error in ELPD computation"):
         compare(compare_dict, ic=ic)
     assert compare(compare_dict, ic=ic, var_name="obs") is not None
 
@@ -248,8 +248,9 @@ def test_calculate_ics(centered_eight, non_centered_eight, ic):
     elpddata_out, _, _ = _calculate_ics(elpddata_dict, ic=ic)
     mixed_out, _, _ = _calculate_ics(mixed_dict, ic=ic)
     for model in idata_dict:
-        assert idata_out[model][ic] == elpddata_out[model][ic]
-        assert idata_out[model][ic] == mixed_out[model][ic]
+        ic_ = f"elpd_{ic}"
+        assert idata_out[model][ic_] == elpddata_out[model][ic_]
+        assert idata_out[model][ic_] == mixed_out[model][ic_]
         assert idata_out[model][f"p_{ic}"] == elpddata_out[model][f"p_{ic}"]
         assert idata_out[model][f"p_{ic}"] == mixed_out[model][f"p_{ic}"]
 
@@ -265,7 +266,7 @@ def test_calculate_ics_ic_override(centered_eight, non_centered_eight):
     with pytest.warns(UserWarning, match="precomputed elpddata: waic"):
         out_dict, _, ic = _calculate_ics(in_dict, ic="loo")
     assert ic == "waic"
-    assert out_dict["centered"]["waic"] == waic(centered_eight)["waic"]
+    assert out_dict["centered"]["elpd_waic"] == waic(centered_eight)["elpd_waic"]
 
 
 def test_summary_ndarray():
@@ -489,7 +490,7 @@ def test_loo(centered_eight, multidim_models, scale, multidim):
     assert loo_pointwise is not None
     assert "loo_i" in loo_pointwise
     assert "pareto_k" in loo_pointwise
-    assert "loo_scale" in loo_pointwise
+    assert "scale" in loo_pointwise
 
 
 def test_loo_one_chain(centered_eight):
diff --git a/arviz/tests/base_tests/test_stats_utils.py b/arviz/tests/base_tests/test_stats_utils.py
index 2967941..ec18484 100644
--- a/arviz/tests/base_tests/test_stats_utils.py
+++ b/arviz/tests/base_tests/test_stats_utils.py
@@ -305,7 +305,7 @@ def test_get_log_likelihood_no_group():
 
 
 def test_elpd_data_error():
-    with pytest.raises(ValueError):
+    with pytest.raises(IndexError):
         repr(ELPDData(data=[0, 1, 2], index=["not IC", "se", "p"]))