3
$\begingroup$

I am trying to model vegetation spectral signatures (grey lines) using a two-part piecewise function (black dotted lines). In it, I am trying to use only a few points (red dots) to fit a linear (first part) and logistic function (second part).

Basically, the linear part of the function stretches a bit further than the second point (how much it stretches, it will depend on "unknown" parameters, but as a rule of thumb I use 35 "X" units).

I currently defined an ifelse function and I am applying a constrained optimization (optim, method="L-BFGS-B"), to find the best parameter values. This has a few handicaps, as the parameters are not normalized/scaled (which makes the search routine less efficient)

logistic.fun<- function(K, C, ro, b, Z, a, lambda){
      ifelse(test = lambda <= Seq_Bands[2]+25, 
           yes = a + b*((Seq_Bands[2]+25)- lambda),
           no = ((ro*K)/((ro + (K-ro)*exp(-C*(lambda-Seq_Bands[2]+Z))))))}

I would be keen to use nls and fit a 3 parameters logistic regression and also use a fourth parameter to estimate the breakpoint of the linear part (Z parameter code above). Also, I would be keen to avoid jumps such as those seen on left plot below. This means that the functions would have to be differentiable (?) at this breakpoint.

I am unsure how to code it.
Cheers and thanks!

enter image description here


UPDATE: It was correctly pointed out that the number of variables is higher than the data points presented; which would render the problem under-determined. Consequently, the logistic equation can only be parameterized using 3 variables (Fig 2: R_0_ , K and r. My understanding is that, necessarily, the breakpoint (Z) cannot be estimated, and should be set prior to the optimization process.
In context, it is also important to state that the R_0_, and K are not "truly" unknowns, as they are expressed by the measurements (second and fourth red datapoints). enter image description here

$\endgroup$

1 Answer 1

3
$\begingroup$

Let $\phi$ be the logistic function

$$\phi(z) = \frac{1}{1 + \exp(-z)}.$$

Your model shifts and scales the argument of $\phi$ and scales its values for arguments exceeding the breakpoint $\zeta,$ thereby requiring three parameters for $x\ge \zeta,$ which we could parameterize as

$$f_{+}(x;\mu,\sigma,\gamma) = \gamma\, \phi\left(\frac{x-\mu}{\sigma}\right).$$

For arguments less than the breakpoint you want a linear function

$$f_{-}(x;\alpha,\beta) = \alpha + \beta x.$$

Assure continuity by matching the values at the breakpoint. Mathematically this means

$$f_{-}(\zeta;\alpha,\beta) = f_{+}(\zeta;\mu,\sigma,\gamma),$$

enabling us to express one of the six parameters in terms of the other five. The simplest choice is the solution

$$\alpha = \gamma\, \phi\left(\frac{\zeta-\mu}{\sigma}\right) - \beta\,\zeta.$$

The resulting model will almost never be differentiable at $\zeta,$ but that doesn't matter.


The illustration in the question shows only four data points, which won't be enough to fit five parameters. But with more data points, measured with a little bit of mean-zero, iid error, a nonlinear least squares algorithm can succeed, especially if provided good starting values (which is an art in itself) and if some care is taken to re-express the parameters that must be positive ($\gamma$ and $\sigma$). Here is a comparable dataset with ten datapoints, obviously measured with substantial error:

Figure

It illustrates what the model looks like, how well it might fit even with such a small dataset, and what a likely 95% confidence interval for the breakpoint $\zeta$ might be (it is spanned by the red band). To find this fit I used $(\zeta,\beta,\mu,\log(\sigma),\log(\gamma))$ for the parameterization, requiring no constraints at all: see the call to nls in the code example below.

You can find effective starting values by eyeballing the data plot, which will clearly indicate reasonable values of $\beta,$ $\zeta,$ and possibly $\gamma.$ You might have to experiment with the other parameters. The model is a little dicey because there may be very strong correlations among $\zeta,$ $\sigma,$ $\gamma,$ and $\mu:$ this is characteristic of the logistic function, especially when only part of that function is reflected in the data.


To give you a leg up on experimenting and developing a solution, here is the R code used to create examples like this, fit the data, and plot the results. For experimentation, comment out the call to set.seed.

#
# The model.
#
f <- function(z, beta=0, mu=0, sigma=1, gamma=1, zeta=0) {
  logistic <- function(z) 1 / (1 + exp(-z))
  alpha <- gamma * logistic((zeta - mu)/sigma) - beta * zeta
  ifelse(z <= zeta,
         alpha + beta * z,
         gamma * logistic((z - mu) / sigma))
}
#
# Create a true model.
#
parameters <- list(beta=-0.0004, mu=705, sigma=20, gamma=0.65, zeta=675)
#
# Simulate from the model.
#
X <- data.frame(x = seq(540, 770, by=25))
X$y0 <- do.call(f, c(list(z=X$x), parameters))
#
# Add iid error, as appropriate for `nls`.
#
set.seed(17)
X$y <- X$y0 + rnorm(nrow(X), 0, 0.05)
#
# Plot the data and true model.
#
with(X, plot(x, y, main="Data with True and Fitted Models", cex.main=1, pch=21, bg="Gray"))
mtext(expression(paste("Black: true; Red: fit.  Vertical lines show ", zeta, ".")),
      side=3, line=0.25, cex=0.9)
curve(do.call(f, c(list(z=x), parameters)), add=TRUE, lwd=2, lty=3)
abline(v = parameters$zeta, col="Gray", lty=3, lwd=2)
#
# Fit the data.
#
fit <- nls(y ~ f(x, beta=beta, mu=mu, sigma=exp(sigma), gamma=exp(gamma), zeta=zeta), 
           data = X, 
           start = list(beta=-0.0004, mu=705, sigma=log(20), zeta=675, gamma=log(0.65)),
           control=list(minFactor=1e-8), trace=TRUE)
summary(fit)
#
# Plot the fit.
#
red <- "#d01010a0"
x <- seq(min(X$x), max(X$x), by=1)
y.hat <- predict(fit, newdata=data.frame(x=x))
lines(x, y.hat, col=red, lwd=2)
#
# Display a confidence band for `zeta`.
#
zeta.hat <- coefficients(fit)["zeta"]
se <- sqrt(vcov(fit)["zeta", "zeta"])
invisible(lapply(seq(zeta.hat - 1.645*se, zeta.hat + 1.645*se, length.out=201), 
      function(z) abline(v = z, col="#d0101008")))
abline(v = zeta.hat, col=red, lwd=2)
$\endgroup$
6
  • $\begingroup$ Thank you very for your elaborate solution. It all makes a lot of sense; however, I skipped some background information: - I am constrained to four data points. Basically, I am trying to reconstruct the original signal (i.e. hyperspectral signature) from these four measurements (multispectral). - You are correct about the need to be continuous (yet not differentiable - it doesn't matter). - Just to confirm our understanding: $\mu{}$ = starting point of the logistic equation. $\gamma{}$ = asymptote; $\zeta{}$ = breakpoint; $\sigma{}$ = scaling factor? thanks! $\endgroup$
    – Gustavo
    Jul 22, 2020 at 13:45
  • $\begingroup$ As a final remark, I am constrained to fit the logistic regression with only 3 data points. Does this pose a limitation to the code/modelling you shared? $\endgroup$
    – Gustavo
    Jul 22, 2020 at 13:48
  • 1
    $\begingroup$ Yes, those interpretations are correct. But your problem even with four points looks insoluble, since you have five parameters to fit: in general there will be a one-parameter family of solutions. It's even worse with three points. You will have to specify at least two of the parameters. $\endgroup$
    – whuber
    Jul 22, 2020 at 13:48
  • $\begingroup$ I have pondered about the under-determination (and consequently made some editing on my original question). Although $R_\text{0}$ and K (even $\beta{}$) are parameters, they are not unknowns: $\beta{}$ is the slope from the first to second point; and K and R0 are inputs (second and fourth red points). Basically, the only true unknown is the logistic growth rate. Not sure if this changes much, though. Finally, thanks a lot for the insights. $\endgroup$
    – Gustavo
    Jul 22, 2020 at 14:12
  • 1
    $\begingroup$ I see where you're going, but I don't fully believe your description. If $\sigma$ (my parameter related to the logistic growth rate) is the only unknown, then you cannot generally achieve a continuous fit. At the very least, the breakpoint $\zeta$ (or its equivalent) will be unknown, too. Regardless, my solution can be used as-is: in the call to nls, simply replace any parameters with their stipulated values in the arguments of f and omit those parameters from the list of starting values. $\endgroup$
    – whuber
    Jul 22, 2020 at 14:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.