Credible intervals incorrect for values below threshold?

[LynnSchmittwilken]

Context

Great work on psignifit. While reproducing some of my results, however, I unfortunately stumbled across a potential issue with the credible intervals. Below a code example (three versions) which shows that credible intervals for performance values below the threshold seem to be incorrect.

Minimal example

import numpy as np
import psignifit as ps
from psignifit import psigniplot
import matplotlib.pyplot as plt


plotVersion = 3  # select which of the proposed versions to plot (1-3)

# Example from documentation:
data = np.array([[0.0010, 45.0000, 90.0000], [0.0015, 50.0000, 90.0000],
                 [0.0020, 44.0000, 90.0000], [0.0025, 44.0000, 90.0000],
                 [0.0030, 52.0000, 90.0000], [0.0035, 53.0000, 90.0000],
                 [0.0040, 62.0000, 90.0000], [0.0045, 64.0000, 90.0000],
                 [0.0050, 76.0000, 90.0000], [0.0060, 79.0000, 90.0000],
                 [0.0070, 88.0000, 90.0000], [0.0080, 90.0000, 90.0000],
                 [0.0100, 90.0000, 90.0000]])

options = {'sigmoid': 'norm',
           'experiment_type' : '2AFC'}

res = ps.psignifit(data, **options)

plt.figure()
psigniplot.plot_psychometric_function(res)


# Plot credible intervals across the whole curve
if plotVersion == 1: # produces an error
    y = np.linspace(0.5, 1, 50)  # get credible intervals between 50-100% performance
    CIlow = []; CIhigh = []
    for i in y:
        threshold, CI = res.threshold(i)  # use scaled sigmoid
        CIlow.append(CI[0][0]); CIhigh.append(CI[0][1])
    plt.plot(CIlow, y, label="lower bound"); plt.plot(CIhigh, y, label="upper bound")

elif plotVersion == 2:
    y = np.linspace(0.02, .5, 50) # get credible intervals according to error message of version 1
    CIlow = []; CIhigh = []
    for i in y:
        threshold, CI = res.threshold(i) # use scaled sigmoid
        CIlow.append(CI[0][0]); CIhigh.append(CI[0][1])
    plt.plot(CIlow, y+0.5, label="lower bound"); plt.plot(CIhigh, y+0.5, label="upper bound")

elif plotVersion == 3:
    y = np.linspace(0., 1, 50)  # get credible intervals over whole range for unscaled sigmoid
    CIlow = []; CIhigh = []
    for i in y:
        threshold, CI = res.threshold(i, unscaled=True) # use unscaled sigmoid instead
        CIlow.append(CI[0][0]); CIhigh.append(CI[0][1])
    
    y = np.linspace(0.5, 1, 50) # redefine to 50-100%
    plt.plot(CIlow, y, label="lower bound"); plt.plot(CIhigh, y, label="upper bound")

plt.legend()
plt.show()

Expected behavior

I would have expected that the psychometric function is always inbetween the credible intervals across the whole performance-range.

Observed behavior

The credible intervals for performance values below threshold (75% in the example) dont include the fitted function values.

Version 2 of code

Version 3 of code

[pberkes]

Hi Lynn, thanks for the report. We'll investigate quickly.

Two things:

1. which version of psignifit are you using? We've been working on the library; ideally, you are using 4.2 (we will soon have a 4.3 release)
2. the confidence intervals are over all the parameters, not only the threshold; the two curves you are generating are the two sigmoids using the parameters at the extreme of the confidence intervals. In that sense, you won't necessarily find the optimal sigmoid inside the two other sigmoids. For example, when the CI of the threshold is tiny, and the CI of the width is large, then you would get two sigmoids that intersect in the same point but with very different slopes

In any case, we'll have a close look. In the latest version of the code, I'm getting to identical plots for versions 1 and 3, and an error for version 2, which I think it's correct because 2AFC should have proportion correct above 0.5

[LynnSchmittwilken]

Hi Pietro! Sounds great!

1. I've been using version 4.2; and I just reproduced the exact same plots with a fresh virtual environment when only installing psignifit 4.2
2. I see, that makes sense! However, I dont think that the upper and lower bound should ever reverse order regardless.

Since in your version of the code, the error does not exist for plotVersion==1 but rather plotVersion==2, it sounds like psignifit 4.3 might have already have fixed this issue (in part?). I agree, that in theory plotVersion==2 should throw an error. I'll just wait for the update then!

[pberkes]

1. I'll check if there have been fixes to the CIs between 4.2 and 4.3dev
2. There might be a misunderstanding about the CIs here. If you think about having only a CI over the slope of the sigmoid, for example, then of course the curves would intersect.

For comparison, this is what happens if I fix all parameters except for the threshold. In that case, the optimal sigmoid lies between the two curves, just as you expect.

