request for Bayes Factor plot

[drbenvincent]

I wonder whether it would be useful to have a Bayes Factor type plot, like this?

[ColCarroll]

It sounds like there is at least one person who would find it useful, and my impression is that they are standard (if somewhat controversial?)

We can leave this issue open for discussion, or for potential contributors to mention that they are working on it?

As an aside, this is the first time I've seen a sketch as a PR, and I love it.

[drbenvincent]

It's not impossible that I might attempt it. But only after I've actually submitted a bunch of papers, late in the summer.

There's a bit of debate between estimation and CI's vs Bayesian hypothesis testing. I've not really dug into it, but good examples of the pro-estimation approach are:

• Kruschke, J. K. (2011). Bayesian Assessment of Null Values Via Parameter Estimation and Model Comparison. Perspectives on Psychological Science, 6(3), 299–312. http://doi.org/10.1177/1745691611406925
• Kruschke, J. K., & Liddell, T. M. (2015). The Bayesian New Statistics: Two historical trends converge, 1–21.

And some examples on the pro hypothesis testing side...

• Morey, R. D., Rouder, J. N., Verhagen, J., & Wagenmakers, E.-J. (2014). Why hypothesis tests are essential for psychological science: a comment on Cumming (2014). Psychological Science, 25(6), 1289–1290. http://doi.org/10.1177/0956797614525969
• Wagenmakers, E.-J., Verhagen, J., Ly, A., Matzke, D., Steringoever, H., Rouder, J. N., & Morey, R. D. (2015). The Need for Bayesian Hypothesis Testing in Psychological Science, 1–29.

I'm not sure, but I think the Bayes Factor camp might have the slight edge, certainly in my field because of the legacy of frequentist hypothesis testing and because of it's use in tools like JASP which makes great plots like this... Personally I always report CI's or HDI's and optionally BF's when practical for the modelling context.

[treszkai]

We should keep in mind that this ArviZ not only allows people to plot the things they want to, but it also shapes how people think about statistical plotting (or it will in the future). So the question is if we think Bayesian hypothesis testing against a point parameter value is useful. Some arguments against this method are:

• BF is quite sensitive to the prior on the parameter.
• Testing against a single parameter value should be replaced with parameter estimation, so suggested statistics are the HDI, posterior median/mean/mode, and the posterior probability of a ROPE'd value are more meaningful statistics. (See e.g. Kruschke: The Bayesian New Statistics.)

Furthermore, the suggested graph looks like the prior and posterior plotted on the same Axes, with a particular value singled out (which is incompatible with Bayesian thinking), and the BF plotted on it. This could be done without new plotting functions, as described in #749.

(Nonetheless, +1 for the sketch!)

[drbenvincent]

Well I guess this just highlights that there is a debate between parameter estimation vs treating specific values as special. I was previously very much in the parameter estimation camp, but then got convinced that there is some usefulness in evaluating point hypotheses with Bayes Factors sometimes. I am not a statistician, so I can't really contribute to that debate as such. But Bayes Factors do seem to have some influence in my field - may that's a hang up from people who've not fully gotten away from NHST yet. Although I don't think that's really it... EJ Wagenmakers is a strong advocate of using Bayes Factors for example.

In terms of sensitivity to priors... yes! Although there are some classes of model (eg linear modelling, but not just that) where it is possible to define sensible priors which you can defend and justify.

In terms of people missing BF's and not knowing what they are doing, all I can say is that that is a danger in many ways when doing Bayesian modelling. I wasn't advocating for a pre-defined workflow, pushing people into BF's, more like just requesting the feature for those who want it and look for it in documentation etc.

[orduek]

I find that many people in my field request this as well. Does anyone currently work on it? I'll be happy to help but I would probably need some assistance on the way.

[drbenvincent]

Hi @orduek Feel free to work on this. I am not a core dev, but I am sure that you'll get some support if you need it :) Ben

[orduek]

Ok, what do you think about this function:

