How to Fit a Psychometric Function in R (Visualizing MLE using lattice package)

Today we will be examining:

1. What is a psychometric function?
2. How to fit a psychometric function in R?
3. How to visualize maximum likelihood estimation (MLE) in R?
4. How to extract a threshold from a psychometric function?

Before we begin, I would like to mention that much has been written on this topic (gross understatement), and there is absolutely no way that I can cover everything in one post. However, I hope by the end of this post you will have gained some knowledge about this topic which is dear to my heart.

With that said, let us begin.

 

1. WHAT IS A PSYCHOMETRIC FUNCTION?

A psychometric function is generally defined as a mathematical model used to represent the detection or discrimination of a stimulus, whereby responses are probabilistic. The probability of a correct or incorrect response is based upon the stimulus intensity shown, where typically correct responses are made more frequently as the stimulus intensity increases. Note, the inverse relationship can also hold such that correct responses decrease with an increase in stimulus intensity. Therefore, psychometric functions model the probability of correct detection or discrimination of a stimulus between 0 and 1 in a manner that approximates human or animal behaviour.

To demonstrate the types of experiments that give rise to behavioural responses in which a psychophysical function (model) can be fit, here we describe a simple experimental procedure.

Say that an observer views a computer monitor, and on every trial the observer must respond whether or not they detect a flash of light onscreen. Throughout the experiment, we vary the intensity of the flash such that at high intensities, the observer almost always correctly detects the flash of light, at low intensities the observer responds at chance, and at intermediate intensities their detection performance is somewhere between these two extremes.

Before we plot these results in R, we load the appropriate libraries into our R workspace that will be required to run the rest of the code shown in this post.

We will now simulate responses for an observer that was run in this experiment. This will also allow us to compare the accuracy of our model (i.e., psychometric function) given the data to the true probability of detection for this observer.

Let us now plot the probability of correct detection of the flash of light on the y-axis as a function of the stimulus intensity on the x-axis:

And from this we get:

There are several things to note about this figure. First, because each trial contains two possible responses that an observer can make, the guess rate (i.e., chance performance) is defined as 1/ m-alternatives (m=2), which results in a probability of 0.5. The guess rate is represented by the dotted red line in the figure.

Next, notice that as the intensity of the flash of light increases, there is an increase in the probability of the observer correctly detecting the flash.  The relationship between the probability of correct detection as a function of the stimulus intensity also seems to be non-linear (sigmoidal).

  Therefore, the function that we use to model the probability of correct detection for this observer should also be sigmoidal in shape. However, it is important in practice to compare models to ensure a non-linear function is necessary.  Often times, a linear function can fit the data much better than a non-linear model.

Here, we will skip this step for the sake of time, as we generated the probability of correct detection from a non-linear probability distribution (i.e., Weibull).

Time to fit a psychometric function to this data.

 

2. HOW TO FIT A PSYCHOMETRIC FUNCTION IN R?

Psychometric functions are typically generated by fitting a sigmoidal function to the probability of correct detection or discrimination of a stimulus using a maximum likelihood fitting procedure. The generic definition of a psychometric function takes the form of:

\psi(x;\alpha,\beta,\gamma,\lambda) = \gamma + (1 - \gamma - \lambda) \text{F}(x;\alpha,\beta)

Where:

  • \alpha – Scale parameter.  Typically determines where along the abscissa the function is located.   This is not entirely true for a Weibull function, as \alpha and \beta jointly determine the shape and scale of the function.
  • \beta – Shape parameter. Typically determines the slope of the function.  Again, this is not entirely true for a Weibull function.
  • \gamma – Guess rate.  If the observer is guessing on every trial, what is the probability of a correct response.  Because an observer can respond with one of two possible responses on every trial in our experiment, this value will be fixed at 0.5 when we fit our function.
  • \lambda – Lapse rate.  The observer will likely get some responses wrong at “easy” intensities as a result of pressing the wrong button by accident, inattention, etc. This is taken into consideration in modelling an observer’s response.  We will set this value at 0.025 which is the known lapse rate that we used when generating responses for our fictitious obserber.  In practice, this value will not be known.  Much discussion has surrounded what value to assign this parameter. Therefore, I urge the reader to read some of the references at the end of this article for more help as to what to assign this value when modelling responses.

\text{F}(x;\alpha,\beta) can be represented by a number of sigmoidal functions.  For this post, we will represent the probability of a correct response using a Weibull function, which is defined by:

\text{F}_W = 1 - \text{exp}(- (\dfrac{x}{\alpha})^\beta )

For simplicity in describing the likelihoods, we will now define the function as:

W(x;\alpha,\beta,\gamma,\lambda) = \gamma + (1 - \gamma - \lambda) \text{F}_W(x;\alpha,\beta)

Now that we have mathematically defined what our psychometric function is, how do we determine what the best parameters for this function are to best fit our data?  Because the response outcome on every trial in this experiment is binary (0 for incorrect or 1 for correct), we turn to the use of “likelihood” as a measure of fit.

We define our likelihood function as:

