How to Easily Perform Piecewise Regression in R: A Step-by-Step Guide

Prerequisites: Installing and Loading the `segmented` Package

Before fitting any segmented model, it is essential to ensure that the necessary statistical infrastructure is in place. In the R statistical environment, the primary tool for this analysis is the segmented package. This package extends standard linear model capabilities by introducing the ability to estimate breakpoints, or “change points,” within the model structure.

If you have not used this package before, you must install it from the Comprehensive R Archive Network (CRAN). The installation process is straightforward, followed by loading the package into your current R session using the library() function. This ensures that the core function, segmented(), becomes accessible for use.

# Install the segmented package (if needed)
# install.packages("segmented")

# Load the segmented library
library(segmented)

Step 1: Defining the Sample Data

The initial step in any regression analysis is establishing a suitable dataset. For this demonstration, we will construct a simple data frame in R that clearly exhibits a change in slope. This artificial dataset consists of two vectors, x (the independent variable) and y (the dependent variable), designed so that the relationship between them accelerates noticeably roughly halfway through the sequence of observations.

The creation of the data frame is crucial for organizing the variables for subsequent modeling. We use the data.frame() function in R to combine the vectors into a structured format, enabling easy referencing of columns (df$x and df$y). Viewing the initial rows confirms the structure and provides a quick verification of the data input, ensuring accuracy before moving to visualization.

# Create the DataFrame
df <- data.frame(x=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16),
                 y=c(2, 4, 5, 6, 8, 10, 12, 13, 15, 19, 24, 28, 31, 34, 39, 44))

# View first six rows of the data frame to confirm structure
head(df)

  x  y
1 1  2
2 2  4
3 3  5
4 4  6
5 5  8
6 6 10

Step 2: Initial Data Visualization and Breakpoint Identification

Visualizing the data is perhaps the most informative initial step in piecewise regression analysis. By generating a scatterplot of the independent variable x against the dependent variable y, we can visually inspect the overall trend and attempt to identify the approximate location of the potential breakpoint. This initial visual assessment guides the parameter settings we use later in the segmented() function.

The scatterplot below clearly illustrates the changing relationship. For the initial values of x, the slope is relatively shallow, indicating a moderate increase in y for every unit increase in x. However, this trend abruptly steepens around the value x=9. This sharp shift in trajectory strongly suggests that a standard single linear regression model would provide a poor fit, thus validating the need for a segmented approach.

# Create scatterplot of x vs. y
plot(df$x, df$y, pch=16, col='steelblue')

From this visual inspection, it is evident that the relationship between x and y experiences an abrupt change around the value of x = 9. This serves as our initial estimate for the breakpoint parameter (psi) when we fit the segmented model.

Step 3: Fitting the Initial Linear Model

The segmented() function requires a pre-existing standard linear regression model object (an lm object) as its primary input. This initial model, typically fitted using the lm() function, acts as the starting point for the iterative estimation process used by the segmented package. Although this simple model will likely provide a poor fit (due to the observed breakpoint), it establishes the basic structure and parameters necessary for the segmentation algorithm to operate.

We fit the simple linear regression model using the syntax y ~ x, indicating that y is predicted by x based on the data in our data frame, df. We store the result of this operation in an object named fit. This step is purely preparatory but crucial for the functionality of the segmented analysis.

# Fit simple linear regression model (required precursor for segmented analysis)
fit <- lm(y ~ x, data=df)

Step 4: Applying Piecewise Regression using `segmented()`

With the initial linear model established, we can now apply the core function of the package: segmented(). This function takes the initial lm object (fit) and refines it by searching for the optimal location of the segmentation point. We must specify two main arguments here: seg.Z and psi.

The argument seg.Z = ~x specifies which independent variable contains the variable where the breakpoint is expected to occur (in this case, x). The argument psi = 9 provides the starting guess for the location of the breakpoint. While the segmented() function will iteratively search for the best fit, providing a reasonable starting value (like 9, which we visually identified in Step 2) helps ensure fast and reliable convergence. The resulting object, segmented.fit, contains the sophisticated segmented model structure.

