[R] Trying to use segmented in a function

From: Rob Knell <R.Knell_at_qmul.ac.uk>
Date: Wed 02 Aug 2006 - 21:04:45 EST


Hi folks

I wonder if anyone can help me. I want to run some simulations to see how big a sample size might be necessary to distinguish a curved bivariate relationship (e.g. something that might be best described by a quadratic model) from a relationship that is two straight lines with a sudden change in slope (e.g. something best described by a breakpoint regression). I am using segmented to do the breakpoint regression: this package seems to be the one that most people use for this, as far as I can see.

Since I want to run some simulations, I'm trying to write functions that use segmented, and it's driving me mad. Here's a simple example:

simdata<-function
(Ns=200,Xmean=20,Xsd=5,SdYerr=0.5,Yint=0,threshold=20,slopebelow=0.5,slo peabove=1)
{

	Xs<-rnorm(Ns,Xmean,Xsd)
	Yerr<-rnorm(Ns,0,SdYerr)
	D<-ifelse(Xs<=threshold,0,1)
	XminusX0<-Xs-threshold
	Ys<-Yint+slopebelow*Xs+slopeabove*XminusX0*D+Yerr

	plot(Xs,Ys)
	
	linmod<-lm(Ys~Xs)
	segment<-segmented(linmod,Z=Xs,psi=threshold)

	segment


}

This code should simply simulate some "breakpoint" data, with the change in slope at "threshold" and then fit a model with segmented. If I just use the code for simulating the data, and run that, and then run segmented as normal in R, then I occasionally get an error when it exceeds the maximum iterations, but 99% of the time it will fit a model happily. When I incorporate it into the function, however, it will sometimes fit a model (about 20% of the time) but most of the time I get this:

> test<-simdata()

Error in segmented.lm(linmod, Z = Xs, psi = threshold) :

        (Some) estimated psi out of its range
>

I emphasise that this is using exactly the same code to simulate the data that gives good results when used without segmented in the function. I'm even giving it the exact right value of the breakpoint to start with in its estimation.

If anyone could give me some advice on where I'm going wrong, I would be very pleased to hear it.

Thanks everyone

Rob Knell

School of Biological Sciences
Queen Mary, University of London

'Phone +44 (0)20 7882 7720
Skype Rob Knell
http://www.qmw.ac.uk/~ugbt794
http://www.mopane.org

"The truth is that they have no clue why the beetles had horns, it's the researchers who have sex on the brain and everything has to have a sexual explanation. And this is reasearch?!" Correspondent known as FairOpinion on Neo-Con American website discussing my research.



R-help@stat.math.ethz.ch mailing list
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. Received on Wed Aug 02 22:18:30 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Thu 03 Aug 2006 - 00:18:07 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.