# Binomial regression asymptotes

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")  • Why is there an issue here? Logistic regression posits that the logit (log odds) of the probability has a linear relationship with the explanatory variables. The valid range of log odds is the entire set of real numbers; there is no possibility of going beyond them! – whuber Apr 28 '12 at 3:58 • Say for instance there is an upper asymptote of probability correct of 0.15. The regression is then poorly fitted to the data. I'll put up an example. – Matt Albrecht Apr 28 '12 at 4:27 • +1 great question. My instinct would be to use 16 as the maximum rather than 100 (cbind(accuracy, 16-accuracy)), but I'm concerned about whether it's mathematically justified. – David Robinson Apr 28 '12 at 4:45 ## 3 Answers 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},...,x_{jk})^T$is a vector of explanatory variables,$\beta=(\beta_0,...,\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")

• This is an interesting approach. What would be the advantages of using this method over a three parameter non-linear regression function? – Matt Albrecht Apr 28 '12 at 14:41
• @MattAlbrecht Thanks for the interest. I can see pros and cons of this approach. One of the pros is the interpretability of the approach, which is similar to the logit regression. On the other hand, a nonlinear regression function might be more flexible. In order to get a good estimation of $p$, it seems necessary to have a good experimental design which is not concentrated on the tails of the link function. I do not know if the model has been studied before. – user10525 Apr 28 '12 at 16:04
• The benefit would be the correct incorporation of the binomial variability. – Aniko Apr 30 '12 at 13:46
• @MattAlbrecht Note that this method restrict the shape of the fitted function to be sigmoidal and the parameter $p$ controls the height while the nonparametric method you are considering does not. BTW the estimated parameters with this model are $(\hat\mu,\hat\sigma,\hat{p})=(8.5121, 0.8987, 0.1483)$. – user10525 Oct 8 '12 at 18:54

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: 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}\$.

• Thanks for the response. This is made up data to simulate a T|F series totalling 100 dichotomous choices for each x point. So the limits are 0 correct or 100 correct but this particular case gets approx 15 correct. Using the betareg package... pacc <- accuracy/100 + 0.00001; b1 <- betareg(pacc ~ x) ... gives me the same looking regression as the binomial. Is this what you meant? Or are you suggesting imposing a post-hoc data based limit? In which case what distinguishes the beta from the binomial when both have been given post-hoc limits? – Matt Albrecht Apr 28 '12 at 14:36