[Rd] nls problems

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Mon, 04 Jun 2007 09:21:13 -0700


Hi, All:

      The following identifies apparent inconsistencies and proposed fixes between the documentation and the behavior of 'nls'.

      Specifically, the help file says 'data' 'Can also be a list'. However, I got an error when I tried it with a list that could not be coerced to a data.frame. This is documented in the first example below.

      Also, the help file says 'data' is 'optional'. Unfortunately, the second example below bombs, apparently while checking the wrong thing for parameters to estimate.

      Attached is a version of 'nls' that fixes these problems and otherwise produces the same answers (per 'all.equal') as the current official version on the "Examples" in 'nls' help file. What do you think about adopting these code changes and adding these examples to the 'nls' help file (see the attached 'nls-new examples.R')?

       I am using this in preparing for a presentation on the 'fda' package at UseR! in Ames (Aug. 8-10). I am temporarily adding the attached version of 'nls' to the 'fda' package with a help file primarily consisting of a pointer to help(nls, package=stats) and the attached examples. However, I'm hoping that someone will find the time to review this and either disabuse me of some of my deficiencies in this area or announce a plan to have something like this incorporated in a future version of R.

      Thanks,
      Best Wishes,
      Spencer Graves

p.s. How can I find the official source code for this? I tried "https://svn.r-project.org/R-packages/trunk" and found 'nlme' and other packages but not core R like 'stats'.
#####################
##
## SLIGHT MODIFICATIONS TO CURRENT EXAMPLES
## THAT DON'T WORK (except with my modifications)
##
# Weighted Michaelis-Menten model
# with data = a list that can not be coerced to a data.frame
TreatIrreg <- with(Treated, list(conc1=conc[1], conc.1=conc[-1], rate=rate))

# Passing arguments using a list that can not be coerced to a data.frame
weighted.MM1 <- function(resp, conc1, conc.1, Vm, K){   conc <- c(conc1, conc.1)
  resp <- TI$rate
#

  pred <- (Vm * conc)/(K + conc)
  (resp - pred) / sqrt(pred)
}
Pur.wt1 <- nls( ~ weighted.MM1(rate, conc1, conc.1, Vm, K), data=TreatIrreg,

               start = list(Vm = 200, K = 0.1), trace = TRUE)
##
## ERROR MESSAGE:
## Error in data.frame(conc1 = 0.02, conc.1 = c(0.02, 0.06, 0.06, 0.11,
0.11, :
## arguments imply differing number of rows: 1, 11, 12
##

# Passing arguments via 'get'

weignted.MM0 <- function(Vm, K){
  TI <- get("TreatIrreg")
  conc <- with(TI, c(conc1, conc.1))
  resp <- TI$rate
#

  pred <- (Vm * conc)/(K + conc)
  (resp - pred) / sqrt(pred)
}
Pur.wt0 <- nls( ~ weighted.MM0(Vm, K), start = list(Vm = 200, K = 0.1),

                   trace = TRUE)

##
## ERROR MESSAGE:
## Error in nls(~weighted.MM0(Vm, K), start = list(Vm = 200, K = 0.1),
trace = TRUE) :
## no parameters to fit
##

############################

