[R] Piecewise regression using segmented package plotted in xyplot

[S Ellison]

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}}