Piecewise regression with constraints

Question

I want to make a piecewise linear regression in R. I have a large dataset with 3 segments where I want the first and third segment to be without slope, i.e. parallel to x-axis and I also want the regression to be continuous. I have a small example dataset and example code below. I'm rather new to R and have tried to write the model formula without success. Can anyone help with the formula and how to extract intercept and slope for the different segments? Three specific questions are inserted as comments in the code. I have looked at package segmented but are unable to understand how to constrain segment 1 and 3 to be parallel to x-axis.

#Example 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 <- c(320,419,650,340,400,800,300,570,720,480,425,460,675,600,850,920,975,1022,450,520,780)  
plot(x, y, col="black",pch=16)

#model 1 not continuous, Q1: how to get that?
fit1<-lm(y ~ I(x<500)+I((x>=500&x<800)*x) + I(x>=800)) 
summary(fit1) #intercepts and slopes extracted from summary(fit1)
lines(c(min(x),500),c(round(summary(fit1)[[4]][1]+summary(fit1)[[4]][2],2),round(summary(fit1)[[4]][1]+summary(fit1)[[4]][2],2)),col="red")
lines(c(800,max(x)),c(round(summary(fit1)[[4]][1]+summary(fit1)[[4]][4],2),round(summary(fit1)[[4]][1]+summary(fit1)[[4]][4],2)),col="red")
lines(c(500,800),c((summary(fit1)[[4]][1])+(500*(summary(fit1)[[4]][3])),
               (summary(fit1)[[4]][1])+(800*(summary(fit1)[[4]][3]))),col="red")

#model 2 continuous but first and third segment not parallell to x-axis, Q2: how to get that?
fit2<-lm(y ~ x + I(pmax(x-500,0)) + I(pmax(x-800,0)))                                                  
summary(fit2) # Q3: how to get intercept and slope from summary(fit2) for model2?
mg <- seq(min(x),max(x),length=400)
mpv <- predict(fit2, newdata=data.frame(x = mg,t = (mg - 800)*as.numeric(mg > 800)))
lines(mg, mpv,col="blue")

Also added a nonlinear least squares solution, with calculation of 95% confidence & prediction intervals... — Tom Wenseleers, Commented Mar 18, 2018 at 8:40

Roland · Accepted Answer · 2015-05-04 15:51:35Z

If the goal is simply to fit a function, you could treat this as an optimization problem:

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 <- c(320,419,650,340,400,800,300,570,720,480,425,460,675,600,850,920,975,1022,450,520,780)  
plot(x, y, col="black",pch=16)


#we need four parameters: the two breakpoints and the starting and ending intercepts
fun <- function(par, x) {
  #set all y values to starting intercept
  y1 <- x^0 * par["i1"]
  #set values after second breakpoint to ending intercept
  y1[x >= par["x2"]] <- par["i2"]
  #which values are between breakpoints?
  r <- x > par["x1"] & x < par["x2"]
  #interpolate between breakpoints
  y1[r] <- par["i1"] + (par["i2"] - par["i1"]) / (par["x2"] - par["x1"]) * (x[r] - par["x1"])
  y1
}

#sum of squared residuals
SSR <- function(par) {
     sum((y - fun(par, x))^2)
   }


library(optimx)
optimx(par = c(x1 = 500, x2 = 820, i1 = 5, i2 = 1), 
       fn = SSR, 
       method = "Nelder-Mead")

#                  x1       x2       i1       i2     value fevals gevals niter convcode kkt1 kkt2 xtimes
#Nelder-Mead 449.8546 800.0002 4.381454 1.512305 0.6404728    373     NA    NA        0 TRUE TRUE   0.06

lines(300:1100, 
      fun(c(x1 = 449.8546, x2 = 800.0002, i1 = 4.381454, i2 = 1.512305), 300:1100))

resulting plot

Tom Wenseleers · Accepted Answer · 2018-03-19 08:34:13Z

If you also want confidence and prediction intervals you can first approximate your three-phase piecewise linear function by a smooth function, do an nls fit and then using the investr package (this helps for the fitting as the function is then continuously differentiable).

In your case:

