Plotting a piecewise regression line Is there a way of plotting the regression line of a piecewise model like this, other than using lines to plot each segment separately, or using geom_smooth(aes(group=Ind), method="lm", fill=FALSE) ?
m.sqft <- mean(sqft)
model <- lm(price~sqft+I((sqft-m.sqft)*Ind))
# sqft, price: continuous variables, Ind: if sqft>mean(sqft) then 1 else 0

plot(sqft,price)
abline(reg = model)
Warning message:
In abline(reg = model) :
  only using the first two of 3regression coefficients

Thank you.
 A: The only way I know how to do this easily is to predict from the model across the range of sqft and plot the predictions. There isn't a general way with abline or similar. You might also take a look at the segmented package which will fit these models and provide the plotting infrastructure for you.
Doing this via predictions and base graphics. First, some dummy data:
set.seed(1)
sqft <- runif(100)
sqft <- ifelse((tmp <- sqft > mean(sqft)), 1, 0) + rnorm(100, sd = 0.5)
price <- 2 + 2.5 * sqft
price <- ifelse(tmp, price, 0) + rnorm(100, sd = 0.6)
DF <- data.frame(sqft = sqft, price = price,
                 Ind = ifelse(sqft > mean(sqft), 1, 0))
rm(price, sqft)
plot(price ~ sqft, data = DF)

Fit the model:
mod <- lm(price~sqft+I((sqft-mean(sqft))*Ind), data = DF)

Generate some data to predict for and predict:
m.sqft <- with(DF, mean(sqft))
pDF <- with(DF, data.frame(sqft = seq(min(sqft), max(sqft), length = 200)))
pDF <- within(pDF, Ind <- ifelse(sqft > m.sqft, 1, 0))
pDF <- within(pDF, price <- predict(mod, newdata = pDF))

Plot the regression lines:
ylim <- range(pDF$price, DF$price)
xlim <- range(pDF$sqft, DF$sqft)
plot(price ~ sqft, data = DF, ylim = ylim, xlim = xlim)
lines(price ~ sqft, data = pDF, subset = Ind > 0, col = "red", lwd = 2)
lines(price ~ sqft, data = pDF, subset = Ind < 1, col = "red", lwd = 2)

You could code this up into a simple function - you only need the steps in the two preceding code chunks - which you can use in place of abline:
myabline <- function(model, data, ...) {
    m.sqft <- with(data, mean(sqft))
    pDF <- with(data, data.frame(sqft = seq(min(sqft), max(sqft),
                                            length = 200)))
    pDF <- within(pDF, Ind <- ifelse(sqft > m.sqft, 1, 0))
    pDF <- within(pDF, price <- predict(mod, newdata = pDF))
    lines(price ~ sqft, data = pDF, subset = Ind > 0, ...)
    lines(price ~ sqft, data = pDF, subset = Ind < 1, ...)
    invisible(model)
}

Then:
ylim <- range(pDF$price, DF$price)
xlim <- range(pDF$sqft, DF$sqft)
plot(price ~ sqft, data = DF, ylim = ylim, xlim = xlim)
myabline(mod, DF, col = "red", lwd = 2)

Via the segmented package
require(segmented)
mod2 <- lm(price ~ sqft, data = DF)
mod.s <- segmented(mod2, seg.Z = ~ sqft, psi = 0.5,
                   control = seg.control(stop.if.error = FALSE))
plot(price ~ sqft, data = DF)
plot(mod.s, add = TRUE)
lines(mod.s, col = "red")

With these data it doesn't estimate the breakpoint at mean(sqft), but the plot and lines methods in that package might help you implement something more generic than myabline to do this job for you diretcly from the fitted lm() model.
Edit: If you want segmented to estimate the location of the breakpoint, then set the 'psi' argument to NA:
mod.s <- segmented(mod2, seg.Z = ~ sqft, psi = NA,
                   control = seg.control(stop.if.error = FALSE))

Then segmented will try K = 10 quantiles of sqft, with K being set in seg.control() and which defaults to 10. See ?seg.control for more.
