R help - Aug 2015 - Piecewise regression using segmented package plotted in xyplot

Hi,

xyplot(threshold ~ age |frequency.a, data=rage,
groups=HL,
  cex=0.5,
  layout=c(7,4),
par.strip.tex=list(cex=0.8),
  xlab="Age (years)",
  ylab="Threshold (dB SPL)",
  na.rm="TRUE",
  panel=function(x,y,groups,...) {
panel.superpose(x,y,groups=HL,...)
# panel.abline(segmented(lm(threshold~age),seg.Z = ~age, psi = NA, control =
seg.control(K=1)))
panel.abline(lm(threshold~age))
  },
  )

Is there anyway to make the commented line work in lattice? I need to fit my
data in each panel using piecewise regression. Being able to use segmented would
make it easy.

The code above works to give me a linear fit.

Thanks for your help in advance.

Regards,
Sumit



	[[alternative HTML version deleted]]

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

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

Hi

Without a reproducible example I am only guessing 

Try this 

xyplot(threshold ~ age |frequency.a, data=subset(rage !is.na(threshold} &
!is.na(age)),
       groups = HL,
       cex=0.5,
       layout=c(7,4),
       par.strip.tex=list(cex=0.8),
       xlab="Age (years)",
       ylab="Threshold (dB SPL)",
       # na.rm="TRUE", # see above
       type = c("p","l"),  # re-read ?xyplot: from the panel
section:
panel.abline(lm(threshold~age))
       panel = panel.superpose,
       panel.groups =function(x,y,groups,...) {

                        panel.xyplot(x,y,...)
                        # panel.abline(segmented(lm(threshold~age),seg.Z ~age,
psi = NA, control = seg.control(K=1)))

                      },
)

I do not know what the commented line means-- you may need to look at the
subscripts section as well as panel.groups of ?xyplot
and ?panel.xyplot

Regards

Duncan

Duncan
Duncan Mackay
Department of Agronomy and Soil Science
University of New England
Armidale NSW 2351
Email: home: mackay at northnet.com.au


-----Original Message-----
From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Sumitrajit
Dhar
Sent: Friday, 28 August 2015 14:12
To: r-help at r-project.org
Subject: [R] Piecewise regression using segmented package plotted in xyplot

Hi,

xyplot(threshold ~ age |frequency.a, data=rage,
groups=HL,
  cex=0.5,
  layout=c(7,4),
par.strip.tex=list(cex=0.8),
  xlab="Age (years)",
  ylab="Threshold (dB SPL)",
  na.rm="TRUE",
  panel=function(x,y,groups,...) {
panel.superpose(x,y,groups=HL,...)
# panel.abline(segmented(lm(threshold~age),seg.Z = ~age, psi = NA, control
seg.control(K=1)))
panel.abline(lm(threshold~age))
  },
  )

Is there anyway to make the commented line work in lattice? I need to fit my
data in each panel using piecewise regression. Being able to use segmented
would make it easy.

The code above works to give me a linear fit.

Thanks for your help in advance.

Regards,
Sumit



	[[alternative HTML version deleted]]

______________________________________________
R-help at 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.