def display_delta(trace, x, varName, prior, null):
    # grab trace, x (linspace to plot) a variable name to compute difference and prior. null is the parameter we want to compare
    from scipy import stats
    # BFs based on density estimation (using kernel smoothing instead of spline)
    # stack trace posterior
    tr = trace.posterior.stack(draws=("chain", "draw"))
    
    tmp = az.summary(trace, var_names=[str(varName)])
    # 94% confidence interval:
    x0 = az.hdi(trace)[str(varName)][0].values
    x1 = az.hdi(trace)[str(varName)][1].values

    t_delt = tr[str(varName)]
    my_pdf = stats.gaussian_kde(t_delt)
    prior_pdf = stats.gaussian_kde(prior)
    plt.plot(
        x, my_pdf(x), "--", lw=2.5, alpha=0.6, label="Posterior"
    )  # distribution function
    plt.plot(x, prior_pdf(x), "r-", lw=2.5, alpha=0.6, label="Prior")
    
    posterior = my_pdf(0)  # this gives the pdf at point delta = 0
    prior = prior_pdf(0)#stats.cauchy.pdf(-.5)  # height of order-restricted prior at delta = 0
    BF10 = posterior / prior
    print("the Bayes Factor is %.3f" % (BF10))
    
    plt.plot(
        [null, null],
        [posterior, prior],
        "k-",
        [null,null],
        [posterior, prior],
        "ko",
        lw=1.5,
        alpha=1,
    )
    plt.xlabel("Delta")
    plt.ylabel("Density")
    plt.legend(loc="upper left")
    plt.show()

we can either add the prior based on prior predictive checks, or add some brewed one (like uniform of normal/cauchy)

Something like that:

x = np.linspace(-1, 1, 4000)
# set prior (taken from the interaction model)
prior = ppChecks['interaction'][:,0,1] - ppChecks['interaction'][:,0,0]
# using uniform distribution 
r = np.random.uniform(low=-10, high=10, size = 4000)

display_delta(trace, x, 'differenceGroups end of treatment', prior, null=0)

Results looks like this:

Very preliminary, and not sure its the best calculation, but that's one way I use to do it

[drbenvincent]

Cool.

So just a few disjointed comments/ideas:

• tmp, x0, x1 seem to be unused
• allow user to specify the x value (e.g. value) where we evaluate the Bayes Factor via a kwarg that defaults to zero
• don't think you need to do str(varName). varName will be a string anyway.
• consider just passing in idata (and not prior samples separately) and just grab the idata.posterior and idata.prior groups.
• change t_delt name to not have any assumption about deltas
• consider adding checks that the target variable is univariate, by checking shapes
• I think the main thing which might get tricky is calculating the density (via kde) when the variable is bounded. For example, if we are considering an absolute difference where values are always >0, or a Beta distribution where all samples will be within 0-1.
• I don't quite understand the part about posterior predictives... this will relate to data space, but we are interested in the parameter space.
• Consider returning BF10 and BF01 values.

[drbenvincent]

Separate issue... if your particular situation evaluates a difference, then I'd strongly recommend avoiding a uniform prior over the difference. It might be better to consider a Normal or a Cauchy, centred on zero. Bonus points if the delta can be expressed in the form of an effect size. But that might be best addressed here https://discourse.pymc.io

[orduek]

Thanks for the feedback. I'll work on that and post the code here one more time before actually implementing it within the arviz env. Indeed - I have expressed it in effect size and used prior predictive for the BF. Although I never feel fully confident in those decisions.

[orduek]

Ok, I have tried to fix most of the comments. I kept the prior variable, as in some cases (like the one presented here), I compare a deterministic posterior to a prior. In this case, there is a prior shaped (2,2), and the deterministic posterior takes prior[0,:] - prior[1,:], so the varName will not be the same. I agree that using the xarray might be better, but I need to think of a good solution for that.

