

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((ExponValueseY0)[1][1]/(ExponValueseY0)[1][(length(ExponValues)1)],(1abs(seq(1,1,length=(length(ExponValues)2)))^3)^3)
alpha < weighted.mean((ExponValues[1]ExponValues[length(ExponValues)])/((E^ExponCycles[1])(E^ExponCycles[length(ExponCycles)])),(1abs(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)


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((ExponValuesmod(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((ExponValueseY0)[1][1]/(ExponValueseY0)[1][(length(ExponValues)1)],(1abs(seq(1,1,length=(length(ExponValues)2)))^3)^3)
> alpha <
> weighted.mean((ExponValues[1]ExponValues[length(ExponValues)])/((E^ExponCycles[1])(E^ExponCycles[length(ExponCycles)])),(1abs(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/StartingestimatesfornlsExponentialFittp932230p932230.html> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.
>
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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: wellridden
hobbyhorse!).
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((ExponValuesmod(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")
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


You used starting values:
pa < c(1,2,3)
but both algorithms (port and GaussNewton) 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: wellridden
> hobbyhorse!).
>
> 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((ExponValuesmod(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")
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.
>
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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 rforge at
http://rforge.rproject.org/R/?group_id=395) gives
par fvalues method fns grs itns conv
4 2.194931, 1.000001, 1.000027 566407.6 nlm NA NA 3 0
1 2.1949335, 9.0413113, 0.7516435 566407.6 BFGS 37 6 NULL 0
2 2.1950009, 0.2548318, 1.1163498 566404.6 Nelder 82 NA NULL 0
3 1.974226, 1.829957, 1.681338 2409.771 SANN 10000 NA NULL 0
6 1.9904045, 0.6151421, 1.7559401 1696.497 ucminf 51 51 NULL 0
5 1.9906488, 0.6012996, 1.7575365 1696.248 nlminb 80 166 51 0
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 GaussNewton) 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: wellridden
>> hobbyhorse!).
>>
>> 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((ExponValuesmod(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")
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/rhelp>> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html>> and provide commented, minimal, selfcontained, reproducible code.
>>
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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_ny_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.9177e01
2 y_n 2.7881e+03 6.7557e01
3 log_E 5.6391e01 1.2876e03
4 alpha 6.0130e04 1.9398e05
5 y_0 1.9906e+03 4.5907e01
6 E 1.7575e+00 4.3935e02
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 2506553364
http://otterrsch.com______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


In reply to this post by Prof J C Nash (U30A)
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 rforge at
http://rforge.rproject.org/R/?group_id=395) gives
par fvalues method fns grs itns conv
4 2.194931, 1.000001, 1.000027 566407.6 nlm NA NA 3 0
1 2.1949335, 9.0413113, 0.7516435 566407.6 BFGS 37 6 NULL 0
2 2.1950009, 0.2548318, 1.1163498 566404.6 Nelder 82 NA NULL 0
3 1.974226, 1.829957, 1.681338 2409.771 SANN 10000 NA NULL 0
6 1.9904045, 0.6151421, 1.7559401 1696.497 ucminf 51 51 NULL 0
5 1.9906488, 0.6012996, 1.7575365 1696.248 nlminb 80 166 51 0
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 GaussNewton) 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: wellridden
>> hobbyhorse!).
>>
>> 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((ExponValuesmod(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")
>>
>> ______________________________________________
>> Rhelp@rproject.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/rhelp>> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html>> and provide commented, minimal, selfcontained, reproducible code.
>>
______________________________________________
Rhelp@rproject.org mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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=(yny1)/(E^Expon[nobs]E^Expon[1])
Y0=y1a*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=(y2y1)/(E^Expon[nobs]E^Expon[1])
Y0=y1a*E^Expon[1];

David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 2506553364
http://otterrsch.com______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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=(yny1)/(E^Expon[nobs]E^Expon[1])
Y0=y1a*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=(y2y1)/(E^Expon[nobs]E^Expon[1])
Y0=y1a*E^Expon[1];

David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 2506553364
http://otterrsch.com______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


I have tried the method proposed by Dave, and I must say it works very well. Not to yield starting estimates for an nlsfit, 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.40e06 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 nlsfit 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=(yny1)/(E^Expon[nobs]E^Expon[1])
Y0=y1a*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=(y2y1)/(E^Expon[nobs]E^Expon[1])
Y0=y1a*E^Expon[1];

David A. Fournier
P.O. Box 2040,
Sidney, B.C. V8l 3S3
Canada
Phone/FAX 2506553364
http://otterrsch.com______________________________________________
Rhelp@rproject.org mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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
Rhelp.
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
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.

