Piecewise regression with a straight line and a horizontal line connecting at the break point

Question

I want to do a piecewise linear regression with one break point, where the second half of the regression line has slope = 0 . There are examples of how to do piecewise linear regression, for example here . The problem I am facing is that I do not know how to correct the tilt of half of the model to 0.

I tried

 lhs <- function(x) ifelse(x < k, kx, 0) rhs <- function(x) ifelse(x < k, 0, xk) fit <- lm(y ~ lhs(x) + rhs(x)) 

where k is the break point, but the segment on the right is not flat / horizontal.

I want to limit the slope of the second segment to 0. I tried:

 fit <- lm(y ~ x * (x < k) + x * (x > k)) 

but then again, I'm not sure how to get the other half to have a zero slope.

Any help is appreciated.

---

My own decision

I have a solution thanks to the comment below. Here is the code that I use to optimize and then plot the graph:

 x <- c(1, 2, 3, 1, 2, 1, 6, 1, 2, 3, 2, 1, 4, 3, 1) y <- c(0.041754212, 0.083491254, 0.193129615, 0.104249201, 0.17280516, 0.154342335, 0.303370501, 0.025503008, 0.123934121, 0.191486527, 0.183958737, 0.156707866, 0.31019215, 0.281890206, 0.25414608) range_x <- max(x) - min(x) intervals <- 1000 coef1 <- c() coef2 <- c() r2 <- c() for (i in 1:intervals) { k <- min(x) + (i-1) * (range_x / intervals) x2 = (x - k) * (x < k) fit <- lm(y ~ x2) coef1[i] <- summary(fit)$coef[1] coef2[i] <- summary(fit)$coef[2] r2[i] <- summary(fit)$r.squared } best_r2 <- max(r2) # get best r squared pos <- which.max(r2) best_k <- min(x) + (pos - 1) * (range_x / intervals) plot(x, y) curve(coef1[pos] - best_k * coef2[pos] + coef2[pos] * x, from=min(x), to=best_k, add = TRUE) segments(best_k, coef1[pos], max(x), coef1[pos]) 

Answer 1

There is a very similar flow in the stack overflow: Piecewise regression with a quadratic polynomial and a straight line smoothly connecting at the break point . The only difference is that we are now considering:

It turns out that the functions est , choose.c and pred defined in my answer do not need to be changed at all; we only need to modify getX to return the design matrix for your piecewise regression:

 getX <- function (x, c) cbind("beta0" = 1, "beta1" = pmin(x - c, 0)) 

Now we follow the example of a toy example to fit the model to your data:

 x <- c(1, 2, 3, 1, 2, 1, 6, 1, 2, 3, 2, 1, 4, 3, 1) y <- c(0.041754212, 0.083491254, 0.193129615, 0.104249201, 0.17280516, 0.154342335, 0.303370501, 0.025503008, 0.123934121, 0.191486527, 0.183958737, 0.156707866, 0.31019215, 0.281890206, 0.25414608) 

x varies from 1 to 6, so consider

 c.grid <- seq(1.1, 5.9, 0.05) fit <- choose.c(x, y, c.grid) fit$c # 4.5 

Finally, we make a prediction:

 x.new <- seq(1, 6, by = 0.1) p <- pred(fit, x.new) plot(x, y, ylim = c(0, 0.4)) matlines(x.new, p[,-2], col = c(1,2,2), lty = c(1,2,2), lwd = 2) 

We have rich information in the established model:

 str(fit) #List of 12 # $ coefficients : num [1:2] 0.304 0.055 # $ residuals : num [1:15] -0.06981 -0.08307 -0.02844 -0.00731 0.00624 ... # $ fitted.values: num [1:15] 0.112 0.167 0.222 0.112 0.167 ... # $ R : num [1:2, 1:2] -3.873 0.258 9.295 -4.37 # $ sig2 : num 0.00401 # $ coef.table : num [1:2, 1:4] 0.3041 0.055 0.0384 0.0145 7.917 ... # ..- attr(*, "dimnames")=List of 2 # .. ..$ : chr [1:2] "beta0" "beta1" # .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)" # $ aic : num -34.2 # $ bic : num -39.5 # $ c : num 4.5 # $ RSS : num 0.0521 # $ r.squared : num 0.526 # $ adj.r.squared: num 0.49 

For example, we can check the table of coefficient summaries:

 fit$coef.table # Estimate Std. Error t value Pr(>|t|) #beta0 0.30406634 0.03840657 7.917039 2.506043e-06 #beta1 0.05500095 0.01448188 3.797915 2.216095e-03 

Answer 2

Try making variables out of expression.

 x2 = (xk)*(x>k) lm( y ~ x2) 

Alternatively you can use I()

 lm(y~ I((xk)*(x>k))) 

I() takes everything that is literally inside and ignores other possible (incorrect) interpretations with any function inside.

If you do not have a well-defined k , you will have to optimize something like a deviation from different values ​​of k .