Regarding the bounded distributions - I think the solution will go forward by extracting the distribution family (i.e., beta, normal, etc) and then using the correct distribution (although I don't know exactly how to generate the same for beta distribution). Although, might be a better way to calculate the BF using loglikelihoods (that are already a product of the data?)

def display_delta(trace, varName, prior, family = 'normal',  ref_val=0):
   
    # grab trace, a variable name to compute difference and prior. ref_val is the parameter we want to compare
    from scipy import stats 
    
    # test some elemtns
    # varName should be string
    if isinstance(varName, str):
        None
    else:
        print('varName is not a string')
    
    # BFs based on density estimation (using kernel smoothing instead of spline)
    # stack trace posterior
    tr = trace.posterior.stack(draws=("chain", "draw"))
    
    # grab prior
    try:
        trPrior = trace.prior.stack(draws=("chain", "draw"))
    except:
        print("No prior in the idata xarray, please use add_group to add prior")
    
    post = tr[varName]
    if post.ndim > 1:
        print("Posterior distribution has {post.ndim} dimensions")
    else:
        None
        
    prior = trPrior[varName]
    if family=='normal':
        # generate vector
        x = np.linspace(-4, 4, 8000)
        my_pdf = stats.gaussian_kde(post)
        prior_pdf = stats.gaussian_kde(prior)
    elif family=='beta':
        # This should be fixed
        x = np.linspace(0, 1, 8000)
        my_pdf = stats.beta(post)
        prior_pdf = stats.beta(prior)
        
    plt.plot(
        x, my_pdf(x), "--", lw=2.5, alpha=0.6, label="Posterior"
    )  # distribution function
    plt.plot(x, prior_pdf(x), "r-", lw=2.5, alpha=0.6, label="Prior")
    
    posterior = my_pdf(ref_val)  # this gives the pdf at point delta = 0
    prior = prior_pdf(ref_val)#stats.cauchy.pdf(-.5)  # height of order-restricted prior at delta = 0
    BF10 = posterior / prior
    BF01 = prior / posterior
    print("the Bayes Factor 10 is %.3f" % (BF10))
    print("the Bayes Factor 01 is %.3f" % (BF01))
    
    plt.plot(
        [ref_val, ref_val],
        [posterior, prior],
        "k-",
        [ref_val,ref_val],
        [posterior, prior],
        "ko",
        lw=1.5,
        alpha=1,
    )
    plt.xlabel("Delta")
    plt.ylabel("Density")
    plt.legend(loc="upper left")
    plt.show()

[drbenvincent]

Tricky to know at what point it makes sense to do further edits in a pull request. But just a few more random thoughts:

• Maybe just me, but I find the plt.plot at the end hard to interpret with the vectors of inputs. Maybe just two separate plots here?
• I think dealing with bounded variables is still going to need kde, but with bounds. Or at least, I think it is too much to assume a parametric form... there's no guarantee that the posterior will behave that nicely. This is probably going to be the most tricky part of making this a generalisable function.
• Try to eliminate the hard coded values (e.g. x = np.linspace(-4, 4, 8000)) and look up ranges of the samples perhaps.

Looks to be taking shape

[orduek]

Thank you again :) So, I've used the prior's min and max (and shape) arguments to generate the x, this is a bit tricky, as now the graphs looks less "pretty" (attached). I have added another example (commented our) which uses the posterior min and max for the x, what do you think is best?

Here's the code.

