[R] Trying to use segmented in a function

[Rob Knell]

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.