Finding the change point in data from a piecewise linear function Greetings,
I'm performing research that will help determine the size of observed space and the time elapsed since the big bang. Hopefully you can help!
I have data conforming to a piecewise linear function on which I want to perform two linear regressions. There is a point at which the slope and intercept change, and I need to (write a program to) find this point.
Thoughts?
 A: R package strucchange might help you. Look at the vignette, it has a nice overview how to solve similar problems. 
A: If the number of points is not too big, you may try all possibilities. Let's assume that the points are $X_i=(x_i,y_i)$ where $i=1,..,N$. Than, you may loop with $j$ from $2$ to $N-2$ and fit two lines to both $\{X_1,...,X_j\}$ and $\{X_{(j+1)},...,X_N\}$. Finally, you pick $j$ for which the sum of sum of squared residuals for both lines is minimal.
A: This is an (offline) changepoint detection problem.  Our previous discussion provides references to journal articles and R code.  Look first at the Barry and Hartigan "product partition model," because it handles changes in slope and has efficient implementations.
A: Also the segmented package has helped me with similar problems in the past.
A: I built on mbq's answer that searching for all the possibilities. Furthermore, I do this:


*

*Check for the significance of the two piecewise models to make sure the coefficients are significant

*Check the difference to the sum of squared residuals for the full model

*Confirm my model visually (make sure it's not something nonsense)


Why check for the significance? That's because the point with the minimum SSE is meaningless if either of piecewise model fits the data very poorly. This can happen for two highly correlated variables without a clear breakpoint where slopes change.
Let's check this simple approach with an easy test case:
x <- c(-50:50)
y <- abs(x)
plot(x,y,pch=19)


The breakpoint is obviously zero. Use the following R script:
f <- function(x, y)
{
    d <- data.frame(x=x, y=y)
    d <- d[order(x),]
    r <- data.frame(k=rep(0,length(x)-4), sums=rep(0,length(x)-4))

    plm <- function(i)
    {
        d1 <- head(d,i)
        d2 <- tail(d,-i)

        # Make sure we've divided the region perfectly        
        stopifnot(nrow(d1)+nrow(d2) == nrow(d))

        m1 <- lm(y~x, data=d1)
        m2 <- lm(y~x, data=d2)

        r <- list(m1, m2)
        r
    }

    lapply(2:(nrow(d)-3), function(i)
    {
        r$k[i-2] <<- d[i,]$x

        # Fit two piecewise linear models
        m <- plm(i)

        # Add up the sum of squares for residuals
        r$sums[i-2] <<- sum((m[[1]]$residuals)^2) + sum((m[[2]]$residuals)^2)
    })

    b <- r[which.min(r$sums),]    
    b
}

Fit piecewise linear models for all possible combinations:
f(x,y)
   k sums
   0    0

If we check the coefficients for the two optimal models, they will be highly significant. Their R2 will be also very high.
A: The mcp package can do this. Say your data is
First, let's simulate some data:
df = data.frame(x = 1:100,
                y = c(rnorm(40, 10 + (1:40)*0.5),
                      rnorm(60, 10 + 40*0.5 -8 + (1:60)*0.2)))

Now let's see if we can recover the change point at 40 (and the parameter values) using mcp:
model = list(
  y ~ 1 + x,  # linear segment
  ~ 1 + x  # another linear segment
)
library(mcp)
fit = mcp(model, df)

Plot it. The gray lines are random draws from the fit, showing that it captures the trend. The blue curve is the estimated change point location:

Let's see the individual parameter estimates. int_ are intercepts, x_ are slopes on x, and cp_ are change points:
summary(fit)

Population-level parameters:
    name  mean lower upper Rhat n.eff
    cp_1 40.48 40.02 41.00    1  2888
   int_1 11.12  9.11 13.17    1   778
   int_2 21.72 20.09 23.49    1   717
 sigma_1  3.23  2.76  3.69    1  5343
     x_1  0.46  0.36  0.54    1   724
     x_2  0.21  0.16  0.26    1   754

Disclaimer: I am the developer of mcp.