def display_delta(trace, varName, prior, family = 'normal',  ref_val=0):
   
    # grab trace, a variable name to compute difference and prior. ref_val is the parameter we want to compare
    from scipy import stats 
    
    # test some elemtns
    # varName should be string
    if isinstance(varName, str):
        None
    else:
        print('varName is not a string')
    
    # BFs based on density estimation (using kernel smoothing instead of spline)
    # stack trace posterior
    tr = trace.posterior.stack(draws=("chain", "draw"))
       
    post = tr[varName]

    if post.ndim > 1:
        print("Posterior distribution has {post.ndim} dimensions")
    else:
        None
        
    if family=='normal':
        # generate vector
        x = np.linspace(np.min(prior), np.max(prior),prior.shape[0])
        #x = np.linspace(np.min(post), np.max(post),prior.shape[0])
        my_pdf = stats.gaussian_kde(post)
        prior_pdf = stats.gaussian_kde(prior)
    elif family=='beta':
        # This should be fixed
        x = np.linspace(0, 1, 8000)
        my_pdf = stats.gaussian_kde(post)
        prior_pdf = stats.gaussian_kde(prior)
        
    plt.plot(
        x, my_pdf(x), "--", lw=2.5, alpha=0.6, label="Posterior"
    )  # distribution function
    plt.plot(x, prior_pdf(x), "r-", lw=2.5, alpha=0.6, label="Prior")
    
    posterior = my_pdf(ref_val)  # this gives the pdf at point delta = 0
    prior = prior_pdf(ref_val)#stats.cauchy.pdf(-.5)  # height of order-restricted prior at delta = 0
    BF10 = posterior / prior
    BF01 = prior / posterior
    print("the Bayes Factor 10 is %.3f" % (BF10))
    print("the Bayes Factor 01 is %.3f" % (BF01))
    
   
    plt.plot(ref_val, posterior, "ko", lw=1.5, alpha=1)
    plt.plot(ref_val, prior, "ko", lw=1.5, alpha=1)

   
    plt.xlabel("Delta")
    plt.ylabel("Density")
    plt.legend(loc="upper left")
    plt.show()

[drbenvincent]

The images aren't working properly. They need to be outside the code block.

