How to do piecewise linear regression with multiple unknown knots? Are there any packages to do piecewise linear regression, which can detect the multiple knots automatically? Thanks.
When I use the strucchange package. I could not detect the change points. I have no idea how it detects the change points. From the plots, I could see there are several points I want it could help me to pick them out.  Could anyone give an example here?
 A: Would MARS be applicable? R has the package earth that implements it.
A: In general, it's a bit odd to want to fit something as piece-wise linear.  However, if you really wish to do so, then the MARS algorithm is the most direct.  It will build up a function one knot at a time; and then usually prunes back the number of knots to combat over-fitting ala decision trees.  You can access the MARS algotithm in R via earth or mda. In general, it's fit with GCV which is not so far removed from the other information criterion (AIC, BIC etc.)
MARS won't really give you an "optimal" fit since the knots are grown one at a time.  It really would be rather difficult to fit a truly "optimal" number of knots since the possible permutations of knot placements would quickly explode.
Generally, this is why people turn towards smoothing splines.  Most smoothing splines are cubic just so you can fool a human eye into missing the discontinuities.  It would be quite possible to do a linear smoothing spline however.  The big advantage of smoothing splines are their single parameter to optimize.  That allows you to quickly reach a truly "optimal" solution without having to search through gobs of permutations.  However, if you really want to seek inflection points, and you have enough data to do so, then something like MARS would probably be your best bet.
Here's some example code for penalized linear smoothing splines in R:
require(mgcv);data(iris);
gam.test <- gam(Sepal.Length ~ s(Petal.Width,k=6,bs='ps',m=0),data=iris)
summary(gam.test);plot(gam.test);

The actual knots chosen would not necessarily correlate with any true inflection points however.
A: There is a pretty nice algorithm described in Tomé and Miranda (1984). 

The proposed methodology uses a least-squares approach to compute the best continuous set of straight lines that fit a given time series, subject to a number of constraints on the minimum distance between breakpoints and on the minimum trend change at each breakpoint.

The code and a GUI are available in both Fortran and IDL from their website: 
http://www.dfisica.ubi.pt/~artome/linearstep.html
A: ... first of all you must to do it by iterations, and under some informative criterion, like AIC AICc BIC Cp; because you can get an "ideal" fit, if number of knots K = number od data  points N, ok. 
... first put K = 0; estimate L = K + 1 regressions, calculate AICc, for instance; 
then assume minimal number of data points at a separate segment, say L = 3 or L = 4, ok
... put K = 1; start from L-th data as the first knot, calculate SS or MLE, ... and step by step the next data point as a knot, SS or MLE, up to the last knot at the N - L data; choose the arrangement with the best fit (SS or MLE) calculate AICc ... 
... put K = 2; ... use all previous regressions (that is their SS or MLE), but step by step divide a single segment  into all possible parts ... choose the arrangement with the best fit (SS or MLE) calculate AICc ... if the last AICc occurs greater then the previous one: stop the iterations ! This is an optimal solution under AICc criterion, ok
A: I once came across a program called Joinpoint. On their website they say it fits a joinpoint model where "several different lines are connected together at the 'joinpoints'". And further: "The user supplies the minimum and maximum number of joinpoints. The program starts with the minimum number of joinpoint (e.g. 0 joinpoints, which is a straight line) and tests whether more joinpoints are statistically significant and must be added to the model (up to that maximum number)."
The NCI uses it for trend modelling of cancer rates, maybe it fits your needs as well.
A: In order to fit to data a piecewise function :

where $a_1 , a_2 , p_1 , q_1, p_2 , q_2 , p_3 , q_3$ are unknown parameters to be approximately computed, there is a very simple method (not iterative, no initial guess, easy to code in any math computer language). The theory given page 29 in paper : https://fr.scribd.com/document/380941024/Regression-par-morceaux-Piecewise-Regression-pdf and from page 30 :

For example, with the exact data provided by Mats Granvik the result is :

Without scattered data, this example is not very signifiant. Other examples with scattered data are shown in the referenced paper.
A: You can use the mcp package if you know the number of change points to infer. It gives you great modeling flexibility and a lot of information about the change points and regression parameters, but at the cost of speed.
The mcp website contains many applied examples, e.g.,
library(mcp)

# Define the model
model = list(
  response ~ 1,  # plateau (int_1)
  ~ 0 + time,    # joined slope (time_2) at cp_1
  ~ 1 + time     # disjoined slope (int_3, time_3) at cp_2
)

# Fit it. The `ex_demo` dataset is included in mcp
fit = mcp(model, data = ex_demo)

Then you can visualize:
plot(fit)


Or summarise:
summary(fit)

Family: gaussian(link = 'identity')
Iterations: 9000 from 3 chains.
Segments:
  1: response ~ 1
  2: response ~ 1 ~ 0 + time
  3: response ~ 1 ~ 1 + time

Population-level parameters:
    name match  sim  mean lower  upper Rhat n.eff
    cp_1    OK 30.0 30.27 23.19 38.760    1   384
    cp_2    OK 70.0 69.78 69.27 70.238    1  5792
   int_1    OK 10.0 10.26  8.82 11.768    1  1480
   int_3    OK  0.0  0.44 -2.49  3.428    1   810
 sigma_1    OK  4.0  4.01  3.43  4.591    1  3852
  time_2    OK  0.5  0.53  0.40  0.662    1   437
  time_3    OK -0.2 -0.22 -0.38 -0.035    1   834

Disclaimer: I am the developer of mcp.
