[Rd] weighted nonlinear fits: `nls' and `eval'

Joerg van den Hoff j.van_den_hoff at fzd.de
Wed Apr 30 14:59:55 CEST 2008


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



More information about the R-devel mailing list