nls function in R within functions and lists

Question

I am working with a package written in R-to fit experimental data to specific models and I want to write my own model and use same package(authors claim it is possible).

So I am digging the source files to find where they define the model functions. I am stuck by the way they use nls. I understood how nls work for the case below

x <- -(1:100)/10
y <- 100 + 10 * exp(x / 2) + rnorm(x)/2
nlmod <- nls(y ~  Const + A * exp(B * x),start=list(Const=5,A=5,B=7))

create a data set and write the model next to it-easy.

However they use nls as shown below:

>  nls(~rescomp(theta = t, d = d, currModel = currModel), 
>     data = list(d = vector(), currModel = currModel))

when I run

>  ~rescomp(theta = t, d = d, currModel = currModel)

to see a formula as in the example above (Const + A * exp(B * x)). I get this

~rescomp(theta = t, d = d, currModel = currModel) 

environment: 0x000000000a7f7530

I want to learn how I can see the resultant formula and how nls work while data is set to list\environment. Any suggestions?

Here is the content of rescomp

  function (theta = vector(), d = vector(), currModel = currModel, 
    currTheta = vector()) 
{
    if (length(currTheta) == 0) 
        currTheta <- getThetaCl(theta, currModel)
    groups <- currModel@groups
    m <- currModel@modellist
    resid <- clpindepX <- list()
    nexp <- length(m)
    for (i in 1:nexp) {
        clpindepX[[i]] <- if (!m[[i]]@clpdep || m[[i]]@getX) 
            getClpindepX(model = m[[i]], theta = currTheta[[i]], 
                multimodel = currModel, returnX = FALSE, rawtheta = theta, 
                dind = 0)
        else matrix()
    }
    for (i in 1:length(groups)) {
        resid[[i]] <- residPart(model = m[[1]], group = groups[[i]], 
            multimodel = currModel, thetalist = currTheta, clpindepX = clpindepX, 
            finished = currModel@finished, returnX = FALSE, rawtheta = theta)
        if (currModel@finished) {
            currModel <- fillResult(group = groups[[i]], multimodel = currModel, 
                thetalist = currTheta, clpindepX = clpindepX, 
                rlist = resid[[i]], rawtheta = theta)
        }
    }
    if (currModel@finished) {
        currModel@fit@nlsres$onls$nclp <- currModel@nclp
        if (currModel@optlist[[1]]@sumnls) {
            if (class(currModel@fit@nlsres$onls) == "nls") 
                class(currModel@fit@nlsres$onls) <- "timp.nls"
            else if (class(currModel@fit@nlsres$onls) == "nls.lm") 
                class(currModel@fit@nlsres$onls) <- "timp.nls.lm"
            else class(currModel@fit@nlsres$onls) <- "timp.optim"
            currModel@fit@nlsres$sumonls <- summary(currModel@fit@nlsres$onls, 
                currModel = currModel, currTheta = currTheta)
        }
        if (currModel@stderrclp) {
            for (i in 1:length(groups)) {
                currModel <- getStdErrClp(group = groups[[i]], 
                  multimodel = currModel, thetalist = currTheta, 
                  clpindepX = clpindepX, rlist = resid[[i]], 
                  rawtheta = theta) 
       }
        if (currModel@stderrclp) {
            for (i in 1:length(groups)) {
                currModel <- getStdErrClp(group = groups[[i]], 
                  multimodel = currModel, thetalist = currTheta, 
                  clpindepX = clpindepX, rlist = resid[[i]], 
                  rawtheta = theta)
            }
        }
    }
    if (currModel@finished && currModel@trilinear) {
        trires <- triResolve(currModel, currTheta)
        currModel <- trires$currModel
        currTheta <- trires$currTheta
    }
    if (currModel@finished && m[[1]]@mod_type == "kin") {
        if (m[[1]]@fullk) {
            for (i in 1:nexp) {
                nocolsums <- length(m[[1]]@lightregimespec) > 
                  0
                eig <- fullKF(currTheta[[i]]@kinpar, currTheta[[i]]@kinscal, 
                  m[[1]]@kmat, currTheta[[i]]@jvec, m[[1]]@fixedkmat, 
                  m[[1]]@kinscalspecial, m[[1]]@kinscalspecialspec, 
                  nocolsums)
                currTheta[[i]]@eigenvaluesK <- eig$values
            }
        }
    }
    if (currModel@finished) {
        return(list(currModel = currModel, currTheta = currTheta))
    }
    if (currModel@algorithm == "optim") 
        retval <- sum(unlist(resid))
    else retval <- unlist(resid)
    retval
}

Putting `rm(list=ls())` in your sample code is not a great idea. — Señor O, Oct 09 '13 at 17:15
@Ben Bolker,I already look into `rescomp` to find the formula but couldn't get it. `rescomp` turns the the content of function.rescomp is too long to write here — behre, Oct 09 '13 at 17:33
`nls` simply evaluates `rescomp` with different starting variables to minimize error. — Señor O, Oct 09 '13 at 18:42

IRTFM · Answer 1 · 2013-10-09T22:08:37.020

1

There is no formula. nls is repeatedly passing a call to a function that is pulling in 'currModel' and parameters (probably theta and d) which are to be minimized based on the scalar that is returned from the rescomp function. Your complaint that the results of typing "rescomp" are "too long to write here" is simply an indication that a) you don't understand that you should be editing your question rather than writing that output in the comments, and b) that your expectations of what is happening are too narrow.

To illustrate writing your first nls problem as a functional form:

myfunc <- function(Const,A,B,y=y,x=x) { abs(y - ( Const + A * exp(B * x)))}

So you will be minimizing the absolute deviation of y from predicted:

nlmod <- nls( ~myfunc(Const,A,B,y=y,x=x) ,start=list(Const=5,A=5,B=7))
nlmod  # same results

edited Oct 09 '13 at 22:08

answered Oct 09 '13 at 18:59

IRTFM

258,963
21
364
487

Thanks for your answer, it is really helpful.I also add the `rescomp` function above – behre Oct 10 '13 at 16:43
If you want to see the formula that was driving the construction of the model then I would suggest you cannot find it in `rescomp` but rather by dissecting the `currModel`-object. My first suggestion would be to try `currModel$call` ... If that is unsuccessful then report the results of `names(currModel)` and `class(currModel)`. – IRTFM Oct 10 '13 at 16:49

nls function in R within functions and lists

1 Answers1