Re: [R] Starting estimates for nls Exponential Fit
I have tried the method proposed by Dave, and I must say it works very well.
Not to yield starting estimates for an nls-fit, but as an independent method
for calculating E (which is, by the way, the only paramater that I am
actually interested in).
The calculated values for E (Esp[3]) are on average 1.40e-06 lower than the
values produced by the nls fit (calculated over 96 reactions). I am using
the found E values to compenstate for differences between reaction
efficiency in genetic quantification assays. An efficiency difference of a
few %'s (+/- 0.05 in absolute value) can cause quantification differences of
several dozen percentages. Since I am in the process of comparing different
efficiency calculation methods (of which exopnential fit is one) I'll
compare the nls-fit results with the results by Dave's method to see if they
yield any significant differnces. I'll probably post the results by the end
of next week (I have an alarming number of reports still due, for which I am
recieving increasingly frequent angry looks by my coworkers. I had the
revision of an article as an excuse for further procrastination, but now
that it has been accepted I'll temporarly have to shift the focus of my
work).
thanks a lot for all your time effort
Antoon
dave fournier wrote:
Thanks to Dennis Murphy who pointed out that ExponCycles
is undefined. It is an R gotcha. I had shortened the name but R still
remembered it so the script worked but only on my computer.
This should fix that.
ExponValues=c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,2468.32,27
78.47)
Expon=c(17,18,19,20,21,22,23,24,25)
# Example starting estimate calculation
E=1000.0
y1=2018
yn=2778.47
nobs=9
#keep y1 and yn fixed and get initial value for E
Esp1 - optim(c(E=E),method =BFGS,
function(x)
{
E=x[1]
a=(yn-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];
diff=ExponValues-(Y0+a*E^Expon)
return(1000*sum(diff*diff))
})$par
E=Esp1[1]
Esp - optim(c(y1=y1,yn=yn,E=E),method =BFGS,
function(x)
{
E=x[3]
a=(x[2]-x[1])/(E^Expon[nobs]-E^Expon[1])
Y0=x[1]-a*E^Expon[1];
diff=ExponValues-(Y0+a*E^Expon)
return(1000*sum(diff*diff))
})$par
y1=Esp[1]
y2=Esp[2]
E=Esp[3]
a=(y2-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];
--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com
__
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.
--
View this message in context:
http://n4.nabble.com/Starting-estimates-for-nls-Exponential-Fit-tp932230p955249.html
Sent from the R help mailing list archive at Nabble.com.
__
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.
Re: [R] Starting estimates for nls Exponential Fit
Hi Antoon, now that you mention trying out different methods, maybe you should consider fitting a sigmoidal curve to the entire dataset and not only the exponential part (which constitutes a very small dataset) as seems to have been the endeavour that initiated the posting to R-help. One options would then be to use the package qpcR and simply fitting a sigmoidal model using for instance drm() in the package drc or nls() combined with SSfpl(). Just a suggestion. Christian __ 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.
Re: [R] Starting estimates for nls Exponential Fit
Thanks everybody, This has been quite helpful, the problem remains tricky but at least now I've got a version of my script that handles all my reactions without error. The DEoptim solution produced good starting values for a lot of reactions, but sadly not for all. I now use scaled parameters and allow more iterations than the standard 50 (1000). I doubt it is completely stable and I think I will run into reactions that will fail the fit, but for now everything works fine so I won't be complaining. Thanks a lot, Antoon Prof. John C Nash wrote: Kate is correct. The parameter scaling helps quite a bit, but not enough to render the problem nice so that many reasonable starting points will give useful results. Indeed, a run using all.methods=TRUE in our optimx package (on r-forge at http://r-forge.r-project.org/R/?group_id=395) gives par fvalues method fns grs itns conv 4 2.194931, 1.01, 1.27 566407.6 nlmNA NA 30 1 2.1949335, -9.0413113, 0.7516435 566407.6 BFGS37 6 NULL0 2 2.1950009, 0.2548318, 1.1163498 566404.6 Nelder82 NA NULL0 3 1.974226, 1.829957, 1.681338 2409.771 SANN 1 NA NULL0 6 1.9904045, 0.6151421, 1.7559401 1696.497 ucminf51 51 NULL0 5 1.9906488, 0.6012996, 1.7575365 1696.248 nlminb80 166 510 KKT1 KKT2 4 FALSE FALSE 1 TRUE FALSE 2 FALSE FALSE 3 FALSE TRUE 6 FALSE TRUE 5 FALSE TRUE A bit of a dog's dinner. On the positive side, this is a simple but nasty problem to illustrate lots of the issues with nonlinear parameter estimation. JN Katharine Mullen wrote: You used starting values: pa - c(1,2,3) but both algorithms (port and Gauss-Newton) fail if you use the slightly different values: pa - c(1,2,3.5) Scaling does not fix the underlying sensitivity to starting values. pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both alg. also fail if there is insufficient spread (e.g., both fail for pa - c(1,1.25,1.5) ). For the record, DE uses (by default at least) a random start and for bad starts will sometimes fail to give good results; decrease the probability this happens by increasing itermax from the default itermax=200, as in: ss - DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10), control=list(trace=FALSE, itermax=1000)) On Wed, 2 Dec 2009, Prof. John C Nash wrote: Kate Mullen showed one approach to this problem by using DEOptim to get some good starting values. However, I believe that the real issue is scaling (Warning: well-ridden hobby-horse!). With appropriate scaling, as below, nls does fine. This isn't saying nls is perfect -- I've been trying to figure out how to do a nice job of helping folk to scale their problems. Ultimately, it would be nice to has an nls version that will do the scaling and also watch for some other situations that give trouble. Cheers, JN ## JN test rm(list=ls()) ExponValues - c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17, 2468.32,2778.47) ExponCycles - c(17,18,19,20,21,22,23,24,25) mod - function(x) x[1] + x[2]*x[3]^ExponCycles modj - function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles) fun - function(x) sum((ExponValues-mod(x))^2) pa-c(1,2,3) lo-c(0,0,0) up-c(20,20,20) names(pa) - c(Y0, a, E) ## fit w/port and GN Emodjn- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, algorithm='port', lower=lo, upper=up, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) Emodjn1 - nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) ## fit matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type=l) __ 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. __ 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. -- View this message in context: http://n4.nabble.com/Starting-estimates-for-nls-Exponential-Fit-tp932230p954292.html Sent from the R help mailing list archive at Nabble.com. __ 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.
Re: [R] Starting estimates for nls Exponential Fit
I thought maybe my suggestion for reparameterizing this simple problem
was ignored because I didn't supply R code for the problem.
Here it is using optim for the optimization. It converges trivially with
an initial value for E of 1000.
As I stated before, there is nothing at all difficult about this
problem. You simply need to parameterize it properly. Of course that is
not to say that rescaling is not useful as well, but the important thing
is to parameterize the model properly.
ExponValues=c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,2468.32,27
78.47)
Expon=c(17,18,19,20,21,22,23,24,25)
# Example starting estimate calculation
E=1000.0
y1=2018
yn=2778.47
nobs=9
#keep y1 and yn fixed and get initial value for E
Esp1 - optim(c(E=E),method =BFGS,
function(x)
{
E=x[1]
a=(yn-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];
diff=ExponValues-(Y0+a*E^ExponCycles)
return(1000*sum(diff*diff))
})$par
E=Esp1[1]
Esp - optim(c(y1=y1,yn=yn,E=E),method =BFGS,
function(x)
{
E=x[3]
a=(x[2]-x[1])/(E^Expon[nobs]-E^Expon[1])
Y0=x[1]-a*E^Expon[1];
diff=ExponValues-(Y0+a*E^ExponCycles)
return(1000*sum(diff*diff))
})$par
y1=Esp[1]
y2=Esp[2]
E=Esp[3]
a=(y2-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];
--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com
__
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.
Re: [R] Starting estimates for nls Exponential Fit
Thanks to Dennis Murphy who pointed out that ExponCycles
is undefined. It is an R gotcha. I had shortened the name but R still
remembered it so the script worked but only on my computer.
This should fix that.
ExponValues=c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,2468.32,27
78.47)
Expon=c(17,18,19,20,21,22,23,24,25)
# Example starting estimate calculation
E=1000.0
y1=2018
yn=2778.47
nobs=9
#keep y1 and yn fixed and get initial value for E
Esp1 - optim(c(E=E),method =BFGS,
function(x)
{
E=x[1]
a=(yn-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];
diff=ExponValues-(Y0+a*E^Expon)
return(1000*sum(diff*diff))
})$par
E=Esp1[1]
Esp - optim(c(y1=y1,yn=yn,E=E),method =BFGS,
function(x)
{
E=x[3]
a=(x[2]-x[1])/(E^Expon[nobs]-E^Expon[1])
Y0=x[1]-a*E^Expon[1];
diff=ExponValues-(Y0+a*E^Expon)
return(1000*sum(diff*diff))
})$par
y1=Esp[1]
y2=Esp[2]
E=Esp[3]
a=(y2-y1)/(E^Expon[nobs]-E^Expon[1])
Y0=y1-a*E^Expon[1];
--
David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 250-655-3364
http://otter-rsch.com
__
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.
Re: [R] Starting estimates for nls Exponential Fit
Kate Mullen showed one approach to this problem by using DEOptim to get some good starting values. However, I believe that the real issue is scaling (Warning: well-ridden hobby-horse!). With appropriate scaling, as below, nls does fine. This isn't saying nls is perfect -- I've been trying to figure out how to do a nice job of helping folk to scale their problems. Ultimately, it would be nice to has an nls version that will do the scaling and also watch for some other situations that give trouble. Cheers, JN ## JN test rm(list=ls()) ExponValues - c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17, 2468.32,2778.47) ExponCycles - c(17,18,19,20,21,22,23,24,25) mod - function(x) x[1] + x[2]*x[3]^ExponCycles modj - function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles) fun - function(x) sum((ExponValues-mod(x))^2) pa-c(1,2,3) lo-c(0,0,0) up-c(20,20,20) names(pa) - c(Y0, a, E) ## fit w/port and GN Emodjn- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, algorithm='port', lower=lo, upper=up, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) Emodjn1 - nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) ## fit matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type=l) __ 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.
Re: [R] Starting estimates for nls Exponential Fit
You used starting values: pa - c(1,2,3) but both algorithms (port and Gauss-Newton) fail if you use the slightly different values: pa - c(1,2,3.5) Scaling does not fix the underlying sensitivity to starting values. pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both alg. also fail if there is insufficient spread (e.g., both fail for pa - c(1,1.25,1.5) ). For the record, DE uses (by default at least) a random start and for bad starts will sometimes fail to give good results; decrease the probability this happens by increasing itermax from the default itermax=200, as in: ss - DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10), control=list(trace=FALSE, itermax=1000)) On Wed, 2 Dec 2009, Prof. John C Nash wrote: Kate Mullen showed one approach to this problem by using DEOptim to get some good starting values. However, I believe that the real issue is scaling (Warning: well-ridden hobby-horse!). With appropriate scaling, as below, nls does fine. This isn't saying nls is perfect -- I've been trying to figure out how to do a nice job of helping folk to scale their problems. Ultimately, it would be nice to has an nls version that will do the scaling and also watch for some other situations that give trouble. Cheers, JN ## JN test rm(list=ls()) ExponValues - c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17, 2468.32,2778.47) ExponCycles - c(17,18,19,20,21,22,23,24,25) mod - function(x) x[1] + x[2]*x[3]^ExponCycles modj - function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles) fun - function(x) sum((ExponValues-mod(x))^2) pa-c(1,2,3) lo-c(0,0,0) up-c(20,20,20) names(pa) - c(Y0, a, E) ## fit w/port and GN Emodjn- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, algorithm='port', lower=lo, upper=up, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) Emodjn1 - nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) ## fit matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type=l) __ 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. __ 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.
Re: [R] Starting estimates for nls Exponential Fit
Kate is correct. The parameter scaling helps quite a bit, but not enough to render the problem nice so that many reasonable starting points will give useful results. Indeed, a run using all.methods=TRUE in our optimx package (on r-forge at http://r-forge.r-project.org/R/?group_id=395) gives par fvalues method fns grs itns conv 4 2.194931, 1.01, 1.27 566407.6 nlmNA NA 30 1 2.1949335, -9.0413113, 0.7516435 566407.6 BFGS37 6 NULL0 2 2.1950009, 0.2548318, 1.1163498 566404.6 Nelder82 NA NULL0 3 1.974226, 1.829957, 1.681338 2409.771 SANN 1 NA NULL0 6 1.9904045, 0.6151421, 1.7559401 1696.497 ucminf51 51 NULL0 5 1.9906488, 0.6012996, 1.7575365 1696.248 nlminb80 166 510 KKT1 KKT2 4 FALSE FALSE 1 TRUE FALSE 2 FALSE FALSE 3 FALSE TRUE 6 FALSE TRUE 5 FALSE TRUE A bit of a dog's dinner. On the positive side, this is a simple but nasty problem to illustrate lots of the issues with nonlinear parameter estimation. JN Katharine Mullen wrote: You used starting values: pa - c(1,2,3) but both algorithms (port and Gauss-Newton) fail if you use the slightly different values: pa - c(1,2,3.5) Scaling does not fix the underlying sensitivity to starting values. pa[3] in particular cannot be above ~3.15 for GN and ~3.3 for port; both alg. also fail if there is insufficient spread (e.g., both fail for pa - c(1,1.25,1.5) ). For the record, DE uses (by default at least) a random start and for bad starts will sometimes fail to give good results; decrease the probability this happens by increasing itermax from the default itermax=200, as in: ss - DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10), control=list(trace=FALSE, itermax=1000)) On Wed, 2 Dec 2009, Prof. John C Nash wrote: Kate Mullen showed one approach to this problem by using DEOptim to get some good starting values. However, I believe that the real issue is scaling (Warning: well-ridden hobby-horse!). With appropriate scaling, as below, nls does fine. This isn't saying nls is perfect -- I've been trying to figure out how to do a nice job of helping folk to scale their problems. Ultimately, it would be nice to has an nls version that will do the scaling and also watch for some other situations that give trouble. Cheers, JN ## JN test rm(list=ls()) ExponValues - c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17, 2468.32,2778.47) ExponCycles - c(17,18,19,20,21,22,23,24,25) mod - function(x) x[1] + x[2]*x[3]^ExponCycles modj - function(x) (1000*x[1] + 0.001*x[2]*x[3]^ExponCycles) fun - function(x) sum((ExponValues-mod(x))^2) pa-c(1,2,3) lo-c(0,0,0) up-c(20,20,20) names(pa) - c(Y0, a, E) ## fit w/port and GN Emodjn- nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, algorithm='port', lower=lo, upper=up, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) Emodjn1 - nls(ExponValues ~ (1000*Y0 + 0.001*a*(E^ExponCycles)), start=pa, control=list(maxiter=1000, trace=TRUE, warnOnly=TRUE)) ## fit matplot(cbind(ExponValues, fitted(Emodjn), fitted(Emodjn1)),type=l) __ 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. __ 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.
Re: [R] Starting estimates for nls Exponential Fit
Figuring out the best parameterization for this kind of model is a bit tricky until you get the hang of it. Let the function be y_t = y_0 + alpha * E^t where uppercase Y_t denotes an observed value and lower case y_t is a predicted value. Index the times by t_1 t_n WLOG assume that t_1 is the smallest time and t_n is the largest time. Now we already have a good idea what the predicted values y_1 and y_n should be as we have observations for them. We have the two equations y_1 = y_0 + alpha * E^t_1 y_n = y_0 + alpha * E^t_n we can solve these for alpha in terms of y_1,y_n,and E giving alpha = (y_n-y_1)/(E^t_n -E^t_1) (1) and solve for y_0 as y_0 = y_1 - alpha * E^t_1 using (1) for alpha Now use the good estimates Y_1 and Y_n as the starting values for y_1 and y_n and try some reasonable value for E (say 0.1 E 100) Do the minimization in two stages first holding y_1 and y_n fixed and then estimate y_1,y_n,and E together. This converges in less than a second using AD Model Builder and the starting values (large value for E. 2018.34 2778.47 exp(10.0) where I deliberately used exp(10) as a large initial value for E. The parameter estimates together with the est std devs are 1 y_1 1.9994e+03 3.9177e-01 2 y_n 2.7881e+03 6.7557e-01 3 log_E 5.6391e-01 1.2876e-03 4 alpha 6.0130e-04 1.9398e-05 5 y_0 1.9906e+03 4.5907e-01 6 E 1.7575e+00 4.3935e-02 There are problems for E near 1 which need to be dealt with if you have to go there, but that is just a technicality. This idea also works well for a logistic say 3 4 or 5 parameter. -- David A. Fournier P.O. Box 2040, Sidney, B.C. V8l 3S3 Canada Phone/FAX 250-655-3364 http://otter-rsch.com __ 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.
Re: [R] Starting estimates for nls Exponential Fit
If you could reformulate your model as alpha * (y0 + E^t) then you could
use nls with alg=plinear (alpha then would be eliminated from the
nonlinear param and treated as conditionally linear) and this would help
with convergence.
Else you can try package DEoptim to get the starting values; the advantage
over optim is that you need then lower and upper bounds on the parameters,
not more starting values; the bounds however should be appropriate and as
tight as possible.
Also: by default nls uses max. 50 iter. Depending on where you start off
you may need more than this.
##
ExponValues - c(2018.34,2012.54,2018.85,2023.52,2054.58,2132.61,2247.17,
2468.32,2778.47)
ExponCycles - c(17,18,19,20,21,22,23,24,25)
library(DEoptim)
mod - function(x) x[1] + x[2]*x[3]^ExponCycles
fun - function(x) sum((ExponValues-mod(x))^2)
ss - DEoptim(fun, lower=rep(0,3), upper=c(10e7, 10, 10),
control=list(trace=FALSE))
pa - ss$optim$bestmem
## now have par est that give OK fit
matplot(cbind(ExponValues, mod(pa)),type=l)
names(pa) - c(Y0, a, E)
## fit w/port and GN
Emodel- nls(ExponValues ~ (Y0 + a*(E^ExponCycles)),
start=pa, algorithm=port,
control=list(maxiter=1000, warnOnly=TRUE))
Emodel1 - nls(ExponValues ~ (Y0 + a*(E^ExponCycles)), start=pa,
control=list(maxiter=1000, warnOnly=TRUE))
## fit
matplot(cbind(ExponValues, fitted(Emodel), fitted(Emodel1)),type=l)
On Tue, 1 Dec 2009, Anto wrote:
Hello everyone,
I have come across a bit of an odd problem: I am currently analysing data
PCR reaction data of which part is behaving exponential. I would like to fit
the exponential part to the following:
y0 + alpha * E^t
In which Y0 is the groundphase, alpha a fluorescence factor, E the
efficiency of the reaction t is time (in cycles)
I can get this to work for most of my reactions, but part fails the nls fit
(Convergence failure: iteration limit reached without convergence). This
mainly has to do with the starting estimates for the nls function. I have
used various 'smart' ways of calculating starting estimates (grid searches,
optim(), etc.) but however close the starting estimates are to the actual
values, nls still fails for many datasets.
The weirdest thing is, I have one set of starting estimates (alpha= 0.5 and
E=1.85) which are totaly arbitray and these DO work for, say 99% of the
datasets. I don't know why this is and I don't why my 'good estimates' do
not work. I am desperatly seeking a way of calculating working starting
estimates for my nls function. Can anybody give me a hand?
thanks,
Anto
R version 2.9.2
example dataset:
ExponValues
[1] 2018.34 2012.54 2018.85 2023.52 2054.58 2132.61 2247.17 2468.32 2778.47
ExponCycles
[1] 17 18 19 20 21 22 23 24 25
Example starting estimate calculation
Y0 - ExponValues[1]
E -
weighted.mean((ExponValues-eY0)[-1][-1]/(ExponValues-eY0)[-1][-(length(ExponValues)-1)],(1-abs(seq(-1,1,length=(length(ExponValues)-2)))^3)^3)
alpha -
weighted.mean((ExponValues[-1]-ExponValues[-length(ExponValues)])/((E^ExponCycles[-1])-(E^ExponCycles[-length(ExponCycles)])),(1-abs(seq(-1,1,length=(length(ExponCycles)-1)))^3)^3)
Optional parameter optimisation:
Esp -
optim(c(Y0=eY0,a=alpha,E=E),function(x){return(sum(abs(ExponValues-(x[1]+x[2]*x[3]^ExponCycles})$par
nls function:
Emodel-try(nls(ExponValues ~ (Y0 + a*(E^ExponCycles)) , start= Esp,
algorithm=port),silent=TRUE)
--
View this message in context:
http://n4.nabble.com/Starting-estimates-for-nls-Exponential-Fit-tp932230p932230.html
Sent from the R help mailing list archive at Nabble.com.
__
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.
__
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.
