How to add raw (frequentist) estimates to forest plot

Let’s say I want to add raw point estimates with 95% confidence intervals to forest plot. What’s the best way to achieve that? Here’s a simple example of creating a forest plot (without frequentist/raw estimates):

import bambi as bmb
import pandas as pd
import arviz as az

df_simple = pd.DataFrame({
    'x': ['A', 'B', 'C'],
    'y': [10, 20, 30],
    'n': [100, 100, 100]
})

m = bmb.Model('p(y, n) ~ 0 + x', data=df_simple, family='binomial')
idata = m.fit(cores=4)

m.predict(idata)

az.plot_forest(idata, var_names='p(y, n)_mean', combined=True)

When using the matplotlib backend, Arviz always returns an array of axis artists. You can use these axes as normal matplotlib axes objects to add more stuff to the plot. In this simple example, I add some dashed vertical lines:

import matplotlib.pyplot as plt

axis = az.plot_forest(idata, var_names='p(y, n)_mean', combined=True)[0]
ymin, ymax = axis.get_ylim()
axis.vlines([0.1, 0.2, 0.3], ymin, ymax, ls='--')
plt.show()

image

If you wanted to add frequentist estimates and confidence intervals, you would just need to axis.scatter the point estimates, then use axis.hlines to draw out the confidence intervals.

As with all matplotlib things, it will take some tweaking to get all the positioning right.

One unsolicited opinion, though: would it not be confusing to present your bayesian CIs along side frequentist CIs? In my mind, this would play to people’s general confusion about what frequentist CIs mean, by showing a measure of parameter uncertainty (the Bayesian CI) along side a statement about the model but not the parameter estimates (the frequentist CI).

3 Likes

Thanks! Do you have any advice on how to tweak y-position? I’d only like to plot raw point estimates.

You can get the tick locations with ax.yaxis.get_majorticklocs() maybe? Something like this:

axes = az.plot_forest(idata, var_names='p(y, n)_mean', combined=True)
axis = axes[0]

# I reverse it so that the y positions are from top to bottom
ylocs = axis.yaxis.get_majorticklocs()[::-1]
xlocs = [0.1, 0.2, 0.3]
axis.scatter(xlocs, ylocs, color='tab:red', marker='^', zorder=100)

image

1 Like