Linear Plateau in R

When working with data in fields such as agriculture, biology, and economics, it’s common to observe a response that increases with an input up to a certain threshold and then levels off. For instance, crop yield might initially increase as more fertilizer is added, but after a certain point (the “breakpoint”), additional fertilizer provides no further yield benefit. In such scenarios, linear plateau models become a powerful tool for identifying the threshold where further increases in the input fail to boost the response.

In this blog post, we’ll walk through:

1. What a linear plateau model is and why it’s useful.
2. How to fit a linear plateau model using R.
3. Best practices for diagnosing the fitted model and interpreting parameter estimates.

Step-by-Step Guide: Fitting a Linear Plateau Model in R

Below is a detailed R example. While we use simulated data here, you can adapt the same approach to your real datasets.

1. Load Libraries

We’ll need a few packages:

• ggplot2 for visualization,
• nls2 for nonlinear least squares (NLS) with more flexible starting values,
• nlstools for diagnostics and confidence intervals.

# Load necessary libraries
library(ggplot2)
library(nls2)       # Provides alternative fitting algorithms for nonlinear models
library(nlstools)   # Useful for diagnostics and confidence intervals in NLS        

2. Generate or Import Your Data

In practice, you will import your dataset using functions like read.csv(), read_excel(), or by connecting to a database. Here, we simulate some data for demonstration. Notice that we define a true breakpoint at x = 10:

# Simulated data for demonstration
set.seed(123)
x <- seq(1, 20, by = 0.5)
y_linear <- ifelse(x <= 10,
                   3 + 2 * x + rnorm(length(x), sd = 2),  # Linear portion
                   23 + rnorm(length(x), sd = 2))         # Plateau portion

data_linear <- data.frame(x = x, y = y_linear)        

3. Fit the Linear Plateau Model

The nls2 function allows us to specify a piecewise model via ifelse(x <= c, a + b*x, a + b*c). We provide initial guesses for parameters ‘start‘`start`‘start‘ based on our knowledge (or best estimates) about the dataset.

# Provide initial guesses for the parameters
start_vals <- list(a = 3, b = 2, c = 10)

# Fit the linear plateau model using nls2
linear_plateau_model <- nls2(
  formula = y ~ ifelse(x <= c, a + b * x, a + b * c),
  data    = data_linear,
  start   = start_vals,
  algorithm = "default"  # Uses Gauss-Newton approach
)

# Check a summary of the model fit
summary(linear_plateau_model)        

Note: If you encounter convergence issues (e.g., the model won’t fit or yields strange parameter values), you can try different fitting algorithms (like "brute-force") or refine your starting values. For instance:

# Example of brute-force search for better initial values:
# linear_plateau_model <- nls2(
#   formula = y ~ ifelse(x <= c, a + b*x, a + b*c),
#   data = data_linear,
#   start = expand.grid(
#     a = seq(2, 4, length = 10),
#     b = seq(1, 3, length = 10),
#     c = seq(5, 15, length = 11)
#   ),
#   algorithm = "brute-force"
# )        

4. Extract Coefficients and Build the Equation

Once the model converges, we can retrieve the parameter estimates ‘a‘, ‘b‘, and ‘c‘:

coef_estimates <- coef(linear_plateau_model)
a_est <- coef_estimates["a"]
b_est <- coef_estimates["b"]
c_est <- coef_estimates["c"]

# Construct a piecewise equation for clarity
piecewise_equation <- paste0(
  "y = ", round(a_est, 2), " + ", round(b_est, 2), " * x,   for x ≤ ", round(c_est, 2), "\n",
  "y = ", round(a_est + b_est * c_est, 2), "         ,   for x > ",  round(c_est, 2)
)

piecewise_equation        

5. Diagnostic Checks and Confidence Intervals

Model diagnostics help ensure we haven’t overlooked poor fit, outliers, or skewed residual distributions. We also want to gauge the uncertainty in parameters (especially the breakpoint ‘c‘):

# Residual diagnostics using nlstools
nls_diag <- nlsResiduals(linear_plateau_model)
plot(nls_diag)  # Residuals vs. fitted, Q-Q plot, etc.

# Approximate 95% confidence intervals for parameters
conf_int <- confint2(linear_plateau_model, level = 0.95)
conf_int        

• Residual plots: Look for random scatter without strong patterns.
• Confidence intervals: Inspect the range of plausible values for a, b, and c. Wide intervals for c might indicate limited data near the plateau transition.

6. Visualize the Fitted Model

A clear plot showcasing the raw data points and the fitted piecewise function helps communicate the model’s performance:

# Predict values from the fitted model
data_linear$y_pred <- predict(linear_plateau_model, newdata = data_linear)

ggplot(data_linear, aes(x = x, y = y)) +
  geom_point(color = "black") +
  geom_line(aes(y = y_pred), color = "blue", size = 1) +
  geom_vline(xintercept = c_est, linetype = "dashed", color = "red") +
  annotate(
    "text", x = c_est + 0.5, y = max(data_linear$y) * 0.9,
    label = paste("Breakpoint =", round(c_est, 2)),
    color = "red", angle = 90, vjust = -0.5
  ) +
  labs(
    title = "Fitted Linear Plateau Model",
    subtitle = "Response vs. Independent Variable with Identified Breakpoint",
    x = "Independent Variable (x)",
    y = "Response Variable (y)"
  ) +
  annotate(
    "text",
    x = mean(range(data_linear$x)),
    y = max(data_linear$y) * 1.05,
    label = piecewise_equation,
    hjust = 0.5,
    color = "blue"
  ) +
  theme_minimal()        

In this visualization:

• Black points represent observed data.
• Blue line is the fitted model’s predictions.
• Dashed red line indicates the estimated breakpoint ‘cest‘`c_est`‘cest‘.
• Annotated text shows the piecewise equation and the breakpoint value for clearer interpretation.

Interpretation and Practical Insights

• Intercept ‘a‘: The baseline or starting level of the response when x = 0.
• Slope ‘b‘: The rate at which the response increases per unit increase in xxx (for x≤cx).
• Breakpoint ‘c‘: The threshold of x beyond which further increases do not raise y.
• Plateau Value: `a + b * c`: This is the maximum value y attains.