From bolker at zoo.ufl.edu  Thu Jun 10 20:10:06 2004
From: bolker at zoo.ufl.edu (Ben Bolker)
Date: Thu Jun 10 19:12:04 2004
Subject: [Rd] overhaul of mle
Message-ID: <40C8A3FE.6030102@zoo.ufl.edu>

   So, I've embarked on my threatened modifications to the mle subset
of the stats4 package.    Most of what I've done so far has *not* been
adding the slick formula interface, but rather making it work properly
and reasonably robustly with real mle problems -- especially ones
involving reasonably complex fixed and default parameter sets.
Some of what I've done breaks backward compatibility (see below), but
there are what I think are some reasonably important design issues -- it
would be nice to get them fixed now, while the mle library is still
in its infancy ...

   I've appended quite a long list of changes and "to do" stuff, which
is also available at http://www.zoo.ufl.edu/bolker/mle/mle-changes.txt
also there are
http://www.zoo.ufl.edu/bolker/mle/mle-diffs.txt  (diffs against 1.9.1) and
http://www.zoo.ufl.edu/bolker/mle/mle.R (R code).

   I'm sure some of the code could be improved, especially in issues 
relating to eval() -- I'm not very good at that stuff.

   I apologize for changing so much at once, I got kind of carried away.
I look forward to comments on the various bits & pieces ...

   Ben Bolker

------------

*** = changes default behavior in a user-visible way that I consider
better, but arguably not just a bug fix: I've tried to justify these

*** added checking code at the beginning of mle() that converts named
numeric vectors into lists [this seems harmless to me -- is there a
reason not to allow the user to specify a named vector rather than a
list?  especially since start gets sapply'd back to a vector before it
gets passed to optim() anyway?]

Added some (probably not quite right) code to mle to make sure that
fullcoef gets evaluated properly when (1) the list contains
expressions and (2) the list contains expressions that depend on other
arguments e.g. fixed = list(a=y,b=x[y])
I used

     fullcoef <- lapply(fullcoef,eval,envir=fullcoef,
                        enclos=parent.frame(100))

but there may well be a better way to do it.

*** Changed type of fullcoef from "numeric" to "list", and return
fullcoef rather than unlist(fullcoef) from mle [couldn't see a rationale 
for this -- it destroys a lot of the information in fullcoef *and* is a
pain, say, when the fixed arguments include a data frame with lots of
information in it]

*** Changed "coef" method to return object@coef, not object@fullcoef
[this really seems to be the better choice to me -- I normally want to
see the *fitted values* of the MLE, not all the other auxiliary
stuff.  Besides, object@fullcoef can be very long, and therefore a
nuisance to see in the default show(object) method]

made a fullcoef accessor *function* to return object@fullcoef --
should this really be a method?

added a cor method for mle objects -- which just normalizes the
variance-covariance matrix to a correlation matrix.  Is this a bad
idea/perversion of the cor method?

changed variable "pi" to "p.i" throughout mle -- less confusing!

changed
call$fixed <- fix
to
call$fixed <- c(fix,eval(call$fixed))
for cases where there are non-trivial fixed arguments

added "follow" argument to profile: this makes profiling use a
continuation method where the starting point for each profile
optimization is the previous best-fit solution, rather than the
overall MLEs of the parameters.  Actually fairly easy to implement (I
think: I haven't really tested that it works on anything hard, just
that it doesn't seem to break profiling) -- requires pfit to be
assigned globally within onestep() and a few lines of code further
down.

added an AIC method for mle objects

collapsed the absVal/!absVal code cases slightly

added a "sqrVal" argument for those who want to see the value of the
log-likelihood, not the square root or signed square root (could be
collapsed into a "scale" argument for the profile plot = "sqrt",
"abssqrt", "lik")

added code and options to plot labels of confidence levels (logical
plot.confstr, character confstr)

added add= argument (to draw profiles on an existing plot)

added arguments for color of minimum values, confidence limits,
profile (col.minval, col.conf, col.prof)

added options for confidence interval: when applied to an mle object,
method "spline" does the previous behavior (profile and apply confint
to the result).  Method "quad" simply presents the quadratic
approximation to the confidence intervals.  Method "exact" uses
uniroot() to find the precise point where the profile crosses the
critical level in each direction.

Added mle.options() command, and .mle.options state variable, to keep
global options (method for optim() and method for confint()): I'm not
at all sure that this is the best way to implement options, this was
just my first crack at it

added a warning to show(mle) if optim() did not converge

Added code that allows (1) default arguments (evaluated
in the frame of the full coefficient list, with fixed values
and starting values substituted and (2) arguments specified in the
start list in arbitrary order (which seems like a reasonable expectation 
since
it *is* specified as a list).  The fundamental problem is that optim() 
loses names
of the parameter vector somewhere.
Example:

x = runif(200)
y = 1+x+x^2+rnorm(200,sd=0.05)
fn <- function(a,b,z=2,c,d) {
    -sum(dnorm(y,mean=a+c*x+d*x^2,sd=exp(b),log=TRUE))
}

m1 = mle(minuslogl=fn,start=list(a=1,b=1,c=1,d=1))
## fails with "missing argument" warning, about wrong argument
m1 = mle(minuslogl=fn,start=list(a=1,b=1,c=1,d=1),fixed=list(z=2))
## works
m2 = mle(minuslogl=fn,start=list(a=1,d=1,c=1,b=1),fixed=list(z=2))
## fails -- coeffs returned in wrong order


TO DO:

torture-test on some real problems!

better documentation?  e.g. ?profile.mle-class doesn't give details on
arguments -- have to look at profile.nls

(allow "which" to be a character vector -- match names)?

HARDER:
fancy formula interface [cf. svymle in survey package] e.g.
mll <-
mLL(type="independent",distrib="normal",resp=y,mean=~a+b*x,sd=~s,
       param=~a+b+s)

allow for fitting of transformed parameters (exp/log, tanh/atanh =
logistic/logit)

2D profiles (quadratic or thin-plate spline???)

EASIER but breaking backward compatibility:
merge absVal/sqrVal into a "scale" argument?

EASIER but ??worthwhile??:
allow spline to be turned off when plotting profiles?
   (method "spline"/"raw")?

code for producing/plotting "slices" (non-optimized transects through
parameter space); other diagnostic tools?

NOT SURE:
change show, show(summary) methods to bring them more in line with
other classes?

add test to confint(profile) that warns if method is supplied?


[Rd] overhaul of mle
Ben Bolker