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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
