Here we use *quickpsy* to illustrate the problem of fitting psychometric functions when participants make lapses.

Wichmann and Hill (2001) showed that when lapses occur and are not taken into account in the fitting process, the estimated threshold could be biased. Consider the following toy example

```
library(quickpsy)
x <- seq(0, 420, 60)
k <- c(0, 0, 4, 18, 20, 20, 19, 20)
dat <- tibble(x, k, n = 20)
fitWithoutLapses <- quickpsy(dat, x, k, n, prob = .75)
curvesWithoutLapses <- fitWithoutLapses$curves %>%
mutate(cond = 'Without Lapses')
pWithout <- ggplot()+
geom_point(data = fitWithoutLapses$averages, aes(x = x, y = prob)) +
geom_line(data = curvesWithoutLapses,
aes(x = x, y = y, color = cond)) +
geom_linerange(data = fitWithoutLapses$thresholds,
aes(x = thre, ymin = 0, ymax = prob), lty =2)
pWithout
```

As lapses do occur, the fit doesn’t seem good and the threshold biased. The problem is not that the likelihood is not correctly maximized; this *weird* fit is the one with the maximum likelihood, given the lapses (Wichmann and Hill, 2001). To reassure that the optimization process is not reaching a local maximum, we repeate the fit using differential evolution algoriths (function `DEoptim`

)

```
fitWithoutLapses2 <- quickpsy(dat, x, k, n, prob = .75, optimization = 'DE',
parini = list(c(1, 500), c(1, 500), c(0,1)),
bootstrap = 'none')
curvesWithoutLapses2 <- fitWithoutLapses2$curves %>% mutate(cond = 'Without Lapses')
pWithout2 <- ggplot()+
geom_point(data = fitWithoutLapses2$averages, aes(x = x, y = prob)) +
geom_line(data = curvesWithoutLapses,
aes(x = x, y = y, color = cond)) +
geom_linerange(data = fitWithoutLapses2$thresholds,
aes(x = thre, ymin = 0, ymax = prob), lty =2)
pWithout2
```

and obtain the the same result.

Wichmann and Hill (2001) reported that allowing the upper assymptote to vary — that is, fitting the lapses — eliminates the bias for the threshold, which seems to work in this example

```
fitWithLapses <- quickpsy(dat, x, k, n, prob = .75, lapses = T)
curvesWithLapses <- fitWithLapses$curves %>% mutate(cond = 'With Lapses')
pWithoutWith <- pWithout +
geom_line(data = curvesWithLapses,
aes(x = x, y = y, color = cond)) +
geom_linerange(data = fitWithLapses$thresholds,
aes(x = thre, ymin = 0, ymax = prob), lty =2)
pWithoutWith
```

```
data(qpdat)
fitWithoutLapses <- quickpsy(qpdat, phase, resp,grouping = c(ecc, interval ,participant))
fitWithoutLapses %>% plot()
```

Due to lapses, some fits, like the one in the upper left panel for the blue dots, do not seem good, so let’s allow the lapses to vary

`fitWithLapses <- quickpsy(qpdat, phase, resp, grouping = .(ecc, interval ,participant), lapses = T) `

`## Error: As y-predictions are not within (0,1), bootstrap should be 'nonparametric'`

We obtain an error in the bootstrap. To understand the error we plot the results without bootstraping

```
fitWithLapses <- quickpsy(qpdat, phase, resp, grouping = c(ecc, interval ,participant), lapses = T, bootstrap = 'none')
ggplot() +
facet_grid(interval~participant) +
geom_point(data = fitWithLapses$averages,
aes(x = phase, y = prob, color = factor(ecc))) +
geom_line(data = fitWithLapses$curves,
aes(x = x, y = y, color = factor(ecc))) +
ylim(0, 1.25) # we include ggplot ylim to be able to visualize what is going on
```

The problem is that, as we did not constraint the values of the upper asymptote, the lapse rate is eventually negative and consequently the fitted curves go through probabilities larger than one, which causes that the parametric bootstrap samples cannot be calculated. For this reason, *quickpsy* indicates that the bootstrap should be non-parametric. We avoid the error by performing a non-parametric bootstrap

```
fitWithLapses <- quickpsy(qpdat, phase, resp,
grouping = c(ecc, interval ,participant), lapses = TRUE,
bootstrap = 'nonparametric')
ggplot() +
facet_grid(interval~participant) +
geom_point(data = fitWithLapses$averages,
aes(x = phase, y = prob, color = factor(ecc))) +
geom_line(data = fitWithLapses$curves,
aes(x = x, y = y, color = factor(ecc))) +
ylim(0, 1.25)
```

Another, possibly better, solution is to perform a parametric bootstrap but constraining the values of the upper asymptote

```
fitWithLapses <- quickpsy(qpdat, phase, resp,
grouping = c(ecc, interval ,participant), lapses = T,
parini = list(c(20, 1000), c(20, 1000), c(0, .1)))
fitWithLapses %>% plot()
```