\text{log} ( \prod_{} W(x_{i})^{(k_{i})} (1 - W(x_{i})^{(1 - k_{i})})

Here, W(x_{i}) is the probability of a correct response as described by our function above evaluated at stimulus intensity x on trial i.  The right hand side of the equation is that probability of incorrect responses evaluated at stimulus intensity x on trial i. We let k_{i} =1 for correct trials and k_{i} = 0 for incorrect trials.

Note, we take the log of this function to ensure that we do not run into computational constraints in computing the product of two very small numbers. Therefore, we can re-write this as:

\sum k_{i}\text{log}(W_{i}) + (1 - k_{i})\text{log}(1 - W_{i})

The last thing we will do is multiply this equation by -1, as optimization routines typically find minima of features spaces (more on this later). So we have:

-\sum k_{i}\text{log}(W_{i}) + (1 - k_{i})\text{log}(1 - W_{i})

Now it is time to define the negative log likelihood function in R, and feed it along with initial parameter estimates for \alpha and \beta into an optimization algorithm:

It is important to choose plausible initial estimates for the parameters being fit to the function. This ensures an optimization routine can find values for the parameters that provide a global minimum in the feature space.

When we run this script, we get the following estimates for the parameters of the function:

\alpha = 0.6614301
\beta = 3.2294249

So how well do these values for the function fit the data for our observer? The first method is to visualize the function by plotting our sigmoidal curve over the data.

Running the code above produces:

In the figure above, the psychometric function that was fit to the data is green, while the true probability function from which the responses were drawn from is displayed in black.  We can see the estimate is slightly off from the true probability, which makes sense given we are sampling from this distribution. However, our function appears to be a good estimate of the true probability of correctly detecting the flash of light for this observer.

Although we are quite confident in this scenario that we chose the best parameters for the function to fit the data, reflected in the value of our log likelihood, how do we know that other functions do not exist that provide an even better fit to the data? Unfortunately, most metrics of fit are quite meaningless when only one function is fit to data, as you have nothing to compare that metric of fit to.  Therefore, several functions should be fit to data and compared on the metric(s) of fit (e.g., AIC, BIC, Mallows’ C, etc.). For the time being, I will leave that as an exercise for the reader.

 

3. HOW TO VISUALIZE MAXIMUM LIKELIHOOD ESTIMATION (MLE) IN R?

It is insightful to visualize how the optimization procedure in R is going about selecting values for the parameters of interest that best fit the data.  Of course, this can be done because we are estimating values for only two parameters in this example.

As was mentioned above, the optimization procedure tries to minimize the function it is passed, in this case being the negative log likelihood.  Therefore, the fitting procedure seeks to find the minima of the parameter space.  This can be visualized through the use of the lattice package in R.

Running the code above produces this beauty:

What exactly is the image above showing us?  Well, it is the parameter space where the optimization procedure is attempting to find values of \alpha and \beta that lead to the lowest values for the negative log likelihod.  We can see that certain combinations of values for \alpha and \beta produce smaller negative log likelihoods.

Therefore, the optimization procedure attempts to find the minimum region of this parameter space, and then selects this combination of values for the parameters being estimated (\alpha and \beta).  This same principle applies to higher n-dimensional optimization problems using MLE.

 

4. HOW TO EXTRACT A THRESHOLD FROM A PSYCHOMETRIC FUNCTION?

Most psychologists are interested in fitting a psychometric function to data from an observer to extract a threshold.  What is a threshold you ask?  Here I provide a very general definition of threshold, which is the minimum intensity of some stimulus x required for that stimulus to be just detectable or discriminable. The value that corresponds to “just” detectable or discriminable is usually chosen by the experimenter to correspond to a value where the observer is able to make a correct judgment about the stimulus just over chance performance (usually a value in the range of 65 – 80%), but is also based upon factors such as the task and function fit.  For example, it is common practice to use the estimated value of \alpha from the function fit to directly estimate the threshold.  The threshold value associated with this level of performance is typically called the Just Noticeable Difference (JND).  Another common metric of performance is known as the Point of Subjective Equality (PSE), and is defined as the stimulus intensity that leads to the probability of a correct response to be 50%.

Here, we extract the 75% threshold from our fitted function above and plot the results:

The code above produces the following plot:

Let us print the threshold value from our code above:

We can see in both the figure and from the printed value above that the luminance level that leads to 75% correct performance is x = 0.45. The threshold is represented in the figure by the dotted black line.  Now it would be time to run additional sessions and observers to get some estimates of variance. In a future post, we will explore how to apply bootstrap procedures to estimate confidence intervals of these psychometric functions.

 

If you found this post useful, please let me know at one of the social links below.

Thank you for reading and until next time, verify first, trust later, and remain skeptical always.

 

References:

  1. Kingdom, F. A., & Prins, N. (2016). Psychophysics: a practical introduction. Academic Press.
  2. Knoblauch, K., & Maloney, L. T. (2012). Modeling psychophysical data in R (Vol. 32). Springer Science & Business Media.
  3. Prins, N. (2012). The psychometric function: The lapse rate revisited. Journal of Vision, 12(6), 25-25.
  4. Weibull, W. (1951). A statistical distribution function of wide applicability. Journal of applied mechanics, 18(3), 293-297.
  5. Wichmann, F. A., & Hill, N. J. (2001). The psychometric function: I. Fitting, sampling, and goodness of fit. Attention, Perception, & Psychophysics, 63(8), 1293-1313.

 

 

Posted in Tutorials and tagged , , , , .