[R] predicting from segmented regression

[Matthieu Dubois]

dear Claudia, 

I was recently in touch with Vito Muggeo (the developer of the 
segmented package) with a similar question. This is an adapted 
version of his answer to your problem. 
In fact, the essential aspect is that predict.segmented has not 
(yet?) been implemented. Nevertheless, you could use the standard 
predict, as long as you provide the predict function with an adequate
data frame (corresponding to the model frame of the segmented 
regression, see head(model.frame(o)) ). 

I thus completed your code according to Vito's suggestions: 

library(segmented)
set.seed(12)

xx <- 1:100
zz <- runif(100)
yy <- 2 + 1.5*pmax(xx-35,0) - 1.5*pmax(xx-70,0) + 
               15*pmax(zz-0.5,0) + rnorm(100,0,2)
dati <- data.frame(x=xx,y=yy,z=zz)

out.lm<-lm(y~x,data=dati)
o<-## S3 method for class 'lm':
     segmented(out.lm,seg.Z=~x,psi=list(x=c(30,60)),
           control=seg.control(display=FALSE))

# prediction
pred <- data.frame(xx = 1:100)
pred$dummy1 <- pmax(pred$xx - o$psi[1,2], 0)
pred$dummy2 <- pmax(pred$xx - o$psi[2,2], 0)
pred$dummy3 <- I(pred$xx > o$psi[1,2]) * coef(o)[3]
pred$dummy4 <- I(pred$xx > o$psi[2,2]) * coef(o)[4]
names(pred)[-1]<- names(model.frame(o))[-c(1,2)]

# compute the prediction, using standard predict function
# computing confidence intervals further 
# suppose that the breakpoints are fixed
pred <- data.frame(pred, predict(o, newdata= pred, 
                       interval="confidence"))

# plot, just for the fun, and using ggplot2
library(ggplot2)
p <- ggplot(aes(x=xx, y=yy), data=dati) + 
		geom_smooth(aes(y=fit, ymin=lwr, ymax=upr), 
                        data=pred, stat="identity") + 
		geom_point()

breakp <- data.frame(confint(o)) # extract breakpoints to 
                                                            # add it to the plot
names(breakp) <- c('est', 'lwr', 'upr')

p + geom_segment(aes(x=lwr, xend=upr, y=rep(-5,2), 
                      yend=rep(-5,2)), data=breakp) + 
		geom_point(aes(x=est, y=-5), data=breakp) + 
		geom_text(aes(x=est, y=-8, 
                      label=paste(rep("95% CI for breakpoint", 2), 1:2)), 
                      data=breakp, size=3)


HTH

Matthieu