By using the same simulations as in Wichmann and Hill (2001), here we reproduce the bias of the estimated thresholds due to response lapses (see the original paper for details).

Wichmann and Hill used 7 sampling squemes defined as follows

```
parweibull <- c(10, 3)
create_xs <- function(i, f) data.frame(squeme = i, y = f, x = inv_weibull_fun(f, parweibull))
s <- list()
s[[1]] <- create_xs(1, c(.3, .4, .48, .52, .6, .7))
s[[2]] <- create_xs(2, c(.1, .3, .4, .6, .7, .9))
s[[3]] <- create_xs(3, c(.3, .44, .7, .8, .9, .98))
s[[4]] <- create_xs(4, c(.1, .2, .3, .4, .5, .6))
s[[5]] <- create_xs(5, c(.08, .18, .28, .7, .85, .99))
s[[6]] <- create_xs(6, c(.3, .4, .5, .6, .7, .99))
s[[7]] <- create_xs(7, c(.34, .44, .54, .8, .9, .98))
s <- do.call('rbind', s)
ggplot(s, aes(x = x, y = squeme, color = factor(squeme))) +
geom_point() + geom_line() +
labs(color = 'Sampling scheme') + theme(legend.position = 'top')
```

They generated parametric bootstrap samples of data coming from a psychometric function with the shape of a weibull function with parameters \(10\) and \(3\), \(guess=0.5\) and variable \(\lambda\) (directly related to the lapses).

For each value of the independent variable \(x\), we simulated 160 trials (in the original paper 20, 40 and 80 trials were also used).

```
create_sim_dat <- function(d) {
psychometric_fun <- create_psy_fun(weibull_fun, .5, d$lambda)
ypred <- psychometric_fun(d$x, parweibull)
k <- rbinom(length(d$x), d$n, ypred)
data.frame(x = d$x, k = k, n = d$n , y = k/d$n)
}
simdat <- merge(s, expand.grid(n = 160, sample = 1:10, lambda = seq(0,.05, .01))) %>%
group_by(squeme, n, sample, lambda) %>% do(create_sim_dat(.))
# we use only 10 samples to increase decrease the computation time
```

To illustrate the simulated data, we plot the first sample for each sampling scheme (columns) and each value of \(\lambda\) (rows)

```
p <- ggplot(simdat %>% filter(sample == 1)) +
facet_grid(lambda ~ squeme) + geom_point(aes(x = x, y = y))
p
```

We use quickpsy to fit a weibull function to each condition using \(\lambda = 0\)

```
fit <- quickpsy(simdat, x, k, n, within = .(squeme, lambda, sample),
fun = weibull_fun, bootstrap = 'none', guess =.5, lapses = 0)
# we use 'none' for the bootstap because we are doing an explicit bootstrap
# (we are fitting by sample), that is, we don't need to bootstrap each of our samples
```

and for the sake of illustration show the fitted psychometric function for the first sample

```
p + geom_line(data = fit$curves %>% filter(sample == 1),
aes(x = x, y =y))
```

```
# we use the function curves instead of fit$curves because we want to
# define the range of the curves manually
```

Replicating Wichmann and Hill, we show that the thresholds are overstimated as \(\lambda\) increases (they showed that the slope is also biased)

```
thre <- fit$thresholds %>% group_by(squeme, lambda) %>%
summarise (threshold = mean(thre))
real_threshold <- inv_weibull_fun((.75 - .5) / (1 - .5 - 0), parweibull)
ggplot(thre) +
geom_point(aes(x = lambda, y = threshold, color = factor(squeme)), size = 4) +
geom_hline(yintercept = real_threshold, lty = 2)
```

Wichmann and Hill (2001) reported that allowing \(\lambda\) to vary within a given small window eliminates the bias for the threshold, but consistently with Prins (2012) we were not able to replicate this finding

```
fit_lapses <- quickpsy(simdat, x, k, n, within = .(squeme, lambda, sample),
fun = weibull_fun, bootstrap = 'none', guess =.5, lapses = T,
parini = list(c(1,30), c(1,10), c(0,.05))) # we introduce the min and max value for each parameter
thre_lapses <- fit_lapses$thresholds %>% group_by(squeme, lambda) %>%
summarise (threshold = mean(thre))
ggplot(thre_lapses) +
geom_point(aes(x = lambda, y = threshold, color = factor(squeme)), size = 4) +
geom_hline(yintercept = real_threshold, lty = 2)
```

So, fitting the lapses does not always eliminate the bias in threshold estimation.

Prins, N. (2012). The psychometric function: The lapse rate revisited. Journal of Vision, 12(6), 25–25.

Wichmann, F. A., & Hill, N. J. (2001). The psychometric function: I. Fitting, sampling, and goodness of fit. Perception and Psychophysics, 63(8), 1293–1313.