# [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, 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[] <- formula[]
formula[] <- 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[]))
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)

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[], 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[], 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[] <- 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 %in% 3:6)
0
else 1)
nls.out\$message <- switch(as.character(iv),
"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, ")", 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, 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()
{"ft_lang_label":"__label__en","ft_lang_prob":0.5244918,"math_prob":0.9895504,"size":10008,"snap":"2020-24-2020-29","text_gpt3_token_len":2999,"char_repetition_ratio":0.1057577,"word_repetition_ratio":0.17225201,"special_character_ratio":0.34372503,"punctuation_ratio":0.20722057,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9972176,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2020-07-16T01:57:45Z\",\"WARC-Record-ID\":\"<urn:uuid:a0128c55-f163-4392-aaf8-381862d7eddd>\",\"Content-Length\":\"23937\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:18b9f39e-ac98-4a8c-b2e0-a2f55d630799>\",\"WARC-Concurrent-To\":\"<urn:uuid:fda6ee1b-366e-43ae-9e30-367d898aa120>\",\"WARC-IP-Address\":\"99.84.191.123\",\"WARC-Target-URI\":\"https://tolstoy.newcastle.edu.au/~rking/R/e2/devel/07/06/3504.html\",\"WARC-Payload-Digest\":\"sha1:TR7PKMIGP5HSRGTKYOVHXPCD4G2P3ENI\",\"WARC-Block-Digest\":\"sha1:ABQABEP7NPVSN2PJ54J3TVPAWI2QSMCN\",\"WARC-Identified-Payload-Type\":\"application/xhtml+xml\",\"warc_filename\":\"/cc_download/warc_2020/CC-MAIN-2020-29/CC-MAIN-2020-29_segments_1593657176116.96_warc_CC-MAIN-20200715230447-20200716020447-00227.warc.gz\"}"}