Breaking the regression line into two pieces My X & Y variables are associated like this below and I am trying to fit a simple linear regression model (y ~ x , data= df) , to estimate β1 .

However I am not confident about extending the linear slope to values beyond 30. I want to restrict the regression slope to only x < 30
For x > 30 I want to slope to be horizontal so that β1 = 0. In other words breaking my linear regression line into two pieces x < 30 and x>30.

How can i do that ?
 A: You can easily extend the solution by @Noah to a non-fixed break point, which you want to estimate from the data. There is the problem, however, that the fitted function $f$ is not differentiable in the break point $x_0$. It is
$$f(x,b_0,b_1,x_0)=\left\{\begin{array}{ll}
b_0 + b_1 x & \mbox{ for } x\leq x_0 \\
b_0+b_1 x_0 & \mbox{ for } x> x_0
\end{array}\right.$$
or, implemented in R
f <- function(x, b0, b1, x0) {
  ifelse(x < x0, b0 + b1*x, b0 + b1*x0)
}

This means that the R function for general least squares optimization nls might report a "singular gradient" error, e.g. for the following example:
set.seed(100)
n <- 50
x <- runif(n, 0, 40)
y <- f(x, 6, 2, 30) + rnorm(n, 0, 30)

fit <- nls(y ~ f(x, b0, b1, x0), start=list(b0=0, b1=1, x0=10))

This can be easily circumvented by using a more robust minimization algorithm that also works for non-differentiable functions, e.g. Nelder-Mead, which is the default algorithm in the R function optim:
# helper function defining the LSQ objective function
ssq <- function(par, x, y) {
  sum((y - f(x, par[1], par[2], par[3]))^2)
}

res <- optim(par=c(0,1,10), ssq, x=x, y=y)
if (res$convergence != 0) {
  stop("optim not converged")
}

# visualize result
# fitted parameters are stored in res$par
plot(x,y)
xr <- range(x)
xx <- seq(xr[1], xr[2], by=0.1)
lines(xx, f(xx, res$par[1], res$par[2], res$par[3]), col="blue")


A: You can split the regression into two parts:
$$y_i = \beta_0^{<30} + \beta_1^{<30} x_i + e_i, \space x < 30 \\
y_i = \beta_0^{\ge30} + e_i, \space x_i \ge 30$$
We can write this as a single equation:
$$
y_i = I(x_i < 30)(\beta_0^{<30} + \beta_1^{<30} x_i) + I(x_i\ge30)(\beta_0^{\ge30}) + e_i
$$
where $I(.)$ is an indicator variable, equal to $1$ if the condition is met and $0$ otherwise, so it acts like a switch, turning on the left side or right side of the model.
We know that $\beta_0^{\ge30}$ needs to be constrained to be whatever value is at the end of the sloped line, i.e., the predicted value of $y$ when $x = 30$ based on the model for the left half. So we can rewrite $\beta_0^{\ge30}$ as $(\beta_0^{<30} + \beta_1^{<30}30)$ and then substitute, so we get
$$
y_i = I(x_i < 30)(\beta_0^{<30} + \beta_1^{<30} x_i) + I(x_i\ge30)(\beta_0^{<30} + \beta_1^{<30}30) + e_i
$$
This is a regression model with two coefficients to be estimated, so we can rearrange the terms to get a model in the form $y_i = \beta_0 + \beta_1 x^*_i + e_i$. First, we'll simplify notation by setting $d_i = I(x_i < 30)$ so that $1 - d_i = I(x_i \ge 30)$. Also, also I'll write $\beta_k^{<30}$ as $\beta_k$:
\begin{align}
y_i &= d_i(\beta_0 + \beta_1 x_i) + (1-d_i)(\beta_0 + \beta_1 30) + e_i \\
    &= \beta_0 d_i + \beta_1 x_i d_i + \beta_0(1-d_i) + \beta_1 30(1-d_i) + e_i \\
&= \beta_0 (d_i + (1-d_i))+ \beta_1 (x_i d_i + 30(1-d_i)) + e_i \\
&= \beta_0 + \beta_1 (x_i d_i + 30(1-d_i)) + e_i \\
&= \beta_0 + \beta_1 x^*_i + e_i
\end{align}
where $x^*_i = x_i d_i + 30(1-d_i)$.
So now we can simply create a new predictor, $x^*$, and fit the regression model. Here's an example of how to do that in R:
set.seed(100)
n <- 50
x <- sort(runif(n, 0, 40))
y <- 6 + 2 * x + rnorm(n, 0, 30)

d <- 1 * (x < 30)
x_star <- x * d + 30*(1-d)

fit <- lm(y ~ x_star)
y_hat <- fitted(fit)

plot(x, y)
lines(x, y_hat, col = "red")



# This lets you generate predictions for any value of x
fit <- lm(y ~ I(x * (x < 30) + 30 * (x >= 30)))

y_hat <- predict(fit, newdata = data.frame(x = 0:60))
plot(x, y, xlim = c(0, 60))
lines(0:60, y_hat, col = "red")


Created on 2023-01-25 with reprex v2.0.2