nls <-
function (formula, data = parent.frame(), start, control = nls.control(),

    algorithm = c("default", "plinear", "port"), trace = FALSE,     subset, weights, na.action, model = FALSE, lower = -Inf,     upper = Inf, ...)
{
##
## 1. Basic set up
##

    formula <- as.formula(formula)
    algorithm <- match.arg(algorithm)
    if (!is.list(data) && !is.environment(data))

        stop("'data' must be a list or an environment")     mf <- match.call()
##
## 2. Names of variables in model and parameters to be estimated
##
# varNames = names of variables in the model

    varNames <- all.vars(formula)
    mWeights <- missing(weights)
    if (length(formula) == 2) {

        formula[[3]] <- formula[[2]]
        formula[[2]] <- 0
    }

# pnames = names of parameters to be estimated
pnames <- { if (missing(start)) { if (!is.null(attr(data, "parameters"))) { names(attr(data, "parameters")) } else { cll <- formula[[length(formula)]] func <- get(as.character(cll[[1]])) if (!is.null(pn <- attr(func, "pnames"))) as.character(as.list(match.call(func, call = cll))[-1][pn]) } } else names(start)

    }
# Confirm we have parameters to estimate
# (This is intended to correct an apparent error below;
# s. Graves 2007.06.02)

    if (!length(pnames))

        stop("no parameters to fit")
# Delete parameters from varNames

    env <- environment(formula)
    if (is.null(env))

        env <- parent.frame()
    if (length(pnames))

        varNames <- varNames[is.na(match(varNames, pnames))]
# if (!length(varNames))stop("no parameters to fit")
##
## 3. Get varNames
##
# 3.1. get lengths of varNames

    lenVar <- function(var) tryCatch(length(eval(as.name(var),

        data, env)), error = function(e) -1)     n <- sapply(varNames, lenVar)
# 3.2. any not found?

    if (any(not.there <- n == -1)) {

        nnn <- names(n[not.there])
        if (missing(start)) {
            if (algorithm == "plinear") 
                stop("No starting values specified")
            warning("No starting values specified for some parameters.\n", 
                "Intializing ", paste(sQuote(nnn), collapse = ", "), 
                " to '1.'.\n",
                "Consider specifying 'start' or using a selfStart model")
            start <- as.list(rep(1, length(nnn)))
            names(start) <- nnn
            varNames <- varNames[i <- is.na(match(varNames, nnn))]

# The following line 'n <- n[i]' seems not to be used before
# being overwritten a few lines below with 'n <- nrow(mf)'
n <- n[i] } else stop("parameters without starting value in 'data': ", paste(nnn, collapse = ", "))
    }
    respLength <- length(eval(formula[[2]], data, env))     {
      if(length(n)<1){
# 3.3. no varNames at all
# proposed addition 2007.06.02 by S. Graves
        varIndex <- logical(0)
        mf <- list(0)
        wts <- numeric(0)
      }
      else{
        varIndex <- n%%respLength == 0
        if(diff(range(n))>0){

# 3.4. 'data' is a list that can not be coerced to a data.frame
# proposed addition 2007.06.02 by S. Graves
mf <- data if (missing(start)) start <- getInitial(formula, mf)
# Store 'start' for initial evaluation of formula
startEnv <- new.env(parent = environment(formula)) for (i in names(start)) assign(i, start[[i]], envir = startEnv) rhs <- eval(formula[[3]], data, startEnv) n <- length(rhs) } else{
# 3.5. Traditional default: 'data' can be coerced to a data.frame
# The following generates a call to 'eval.parent'
# that generates an Error
# if 'data' is a list that can NOT be coerced to a data.frame
mf$formula <- as.formula(paste("~", paste(varNames[varIndex], collapse = "+")), env = environment(formula)) mf$start <- mf$control <- mf$algorithm <- mf$trace <- mf$model <- NULL mf$lower <- mf$upper <- NULL mf[[1]] <- as.name("model.frame") mf <- eval.parent(mf) n <- nrow(mf) mf <- as.list(mf) } wts <- { if (!mWeights) model.weights(mf) else rep(1, n) } if (any(wts < 0 | is.na(wts))) stop("missing or negative weights not allowed") }

    }
##
## 4. Set up iteration
##

    if (missing(start))

        start <- getInitial(formula, mf)     for (var in varNames[!varIndex]) mf[[var]] <- eval(as.name(var),

        data, env)
    m <- switch(algorithm,

                plinear = stats:::nlsModel.plinear(formula, mf, start, wts),
                port = stats:::nlsModel(formula, mf, start, wts, upper),
                stats:::nlsModel(formula, mf, start, wts))
    ctrl <- nls.control()
    if (!missing(control)) {
        control <- as.list(control)
        ctrl[names(control)] <- control

    }
##
## 5. Iterate
##

    if(algorithm != "port") {

        if (!missing(lower) || !missing(upper)) 
            warning("Upper or lower bounds ignored ",
                    "unless algorithm = \"port\"")
        convInfo <- .Call(stats:::R_nls_iter, m, ctrl, trace)
        nls.out <- list(m = m, convInfo = convInfo, data = substitute(data), 
            call = match.call())

    }
    else {
        iv <- stats:::nls_port_fit(m, start, lower, upper, control, trace)
        nls.out <- list(m = m, data = substitute(data), call = match.call())
        nls.out$convergence <- as.integer(if (iv[1] %in% 3:6) 
            0
        else 1)
        nls.out$message <- switch(as.character(iv[1]),
            "3" = "X-convergence (3)", "4" = "relative convergence (4)",
            "5" = "both X-convergence and relative convergence (5)", 
            "6" = "absolute function convergence (6)",
            "7" = "singular convergence (7)", 
            "8" = "false convergence (8)",
            "9" = "function evaluation limit reached without convergence (9)", 
            "10" = "iteration limit reached without convergence (9)", 
            "14" = "storage has been allocated (?) (14)",
            "15" = "LIV too small (15)", "16" = "LV too small (16)",
            "63" = "fn cannot be computed at initial par (63)", 
            "65" = "gr cannot be computed at initial par (65)", 
            "300" = "initial par violates constraints")
        if (is.null(nls.out$message)) 
            nls.out$message <- paste("See PORT documentation.\tCode (", 
                iv[1], ")", sep = "")
        if (nls.out$convergence) {
            msg <- paste("Convergence failure:", nls.out$message)
            if (ctrl$warnOnly) {
                warning(msg)
            }
            else stop(msg)
        }
        nls.out$call$lower <- lower
        nls.out$call$upper <- upper

    }
##
## 6. Done
##
    nls.out$call$algorithm <- algorithm
    nls.out$call$control <- ctrl
    nls.out$call$trace <- trace
    nls.out$na.action <- attr(mf, "na.action")
    nls.out$dataClasses <- attr(attr(mf, "terms"), "dataClasses")
    if (model) 
        nls.out$model <- mf
    if (!mWeights) 
        nls.out$weights <- wts

    nls.out$control <- control
    class(nls.out) <- "nls"
    nls.out
}

