---
title: "Five-parameters logistic regression"
author: "Stéphane Laurent"
date: '2019-11-20'
tags: maths, statistics, R
rbloggers: yes
output:
  md_document:
    variant: markdown
    preserve_yaml: true
  html_document:
    highlight: kate
    keep_md: no
highlighter: pandoc-solarized
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.path = "./figures/5pl-", 
                      fig.width = 3.5, fig.align = "center",
                      attr.source = ".numberLines")
```

The five-parameters logistic curve is commonly defined by
$$
f(x) = A + \frac{D-A}{\Bigl(1+\exp\bigl(B(C-x)\bigr)\Bigr)^S}.
$$
Assuming $B>0$ and $S>0$, 

- $A$ is the value of the horizontal asymptote when $x \to -\infty$;

- $D$ is the value of the horizontal asymptote when $x \to +\infty$;

- $B$ describes how rapidly the curve makes its transition between the two
asymptotes;

- $C$ is a location parameter, which does not have a nice interpretation 
(except if $S=1$);

- $S$ describes the asymmetry of the curve (the curve is symmetric when $S=1$).

In the case when $S=1$, the parameter $C$ is the value of $x$ for which the 
corresponding value $f(x)$ is the midpoint between the two asymptotes; 
moreover, the curve has an inflection point at $x = C$.

In the general case, the value of $x$ for which the 
corresponding value $f(x)$ is the midpoint between the two asymptotes is 
$$
x_{\text{mid}} = C - \frac{\log\Bigl(2^{\frac{1}{S}}-1\Bigr)}{B}.
$$
It is obtained by solving 
$\Bigl(1+\exp\bigl(B(C-x)\bigr)\Bigr)^S = 2$.

```{r plot_5pl, collapse=TRUE, fig.height=2.5}
n <- 100
x <- seq(49, 60, length.out = n)
A <- 30; D <- 100; B <- 1; C <- 50; S <- 10
f <- function(x) A + (D-A) / (1 + exp(B*(C-x)))^S
y0 <- f(x) 
par(mar = c(4, 4, 0.5, 1))
plot(x, y0, type = "l", cex.axis = 0.5, ylab = "f(x)")
abline(v = C, col = "green", lty = "dashed")
( xmid <- C - log(2^(1/S) - 1)/B )
abline(v = xmid, col = "red", lwd = 2) 
abline(h = (A+D)/2, col = "red", lwd = 2)
```

Note that the inflection point of the curve is *not* the point correspoding 
to $x_{\text{mid}}$:

```{r inflection, fig.height=2.5}
library(numDeriv)
df <- grad(f, x)
par(mar = c(4, 4, 0.5, 1))
plot(x, df, type = "l", cex.axis = 0.5, ylab="f'(x)")
abline(v = xmid, col = "red", lwd = 2) 
```


In practice, we are often interested in estimating $x_{\text{mid}}$. 
So it is better to use this other parameterization of the five-parameters 
logistic curve:
$$
g(x) = 
A + \frac{D-A}{{\biggl(1+\exp\Bigl(\log\bigl(2^{\frac{1}{S}}-1\bigr) + B(x_{\text{mid}}-x)\Bigr)\biggr)}^S}
$$
because fitting this curve will directly give the estimate of $x_{\text{mid}}$ 
and its standard error.

Another advantage of this parameterization is that there is a way to get 
a good starting value of $x_{\text{mid}}$ when one wants to fit the 
five-parameters logistic regression model: 

```{r getInitialSSfpl}
getInitial1 <- function(x, y){
  s <- getInitial(y ~ SSfpl(x, A, D, xmid, inverseB),
             data = data.frame(x = x, y = y))
  c(A = s[["A"]], 
    B = 1/s[["inverseB"]], 
    xmid = s[["xmid"]], 
    D = s[["D"]], 
    S = 1)
}
```

I don't know how to get a good starting value for $S$, so I always take $1$.

Sometimes, `SSfpl` can fail. Here is another function which returns some 
starting values:

```{r getInitial2}
getInitial2 <- function(x, y){
  NAs <- union(which(is.na(x)), which(is.na(y)))
  if(length(NAs)){
    x <- x[-NAs]
    y <- y[-NAs]
  }
  low_init <- min(y)
  high_init <- max(y)
  minmax <- c(which(y == low_init), which(y == high_init))
  X <- cbind(1, x[-minmax])
  Y <- log((high_init-y[-minmax])/(y[-minmax]-low_init))
  fit <- lm.fit(x = X, y = Y)
  b_init <- fit$coefficients[[2]]
  xmid_init <- -fit$coefficients[[1]] / b_init
  if(b_init < 0){
    b_init <- -b_init
    A <- low_init
    D <- high_init
  }else{
    A <- high_init
    D <- low_init
  }
  c(A = A, B = b_init, xmid = xmid_init, D = D, S = 1)
}
```

Now we wrap these two functions into a single one:

```{r getInitial5pl}
getInitial5PL <- function(x, y){
  tryCatch({
    getInitial1(x, y)
  }, error = function(e){
    getInitial2(x, y)
  })
}
```

And finally we can write a function for the fitting: 

```{r fit5pl, message=FALSE}
library(minpack.lm)
fit5pl <- function(x, y){
  startingValues <- getInitial5PL(x, y)
  fit <- tryCatch({
    nlsLM(
      y ~ A + (D-A)/(1 + exp(log(2^(1/S)-1) + B*(xmid-x)))^S,
      data = data.frame(x = x, y = y),
      start = startingValues,
      lower = c(-Inf, 0, -Inf, -Inf, 0),
      control = nls.lm.control(maxiter = 1024, maxfev=10000))
  }, error = function(e){
    paste0("Failure of model fitting: ", e$message)
  })
  if(class(fit) == "nls" && fit[["convInfo"]][["isConv"]]){
    fit
  }else if(class(fit) == "nls" && !fit[["convInfo"]][["isConv"]]){
    "Convergence not achieved"
  }else{ # in this case, 'fit' is the error message
    fit
  }
}
```

Let's try it on a couple of simulated samples:

```{r simulations}
set.seed(666)
nsims <- 25
epsilon <- matrix(rnorm(nsims*n, 0, 5), nrow = nsims, ncol = n)
estimates <- matrix(NA_real_, nrow = nsims, ncol = 5)
colnames(estimates) <- c("A", "B", "xmid", "D", "S")
for(i in 1:nsims){
  fit <- fit5pl(x, y0 + epsilon[i,])
  if(class(fit) == "nls"){
    estimates[i, ] <- coef(fit)
  }else{
    estimates[i, ] <- c(NaN, NaN, NaN, NaN, NaN)
  }
}
summary(estimates)
```

The estimate of $x_{\text{mid}}$ is excellent. 
As you can see, the estimate of $S$ is sometimes much larger than the true 
value. Let's have a look at the worst case:

```{r worstcase, collapse=TRUE, fig.height=3}
i0 <- match(max(estimates[, "S"]), estimates[, "S"])
estimates[i0, ]
# sample
par(mar = c(4, 4, 0.5, 1))
plot(x, y0 + epsilon[i0, ], col = "yellow", cex.axis = 0.6)
# true curve
curve(A + (D-A)/(1 + exp(log(2^(1/S)-1) + B*(xmid-x)))^S, 
      add = TRUE, col = "red", lwd = 2)
# fitted curve
with(as.list(estimates[i0, ]), 
     curve(A + (D-A)/(1 + exp(log(2^(1/S)-1) + B*(xmid-x)))^S, 
           add = TRUE, col = "blue", lwd = 2, lty = "dashed")
)
```

Thus, while the estimate of $S$ is very far from the true value of $S$, 
the fitted curve correctly estimates the true curve. 
And in such cases, the standard error of the estimate of $S$ is big:

```{r fitworstcase}
fit <- fit5pl(x, y0 + epsilon[i0,])
summary(fit)
```

Note that `nlsLM` provides a test of the nullity of $S$. 
This is not interesting, whereas the equality $S = 1$ is of interest. So it is 
better to parametrize the logistic function with $L = \log(S)$ instead of $S$:
$$
h(x) = 
A + \frac{D-A}{{\biggl(1+\exp\Bigl(\log\bigl(2^{\exp(-L)}-1\bigr) + B(x_{\text{mid}}-x)\Bigr)\biggr)}^{\exp(L)}}.
$$
In this way we can get a test of $L = 0$, that is $S = 1$.