Skip to content

R-beta: SEs for one-param MLE in R?

3 messages · Martin Maechler, Douglas Bates, Peter Dalgaard

Original
Martin Maechler
# 
[--- diverting from  R-help to R-devel ... MM ---]
Jim> Just back from a week camping in the snow in the English Lakes,
    Jim> and trying to catch up...

I hope the adventurous vacation was also relaxing in some ways..

    MM> ....
    MM> BTW: I am (we are) interested in the functions that you are writing
    MM> for nlm(.)  It certainly is worthwhile to have nlm(.) return a class
    MM> "nlm" result and provide print.nlm(.) and summary.nlm(.) functions
    MM> {{ Jim Lindsey already posted something like this, unfortunately
    MM> using "nls" which we don't want as long as it is not very close to
    MM> S' nls(.) function }}

    Jim> I am afraid that I don't understand the logic of this requirement
    Jim> of closeness to R for such functions. nlm() itself is not close
    Jim> and even hist() has never been very similar.

Yes, there are some cases where we (the R-core team, or the "R-devel" group)
have decided that S-plus is so wrong that we don't want to emulate or stay
close to it.   hist(.) is one such example.

With  nls(.), this is quite different, I think.
nls(.) has several nice features, (your "nls" did too!).
Most notably, the calling syntax of nls(.) using model notation,
is something I would want before I'd call a function "nls".

At the moment, I think it would be better to add an "nlm" class and
corresponding print and summary methods to  nlm(.), rather than calling
such a thing "nls".

Also, a few weeks ago, Ross said that he planned to add code to nlm
which would make use of  gradient and hessian function when provided
(or also, using D(.) ?).

    Jim> On the other hand,
    Jim> remember that the nls() I sent was a very cut-down version of one
    Jim> in one of my libraries. It had the above solution for the
    Jim> inversion with one dimensional parameters. By the way, the
    Jim> original Fortran for nlm that I ported to R printed out a warning
    Jim> that the algorithm is very inefficient for one-dimensional
    Jim> problems.
Yes, your "nls" function was useful!
and I agree that "nlm" should be pushed in the direction of what you had.

    Jim>   With respect to Bill's negative binomial that started another
    Jim> discussion, one of the functions in my nonlinear regression
    Jim> library does negative binomial nonlinear regression for both the
    Jim> mean and the dispersion parameters, along with twenty odd other
    Jim> distributions. I used it for beta-binomial regression in my paper
    Jim> with Pat Altham in the latest Applied Statistics. (Also another
    Jim> similar function for a finite mixture with these distributions,
    Jim> for example for negative binomial with inflated zeroes.) In my
    Jim> repeated measures library, there is a similar function for
    Jim> regressions with the same collection of distributions, but having
    Jim> a random intercept. Once R stabilizes...

well, isn't  R  stabilizing more and more?
Since 0.61, at least the way extension packages (formerly known as "libraries")
should be written is pretty stable.
Jim, I think several people really would be interested in your 
``nonlinear regression library'' and the ``repeated measures library''.

If they don't work in R 0.61, they could be at least put into 
CRAN/src/contrib/devel/.

Best regards!
Martin
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel 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-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Douglas Bates
# 
Martin Maechler <maechler@stat.math.ethz.ch> writes:
I have to admit to some bias here being one of the people who wrote
nls for S.  I find a couple of other features for nls to be very
helpful.  The "plinear" algorithm that profiles over any linear
parameters in the model can be very effective.  It doesn't always
converge in fewer iterations but it starts at values of the parameters
that are much closer to the optimum and it requires fewer starting
estimates to be calculated.

I have just finished preparing a short course on nonlinear regression
to be presented this weekend.  The slides are available at
	     http://franz.stat.wisc.edu/pub/NLME/Kansas/
Section 4 shows the use of some of the current and new facilties in
S for nonlinear regression.  

I refer in those slides to the nonlinear regression test cases at the
U.S. National Institutes of Standards and Technology (NIST) web site.
It is interesting to compare how numerical analysts view nonlinear
least squares with the approach of statisticians.  The NIST folks make
a big deal about using 128 bit arithmetic to produce "certified"
values of the parameter estimates.  In their favor they also produce
the approximate standard errors for the estimates but they don't seem
to notice that the estimates are things like 2.3 +/- 0.7 and they are
worrying about whether the 8th decimal place should be a 6 or a 7.
Also their idea of a "difficult" problem is to take a relatively easy
model fit and multiply all the starting estimates by 100.  In one
case, for the BoxBOD data where there is an increase to an asymptote,
all the responses are between 100 and 225 so they start the asymptote
at 1 and find it is difficult to converge.  With the "plinear"
algorithm you don't even need a starting estimate for the asymptote.

Those slides show the new approach to selfStart models in S.  We have
found that with a minimal effort it is possible to get convergence on
any of the examples on the NIST site simply by considering what the
parameters represent and by exploiting conditional (or "partial")
linearity.

For R it is still not clear to me how to manipulate the expressions
and functions.  I saw some of the discussion about function closures
and why it is not easy to coerce to a function but that is still black 
magic to me.  I regret that I am not likely to be able to take the
time to understand it better in the near future.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel 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-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Peter Dalgaard
# 
Douglas Bates <bates@stat.wisc.edu> writes:
You could always take Bill V's approach of blatantly claiming "this
can't be done in R" and watching the reaction... ;)

<tutorial audience=anyone>

Expressions are fairly easy: Just mentally convert them to LISP to
access the parse tree. Convert everything to function calls and
represent each as a list with the first element being the function
name and the others being it's arguments. (This is not particularly
useful to people who don't know Lisp, I know.)

This little function gets you some of the way, although it won't do
things like tagged arguments (it can certainly be improved in other
ways too):
function (e) 
{
        level <<- level + 1
        cat(substring("||||||||||||||||", 1, level))
        cat(deparse(e), "\n")
        if (is.call(e) || is.expression(e)) 
                for (i in 1:length(e)) dissect(e[[i]])
        level <<- level - 1
}
1 + 2 + 3 * 4 
|1 + 2 + 3 * 4 
||+ 
||1 + 2 
|||+ 
|||1 
|||2 
||3 * 4 
|||* 
|||3 
|||4 

e.g. like LISP's (+(+ 1 2)(* 3 4)). 

This is all pretty much like S. What was the problem in the first place?

</tutorial>

Regarding function coercion, I'm still unconvinced why this shouldn't
be possible. I'd kind of like to keep the "everything is a list"
paradigm. Obviously, one cannot keep perfect S compatibility (because
of the environment/closure issue), but some way of constructing a
function from it's arguments, body, and environment would be useful.

Many things are in fact possible if you are sneaky enough, e.g. here's
a list.to.envir conversion:

environment(eval(expression(function(){}),envir=list(a=1,b=2))))