I’d like to fit a psychometric function to some binomial data. In the past I’ve always used vanilla logistic regression, but if there are random errors this can bias the estimates of the slope and intercept of the curve (see for example Wichman & Hill, 2001). Thus I’d like to include the possibility that there’s a mixture of random/lapse responses (e.g. sampled from a single binomial error distribution) and on-task responses (as in normal logistic regression). The “classical” functional form of the model is that the response probability is a saturating logistic function with asymptotes not exactly equal to 0 or 1.
My question is whether I can express this in brms other than through the mechanism for specifying a “generic” non-linear function, and (if not) whether that comes with a big performance penalty. Of course I could write the stan code by hand, but I’d really like to be able to take advantage of brms’s ability to concisely specify random effects as well.
To be concrete, here’s an example dataset and three brms models: first is vanilla logistic regression, second is the manual nonlinear specification of logistic regression (which on my system is about 50% slower to evaluate, which I can bear), and the third is my first attempt at a lapsing logistic regression.
Simulated data
library(tidyverse)
library(brms)
options(mc.cores = parallel::detectCores())
dat <- cross_df(list(x = seq(0, 10),
rep = seq(1, 5),
subj = seq(1, 10))) %>%
mutate(theta = 0.02 + 0.95*plogis(1 * (x-5)),
y = rbernoulli(n=length(theta), p=theta))
Vanilla logistic regression
fit1 = brm(y ~ 1 + x, data=dat, family=bernoulli())
Nonlinear formula logistic regression
fit2 = brm(bf(y ~ 1 / (1+exp(-(eta))),
eta ~ 1+x,
family=bernoulli(link="identity"),
nl=TRUE),
data=dat,
prior = set_prior("student_t(3,0,10)", class="b", nlpar="eta"))
Lapsing logistic regression
fit3 = brm(bf(y ~ lapselow + (1-lapsehigh-lapselow) * 1 / (1+exp(-(eta))),
eta ~ 1+x,
lapselow ~ 1,
lapsehigh ~ 1,
family=bernoulli(link="identity"),
nl=TRUE),
data=dat,
inits="0",
control=list(adapt_delta=0.99),
prior = c(set_prior("student_t(3,0,10)", class="b", nlpar="eta"),
set_prior("uniform(0,0.1)", class="b", nlpar="lapselow"),
set_prior("uniform(0,0.1)", class="b", nlpar="lapsehigh")))
As my attempt to increase adapt_delta
might suggest, I’m getting a ton of divergent transitions with the last model (nearly 25%). So I’m afraid that it’s going to be hard to make this work without digging into the stan code (or doing some kind of hack to mimic stan’s transformation of bounded parameters)
References
Wichmann, F. A., & Hill, N. J. (2001). The psychometric function: I. Fitting, sampling, and goodness of fit. Perception and Psychophysics, 63(8), 1293–1313.