On 2010-06-10 8:27, Graves, Gregory wrote:
What am I failing to understand here?
Several things; see below.
The script below works fine if the dataset being used is
DNase1<- DNase[ DNase$Run == 1, ] per the example given in
help(nlrob).
Obviously, I am trying to understand how to use nls and nlrob to fit
curves to data using R.
#package=DAAG
attach(codling)
plot(pobs~dose)
#next command returns 'step factor reduced below min factor error'
m.nls<- nls( pobs ~ a/(1 + exp(( b - log(dose) )/c ) ),
data = codling,
start = list( a = 3, b = 0, c = 1 ),
trace = TRUE )
s<-seq(min(pobs), max(pobs), .01)
If pobs is your y-value (response) then what is 's' for?
p.nls<-predict(m.nls,list(pobs=s))
Ah, you want predicted *response* values; so feed
appropriate *predictor* values to predict().
But I don't see how this would not give you an error
anyway if your nls() call resulted in an error.
lines(s,p.nls,col='blue') #generates 'x and y lengths differ' error
This error would result from incorrectly generating p.nls,
(but I don't see how you got R to give you that anyway).
I have 4 suggestions:
1. Define a new variable ldose=log(dose);
It's usually less confusing to work on
the log-scale in these cases.
Then plot pobs vs ldose.
2. From the plot you should notice that the
left asymptote is not likely to be zero
(which is what your model assumes). It's
roughly 0.2.
3. The right asymptote should probably be 1.0
since pobs is a proportion. So a start value
of a=3 in your model makes no sense as you
would also quickly see if you plotted the
curve represented by your model on the
plot of the data (use: curve(...,add=TRUE)).
4. Try this model:
pobs ~ (A + exp((ldose - B)/C)) / (1 + exp((ldose - B)/C))
with starting values: A = 0.2, B = 3, C = 1
-Peter Ehlers
Gregory A. Graves, Lead Scientist
Everglades REstoration COoordination and VERification (RECOVER)
Restoration Sciences Department
South Florida Water Management District
Phones: DESK: 561 / 682 - 2429
CELL: 561 / 719 - 8157
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org 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.