# Fit piecewise regression model to original model, estimating a breakpoint starting at x=9
segmented.fit <- segmented(fit, seg.Z = ~x, psi=9)

# View summary of the segmented model results
summary(segmented.fit)

Call: 
segmented.lm(obj = fit, seg.Z = ~x, psi = 9)

Estimated Break-Point(s):
         Est. St.Err
psi1.x 8.762   0.26

Meaningful coefficients of the linear terms:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.32143    0.48343   0.665    0.519    
x            1.59524    0.09573  16.663 1.16e-09 ***
U1.x         2.40476    0.13539  17.762       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6204 on 12 degrees of freedom
Multiple R-Squared: 0.9983,  Adjusted R-squared: 0.9978 

Convergence attained in 2 iter. (rel. change 0)

Step 5: Interpreting the Model Summary and Coefficients

The output of the summary() function provides crucial details regarding the optimized piecewise regression model. The most critical component is the estimation of the breakpoint itself, listed under Estimated Break-Point(s). The algorithm successfully refined our initial guess of 9 and determined the optimal separation point (psi) to be x = 8.762, with a relatively small standard error, indicating high confidence in this location.

The coefficients define the two distinct regression lines. The coefficients (Intercept) and x define the regression line for the first segment (where X is less than the breakpoint). The coefficient U1.x represents the change in slope after the breakpoint. This is the difference between the slope of the second segment and the slope of the first segment.

Specifically:

• First Segment Slope: The coefficient for x is 1.59524. This is the slope when $X le 8.762$.
• Second Segment Slope: This is calculated as the initial slope (1.59524) plus the change in slope (U1.x, which is 2.40476). Thus, the second slope is $1.59524 + 2.40476 = 4.00000$. This highly significant increase shows the relationship nearly triples in steepness after the breakpoint.
• Model Fit: The high R-squared values (0.9983) and low Residual standard error (0.6204) confirm that the segmented model provides an extremely accurate fit to the data, significantly better than what a standard single linear regression model could achieve.

Step 6: Calculating Estimated Values from the Model

Based on the coefficients derived in the previous step, we can formally write the equation for the fitted piecewise regression model, using the estimated breakpoint at $X_{psi} = 8.762$. The general form of the two equations is:

For the first segment (when $x le 8.762$):

If $x le 8.762$: $y = 0.32143 + 1.59524 times (x)$

For the second segment (when $x > 8.762$): the equation must account for the intersection point to maintain continuity.

If $x > 8.762$: $y = 0.32143 + 1.59524 times (8.762) + (1.59524 + 2.40476) times (x – 8.762)$

We can use these derived equations to predict dependent variable values for specific inputs of x. Suppose we want to estimate y when x = 5, which falls into the first segment ($5 le 8.762$):

• $y = 0.32143 + 1.59524 times (5)$
• $y = 0.32143 + 7.9762$
• $y = 8.29763$

Now, consider a value of x = 12, which falls into the second segment ($12 > 8.762$). The estimated y value requires the second, more complex formula:

• $y = 0.32143 + 1.59524 times (8.762) + (1.59524 + 2.40476) times (12 – 8.762)$
• $y = 0.32143 + 13.9743 + (4.00000) times (3.238)$
• $y = 14.29573 + 12.952$
• $y = 27.24773$

Step 7: Visualizing the Final Piecewise Fit

The final and most satisfying step is to overlay the derived piecewise regression model onto the original scatterplot to visually confirm the goodness of fit. This is accomplished easily within R using the plot() function for the original data, followed by applying plot() to the segmented.fit object with the add=T argument.

The generated visualization confirms that the two distinct linear segments accurately capture the underlying trend of the data. The model effectively handles the sharp change in slope that occurs near $x = 8.762$, demonstrating the superiority of segmented regression over a standard single linear fit for datasets exhibiting such nonlinear characteristics.

# Plot original data points
plot(df$x, df$y, pch=16, col='steelblue')

# Add the segmented regression model line
plot(segmented.fit, add=T)

As seen in the figure, the line segments align extremely well with the observations across the entire range of the independent variable, confirming the successful application of the piecewise regression methodology.

For additional learning and information on advanced regression models in R, consult the following tutorials: