So the data I have is whether a subject has performed a test correctly, or incorrectly. They have to match choose which of a pair of stimuli matches one they have memorized, and this gets harder and harder to figure out (represented by the "angle" parameter approaching zero) and at that point they get the test right 50% of the time.

So when "angle" is larger, they get the test right 100% of the time (so "correct" = 1, with no variance), but as "angle" approaches 0, "correct" is either 1 or 0, with a mean of 0.5, and a variance of 0.25. As "angle" changes, the subjects responding moves from getting the answer correct, to getting it right 50/50. Each subject is only tested at a given angle once.

I then give the subjects a intervention which I believe will make them worse at this task, hence shifting this curve to the right.

Small Edit: The above is what I want to know: Has the intervention made the subjects better or worse at the task. Because of ethical and financial reasons, each subject will experience the stimulus only once. Perhaps a mixed-effect model would be better. But right now I'm binning the data from all the subjects together. Ultimately, I'm not sure that is a huge sin, because this is not a within subjects experiments: some subjects get the intervention, and some do not.

This data is also like a logistic regression, but no matter how small "angle" becomes, the data will never become all zeros. Hence, I believe I need to run a binomial regression with a custom link function.

The curve should take on the form of:

y = 0.5 + (0.5 * exp(Beta*X) ) /( exp(Beta*X) + 1 ) 

The inverse of this (i.e. the link function) is then

log( (1-2*y)/(2*y-2) )

The deriviate of the inverse link function with respect to eta is

 0.5*exp(Beta*X) /( ( exp(Beta*X) + 1 ) * ( exp(Beta*X) + 1  ) )

So I have written my custom link function as

bills <- function() {
  ## link
  linkfun <- function(y) { 
    log( (1-2*y)/(2*y-2) )

  ## inverse link
  linkinv <- function(eta) {
    tmp = exp(eta)
    0.5 + (0.5 * tmp) /( tmp + 1 ) 

  ## derivative of invlink wrt eta
  mu.eta <- function(eta) {
    tmp = exp(eta)
    0.5*tmp /( ( tmp + 1 ) * ( tmp + 1  ) )

  valideta <- function(eta) TRUE
  link <- "bills"
  structure(list(linkfun = linkfun, linkinv = linkinv,
                 mu.eta = mu.eta, valideta = valideta, 
                 name = link),  class = "link-glm")
vv <- bills()

In order to test that this function behaves, I have written the following script to make data like my real data:

n_reps = 10
differences = c(5,10,15,20)

output <- matrix(nrow = n_reps, ncol = length(differences))

d <- 1

for(dif in differences) {
  for(n in 1:n_reps) {
    sample_size <- 5
    dat_c <-data.frame(angle = rep(seq(0, 90, by=10), times = sample_size), treated = rep(0, times=sample_size*10))
    dat_d <-data.frame(angle = rep(seq(0, 90, by=10), times = sample_size), treated = rep(1, times=sample_size*10))

    my_sigmoid <- function(x, v50) {
      return( 0.5+(1-0.5)/(1+exp((v50-x)/5)) ) 

    control_v50 = 45
    drug_v50 = control_v50 + dif
    dat_c$p <- my_sigmoid(dat_c$angle, control_v50)
    dat_d$p <- my_sigmoid(dat_d$angle, drug_v50)

    dat_c$correct <- as.integer(runif(nrow(dat_c))<dat_c$p)
    dat_d$correct <- as.integer(runif(nrow(dat_d))<dat_d$p)

    dat = merge(dat_c, dat_d, all=TRUE)

    mdl_simp <- glm(correct ~ angle, family = binomial(link=vv), data = dat, start=c(-3,.1))
    mdl_full <- glm(correct ~ angle + treated, family = binomial(link=vv), data = dat, start=c(-3,.1, .5))

    a <- anova(mdl_simp, mdl_full, test="Chisq")

    p <- a$"Pr(>Chi)"[2]
    if ( is.na(p) ) {
      p <- 1

    output[n, d] <- p
  d <- d+1

However, when I run that script, I get constant errors like:

step size truncated: out of bounds
glm.fit: algorithm stopped at boundary value
glm.fit: fitted probabilities numerically 0 or 1 occurred
Error: inner loop 1; cannot correct step size

I have double and triple checked my math, but I can't see any problem. The link function passes the typical suggested tests i.e. that running linkinv(linkinv(x)) returns x, and that mu.eta returns values that match a numerical calculation of the gradient of linkfunction.

But perhaps I am doing something else wrong? Any ideas appreciated.

  • $\begingroup$ Wouldn't it be easier to just model what you call angle with a non-negative distribution, like a gamma GLM, or beta regression? I would go with the latter, although I can't tell from your question whether angle is in $[0, 1]$ or $[0, \infty)$... Alternatively, if this 'angle' is hard to model, it might be easier to just model the ratio $\text{correct}:\text{incorrect}$. Can you elaborate why you want to model your data as you propose? $\endgroup$ – Frans Rodenburg Jul 8 at 1:21
  • $\begingroup$ @FransRodenburg angle is in degrees. It's [0, 90]. But I don't understand that rest of your statement. Angle isn't hard to model. It's set by the experimenter. What I'm trying to model is how angle alters the probability of the subject getting the answer correct. Why don't I just look at the ratio? Well because I generally try to avoid anything like that. I'm not actually sure what the correct terminology is, but it masks uncertainty. I.e. if a subject does 2 trials and gets half of them right, is really a very different result that if a subject does 50 trials and gets 25 of them right. $\endgroup$ – Bill Connelly Jul 8 at 2:37
  • $\begingroup$ @Alex I don't think so. Even with sample size <- 100 I get the same errors $\endgroup$ – Bill Connelly Jul 8 at 2:39
  • 2
    $\begingroup$ I'm trying to help you add the necessary details to the question so that I or someone else can write a good answer. I understand that you have thought about this experimental design a lot more than me, but by simply assuring me I don't need to know more, I don't suddenly understand the question better. "What I'm trying to model is how angle alters the probability of the subject getting the answer correct." sounds to me like the outcome is correct/incorrect and angle is an explanatory variable. If that is the case, wouldn't it make more sense to model the problem as a binomial GLM? $\endgroup$ – Frans Rodenburg Jul 8 at 3:18
  • 1
    $\begingroup$ It seems you are talking about English, instead of statistics or math. I am not interesting on English. $\endgroup$ – user158565 Jul 8 at 3:43