Alright, let's work through this three-segment piecewise regression problem in R. You need a model where the first and last segments are flat (slope = 0), with a linear segment connecting them, and we’ll estimate both breakpoints directly from the data. Here’s a complete, step-by-step solution using your sample data, plus tips for scaling this to your 60 datasets.

First, let’s formalize the model to ensure continuity between segments. Let ( t_1 ) and ( t_2 ) be our two breakpoints (( t_1 < t_2 )). The model looks like this:

• For ( x \leq t_1 ): ( y = \beta_0 ) (flat baseline)
• For ( t_1 < x < t_2 ): ( y = \beta_0 + \beta_1(x - t_1) ) (linear transition)
• For ( x \geq t_2 ): ( y = \beta_0 + \beta_1(t_2 - t_1) ) (final flat level)

This structure guarantees the curve is continuous at both breakpoints, which is usually desirable for this type of regression.

Since we’re estimating breakpoints (non-linear parameters), we’ll use R’s nls() function (nonlinear least squares). Let’s start with your sample data:

# Load sample data
x <- c(1.306, 1.566, 1.736, 1.854, 2.082, 2.328, 2.650, 2.886, 3.162, 3.392)
y <- c(176.4, 188.0, 193.8, 179.4, 134.4, 119.0, 66.2, 58.2, 58.2, 41.2)
df <- data.frame(x, y)

Next, define the piecewise model function:

# Custom piecewise model function
piecewise_flat_linear <- function(x, beta0, beta1, t1, t2) {
  ifelse(x <= t1, 
         beta0,  # First flat segment
         ifelse(x < t2, 
                beta0 + beta1*(x - t1),  # Linear transition
                beta0 + beta1*(t2 - t1)  # Final flat segment
                )
         )
}

Critical: Set Initial Parameter Guesses #

nls() requires initial guesses for all parameters. We’ll base these on the data’s trend:

• beta0: Average of the first few y-values (baseline)
• beta1: Rough slope between the early and late y-values
• t1/t2: Guesses for where the slope changes (we’ll use visual intuition from the sample data)

# Initial parameter guesses tailored to the sample data
init_params <- list(
  beta0 = mean(y[x <= 1.8]),  # Avg y where x is in the first segment
  beta1 = (mean(y[x >= 2.8]) - mean(y[x <= 1.8])) / (2.8 - 1.8),  # Rough slope
  t1 = 1.8,  # Guess for first breakpoint
  t2 = 2.8   # Guess for second breakpoint
)

# Fit the model
fit <- nls(
  formula = y ~ piecewise_flat_linear(x, beta0, beta1, t1, t2),
  data = df,
  start = init_params,
  control = nls.control(maxiter = 1000)  # Increase iterations if needed
)

# View results
summary(fit)

The summary will give you estimated values for ( \beta_0 ), ( \beta_1 ), ( t_1 ), and ( t_2 ), plus standard errors and significance tests.

To verify the model fits well, plot the data and the regression line:

# Generate smooth prediction data
x_pred <- seq(min(x), max(x), length.out = 100)
y_pred <- predict(fit, newdata = data.frame(x = x_pred))

# Create plot
plot(df$x, df$y, pch = 16, col = "darkgray",
     main = "Three-Segment Piecewise Regression",
     xlab = "x", ylab = "y")
lines(x_pred, y_pred, col = "red", lwd = 2)
# Add dashed lines for breakpoints
abline(v = coef(fit)[c("t1", "t2")], lty = 2, col = "blue")

If you have 60 groups of data (e.g., stored in a long-format data frame with a group_id column), use a combination of dplyr and purrr to fit models in bulk. Here’s a reusable workflow:

library(dplyr)
library(purrr)
library(tibble)

# Example: Assume df_all has columns group_id, x, y
# Function to fit model for a single group
fit_single_group <- function(group_data) {
  x <- group_data$x
  y <- group_data$y
  
  # Auto-generate initial guesses using data quantiles (adapts to each group)
  q20 <- quantile(x, 0.2)
  q80 <- quantile(x, 0.8)
  init_beta0 <- mean(y[x <= q20])
  init_beta1 <- (mean(y[x >= q80]) - init_beta0) / (q80 - q20)
  
  # Try to fit the model, handle convergence failures gracefully
  tryCatch({
    group_fit <- nls(
      y ~ piecewise_flat_linear(x, beta0, beta1, t1, t2),
      data = group_data,
      start = list(beta0 = init_beta0, beta1 = init_beta1, t1 = q20, t2 = q80),
      control = nls.control(maxiter = 1000)
    )
    
    # Extract key results
    tibble(
      group_id = unique(group_data$group_id),
      beta0 = coef(group_fit)["beta0"],
      beta1 = coef(group_fit)["beta1"],
      t1 = coef(group_fit)["t1"],
      t2 = coef(group_fit)["t2"],
      final_level = coef(group_fit)["beta0"] + coef(group_fit)["beta1"]*(coef(group_fit)["t2"] - coef(group_fit)["t1"]),
      r_squared = cor(y, predict(group_fit))^2  # R-squared for fit quality
    )
  }, error = function(e) {
    # Return NA for groups where model fails to converge
    tibble(
      group_id = unique(group_data$group_id),
      beta0 = NA, beta1 = NA, t1 = NA, t2 = NA, final_level = NA, r_squared = NA
    )
  })
}

# Run bulk fitting
all_results <- df_all %>%
  group_split(group_id) %>%
  map_dfr(fit_single_group)

# View results for all groups
print(all_results)

• Convergence issues: If nls() throws an error, adjust your initial parameter guesses. Plot the group’s data first to get a better sense of where the breakpoints lie.
• Constraining breakpoints: If you need ( t_1 < t_2 ) enforced, ensure your initial guesses follow this rule (the model will maintain it during fitting).
• Alternative packages: If you prefer a more streamlined approach, check out the segmented package, but it requires a bit more customization to enforce the flat first/last segments.

内容的提问来源于stack exchange，提问作者Jana