I was playing with some examples of piecewise regression using lm() and have come across a behavior I'm uncertain about. Below is simple 3-segment dataset. I compare predicted output of a model created by one call to lm() to that of 3 models created by 3 calls to lm(). In case A and B, the results are the same. However, in case C the results differ for the middle segment. Is the output of lm() for case C to be expected or not and if so, why? Thank you, Jill
set.seed(133) y <- c(21:30, rep(15,10), 10:1) + runif(30, -2, 2) x <- 1:30 # Case A: 3 segments, each fit with a constant value plot(x, y) # 3 different lm fits p.c3 <- c(predict(lm(y[1:10]~1)), predict(lm(y[11:20]~1)), predict(lm(y[21:30]~1))) lines(x, p.c3, col="red") # piecewise via lm p.lm1 <- predict(lm(y ~ I(x<=10) + I(x>10 & x<=20) + I(x>20))) lines(x, p.lm1, col="blue") max(abs(p.c3-p.lm1)) # Case B: 3 segments, each fit with a line plot(x, y) # 3 different lm fits p.c3 <- c(predict(lm(y[1:10]~x[1:10])), predict(lm(y[11:20]~x[11:20])), predict(lm(y[21:30]~x[21:30]))) lines(x, p.c3, col="red") # piecewise via lm p.lm1 <- predict(lm(y ~ (x<=10)*x + (x>10 & x<=20)*x + (x>20)*x)) lines(x, p.lm1, col="blue") max(abs(p.c3-p.lm1)) # Cast C: 3 segments - the middle fit with a constant value, the outer by a line plot(x, y) # 3 different lm fits p.c3 <- c(predict(lm(y[1:10]~x[1:10])), predict(lm(y[11:20]~1)), predict(lm(y[21:30]~x[21:30]))) lines(x, p.c3, col="red") # piecewise via lm p.lm1 <- predict(lm(y ~ (x<=10)*x + I(x>10 & x<=20) + (x>20)*x)) lines(x, p.lm1, col="blue") max(abs(p.c3-p.lm1)) ## the single call to lm does not have a constant value fit to the middle section plot(x, p.c3-p.lm1 ) The information contained in this transmission may contain confidential information. If the reader of this message is not the intended recipient, you are hereby notified that any review, dissemination, distribution or duplication of this communication is strictly prohibited. If you are not the intended recipient, please contact the sender by reply email and destroy all copies of the original message. ______________________________________________ [email protected] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