Code (v4.3dev):

import numpy as np
import psignifit as ps
from psignifit import psigniplot
import matplotlib.pyplot as plt

# Example from documentation:
data = np.array([[0.0010, 45.0000, 90.0000], [0.0015, 50.0000, 90.0000],
                 [0.0020, 44.0000, 90.0000], [0.0025, 44.0000, 90.0000],
                 [0.0030, 52.0000, 90.0000], [0.0035, 53.0000, 90.0000],
                 [0.0040, 62.0000, 90.0000], [0.0045, 64.0000, 90.0000],
                 [0.0050, 76.0000, 90.0000], [0.0060, 79.0000, 90.0000],
                 [0.0070, 88.0000, 90.0000], [0.0080, 90.0000, 90.0000],
                 [0.0100, 90.0000, 90.0000]])

options = {
    'sigmoid': 'norm',
    'experiment_type' : '2AFC',
    'fixed_parameters': {
        'gamma': np.float64(0.5),
        'eta': np.float64(0.0),
        'lambda': np.float64(0.0),
        'width': np.float64(0.005)
    }
}

res = ps.psignifit(data, **options)

plt.figure()
psigniplot.plot_psychometric_function(res)

y = np.linspace(0.5, 0.98, 50)  # get credible intervals between 50-100% performance
threshold, CI = res.threshold(y)  # use scaled sigmoid
plt.plot(CI['0.95'][0], y, label="lower bound")
plt.plot(CI['0.95'][1], y, label="upper bound")
plt.legend()

[pberkes]

and if I fix everything but the width

[LynnSchmittwilken]

Great, this behavior looks good to me!

So, to me, it sounds like 4.3 is fixing a (potentially minor) issue in the credible intervals related to why in 4.2 plotVersion==1 in my example code throws an error.
And then, it seems like there might just be a difference in default behavior as compared to the pre-pip version of psignifit. The old psignifit might have kept all parameters except for the threshold constant by default when using the "threshold"-method to compute the credible intervals, but this I am not sure about without checking the code more thoroughly.

Happy to give this another look, once 4.3 is out!

[pberkes]

My investigative report:

• In 4.2 there is a bug where Result.threshold calls inverse with the gamma and lambda parameters inverted. This was fixed in fixes threshold method, order of variables gamma and lambda was reversed #139 . This is why plot version 2 works even though plot version 1 is the correct one
• I believe the output of the method is correct, even though I have been unable to explain clearly why ... :-) the blue and orange curves are not bounds on the optimal sigmoids, they are the two most extreme sigmoids that you could get when you set the parameters at the extreme of their confidence intervals. Maybe @guillermoaguilar can jump in and help me out? :-)
• The big difference in the curve between version 2 and 3 is due to the same fixes threshold method, order of variables gamma and lambda was reversed #139 bug, otherwise they should be identical if you rescale the input variable y

[guillermoaguilar]

I'm happy to see engagement here, @LynnSchmittwilken and @pberkes ! :D

• I believe the output of the method is correct, even though I have been unable to explain clearly why ... :-) the blue and orange curves are not bounds on the optimal sigmoids, they are the two most extreme sigmoids that you could get when you set the parameters at the extreme of their confidence intervals. Maybe @guillermoaguilar can jump in and help me out? :-)

Indeed I think the method is correct, which is maybe counter-intuitive. Maybe another way of understanding is by looking at the posterior of parameters threshold and width. You run psignifit with argument debug=True, and then do the plot with psigniplot.plot_2D_margin(res, 'threshold', 'width').

For the completely unbounded case, no fixed parameters, the 2D posterior has a mass that is correlated, meaning a blob that is not just circular but tilted

The CIs are computed from the projection from the high dimensional posterior (in this case 4-D) into each parameter axes. Here we visualized only 2 of those axes, there are two more. The observed behavior (the crossings in the functions) can happen when the mass is like this, because there will be values at the boundary of the blob which is inside the interval for one parameter but outside for the other.

In other words the range in CI from the posterior is telling you that with your data all those possible shapes of sigmoids are possible within confidence, all that family of sigmoids shifted horizontally but also with higher or lower slope.

As a counter-example when all parameters are fixed but the width (the case of @pberkes), you have the following 2D marginals:

which is trivial as the width is fixed, so the mass is aligned to the axes. And here the CI follows the intuition for the range of possible sigmoids.

More generally, CIs are thought to give you the interval for the fitted parameter, not the whole function.

[guillermoaguilar]

So I'd say this is not an issue per se.

[guillermoaguilar]

