[R] Starting estimates for nls Exponential Fit

[Katharine Mullen]

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 at 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.
>