STAN is a new system for Bayesian inference using MCMC. It's like BUGS or JAGS, but based on a better MCMC algorithm. It's certainly overkill for something as simple as psychometric function fitting, but the code below might provide a good starting point for more complicated models. For an intro to psychometric function fitting in general, see Wichmann and Hill (2011), or this recent book by Ken Knoblauch and Larry Maloney (2012). Psychometric function fitting is essentially fitting a GLM with a different parameterisation, so for more advanced models you might want to use the GLM parameterisation instead. Here I'm using the traditional psychophysics parameterisation. Also, I'm not modelling performance (% correct) but response probability. The data could come from example from an experiment in which you vary the orientation of a pattern and the subject has to judge whether the pattern is tilted to the left or to the right: the point of subjective verticality would be the physical orientation that induces 50% “right-tilted” judgment. The traditional model is \[ p(y=1|x) = \Phi((x-\mu)/\sigma) \] where \( y \) is the response, \( x \) is the stimulus (e.g, orientation), and \( \mu \) and \( \sigma \) model location and (inverse) slope respectively. In the orientation experiment, \( \mu \) is the point of subjective verticality. As Wichmann & Hill (2001) point out, every once in a while, the observer gets distracted and pushes buttons at random. The probability that this happens is called the lapse rate, and including it in the model leads to: \[ p(y=1|x) =(1-\lambda)\Phi((x-\mu)/\sigma)+\lambda/2 \]
In total we have three parameters to estimate. I put in priors that should behave OK if the stimulus range is approximately scaled between -3 and 3 like it is here. \( \lambda \) is assumed to be < 5%. You really should adjust the priors to the actual setting if you ever use this code.
STAN needs a model file that describing the model. Here it is, see the STAN documentation for an explanation of the various parts:
data {
int<lower=1> N;
vector[N] x;
int<lower=0,upper=1> y[N];
}
parameters {
real<lower=0> sigma;
real mu;
real<lower=0,upper=0.05> lambda;
}
model {
vector[N] eta;
mu ~ normal(0,3);
eta <- (x-mu)/sigma;
for (i in 1:N)
{
y[i] ~ bernoulli((1-lambda)*Phi(eta[i])+lambda/2);
}
}
I've stored the model file in “psychfun_lapse.stan”.
We generate some random data in R:
require(plyr)
## Loading required package: plyr
x <- rep(seq(-3, 3, l = 8), each = 60) #8 stim. levels with 60 trials each
# The psych. function
pf <- function(x, mu, sigma, lambda = 0) (1 - lambda) * pnorm((x - mu)/sigma) +
lambda/2
y <- (runif(length(x)) < pf(x, 0.5, 1.6, 0.03)) + 0
df <- data.frame(x = x, y = y)
ggplot(df, aes(x, y)) + stat_summary(fun.data = "mean_cl_boot") + labs(x = "Stim. level",
y = "Response frequency")
The bars show 95% bootstrap intervals on response probability.
Now we call STAN:
require(rstan)
require(reshape2)
res <- stan("psychfun_lapse.stan", data = list(x = x, y = y, N = length(x)),
iter = 1000, chain = 4, thin = 10)
samples <- extract(res)
smpl <- mutate(as.data.frame(samples), index = 1:length(mu))
ggplot(melt(smpl, id.var = "index"), aes(index, value)) + geom_point() + facet_wrap(~variable,
scale = "free_y")
STAN needs to compile the model into C++ code. On my computer this is the really slow bit, but it only needs to be done once. The plot shows the MCMC samples generated by STAN. They look OK and for a model with three parameters STAN shouldn't run into too many problems, but this needs to be checked always. See the STAN manual for more on convergence checking.
We now replot the data, overlaying the psychometric functions corresponding to posterior samples:
xv <- seq(-3, 3, l = 100)
dpost <- adply(smpl, 1, function(s) data.frame(xv = xv, pfun = pf(xv, s$mu,
s$sigma, s$lambda), index = s$index))
ggplot(df, aes(x, y)) + stat_summary(fun.data = "mean_cl_boot") + geom_line(data = dpost,
aes(xv, pfun, group = as.factor(index)), alpha = 0.03, col = "red") + labs(x = "Stim. level",
y = "Response")
And here are posterior distributions for the parameters:
ggplot(subset(melt(smpl, id.var = "index"), variable != "lp__"), aes(value)) +
geom_density() + facet_wrap(~variable, scale = "free")
True values were \( \sigma=1.6 \), \( \mu=0.5 \) and \( \lambda=.03 \) Fitting performance (% correct) data is just a small change to the model.
Knoblauch, K. and Maloney, L. T. (2012). Modeling Psychophysical Data in R (Use R!). Springer, 2012 edition. Wichmann, F. A. and Hill, N. J. (2001). The psychometric function: I. fitting, sampling, and goodness of fit. Perception and Psychophysics, 63(8):1293-1313.