# Weighted Michaelis-Menten model
# with data = a list that can not be coerced to a data.frame
Treated <- Puromycin[Puromycin$state == "treated", ] TreatIrreg <- with(Treated,

        list(conc1=conc[1], conc.1=conc[-1], rate=rate)) sapply(TreatIrreg, length)

# Passing arguments via 'get'

weighted.MM0 <- function(Vm, K){
  TI <- get("TreatIrreg")
  conc <- with(TI, c(conc1, conc.1))
  resp <- TI$rate
#

  pred <- (Vm * conc)/(K + conc)
  (resp - pred) / sqrt(pred)
}
gc()
start.time <- proc.time()
Pur.wt0 <- nls( ~ weighted.MM0(Vm, K), start = list(Vm = 200, K = 0.1),

                   trace = TRUE)

(et0 <- proc.time()-start.time)
# NOTE: 'get' works but is not recommended,
# because it's too subtle and inflexible.

# Passing arguments using a list that can not be coerced to a data.frame
weighted.MM1 <- function(resp, conc1, conc.1, Vm, K){   conc <- c(conc1, conc.1)
#

  pred <- (Vm * conc)/(K + conc)
  (resp - pred) / sqrt(pred)
}
gc()
start.time <- proc.time()
Pur.wt1 <- nls( ~ weighted.MM1(rate, conc1, conc.1, Vm, K),

     data=TreatIrreg, start = list(Vm = 200, K = 0.1), trace = TRUE) (et1 <- proc.time()-start.time)

gc()
start.time <- proc.time()
Pur.wt.1port <- nls( ~ weighted.MM1(rate, conc1, conc.1, Vm, K),

                    data=TreatIrreg, start = list(Vm = 200, K = 0.1),
                    algorithm="port", trace = TRUE)
(et.1port <- proc.time()-start.time)

# Chambers and Hastie (1992) Statistical Models in S (Wadsworth, p. 537)
# say, "If the value of the right side [of formula] has an attribute
# called 'gradient' this should be a matrix with the number of rows
# equal to the length of the response and one column for each
# parameter.

weighted.MM.gradient <- function(resp, conc1, conc.1, Vm, K){   conc <- c(conc1, conc.1)
#

  K.conc <- (K+conc)
  dy.dV <- conc/K.conc
  dy.dK <- (-Vm*dy.dV/K.conc)
  pred <- Vm*dy.dV
# Avoid NAs

  pred.5 <- sqrt(pmax(0, pred))
  dev <- (resp - pred) / pred.5
  Ddev <- (-0.5*(resp+pred)/(pred.5*pred))   attr(dev, "gradient") <- (Ddev * cbind(Vm = dy.dV, K = dy.dK))   dev
}
gc()
start.time <- proc.time()
Pur.wt.grad <- nls(

   ~ weighted.MM.gradient(rate, conc1, conc.1, Vm, K), data=TreatIrreg,    start = list(Vm = 200, K = 0.1), trace = TRUE) (et.grad <- proc.time()-start.time)



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Mon 04 Jun 2007 - 16:30:48 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 04 Jun 2007 - 18:33:30 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.