x <- c(478, 525, 580,  650,  700,  720,  780,  825,  850,  900,  930,  980, 1020, 1040, 1050, 1075, 1081, 1100, 1160, 1180, 1200)
y <- c(1.70, 1.45, 1.50, 1.42, 1.39, 1.90, 2.49, 2.21, 2.57, 2.90, 3.55, 3.80, 4.27, 4.10, 4.60, 4.42, 4.30, 4.52, 4.40, 4.50, 4.15)
# calculate rolling slopes at each point to provide good initial estimates for slope parameter b
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)

# smooth approximation is
# y ~ a + (1/2)*b*(B2-B1) + 
#         (1/2)*sqrt(abs(b*(4*s+b*(B1-x)^2))) -  
#         (1/2)*sqrt(abs(b*(4*s+b*(B2-x)^2)))
# this smooth approximation approaches the piecewise linear model more as s -> 0

require(minpack.lm)
nlslmfit = nlsLM(y ~ a + (1/2)*exp(logb)*(B2-B1) + # we fit exp(logb) to force b > 0, if you don't want this just fit b instead
                   (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")

You can also do a robust nls fit (slightly more robust to outliers) using the nlrob function in the robustbase package, rest is the same as above:

require(robustbase)
nlsrobfit  <- nlrob(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),
                    maxit = 1000,
                    method="M",
                    algorithm="port",
                    doCov=TRUE,
                    start = c(B1=min(x)+1E-10, B2=max(x)-1E-10, a=min(y)+1E-10, logb=log(mean(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(nlsrobfit)
class(nlsrobfit)="nls" # for compatibility with investr

Comparison with model where s parameter is also fitted:

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*exp(logs)+exp(logb)*(B1-x)^2))) - # we now fit exp(logs) 
                     (1/2)*sqrt(abs(exp(logb)*(4*exp(logs)+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(mean(slopes)), logs=-10),
                 control = nls.control(maxiter=1000, warnOnly=TRUE) )
summary(nlslmfit)
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
#   B1    7.000e+02  2.079e+01   33.67 2.78e-16 ***
#   B2    1.051e+03  1.614e+01   65.08  < 2e-16 ***
#   a     1.514e+00  1.000e-01   15.13 6.70e-11 ***
#   logb -4.806e+00  7.131e-02  -67.39  < 2e-16 ***
#   logs -1.805e+01  4.561e+04    0.00        1    
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")

Comparison with smooth 4-parameter logistic model:

M.4pl <- function(x, lower.asymp, upper.asymp, inflec, hill){
  f <- lower.asymp + ((upper.asymp - lower.asymp)/
                        (1 + (x / inflec)^-hill))
  return(f)
}

require(minpack.lm)
nlslmfit = nlsLM(y ~ M.4pl(x, lower.asymp, upper.asymp, inflec, hill),
                 data = data.frame(x=x, y=y),
                 start = c(lower.asymp=min(y)+1E-10, upper.asymp=max(y)-1E-10, inflec=mean(x), hill=1),
                 control = nls.control(maxiter=1000, warnOnly=TRUE) )
summary(nlslmfit)
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
#   lower.asymp   1.5371     0.1080   14.24 7.06e-11 ***
#   upper.asymp   4.5508     0.1497   30.40 2.93e-16 ***
#   inflec      889.1543    14.0924   63.09  < 2e-16 ***
#   hill         13.1717     2.5475    5.17 7.68e-05 ***
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")

JJacquelin · Accepted Answer · 2018-07-07 15:43:53Z

In the general case, the linear piecewise regression for three segments leads to this kind of function :

The parameters where computed according to the direct method (not iterative) given page 30 in the paper https://fr.scribd.com/document/380941024/Regression-par-morceaux-Piecewise-Regression-pdf

But the first and third segment are not parallel to the x-axis, which probably is not what is expected.

The particular case of piecewise regression for thee segments with first and third segments parallel to the x-axis is given page 18 of the paper. The result is :

Stack Exchange Network

Piecewise regression with constraints

3 Answers 3

Your Answer

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

Linked

Hot Network Questions

Piecewise regression with constraints

3 Answers 3

Your Answer

Sign up or log in

Post as a guest

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

Linked

Related

Hot Network Questions

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