What's wrong with the piecewise fitting

I am new to R. I want to ask a question below:

Here is the data:

year <- c(2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022)
score <- c(85,  85,  88,  88,  94,  94,  94,  82,  82,  84,  84,  84,  84,  84, 84)

I want to set 2015 as the breakpoint to do the piecewise linear fitting (2008-2015 and 2015-2022). I have tried the following code, it get the results below. However, I think the result is not correct, especially for stage2, which shoube an increasing trend.

stage1 <- year - 2008
stage2 <- (year - 2015) * (year >= 2015)
fm <- lm(score~ stage1 + stage2)
summary(fm)

library(car)
linearHypothesis(fm, "stage1 + stage2", verbose = TRUE)

plot(score ~ year)
lines(fitted(fm) ~ year, col = "red")
abline(v = 2015, lty = 2)

The piecewise linaer fitting result

The main problem is using an inappropriate model. The model described the data in the SO post it was taken from (https://stackoverflow.com/questions/76480532/how-can-i-set-the-breakpoints-myself-to-do-the-piecewise-linear-fitting-with-man/76482328#76482328) but not this data. In this case a model consisting of discontinuous rather than continuous line segments seems more appropriate.

stage1.slope <- (year < 2015) * (year - 2015)
stage1.icept <- +(year < 2015)
stage2.slope <- (year >= 2015) * (year - 2015)
stage2.icept <- +(year >= 2015)
fm <- lm(score ~ stage1.icept + stage1.slope + stage2.icept + stage2.slope + 0)
summary(fm)

## Call:
## lm(formula = score ~ stage1.icept + stage1.slope + stage2.icept + 
##     stage2.slope + 0)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.71429 -0.64286  0.07143  0.64286  2.46429 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## stage1.icept  97.0000     0.9904  97.936  < 2e-16 ***
## stage1.slope   1.8214     0.2215   8.224 5.02e-06 ***
## stage2.icept  82.5000     0.7565 109.060  < 2e-16 ***
## stage2.slope   0.2857     0.1808   1.580    0.142    
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 1.172 on 11 degrees of freedom
## Multiple R-squared:  0.9999,    Adjusted R-squared:  0.9998 
## F-statistic: 2.043e+04 on 4 and 11 DF,  p-value: < 2.2e-16

or given that stage2.slope is not significant we could consider dropping that term. We can optionally replace the fm2<- line with the equivalent commented out line.

# fm2 <- update(fm, . ~ . - stage2.slope)
fm2 <- lm(score ~ stage1.icept + stage1.slope + stage2.icept + 0)
summary(fm2)

## Call:
## lm(formula = score ~ stage1.icept + stage1.slope + stage2.icept + 0
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.714 -1.125  0.500  0.500  2.464 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## stage1.icept  97.0000     1.0504  92.347  < 2e-16 ***
## stage1.slope   1.8214     0.2349   7.755 5.16e-06 ***
## stage2.icept  83.5000     0.4394 190.028  < 2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 1.243 on 12 degrees of freedom
## Multiple R-squared:  0.9998,    Adjusted R-squared:  0.9998 
## F-statistic: 2.422e+04 on 3 and 12 DF,  p-value: < 2.2e-16


plot(score ~ year)
lines(fitted(fm) ~ year, col = "red")
lines(fitted(fm2) ~ year, col = "blue", lty = 2, lwd = 2)
legend("topright", c("fm", "fm2"), col = c("red", "blue"), lty = 1:2, lwd = 1:2)

given your data as dataframe d:

d <- 
  data.frame(
    year = c(2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022),
    score = c(85,  85,  88,  88,  94,  94,  94,  82,  82,  84,  84,  84,  84,  84, 84)
  )

• with base R, you could plot the whole set, then split the data at the desired point, and Map the resulting list of dataframes to a list of appropriate objects (here: the ablines for the coefficients of separate models). Add the ablines as you go along:

plot(score ~ year, data = d) 

d |>
  split(list(d$year >= 2015)) |>
  Map(f = \(chunk) abline(coef(lm(score ~ year, data = chunk))))

• alternative, you could use the split criterion for groupwise linear smoothing in ggplot:

library(ggplot2)

d |>
  ggplot(aes(year, score, group = year >= 2015)) +
  geom_point() +
  geom_smooth(method = 'lm',
              se = FALSE ## hide confidence bands
  )