How to fit a generalized logistic function?

I'm trying to fit models that predict probabilities: $h(X,B) \rightarrow (0,1)$

It struck me that for for a great many cases logistic regression seems like a bad link function, I expect an S-curve, but I also expect bounds on the predicted probability to be much larger than 0, and much less than 1.

For example, let's say I want to predict the likelihood of a first bankruptcy within 10 years, and I have some feature $x$ (e.g. number of years a loan in is default). Here I might expect, that for all x, $p>>0$ because many people go bankrupt for reasons unrelated to loan defaults, and that $p<<1$ because beyond some point more years in loan default doesn't increase the risk of new bankruptcy, because the number of years a loan in is default doesn't completely determine probability, you'd expect $p<<1$.

Fitting such a probability function with logistic regression leads to a very poor fit:

The target function above is a (special case) of "generalized logistic function". In this case:

$${prob} = p_{min} + (p_{max}-p_{min})*logistic(X,B)$$

with logistic inputs $X$ and coefficients $B$. Are there good ways to optimize $p_{min}$, $p_{max}$ and $B$ for this type of function such that a measure of regression error is minimized (log error, squared error, ...), or a similar 'S'-shaped model that's likely to provide a better fit than a logistic curve?

• The problem you're trying to solve sounds a lot like omitted variable bias.
– Sycorax
Commented Mar 8, 2017 at 18:35
• Of course -- there are always omitted variables in most practical problems. There's what we'd like to measure, and what we can actually measure. Only the latter is available as data to train the model. Commented Mar 8, 2017 at 19:10
• Do you really want a ceiling for estimated probabilities that is less than 1 and a floor greater than 0 or is it possible that you just want more of the estimated probabilities to fall further away from 0 and 1? If the latter it may be that you just need to add additional terms to the linear predictor. Commented Mar 8, 2017 at 19:18
• Yes -- something a little softer than hard limits would be preferable. Mostly I'm trying to find an alternative to LR that could reasonably reduce the error in fit (with the constraint that the output is in 0..1 and that the function is kinda S-shaped). For my data, the insight for me has been that, for the lowest probability portion of the sample, the best prediction is >>0 and for the highest probability portion, the best prediction is <<1. Where this is true, LR looks much more like a linear fit than an S-shape. Commented Mar 8, 2017 at 19:28
• Sycorax is absolutely correct. Commented Mar 8, 2017 at 20:43

Given the binary response $$y_i$$ and the covariate $$x_i$$, $$i=1,2,\dots,n$$, the likelihood for your model is $$L(\beta_0,\beta_1,p_\text{min},p_\text{max})=\prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i}$$ where each $$p_i=p_\text{min} + (p_\text{max} - p_\text{min})\frac1{1+\exp(-(\beta_0 + \beta_1 x_i)}.$$ Just write a function computing the log of this an apply some general purpose optimization algorithm to maximise this numerically with respect to the four parameters. For example, in R do:

# the log likelihood
loglik <- function(par,y,x) {
beta0 <- par[1]
beta1 <- par[2]
pmin <- par[3]
pmax <- par[4]
p <- pmin + (pmax - pmin)*plogis(beta0 + beta1*x)
sum(dbinom(y, size=1, prob=p, log=TRUE))
}
# simulated data
x <- seq(-10,10,len=1000)
y <- rbinom(n=length(x),size=1,prob=.2 + .6*plogis(.5*x))
# fit the model
optim(c(0, 0.5 ,.1, .9), loglik, control=list(fnscale=-1), y=y,x=x, lower=c(-Inf,-Inf,0,0),upper=c(Inf,Inf,1,1))


Note that to test for evidence of a lower plateau at $$p_\text{min}$$ in your data, your $$H_0:p_\text{min}=0$$ is at the boundary of the parameter space and the approximate/asymptotic distribution of $$2(\log L(\hat\theta_1)-\log L(\hat\theta_0))$$ is going to be a mixture of chi-square distributions with 1 and 0 degrees of freedom, see Self, S. G. & Liang, K. Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions J. Amer. Statist. Assoc., 1987, 82, 605-610.

In the simpler case that there is only one plateau (so $$p_\text{max}=1$$ or $$p_\text{min}=0$$) the model is equivalent to a zero-inflated binary regression model that can be fitted with e.g. the glmmTMB R-package.

• This a good solution -- I had a similar idea and implemented (within Python) on squared loss (log loss seems better). One of the optimizers I tried for this (on squared loss) didn't seem to converge on a useful answer. The other worked fine. What would be better is something that does have such hard limits, p_min,p_max, but some kind of squish factor at the tails. For p_min. Are you suggesting to calculate it directly from the data and set it? (e.g. we're 98% sure p_min>0.1, so set p_min=0.1) (and skip the limits in the optimization). Commented Mar 9, 2017 at 19:07
• If you have precise prior information about $p_\text{min}$ you may perhaps treat is as known in the above model. If you only have vague information you could incorporate this into a prior and do Bayesian inference (maybe using MCMC). I also believe using maximum likelihood (or Bayesian inference) here will be more efficient than a machine learning approach minimising some arbitrary loss function. Doesn't using a squared loss correspond to assuming that the variance of the binary response is homoscedastic? This will give you some loss in efficiency because the true variance is $p_i(1-p_i)$. Commented Mar 10, 2017 at 8:33