[R] weighted nonlinear fits: `nls' and `eval'
Joerg van den Hoff
j.van_den_hoff at fzd.de
Mon Apr 28 14:56:52 CEST 2008
my question concerns the use of `eval' in defining the model formula
for `nls' (version 2.6.2.).
consider the following simple example, where the same model and data
are used to perform unweighted and weighted fits. I intentionally
used very uneven weights to guarantee large differences in the results
ln <- 100
x <- 1:ln
y0 <- exp(-.02*x)
y <- rnorm(y0, y0, .01)
wts <- rep(1, ln)
y <- 1.2*y
wts <- 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[]
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)
as 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'?) but not what exactly is the problem
`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 `(weights)` component here, too. but that probably is not
the way to go ...
while it is clear in the case of `r2' that my example is artificial, `r3' is what
I have used up to now for unweighted fits without problems. moreover there _are_
circumstances where `eval' constructs for defining the model formula occur quite
is it actually not intended/supported to use `eval' in defining the
model formula argument of `nls'? if so, why not? or what am I missing (e.g.
necessity to provide explicitely an `envir' argument to `eval': which one?)?
could/should not `nls' be modified to cleanly support the use of `eval'
(remarks along the line of "submit changes" would not help much :-)).
if the policy of the core group is that the current state of affairs is
quite perfect, should not 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: `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.
More information about the R-help