I perhaps should have added a stronger warning here; note that the model
fitting in my previous post (below) uses explicit initial breakpoints for
segmented (specifically, c(30,60) at line 1 of the get.segments() ). if you
know where yours are, substitute them there. Otherwise, you'd need to use the
automated breakpoint routines documented in ?segmented, typically adjusting the
control list. To be general, you'd need a way to get a choice of at least
number of breakpoints into the model. xyplot doesn't (I think) pass extra
parameters to its panel function, but it will pick up environment parameters so
you _could_ just rely on that, or perhaps on setting something in options(), as
a quick and dirty work-round.
S Ellisaon
--------------------------
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.
