Thanks for responding!
I was originally constructing the call and then evaluating it because in an
earlier version of this problem, I needed to specify different arguments,
depending upon how the function was called; this is no longer a problem, so
I've just changed to calling gnls() directly. Unfortunately, that does
not solve my problem. Below I've copied a full script that causes the
error, along with the results of a sample run, with traceback. When I run
this (on rw1023, Windows 98, with nlme version 3.1-10 dated 2001/01/10.
First, the results:
> source("C:/home/tmp/Rfixes/test3.R")
Error in eval(expr, envir, enclos) : couldn't find function
"Cexp"> traceback()
9: eval(expr, envir, enclos)
8: eval(model, data.frame(data, getParsGnls(plist, pmap, beta, N)))
7: eval(expr, envir, enclos)
6: eval(modelResid[[2]], envir = nlEnv)
5: gnls(Model, data = dta, weights = wtFunc, start = start)
4: BMDS.Cexp(T4 ~ Dose, data = testdata, start = c(A = 3.9365902,
m = 0.4666282), fixed = c(B = 0, g = 1))
3: eval.with.vis(expr, envir, enclos)
2: eval.with.vis(ei, envir)
1: source("C:/home/tmp/Rfixes/test3.R")
This is test3.R:
if (!require(nlme)) stop("nlme must be installed to try this")
CexpB <-function (dose, A, B, m, g) {
### exp(B) + (exp(A)-exp(B))*exp(-(exp(m)*dose)^(1 + g^2))
mCall <- match.call()
fixed <- attr(eval(mCall[[1]]),"fixed")
.expr1 <- if (!("B" %in% names(fixed))) exp(B) else
fixed["B"]
.expr2 <- if (!("A" %in% names(fixed))) exp(A) else
fixed["A"]
.expr5 <- (if (!("m" %in% names(fixed))) exp(m) else
fixed["m"]) * dose
.expr7 <- if (!("g" %in% names(fixed))) 1 + g^2 else
fixed["g"]
if ("g" %in% names(fixed)) g <- sqrt(fixed["g"] - 1)
.expr3 <- .expr2 - .expr1
.expr8 <- .expr5^.expr7
.expr10 <- exp(-.expr8)
.value <- .expr1 + .expr3 * .expr10
.grad <- array(0, c(length(.value), 4), list(NULL, c("A",
"B",
"m", "g")))
.grad[, "A"] <- .expr2 * .expr10
.grad[, "B"] <- .expr1 - .expr1 * .expr10
.grad[, "m"] <- ifelse(dose > 0,
-.expr3 * (.expr10 * (.expr5^(.expr7 - 1) *
(.expr7 * .expr5))),
0)
.grad[, "g"] <- ifelse(dose > 0,
-.expr3 * (.expr10 * (.expr8 * (log(.expr5) *
(2 * g)))),
0)
if (!is.null(fixed))
.grad <-
.grad[,-match(names(fixed),c("A","B","m","g"))]
attr(.value, "gradient") <- .grad
.value
}
BMDS.Cexp <- function(Model,data=sys.frame(sys.parent()),
start=NULL,fixed=NULL) {
DoseName <- deparse(Model[[3]])
Dose <- eval(Model[[3]],data)
MName <- deparse(Model[[2]])
M <- eval(Model[[2]],data)
dta <- data.frame(Dose,M)
names(dta) <- c(DoseName,MName)
idx <- order(dta[,DoseName])
dta <- dta[idx,]
wtFunc <- varPower()
### Reconstruct the arguments for the model function
Cexp <- CexpB
formals(Cexp) <- formals(Cexp)[-match(names(fixed),names(formals(Cexp)))]
attr(Cexp,"fixed") <- fixed
### Construct the expression to fit the model
newforms <- names(formals(Cexp))[-1]
Model <- as.formula(paste(MName," ~ Cexp(", DoseName,",
",
paste(newforms,collapse=", "),
")",
sep=""))
### Fit the model
gnls(Model,data=dta, weights=wtFunc, start=start)
}
testdata <- data.frame(Dose=rep(c(0,0.003,0.03,0.3),c(14,12,6,6)),
T4=c(51.61, 54.37, 51.32, 62.17, 30.54, 43.56, 53.56,
53.01, 51.23, 44.47, 54.63, 63.14, 43.73, 60.07,
55.73, 61.13, 41.79, 48.91, 40.50, 50.29, 49.45,
52.50, 44.59, 37.31, 45.98, 46.37, 44.15, 46.77,
48.45, 54.64, 51.59, 41.61, 31.63, 29.85, 26.34,
31.72, 40.54, 30.48))
out <- BMDS.Cexp(T4 ~ Dose, data=testdata, start=c(A=3.9365902, m=0.4666282),
fixed=c(B=0,g=1))
R. Woodrow Setzer, Jr. Phone:
(919) 541-0128
Experimental Toxicology Division Fax: (919) 541-5394
Pharmacokinetics Branch
NHEERL MD-74; US EPA; RTP, NC 27711
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._