Binomial regression asymptotes

Question

Binomial logistic regression has upper and lower asymptotes of 1 and 0 respectively. However, accuracy data (just as an example) may have upper and lower asymptotes vastly different to 1 and/or 0. I can see three potential solutions to this:

1. Don't worry about it if you are getting good fits within the area of interest. If you aren't getting good fits then:
2. Transform the data so that the minimum and maximum number of correct responses in the sample give proportions of 0 and 1 (instead of say 0 and 0.15).
   or
3. Use non-linear regression so that you can either specify the asymptotes or have the fitter do it for you.

It seems to me that options 1 & 2 would be preferred over option 3 largely for simplicity reasons, in which case option 3 is perhaps the better option because it can yield more information?

edit
Here's an example. Total possible correct for accuracy is 100, but the maximum accuracy in this case is ~ 15.

accuracy <- c(0, 0, 0, 0, 0, 1, 3, 5, 9, 13, 14, 15, 14, 15, 16, 
              15, 14, 14, 15)
x <- 1:length(accuracy)
glmx <- glm(cbind(accuracy, 100-accuracy) ~ x, family=binomial)
ndf <- data.frame(x=x)
ndf$fit <- predict(glmx, newdata=ndf, type="response")
plot(accuracy/100 ~ x)
with(ndf, lines(fit ~ x))  

Option 2 (as per comments and to clarify my meaning) would then be the model

glmx2 <- glm(cbind(accuracy, 16-accuracy) ~ x, family=binomial)

Option 3 (for completeness) would be something akin to:

fitnls <- nls(accuracy ~ upAsym + (y0 - upAsym)/(1 + 
           (x/midPoint)^slope), 
  start = list("upAsym" = max(accuracy), "y0" = 0, 
  "midPoint" = 10, "slope" = 5), 
  lower = list("upAsym" = 0, "y0" = 0, "midPoint" = 1, 
               "slope" = 0), 
  upper = list("upAsym" = 100, "y0" = 0, "midPoint" = 19, 
               hillslope = Inf), 
  control = nls.control(warnOnly = TRUE, maxiter=1000),
  algorithm = "port")

Answer 1

Interesting question. A possibility that comes to my mind is including an additional parameter $p\in[0,1]$ in order to control the upper bound of the 'link' function.

Let $\{{\bf x}_j,y_j,n_j\}$, $j=1,...,n$ be independent observations, where $y_j\sim \text{Binomial}\{n_i,pF({\bf x}_j^T\beta)\}$, $p\in[0,1]$, ${\bf x}_j=(1,x_{j1}, \dotsc ,x_{jk})^T$ is a vector of explanatory variables, $\beta=(\beta_0,\dotsc,\beta_k)$ is a vector of regression coefficients and $F^{-1}$ is the link function. Then the likelihood function is given by

$${\mathcal L}(\beta,p) \propto \prod_{j=1}^n p^{y_j}F({\bf x}_j^T\beta)^{y_j}[1-pF({\bf x}_j^T\beta)]^{n_j-y_j}$$

The next step is to choose a link, say the logistic distribution and find the corresponding MLE of $(\beta,p)$.

Consider the following simulated toy example using a dose-response model with $(\beta_0,\beta_1,p)=(0.5,0.5,0.25)$ and $n=31$

dose = seq(-15, 15, 1)
a = 0.5
b = 0.5
n=length(dose)
sim = rep(0, n)
for(i in 1:n) sim[i] = rbinom(1, 100, 0.25*plogis(a+b*dose[i]))

plot(dose, sim/100)

lp = function(par){
if(par[3]>0& par[3]<1) return(-(n*mean(sim)*log(par[3]) +  
   sum(sim*log(plogis(par[1]+par[2]*dose)))  + 
    sum((100-sim)*log(1-par[3]*plogis(par[1]+par[2]*dose))) ))
else return(-Inf)
}

optim(c(0.5, 0.5, 0.25), lp)

One of the outcomes I got is $(\hat\beta_0,\hat\beta_1,\hat p)=( 0.4526650, 0.4589112, 0.2395564)$. Therefore it seems to be accurate. Of course, a more detailed exploration of this model would be necessary because including parameters in a binary regression model can be tricky and problems of identifiability or existence of the MLE may jump on the stage 1 2.

Edit

Given the edit (which changes the problem significantly), the method I proposed previously can be modified for fitting the data you have provided. Consider the model

$$\mbox{accuracy} = pF(x;\mu,\sigma),$$

where $F$ is the logistic CDF, $\mu$ is a location parameter, $\sigma$ is a scale parameter, and the parameter $p$ controls the height of the curve similarly as in the former model. This model can be fitted using Nonlinear Least Squares. The following R code shows how to do this for your data.

rm(list=ls())
y = c(0, 0, 0, 0, 0, 1, 3, 5, 9, 13, 14, 15, 14, 15, 16, 15, 14, 
      14, 15)/100
x = 1:length(y)
N = length(y)

plot(y ~ x)

Data = data.frame(x,y)

nls_fit = nls(y ~ p*plogis(x,m,s), Data, start = 
              list(m = 10, s = 1,  p = 0.2) )

lines(Data$x, predict(nls_fit), col = "red")

Answer 2

I would use the maximum of the X vector as the total possible number of successes. (This is a biased estimate of the true maximum number of successes, but it should work fairly well if you have enough data).

accuracy <- c(0, 0, 0, 0, 0, 1, 3, 5, 9, 13, 14, 15, 14, 15, 16, 
              15, 14, 14, 15)
x <- 1:length(accuracy)
glmx <- glm(cbind(accuracy, max(accuracy)-accuracy) ~ x, 
        family=binomial)
ndf <- data.frame(x=x)
ndf$fit <- predict(glmx, newdata=ndf, type="response")
plot(accuracy/max(accuracy) ~ x)
with(ndf, lines(fit ~ x))

This creates a plot that looks like:

Answer 3

Note that binomial regression is based on having a binary response for each individual case. each individual response has to be able to take one of two values. If there is some limit to the proportion then there must also have been some cases which could only take one value.

It sounds like you are not dealing with binary data but with data over a finite range. if this is the case, then beta regression sounds more appropriate. We can write the beta distribution as:

$$p(d_i|LU\mu_i\phi)=\frac{(d_i-L)^{\mu_i\phi-1}(U-d_i)^{(1-\mu_i)\phi-1}}{B(\mu_i\phi,(1-\mu_i)\phi)(U-L)^{\phi-1}}$$

You then set $g(\mu_i)=x_i^T\beta$ same as any link function which maps the interval $[L,U]$ into the reals. There is an R package which can be used to fit these models, though i think you need to know the bounds. If you do, then redefine the new variable $y_i=\frac{d_i-L}{U-L}$.