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
______________________________________________
R-help@r-project.org mailing list
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.
