# How to write a piecewise regression model as a linear model?

Let's write the following piecewise regression model $$y= \alpha_0 + \alpha_1 x +\epsilon ;\ \ x\le x_0$$

$$y=\beta_0 +\beta_1 x + \epsilon;\ \ x\gt x_0$$

When the value $$x_0$$ is known, this regression model is a linear model. How to write this model in the form $$Y= X\beta +\epsilon$$ where $$X$$ is a matrix and $$\beta$$ is a parameter vector?

• Are you thinking of something like so: r-bloggers.com/… Commented Nov 7, 2015 at 19:34

The overall model has four parameters: $$\alpha_0,$$ $$\alpha_1,$$ $$\beta_0,$$ and $$\beta_1.$$ Therefore, if a solution is at all possible, we must be able to construct four corresponding variables $$z_1, z_2,$$ $$z_3,$$ and $$z_4.$$

One solution dedicates the first two parameters to the first model and the second two parameters to the second model. Thus, when $$x\le x_0,$$ we may set $$\mathbf z = (z_1,z_2,z_3,z_4)=(1,x,0,0)$$ and when $$x \gt x_0$$ set $$\mathbf z = (z_1,z_2,z_3,z_4)=(0,0,1,x).$$ The zeros ensure the parameters for one part of the model do not affect the other part of the model.

To put it another way, for any parameter vector $$\gamma = (\gamma_1,\gamma_2,\gamma_3,\gamma_4)^\prime$$ notice that

\begin{aligned} \mathbf z \gamma &= (1,x,0,0)(\gamma_1,\gamma_2,\gamma_3,\gamma_4)^\prime = \gamma_1 + \gamma_2 x,& \ x \le x_0\\ \mathbf z \gamma &= (0,0,1,x)(\gamma_1,\gamma_2,\gamma_3,\gamma_4)^\prime = \gamma_3 + \gamma_4 x,& \ x \gt x_0. \end{aligned}

This is precisely the desired model with $$\gamma_1=\alpha_0, \gamma_2=\alpha_1,$$ $$\gamma_3=\beta_0,$$ and $$\gamma_4=\beta_1.$$ As is usual, when you assemble the observations as rows in a "model matrix" $$Z,$$ the model can be written $$y = Z\gamma + E$$ where $$E$$ is the vector of errors.

The $$x$$ values in the plot are $$0, 1, 2, \ldots, 20$$ and $$x_0=11.$$ The model matrix $$Z$$ is

      [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    1    1    0    0
[3,]    1    2    0    0
[4,]    1    3    0    0
[5,]    1    4    0    0
[6,]    1    5    0    0
[7,]    1    6    0    0
[8,]    1    7    0    0
[9,]    1    8    0    0
[10,]    1    9    0    0
[11,]    1   10    0    0
[12,]    1   11    0    0
[13,]    0    0    1   12
[14,]    0    0    1   13
[15,]    0    0    1   14
[16,]    0    0    1   15
[17,]    0    0    1   16
[18,]    0    0    1   17
[19,]    0    0    1   18
[20,]    0    0    1   19
[21,]    0    0    1   20


The pattern is clear.

These data and fits were produced by the following R implementation of this procedure.

x <- seq(0, 20)
gamma <- c(12, -1, 0, 1/2)
x0 <- 11
sigma <- 1
#
# Create the model matrix Z and true responses y.
#
Z <- cbind(x <= x0, x*(x <= x0), x > x0, x*(x > x0))
y. <- Z %*% gamma
#
# Create random responses according to the model.
#
set.seed(17)
y <- y. + rnorm(length(y.), 0, sigma)
#
# Find the least-squares fit.
#
fit <- lm.fit(Z, y)
#
# Plot the data, model, and fit.
#
plot(x, y, type="n", ylim=c(0,15), main="Data, True Model (Red), and Fit (Dashed Blue)")

f <- function(x, gamma) cbind(x <= x0, x*(x <= x0), x > x0, x*(x > x0)) %*% gamma
curve(f(x, gamma), add=TRUE, n=2001, col="Red", lwd=2)
curve(f(x,  fit\$coefficients), add=TRUE, n=2001, lty=3, col="Blue", lwd=2)
abline(v = x0, lwd=2, col="Gray")
text(x0, 14, expression(x[0]), pos=4)
text(1,1, expression(y==gamma[1]+gamma[2]*x+epsilon), pos=4)
text(19,1, expression(y==gamma[3]+gamma[4]*x+epsilon), pos=2)
points(x, y, pch=19)


For a known discontinuous break point $$x_0$$, the following piecewise regression model:

$$y= \alpha_0 + \alpha_1 x +\epsilon ;\ \ x\le x_0$$

$$y=\beta_0 +\beta_1 x + \epsilon;\ \ x\gt x_0$$

can instead be expressed using an indicator function of $$x$$ and $$x_0$$ for the change in intercept, and a hinge function of $$x_i$$ and $$x_0$$ for the change in slope:

$$y_{i} = \beta_{0} + \beta_{x}x_{i} + \beta_{0\text{c}}I(x_{i} \ge x_0) + \beta_{x\text{c}} \max(x_i - x_0,0) + \varepsilon_{i}$$

Explanation

Change in intercept: If $$x_i < x_0$$, then $$I(x_i \ge x_0) = 0$$, so the $$\beta_{0\text{c}}I(x_{i} \ge x_0)$$ term also equals $$0$$. However, if $$x_i \ge x_0$$, then $$I(x_i \ge x_0) = 1$$, and the $$\beta_{0\text{c}}I(x_{i} \ge x_0)$$ term equals $$\beta_{0\text{c}}$$… which is just a constant, so the intercept for values at $$x_0$$ and higher equals $$\beta_0 + \beta_{0\text{c}}$$.

Change in slope: If $$x_i < x_0$$, then $$\max(x_i - x_0,0) = 0$$, so the $$\beta_{x\text{c}} \max(x_i - x_0,0)$$ term also equals $$0$$. However, if $$x_i \ge x_0$$, then the $$\max(x_i - x_0,0)$$ term increases at exactly the same rate as $$x_i$$: a one-unit increase in $$x_i$$ corresponds to a one-unit increase in $$\max(x_i-x_0,0)$$. For example, if $$x_i = x_0 + 2$$, then $$\max(x_i - x_0,0)=2$$. The effect of $$x$$ on $$y$$ when $$x_i < x_0$$ is just $$\beta_x$$, but the effect of $$x$$ on $$y$$ changes from $$\beta_{x}$$ to $$\beta_x + \beta_{x\text{c}}$$ for values of $$x_i \ge x_0$$.

A single linear model

You can thus create a single linear model with two new variables, where $$x_{1i} = x_i$$ and, say, $$x_{2i} = I(x_{1i} \ge x_0)$$, and $$x_{3i} = \max(x_{1i}-x_0,0)$$ and estimate (e.g., using OLS):

$$\boldsymbol{y_i = \beta_0 + \beta_x x_{1i} + \beta_{0\text{c}}x_{2i} + \beta_{x\text{c}}x_{3i} + \varepsilon_{i}} = BX_{i} + \varepsilon_{i}$$