10
$\begingroup$

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.

$\endgroup$
6
$\begingroup$

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.

$\endgroup$
11
  • $\begingroup$ @Gavin (+1) Far more complete response than mine; I just like it. $\endgroup$ – chl Jan 4 '11 at 10:25
  • $\begingroup$ @Gavin The "Via the segmented package" section didn't work for my data. I got a "No breakpoint estimated" after running the segmented command. $\endgroup$ – George Dontas Jan 4 '11 at 10:35
  • $\begingroup$ @gd047: Apologies, there was an error in the code I showed. You need to supply argument seq.Z with a one sided formula of the variable(s) that have a segmented relationship with the response. I've edited my answer to include seq.Z = ~ sqft and added a note about having segmented choose values of psi for you. $\endgroup$ – Gavin Simpson Jan 4 '11 at 11:52
  • $\begingroup$ @gd047 I would like to remove my answer as this one addresses your original question in a more better way. Would mind accepting this one instead of mine? $\endgroup$ – chl Jan 4 '11 at 12:00
  • $\begingroup$ @chl Of course, even though I still get an error : Error in if (model) objF$model <- mf : argument is not interpretable as logical In addition: Warning message: In if (model) objF$model <- mf : the condition has length > 1 and only the first element will be used $\endgroup$ – George Dontas Jan 4 '11 at 12:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.