# How can I make my R nls model more stable?

I am new to R, and I am trying to automate some nonlinear modeling of several thousand sets of x,y data. These data are experimental results representing exponential decay. I want to model my data to this formula:

y = (y0 - plateau)*exp(-k*x) + plateau


I have created functions for the model and the initial value estimation:

expDecayModel <- function(predictor, k, y0, plateau) {
(y0 - plateau)*exp(-k*predictor) + plateau
}

expDecayModelInit <- function (mCall, LHS, data)
{
xyframe <- sortedXyData(mCall[["predictor"]], LHS, data)
x <- xyframe$x y <- xyframe$y

y0 <- max(y)
plateau <- min(y)
ymid <- (y0 + plateau)/2

closesty <- y[1]
closestyindex <- 1
for (i in 2:length(y)) {
if (abs(y[i] - ymid) < abs(closesty - ymid)) {
closesty <- y[i]
closestyindex <- i
}
}
xatclosesty = x[closestyindex]
k = 1 / xatclosesty
value = c(k, y0, plateau)
names(value) <- mCall[c("k", "y0", "plateau")]
value
}


I created the self-starter like so:

SSexpDecay <- selfStart(expDecayModel, expDecayModelInit, c("k", "y0", "plateau"))


I then try to execute the model like this:

fitmodel <- nls(y ~ SSexpDecay(x, k, y0, plateau), data = AllData)


For some data sets, this seems to work nicely. For others, it does not. I'd like to know either what I'm doing wrong, or things that I can do to make my automated task more robust. My initial value selection is (intentionally) intended to mimic those in GraphPad Prism, and at least for the data set below, the initial values for my three parameters are identical.

So, what can I do to improve this?

An example data set that fails is at this fiddle.

-

## migrated from stackoverflow.comOct 17 '11 at 16:00

This question came from our site for professional and enthusiast programmers.

the built-in self-starting models are usually sufficient. I just got it to work for your data with SSasymp, so perhaps use that as a base if you need some option that it doesn't provide.

AllData = read.table('data.txt', header=T)

fit = nls(y ~ SSasymp(x, Asym, r0, lrc), data=AllData)

x.new = data.frame(x=with(AllData, seq(min(x), max(x), len=100)))
pred = predict(fit, x.new)

dev.new(width=5, height=5)
plot(AllData)
lines(x.new\$x, pred, col='red')


-