There isn't an abline method for segmented, and even if there were you'd need
segments() for a segmented line plot. You're going to have to roll your own.
That will need a function to extract the break locations and predicted values
at those points
I don't have your data, so I can't do one specifically for you. But here's a
version that works on the data in the first example in ?segmented. I've
deliberately separated segmented model fitting and line segment extraction
(get.segments) from the panel function that plots them, as that would then work
with the base graphics segments function
#Construct the data set (from ?segmented, though the z covariate is not used
here)
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)
# Function to fit a simple model using segmented
# and then get line segment data from the model
# For the latter, I've used code from plot. segmented to locate extremes
# and break points (in psi) and then just used
# predict() rather crudely to get the corresponding ordinate values
# you would have to do something more clever
# if your model is not just y~x
get.segments <- function(x, y, term, ...) {
#Returns a data frame of line segment coordinates
#from a simple segmented() model
S <- segmented(lm(y~x), seg.Z=~x,psi=list(x=c(30,60)),
control=seg.control(display=FALSE))
if(missing(term)) term <- S$nameUV$Z[1]
x<-with(S, sort(c(psi[,'Est.'], rangeZ[,term])))
#Then cheat a bit
new <- data.frame(x=unname(x))
names(new) <- term
y <- predict(S, new=new)
L <- length(x) -1
#return the line segment data in an easy-to-use format
data.frame(x0=x[1:L], y0=y[1:L], x1=x[1:L + 1], y1=y[1:L + 1])
}
#Write a panel function using that...
panel.segmented <- function(x, y, ...) {
## Fit the simple model
m <- get.segments( x, y )
#Plot the segments
with(m, panel.segments(x0, y0, x1, y1, ... ))
}
library(lattice)
xyplot(y~x, data=dati,
panel = function(x, y, ...) {
panel.xyplot(x, y, ...)
panel.segmented(x, y, ...)
}
)
#And just to see if it works for panel-grouped data:
set.seed(1023)
dati$parts <- sample(gl(2, 50))
xyplot(y~x|parts, data=dati,
panel = function(x, y, ...) {
panel.xyplot(x, y, ...)
panel.segmented(x, y, ...)
}
)
S Ellison
*******************************************************************
This email and any attachments are confidential. Any use...{{dropped:8}}
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.
