R : Robust nonlinear least squares fitting of three-phase linear model with confidence & prediction intervals

Question

I would like to fit a monotonically increasing three-phase linear model using nls in R. Say I have data

y <- c(4.5,4.3,2.57,4.40,4.52,1.39,4.15,3.55,2.49,4.27,4.42,4.10,2.21,2.90,1.42,1.50,1.45,1.7,4.6,3.8,1.9)
x <- 1500-c(320,419,650,340,400,800,300,570,720,480,425,460,675,600,850,920,975,1022,450,520,780)

I would like to get something like i.e. with x breakpoints at x=B1 and B2, plus 95% confidence & prediction intervals, which I would like to calculate based on the nls fit using the predFit function in the investr package.

The model here would be

y=(x < B1)*a +
  (x >= B1)*(x <= B2)*(a + b*(x - B1)) +
  (x > B2)*(a + b*(B2 - B1))

To take into account the desired constraints that B1 should be >min(x), I then set

B1 = min(x)+exp(logB1minminx)

To make sure that B2 > B1 I set

B2 = B1+exp(logB2minB1)

And to make sure that the slope of the middle section line b > 0 I set

b = exp(logb)

[I didn't quite know how to put in the remaining constraint that B2<max(x)]

To get an idea of sensible starting values for the slope parameter b I first calculated

f <- function (d) {
  m <- lm(y~x, as.data.frame(d))
  return(coef(m)[2])
}
require(zoo)
slopes <- rollapply(data.frame(x=x,y=y), 3, f, by.column=F)

Optimizing the parameters then kind of works using optimx with method="nlminb" (=port algo):

    preds = function (par) { 
  B1 = min(x)+exp(par[["logB1minminx"]]) # to make sure that B1 > min(x)
  B2 = B1+exp(par[["logB2minB1"]]) # to make sure that B2 > B1
  b = exp(par[["logb"]]) # to make sure that slope b > 0
  a = par[["a"]]
  pred = (x < B1)*a +
    (x >= B1)*(x <= B2)*(a + b*(x - B1)) +
    (x > B2)*(a + b*(B2 - B1))
  return(pred)
  }
SSR <- function (par, x=x, y=y) { # sums of squares
  fitted = preds(par)
  SS = sum((y - fitted)^2)
  return(SS) }
library(optimx)
fits = optimx(par = c(logB1minminx=log(650-min(x)), logB2minB1=log(1000-650), a=1.5, logb=log(mean(slopes))),
       lower = c(logB1minminx=-100, logB2minB1=-100, a=min(y), logb=-100),
       upper = c(logB1minminx=log(mean(x)-min(x)), logB2minB1=log(max(x)-min(x)), a=max(y), logb=log(max(slopes))),
       fn = SSR, 
       x = x,
       y = y,
       method = "nlminb",
       hessian=TRUE,
       control=list(all.methods=TRUE, maxit=1000, starttests=FALSE))
fits
#          logB1minminx logB2minB1        a      logb         value fevals gevals niter convcode kkt1 kkt2 xtimes
# L-BFGS-B     5.402100   5.859305 1.511979 -4.804957  6.405210e-01     41     41    NA        0   NA   NA   0.01
# nlminb       5.402677   5.858434 1.512409 -4.804421  6.404725e-01     65    155    31        1   NA   NA   0.00
# spg          5.402677   5.858560 1.512154 -4.804395  6.404726e-01    349     NA   195        0   NA   NA   0.11
# Rcgmin             NA         NA       NA        NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
# Rvmmin             NA         NA       NA        NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
# bobyqa       5.402677   5.859331 1.511529 -4.804637  6.404949e-01    148     NA    NA        0   NA   NA   0.00
# nmkb               NA         NA       NA        NA 8.988466e+307     NA     NA    NA     9999   NA   NA   0.00
# hjkb         5.147494   5.857933 1.500000 -5.218677  9.533185e+00      1     NA     0     9999   NA   NA   0.00
xvals=seq(min(x),max(x),length.out=1000)
plot(x, y, col="black",pch=16)
lines(xvals, 
      preds(coef(fits)["nlminb",], xvals), col="blue")

So this gives the fit as shown above. Most of the algorithms don't seem to converge though. And calculation of the Hessian fails, which is a problem since I need that to be able to calculate standard errors on the coefficients & confidence intervals & prediction intervals on the overall fit.

Likewise, nls and nlsLM both return errors when I ask the coefficient summaries, and this prevents me from calculating the confidence & prediction intervals using the investr package :

nlsfit = nls(y ~ (x < (min(x)+logB1minminx))*
               a + 
               (x >= (min(x)+logB1minminx))*
               (x <= (min(x)+exp(logB1minminx)+exp(logB2minB1)))*(a + exp(logb)*(x - (min(x)+logB1minminx))) + 
               (x > (min(x)+exp(logB1minminx)+exp(logB2minB1)))*
               (a + exp(logb)*((min(x)+exp(logB1minminx)+exp(logB2minB1)) - (min(x)+logB1minminx))),
             data = data.frame(x=x, y=y),
             algorithm = "port",
             start = c(logB1minminx=log(650-min(x)), logB2minB1=log(1000-650), a=1.6, logb=log(mean(slopes))),
             control = nls.control(maxiter=1000, warnOnly=TRUE) )
summary(nlsfit)
# Error in chol2inv(object$m$Rmat()) : 
#  element (4, 4) is zero, so the inverse cannot be computed
library(investr)
predFit(nlsfit, newdata=data.frame(x=xvals), interval="prediction")
# Error in solve.default(crossprod(R1)) : 
#   Lapack routine dgesv: system is exactly singular: U[4,4] = 0

require(minpack.lm)
nlslmfit = nlsLM(y ~ (x < (min(x)+logB1minminx))*
                   a + 
                   (x >= (min(x)+logB1minminx))*
                   (x <= (min(x)+exp(logB1minminx)+exp(logB2minB1)))*(a + exp(logb)*(x - (min(x)+logB1minminx))) + 
                   (x > (min(x)+exp(logB1minminx)+exp(logB2minB1)))*
                   (a + exp(logb)*((min(x)+exp(logB1minminx)+exp(logB2minB1)) - (min(x)+logB1minminx))),
                 data = data.frame(x=x, y=y),
      start = c(logB1minminx=log(650-min(x)), logB2minB1=log(1000-650), a=1.6, logb=log(mean(slopes))) )
# Error in nlsModel(formula, mf, start, wts) : 
#   singular gradient matrix at initial parameter estimates

Does anybody know how I could fit this kind of model robustly using either nls or nlsLM, perhaps through the use of a smooth continuous differentiable function that approaches the 3-phase linear model above, thereby allowing the first derivative to be passed to the optimizer? I tried myself with the 4 parameter logistic model, but couldn't quite find a good smooth centrosymmetric function that was close enough to the 3-phase linear model. If there is no clear breakpoint in the data I would like B1 to be estimated at min(x) and B2 at max(x), if there is no lower breakpoint I would like B1 to be estimated at min(x) and if there is no upper breakpoint I would like B2 to be estimated at max(x). In other words, the fit should ideally also work for data where the points just follow a linear model. Any thoughts?

EDIT: made a little bit of progress - I found a good smooth approximation and that fits OK with nlsLM. If I try it on the data without the upper breakpoint it still doesn't work though - I guess I have to try to fit several models - with 2 breakpoints, just a single breakpoint at either the lower or upper end or no breakpoints and see which one has the best AIC or BIC...

require(minpack.lm)
nlslmfit = nlsLM(y ~ a + (1/2)*exp(logb)*(B2-B1) + # we fit exp(logb) to force b > 0
                   (1/2)*sqrt(abs(exp(logb)*(4*1E-10+exp(logb)*(B1-x)^2))) - # now set s to 1E-10, we could also fit exp(logs) 
                   (1/2)*sqrt(abs(exp(logb)*(4*1E-10+exp(logb)*(B2-x)^2))),
                 data = data.frame(x=x, y=y),
                 start = c(B1=min(x)+1E-10, B2=max(x)-1E-10, a=min(y)+1E-10, logb=log(max(slopes))),
                 # lower = c(B1=min(x), B2=mean(x), a=min(y), logb=log(min(slopes[slopes>0]))),
                 # upper = c(B1=mean(x), B2=max(x), a=mean(y), logb=log(max(slopes))),
                 control = nls.control(maxiter=1000, warnOnly=TRUE) )
# as s->0 this smooth model approximates more closely the piecewise linear one
summary(nlslmfit)
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
#   B1    699.99988   19.23569   36.39  < 2e-16 ***
#   B2   1050.00069   15.49283   67.77  < 2e-16 ***
#   a       1.50817    0.09636   15.65 1.57e-11 ***
#   logb   -4.80172    0.06347  -75.65  < 2e-16 ***
require(investr)
xvals=seq(min(x),max(x),length.out=100)
predintervals = data.frame(x=xvals,predFit(nlslmfit, newdata=data.frame(x=xvals), interval="prediction"))
confintervals = data.frame(x=xvals,predFit(nlslmfit, newdata=data.frame(x=xvals), interval="confidence"))
require(ggplot2)
qplot(data=predintervals, x=x, y=fit, ymin=lwr, ymax=upr, geom="ribbon", fill=I("red"), alpha=I(0.2)) +
  geom_ribbon(data=confintervals, aes(x=x, ymin=lwr, ymax=upr), fill=I("blue"), alpha=I(0.2)) +
  geom_line(data=confintervals, aes(x=x, y=fit), colour=I("blue"), lwd=2) +
  geom_point(data=data.frame(x=x,y=y), aes(x=x, y=y, ymin=NULL, ymax=NULL), size=5, col="blue") +
  ylab("y")

# on subset of data without lower breakpoint:
nlslmfit = nlsLM(y ~ a + (1/2)*exp(logb)*(B2-B1) + # we fit exp(logb) to force b > 0
                   (1/2)*sqrt(abs(exp(logb)*(4*1E-10+exp(logb)*(B1-x)^2))) - # now set s to 1E-10, we could also fit exp(logs) 
                   (1/2)*sqrt(abs(exp(logb)*(4*1E-10+exp(logb)*(B2-x)^2))),
                 data = data.frame(x=x, y=y),
                 subset = x>760,
                 start = c(B1=min(x[x>760])+1E-10, B2=max(x)-1E-10, a=min(y)+1E-10, logb=log(max(slopes))),
                 # lower = c(B1=min(x), B2=mean(x), a=min(y), logb=log(min(slopes[slopes>0]))),
                 # upper = c(B1=mean(x), B2=max(x), a=mean(y), logb=log(max(slopes))),
                 control = nls.control(maxiter=1000, warnOnly=TRUE) )
summary(nlslmfit)
require(investr)
xvals=seq(760,max(x),length.out=100)
predintervals = data.frame(x=xvals,predFit(nlslmfit, newdata=data.frame(x=xvals), interval="prediction"))
confintervals = data.frame(x=xvals,predFit(nlslmfit, newdata=data.frame(x=xvals), interval="confidence"))
require(ggplot2)
qplot(data=predintervals, x=x, y=fit, ymin=lwr, ymax=upr, geom="ribbon", fill=I("red"), alpha=I(0.2)) +
  geom_ribbon(data=confintervals, aes(x=x, ymin=lwr, ymax=upr), fill=I("blue"), alpha=I(0.2)) +
  geom_line(data=confintervals, aes(x=x, y=fit), colour=I("blue"), lwd=2) +
  geom_point(data=data.frame(x=x,y=y)[x>760,], aes(x=x, y=y, ymin=NULL, ymax=NULL), size=5, col="blue") +
  ylab("y")

# on subset of data without upper breakpoint - here I still get an error:
nlslmfit = nlsLM(y ~ a + (1/2)*exp(logb)*(B2-B1) + # we fit exp(logb) to force b > 0
                   (1/2)*sqrt(abs(exp(logb)*(4*1E-10+exp(logb)*(B1-x)^2))) - # now set s to 1E-10, we could also fit exp(logs) 
                   (1/2)*sqrt(abs(exp(logb)*(4*1E-10+exp(logb)*(B2-x)^2))),
                 data = data.frame(x=x, y=y),
                 subset = x<1040,
                 start = c(B1=min(x)+1E-10, B2=max(x[x<1040])-1E-10, a=min(y)+1E-10, logb=log(max(slopes))),
                 # lower = c(B1=min(x), B2=mean(x), a=min(y), logb=log(min(slopes[slopes>0]))),
                 # upper = c(B1=mean(x), B2=max(x), a=mean(y), logb=log(max(slopes))),
                 control = nls.control(maxiter=1000, warnOnly=TRUE) )
summary(nlslmfit)
require(investr)
xvals=seq(min(x),1040,length.out=100)
# here prediction & confidence intervals still fail though:
predintervals = data.frame(x=xvals,predFit(nlslmfit, newdata=data.frame(x=xvals), interval="prediction"))
# Error in solve.default(crossprod(R1)) : 
# system is computationally singular: reciprocal condition number = 2.65525e-23
confintervals = data.frame(x=xvals,predFit(nlslmfit, newdata=data.frame(x=xvals), interval="confidence"))
require(ggplot2)
qplot(data=predintervals, x=x, y=fit, ymin=lwr, ymax=upr, geom="ribbon", fill=I("red"), alpha=I(0.2)) +
  geom_ribbon(data=confintervals, aes(x=x, ymin=lwr, ymax=upr), fill=I("blue"), alpha=I(0.2)) +
  geom_line(data=confintervals, aes(x=x, y=fit), colour=I("blue"), lwd=2) +
  geom_point(data=data.frame(x=x,y=y)[x<1040,], aes(x=x, y=y, ymin=NULL, ymax=NULL), size=5, col="blue") +
  ylab("y")

I tried it before but didn't have much luck with it unfortunately... The smooth approximation of a segmented model above worked best for me so far, but still it doesn't seem to work in all cases... — Tom Wenseleers, Mar 18 '18 at 07:36
Plus I need 95% confidence & prediction intervals on the model predictions and I want the first and last segment to be flat, which the segmented package doesn't provide or allow for... — Tom Wenseleers, Mar 18 '18 at 08:47

G. Grothendieck · Answer 1 · 2018-03-18T21:52:36.100

1

library(minpack.lm)
fo <- y ~ pmax(a1, pmin(a2 + b * x, a3)) 
co <- coef(lm(y ~ x))
fm <- nlsLM(fo, start = list(a1 = min(y), a2 = co[[1]], b = co[[2]], a3 = max(y)))

o <- order(x)
plot(y ~ x, subset = o)
lines(fitted(fm) ~ x, subset = o, col = "red")

summary(fm)

library(investr)
predFit(fm, data.frame(x), se = TRUE)

edited Mar 18 '18 at 21:52

answered Mar 18 '18 at 13:49

G. Grothendieck

254,981
17
203
341

Thanks - I tried things similar to this too, but it would still frequently fail on real data right? E.g. with data that were not really trilinear etc. Also I get an error Error in pmin(a2 + b * x, a3) : object 'a2' not found when I evaluate your predFit line. Any thoughts? – Tom Wenseleers Mar 18 '18 at 15:07
E.g. for me calculation of standard errors on the parameter fails if you leave out the data with x > 1040 (i.e. if there is no upper plateau in the data and you try to fit this 3-phase linear model). – Tom Wenseleers Mar 18 '18 at 21:23
`co[1]` and `co[2]` should have been `co[[1]]` and `co[[2]]`. Have fixed now. In general, nonlinear optimization may fail if the data is far from the assumed model. You can't just blindly apply it. What you can do is come up with a set of models and then check if the first one fails and try the second if it does and so on for as many models as you need to cover all your cases. – G. Grothendieck Mar 18 '18 at 21:56
Ha OK that's a good idea - I'll do that, and then select the model that ran OK and that had the best AIC... Your model above is still discontinuous though - wouldn't one expect a smooth continuously differentiable model like the approximation I provided above to fit more easily, e.g. using nlmrt::nlxb, which takes advantage of the analytically calculated Jacobian? – Tom Wenseleers Mar 18 '18 at 22:23
y is a continuous function of x, just not differentiable at the transition points but from a practical matter it will likely work if the data comes from the assumed model and if it does not then you should be using some other model. – G. Grothendieck Mar 18 '18 at 22:34

R : Robust nonlinear least squares fitting of three-phase linear model with confidence & prediction intervals

1 Answers1