R help - Sep 2009 - predicting from segmented regression

Hello

I'm having trouble figuring out how to use the output of
"segmented()"
with a new set of predictor values.
Using the example of the help file:

??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-.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))

Now I would like to create a new set of values between 0 and 100 and
predict the value of the piecewise linear function at those. Would you
please help me figuring out how I would do that?

Thank you so much in advance

claudia




-- 
Claudia Tebaldi
Research Scientist, Climate Central
http://www.climatecentral.org
& Adjunct Professor
Department of Statistics - UBC Vancouver
office 604 822 3595 (Canadian area code)
cell 303 775 5365 (US area code)

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