On second thought, maybe I'm wrong and it doesn't have to do with the correlation in the posterior but just from the fact that here you are estimating threshold and width, both are free, and of course a sigmoid on the lower 'bound' for the width parameter will be shallower than the fitted, and then the function will be above the fitted for values less than the threshold. And the other way around, like here:

So maybe my explanation above is not the reason why this happens, it's just what happens when you have both threshold and width free to vary.
In any case the method seem to be correct.

@LynnSchmittwilken : if you want to visualize the range of functions that are 'compatible' with your data, I'd generate ones corresponding to varying both width and threshold simultaneously.

[pberkes]

Reporting @FelixWichmann 's reply by email:
"""
Well … I certainly sympathise with Lynn here ;-)

The plot of the psychometric function should really show the CIs, not two samples from the posterior -- thus I'd expect as a plot:
a. PF and CIs around the PF. Just a single PF and "error bar" (i.e. without blue & orange PFs)
b. PF and the CI as "sausage" around the PF. This can be obtained by sampling a few hundred PFs from the posterior: The fitted PF will never be outside the "sausage" of the random samples from the posterior (once you have enough samples).

In Lynn's plot the light blue and orange are not the lower and upper bound but two specific samples from the posterior with the threshold at the lower/upper and the slope at the upper/lower bound. To me this is misleading---as it was for Lynn---because I'd initially interpreted this as the bounds of the "sausage" (see b. above) ... and then it would be wrong!
If you really want to show a "quick and dirty" visualisation of the "sausage" you should draw four (!) half-PFs:

1. lower bound for threshold and lower bound for slope for proportion correct < 0.5 (nominal)
2. lower bound for threshold and upper bound for slope for proportion correct > 0.5 (nominal)
3. upper bound for threshold and lower bound for slope for proportion correct > 0.5 (nominal)
4. upper bound for threshold and upper bound for slope for proportion correct < 0.5 (nominal)

In the figure only the half-PFs for p > 0.5 are correct, and then you have to swap the slopes below the threshold ("fuse" the orange lower half to the upper blue half and vice versa).

BUT: I'd stick to either a.) or b.) above.
"""

[pberkes]

@FelixWichmann I can't really disagree... but the function works as advertized, don't you agree ?

This code implements Felix' proposal b. :

options = {
    'sigmoid': 'norm',
    'experiment_type' : '2AFC',
    'fixed_parameters': {
        'eta': np.float64(0.0),
        'lambda': np.float64(0.0),
        'threshold': np.float64(0.005)
    }
}

res = ps.psignifit(data, debug=True, **options)


params_grid = np.meshgrid(*(res.parameter_values.values()), indexing='ij')
params_stack = np.dstack([p.flatten() for p in params_grid]).squeeze()

# Sample from the posterior
samples_idx = np.random.choice(np.arange(params_stack.shape[0]), size=(100,), replace=True, p=res.debug['posteriors'].flatten())
params_samples = params_stack[samples_idx]

# --- Plot
plt.figure()

# Plot the samples from the posterior
sigmoid = res.configuration.make_sigmoid()
x = np.linspace(0.001, 0.01, num=1000)
for idx in range(params_samples.shape[0]):
    y = sigmoid(x, threshold=params_samples[idx, 3], width=params_samples[idx, 4], gamma=params_samples[idx, 1], lambd=params_samples[idx, 2])
    plt.plot(x, y, alpha=0.1, color='gray')

# Plot the MAP estimate
psigniplot.plot_psychometric_function(res, line_width=3)

Which gives this result (100 samples from the posterior):

[pberkes]

For the case with no fixed parameters, the plot looks like this

[guillermoaguilar]

thanks for clarifying Felix, @FelixWichmann , and thanks @pberkes for implementing the plots, they're really good and beautiful as well.

I think we should definitely include this type of plot to the package, in psigniplot.py, and add it to the plotting demos in the documentation. It is really clear to show what is the range of functions for the given confidence level.

[LynnSchmittwilken]

At this point, I am really not sure whether my original issue "just" stems from the bug in version 4.2. or whether there is also a more conceptual misunderstanding / disagreement.
I'd be happy to wait until 4.3. and give this another look, maybe with my reproduction attempts of what the old "threshold-method" did (pre-pip psignifit) in comparison to the new version. If there is actually a conceptual difference between the two, maybe having a real example helps to make a more informed decision about what the best default behavior of that method should be, and how to get from one to the other (ie which variables need to be changed)

Thanks already for all the effort and quick reactions!

[pberkes]

There is a consensus that the method is working correctly (in v4.3), with a doubt if the two extreme samples returned by the current method are useful (the current method being picking the parameter values at the extreme of the confidence interval, for all parameters independently).

My suggestion:

1. Remove the return_ci option from Result.threshold, it's not very useful
2. Add a Result.posterior_samples method to sample parameters from the posterior
3. Add a method to plot posterior samples