4
$\begingroup$

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.

$\endgroup$
0

1 Answer 1

9
$\begingroup$

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')

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.