Custom binomial link functions for rstanarm

Summary:

I would like to add some custom binomial link functions for rstanarm; I'm seeking guidance to do that.

Description:

In psychophysics and neuroscience, it is common to use forced-choice experimental designs that result in a nonzero lower bound for binomially-distributed data. That is, the data will asymptote to something like 0.5 or 0.25 rather than 0 (see further example below). It would be useful if one could include custom versions of the binomial family link functions that show this behaviour. Specifically, the custom inverse link functions take the standard binomial inverse link functions (e.g. the logit) and squish them into a new range:

yhat = lower_bound + (1 - lower_bound - lapse_rate) * inv_logit(eta)

In a task with two possible responses (see below), the lower bound would be fixed at 1/2. In a task with four possible responses, the lower bound would be 1/4. Sometimes it is even desirable to have the lower and upper bound be free parameters (i.e. fitted to data), but I expect this would require a major recoding of the relevant rstanarm functions (because these free parameters would be outside the linear predictor of the GLM / GLMER). So I being able to set the lower bound (and maybe the lapse rate) to nonzero values would fulfil the bulk of use cases here. Obviously this changes how one can interpret the linear predictor coefficients (e.g. they won't be log odds any more for the logit link). Nevertheless, this would be a very useful addition to rstanarm for people like me. I would envisage it being added to stan_glm and stan_glmer at this stage.

Steps to add new link functions?

In this related issue, @jgabry mentioned he could outline the steps needed to add a new link function. I've forked the repo and started browsing bernoulli.stan, but I'm a little lost in the internals of the package. Some high level guidance would be useful. Also, what tests / docs should I aim to provide for a successful pull request?

Envisaged user interface

fit <- stan_glm(y ~ x1 + x2, data = dat, 
                 family = binomial(link = "mafc_logit", lower = 0.5, upper = 0.95))

The above example would use an mafc_logit link function (mAFC stands for "m alternative forced-choice") with lower set to 0.5 and upper set to 0.95, resulting in a logistic link function whose inverse link spans from 0.5 to 0.95.

I'm happy to take changes on this interface. I'm also not sure whether this envisaged use would require adjustment to the internal validate_family function.

More intuitive example:

Imagine I show you a screen that contains an image on one side and nothing (grey screen) on the other. Your task is to tell me which side the image is on ("left" or "right"). Over consecutive trials the image is presented randomly to the left or right. We can recode each response to be "correct" or "incorrect", and we usually assume that responses over many trials are independent and binomially distributed.

If I make the image very faint you won't be able to see it. On average, your performance (proportion correct) will be 0.5. These link functions allow one to fit a sigmoid that asymptotes to this theoretical chance performance level (0.5 in this case) rather than to zero (like a standard binomial link function). Though real data could fall below 0.5 (due to measurement noise, or the observer mistakenly flipping responses), new data would be expected to fall no lower than 0.5 on average. There is a useful book that describes approaches to fitting this data in R by Knoblauch and Maloney.

Start with exec/bernoulli.stan, then trace how stan_glm and stan_glmer handle binomial families and whether validate_family needs adjustment. Review the proposed interface for fixed lower and upper bounds, and identify the tests and documentation needed to cover the new link functions. Done should include guidance-supported behavior for the stated forced-choice use cases.