Solved – Piecewise regression with constraints

piecewise linearr

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")

Best Answer

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

Related Question