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.