PS If you do ```python ... then you get syntax highlighting :)

[drbenvincent]

Thank you again :) So, I've used the prior's min and max (and shape) arguments to generate the x, this is a bit tricky, as now the graphs looks less "pretty" (attached). I have added another example (commented our) which uses the posterior min and max for the x, what do you think is best?

Can always use percentiles rather than min/max. But can always focus on the aesthetic and Arviz style issues at a later point.

[orduek]

ok, updated using the python magic and figures

[drbenvincent]

Looks good. Could consider adding a kwarg xlim=None where a user can optionally pass a tuple or list of x axis limits. So you can apply that zoom at the end, after plotting. Or even HDI value (e.g. 0.95) which automatically zooms in to that range in the posterior?

[orduek]

ok. I've checked both options for xlim, using the hdi might create a situation where zero is not present (i.e. the ref_val), so I think for now it might be better to keep the xlim, which will enable the user to "zoom in" when needed.

What do you think?

Here's the latest code. The big issue is still the bounded distributions

def display_delta(trace, varName, prior, family = 'normal',  ref_val=0, xlim=None):
   
    # grab trace, a variable name to compute difference and prior. ref_val is the parameter we want to compare
    from scipy import stats 
    
    # test some elemtns
    # varName should be string
    if isinstance(varName, str):
        None
    else:
        print('varName is not a string')
    
    # BFs based on density estimation (using kernel smoothing instead of spline)
    # stack trace posterior
    tr = trace.posterior.stack(draws=("chain", "draw"))
       
    post = tr[varName]
    if post.ndim > 1:
        print("Posterior distribution has {post.ndim} dimensions")
    else:
        None
        
    if family=='normal':
        # generate vector
        if xlim is None:
            x = np.linspace(np.min(prior), np.max(prior),prior.shape[0])
        else:
            x = np.linspace(xlim[0], xlim[1],prior.shape[0])
        #x = np.linspace(np.min(post), np.max(post),prior.shape[0])
        my_pdf = stats.gaussian_kde(post)
        prior_pdf = stats.gaussian_kde(prior)
    elif family=='beta':
        # This should be fixed
        x = np.linspace(0, 1, 8000)
        my_pdf = stats.gaussian_kde(post)
        prior_pdf = stats.gaussian_kde(prior)
        
    plt.plot(
        x, my_pdf(x), "--", lw=2.5, alpha=0.6, label="Posterior"
    )  # distribution function
    plt.plot(x, prior_pdf(x), "r-", lw=2.5, alpha=0.6, label="Prior")
    
    posterior = my_pdf(ref_val)  # this gives the pdf at point delta = 0
    prior = prior_pdf(ref_val)#stats.cauchy.pdf(-.5)  # height of order-restricted prior at delta = 0
    BF10 = posterior / prior
    BF01 = prior / posterior
    print("the Bayes Factor 10 is %.3f" % (BF10))
    print("the Bayes Factor 01 is %.3f" % (BF01))
    
   
    plt.plot(ref_val, posterior, "ko", lw=1.5, alpha=1)
    plt.plot(ref_val, prior, "ko", lw=1.5, alpha=1)

   
    plt.xlabel("Delta")
    plt.ylabel("Density")
    plt.legend(loc="upper left")
    plt.show()

[drbenvincent]

I think my strategy would be to include kwargs of lower_bound=None and upper_bound=None and allow users to optionally specify lower and/or upper bounds. But for now, you could just throw a NotImplementedError perhaps.

I'd also perhaps start to separate a function that you are using for your own research purposes to get the job done, versus a more general function that could be used in aviz. So the latter would avoid having situation specific variable names (like delta).

Could also start to use the more OOP way of Matplotlib with fig, ax = plt.subplots() then ax.plot(...

Maybe also return ax, BF10, BF01

Think that's most of my ideas at this point :)

[orduek]

Thank you again, I learn a lot from this conversation! Here's the latest version.

def plot_BF(trace, varName, prior, family = 'normal',  ref_val=0, xlim=None):
   
    # grab trace, a variable name to compute difference and prior. ref_val is the parameter we want to compare
    from scipy import stats 
    
    # test some elemtns
    # varName should be string
    if isinstance(varName, str):
        None
    else:
        print('varName is not a string')
    
    # BFs based on density estimation (using kernel smoothing instead of spline)
    # stack trace posterior
    tr = trace.posterior.stack(draws=("chain", "draw"))
       
    post = tr[varName]
    if post.ndim > 1:
        print("Posterior distribution has {post.ndim} dimensions")
    else:
        None
        
    if family=='normal':
        # generate vector
        if xlim is None:
            x = np.linspace(np.min(prior), np.max(prior),prior.shape[0])
        else:
            x = np.linspace(xlim[0], xlim[1],prior.shape[0])
        #x = np.linspace(np.min(post), np.max(post),prior.shape[0])
        my_pdf = stats.gaussian_kde(post)
        prior_pdf = stats.gaussian_kde(prior)
    elif family!='normal':
        # for now and error of notImplemented
        print("The method for {family} distribution is not implemented yet")
    
    fig, ax = plt.subplots()    
    ax.plot(
        x, my_pdf(x), "--", lw=2.5, alpha=0.6, label="Posterior"
    )  # distribution function
    ax.plot(x, prior_pdf(x), "r-", lw=2.5, alpha=0.6, label="Prior")
    
    posterior = my_pdf(ref_val)  # this gives the pdf at point delta = 0
    prior = prior_pdf(ref_val)#stats.cauchy.pdf(-.5)  # height of order-restricted prior at delta = 0
    BF10 = posterior / prior
    BF01 = prior / posterior
    print("the Bayes Factor 10 is %.3f" % (BF10))
    print("the Bayes Factor 01 is %.3f" % (BF01))
    
   
    ax.plot(ref_val, posterior, "ko", lw=1.5, alpha=1)
    ax.plot(ref_val, prior, "ko", lw=1.5, alpha=1)

   
    plt.xlabel("Delta")
    plt.ylabel("Density")
    plt.legend(loc="upper left")
    plt.show()
    return {'BF10': BF10, 'BF01':BF01}, ax

[drbenvincent]

Oh... with the xlim that was meant to apply purely to setting the axis limits, so a plotting thing, not a calculation thing. That will likely change the results.

Stick with x = np.linspace(np.min(prior), np.max(prior),prior.shape[0]) but then something like

if xlim is not None:
    ax.set(xlim=xlim)

after your plan plot functions

[drbenvincent]

You could also have a kwarg ax=None and in the function

if ax is None:
    fig, ax = plt.subplots()

That would allow users to get the plot targeted at an axis of their choice (eg part of a figure with multiple subplots), but if they don't specify an axis, a new fig, ax will be generated.

[orduek]

I'm not sure I follow the xlim comment. I've checked, and in the current code, it just affects the vector we send to the plots, no? The BFs are the same whether I don't use xlim (i.e. None) or I use [-100,100] or [-1,.5].

Here's the latest code, adding the ax option:

def plot_BF(trace, varName, prior, family = 'normal',  ref_val=0, xlim=None, ax=None):
   
    # grab trace, a variable name to compute difference and prior. ref_val is the parameter we want to compare
    from scipy import stats 
    
    # test some elemtns
    # varName should be string
    if isinstance(varName, str):
        None
    else:
        print('varName is not a string')
    
    # BFs based on density estimation (using kernel smoothing instead of spline)
    # stack trace posterior
    tr = trace.posterior.stack(draws=("chain", "draw"))
       
    post = tr[varName]
    if post.ndim > 1:
        print("Posterior distribution has {post.ndim} dimensions")
    else:
        None
        
    if family=='normal':
        # generate vector
        if xlim is None:
            x = np.linspace(np.min(prior), np.max(prior),prior.shape[0])
        else:
            x = np.linspace(xlim[0], xlim[1],prior.shape[0])
        #x = np.linspace(np.min(post), np.max(post),prior.shape[0])
        my_pdf = stats.gaussian_kde(post)
        prior_pdf = stats.gaussian_kde(prior)
    elif family!='normal':
        # for now and error of notImplemented
        print("The method for {family} distribution is not implemented yet")
    
    if ax is None:
        fig, ax = plt.subplots()
    
    ax.plot(
        x, my_pdf(x), "--", lw=2.5, alpha=0.6, label="Posterior"
    )  # distribution function
    ax.plot(x, prior_pdf(x), "r-", lw=2.5, alpha=0.6, label="Prior")
    
    posterior = my_pdf(ref_val)  # this gives the pdf at point delta = 0
    prior = prior_pdf(ref_val)#stats.cauchy.pdf(-.5)  # height of order-restricted prior at delta = 0
    BF10 = posterior / prior
    BF01 = prior / posterior
    print("the Bayes Factor 10 is %.3f" % (BF10))
    print("the Bayes Factor 01 is %.3f" % (BF01))
    
   
    ax.plot(ref_val, posterior, "ko", lw=1.5, alpha=1)
    ax.plot(ref_val, prior, "ko", lw=1.5, alpha=1)

   
    plt.xlabel("Delta")
    plt.ylabel("Density")
    plt.legend(loc="upper left")
    plt.show()
    return {'BF10': BF10, 'BF01':BF01}, ax

[drbenvincent]

Sorry, I read the code wrongly. The concern was: if your actual prior was Uniform(-1, 1) then the pdf at zero is going to be different then if your prior was Uniform(-20, 20). But I don't think that's actually a concern here... all you changed was the range over which you plot, so I believe it would be the same as if you plot the whole range then zoom in.

[orduek]

Yes, the prior decision is affecting BF dramatically. I can easily move this one to Uniform(-100,100) and get a BF10 >50. I think this is one of the main criticism regarding BF

[drbenvincent]

Agreed. Which is why it's wise to place a prior on effect sizes because then the scale is at least interpretable and comparable with effect sizes seen in your area of the literature.

[orduek]

How do you suggest to proceed from this stage?

[drbenvincent]

I've not yet properly contributed to this repo, so I'm not 100% sure. Could either put together a pull request and then go through the review process to get things polished. I'd expect it to involve a reasonable amount of work to get it up to scratch to the reviewers. But doing all this via a pull request is probably best in terms of people making comments. Other option would be to try to deal with the kde for bounded variables before the pull request. Personally I'd be tempted to go for the pull request route to get more eyes and feedback from the core devs.

[orduek]

Ok. I've cloned and branched (plot_bf), but can't push the changes in order to create a PR. I probably need some guidance

[drbenvincent]

I'm only half way through my own first PR, so not an expert here. But I believe you have to fork (not just clone) the repo. I'd check out the Contributing guide https://python.arviz.org/en/latest/contributing/index.html