2 days ago I asked this on r-help, but no luck... since
this is actually a programming question, I post it here
again:
my question concerns the use of `eval' in defining the model
formula for `nls' when performing weighted fits. (I use
version 2.6.2., but according to NEWS there were no changes
to `nls' in 2.7.0, so the problem is still present). in this
scenario their exists a serious problem.
consider the following simple example, where some fixed
model plus data are used to perform both unweighted and
weighted fits. (I intentionally used very uneven weights to
guarantee large differences in the results.)
#================================CUT==========================ln <- 100
x <- 1:ln
y0 <- exp(-.02*x)
y <- rnorm(y0, y0, .01)
wts <- rep(1, ln)
y[30] <- 1.2*y[30]
wts[30] <- 1e3
model <- y ~ exp(-k*x)
xydata <- list(x=x, y=y)
#simple unweighted fit works as expected:
r0 <- nls(model, start = c(k=.01), data = list(x=x, y=y))
#simple weighted fit works as expected:
r1 <- nls(model, start = c(k=.01), data = xydata, weights = wts)
#this actually performs an unweighted fit (issuing a warning):
mdl <- model[[3]]
r2 <- nls(y ~ eval(mdl), start = c(k=.01), data = xydata, weights = wts)
#this, too, actually performs an unweighted fit (issuing a warning):
dv1 <- deriv(model, "k")
r3 <- nls(y ~ eval(dv1), start = c(k=.01), data = xydata, weights = wts)
#weighted fit, works as expected
dv2 <- deriv(model, "k", c("k", "x"))
r4 <- nls(y ~ dv2(k, x), start = c(k=.01), data = xydata, weights = wts)
#================================CUT==========================
if you copy this to the R prompt you can see the fits
producing `r2' and `r3' do _not_ do what's intended. looking
at the source code of `nls' I get some ideas where it is
going wrong (in setting up the model frame in `mf') what
exactly is the problem.
moreover, I found out that `r2' and `r3' can be convinced
to do what's intended by modifying `xydata' to
xydata <- list(x=x, y=y, `(weights)` = wts)
i.e. by adding a `(weights)` component to the `data' field.
but that probably is not the way to go ...
while it is clear that my example is artificial, there _are_
circumstances where `eval' constructs for defining the model
formula occur quite naturally when writing some
wrappers/drivers for `nls' usage.
questions:
=========
1.) is it actually not intended/supported to use `eval' in
defining the model formula argument of `nls'? if so, why
not (it's a legitimate and not that esoteric language
feature)? or what am I missing (e.g. necessity to provide
explicitely an `envir' argument to `eval'. if so, which
one)?
2.) could/should not `nls' be modified to cleanly support
the use of `eval' (remarks along the line of "submit
changes" would not help much, I'm afraid :-))? apart from
the problem described here there are others where `nls'
simply parses `eval' based formulas not as intended
(erroneously) and bails out.
3.) if the policy of the core group is that the current
state of affairs is quite perfect (an assessement which I
would'nt share), should not at least the manpage of `nls'
be augmented to warn people very clearly that any use of
`eval' in the model formula is calling for trouble.
I find the example above especially annoying: the fits
producing `r2' and `r3' run apparently fine, only issuing a
warning, which sure is cryptic for joe user. so, if this
happens in some wrapper function provided to others (which I
do) they very well might believe having performed a weighted
fit when they did not.
thanks,
joerg
