[R] Piecewise regression using segmented package plotted in xyplot

[S Ellison]

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