Using a custom link function with rstanarm models

[tomwallis]

Summary:

I would like to use a custom link function with an rstanarm glm, but this throws an error. The link function works fine for base R glm.

Description:

The psyphy package contains custom link functions for psychophysical data. Rstanarm throws an error if I try to fit this model. The documentation of rstanarm says that family should work as in glm.

Reproducible Steps:

library(rstanarm)
library(psyphy)

dat <- mtcars

prior <- student_t(df = 3, location = 0, scale = 2.5)

fit1 <- stan_glm(am ~ 1, data = dat, family = binomial(link=mafc.logit(2)), prior_intercept = prior)

Error in stan_glm.fit(x = X, y = Y, weights = weights, offset = offset,  : 
  'link' must be one of logit, probit, cauchit, log, cloglog

fit2 <- glm(am ~ 1, data = dat, family = binomial(link=mafc.logit(2)))

fit2 works.

RStanARM Version:

‘2.9.0.3’

R Version:

‘3.2.3’

Operating System:

OS X 10.11.3

[jgabry]

In order for new link functions to work the Stan programs used to estimate rstanarm models would need to be amended. In other words, right now if you use link=mafc.logit(...) there's no Stan code to do the necessary computations for this link function when estimating the model. It's possible to add new link functions, but not on-the-fly.

For an example of what's required to add a new link function, take a look at the Stan code in exec/bernoulli.stan. You'll see functions like linkinv_bern and ll_bern_lp, for the inverse-link function and log-likelihood of a Bernoulli model, respectively. Functions like these would need to be updated to accommodate the new link function.

[tomwallis]

Ah yes, of course. How silly of me not to realise that it would require a stan code backend. Thanks.

[jgabry]

If those link functions are commonly used and you think it would be good to have any of them in rstanarm then we would certainly consider implementing them. Feel free to open a new issue if you want to suggest that we implement any particular link functions. Also, if you wanted to contribute the code yourself then we can also walk you through the process of adding something like that. It's not too difficult once you know where everything goes.

[tomwallis]

Hi @jgabry: those link functions would be very useful for people fitting data from forced-choice experiments, as often used in psychophysics and neuroscience.

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). This is desirable from theoretical considerations about the task, even though real data could fall below 0.5 (due to measurement noise, or the observer mistakenly flipping responses). There is a useful book that describes approaches to fitting this data in R by Knoblauch and Maloney.

I would be very happy to try contributing the code myself, but my time is quite limited right now (on part-time parental leave) so I'm not sure when I will get to it. Nevertheless, please do outline where things would need to go, and I will give it a shot.