[Rd] overhaul of mle

[Peter Dalgaard]

Ben Bolker <bolker at zoo.ufl.edu> writes:

> *** 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]

Wait a minute. How can a likelihood function have an argument that is
a data frame? I think you're abusing the fixed arguments if you use it
to pass in data. The natural paradigm for that would be to pass data
via a closure, i.e. 

mll <- with(data,
    function(lambda=1,theta=0)sum(dpois(y, lambda+theta*x, log=TRUE))
)

 
> *** Changed "coef" method to return object at coef, not object at 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 at fullcoef can be very long, and therefore a
> nuisance to see in the default show(object) method]

See above. This was never intended to contain auxiliary stuff (and
AFAIR has already been changed once in the opposite direction, by Brian)
 
> made a fullcoef accessor *function* to return object at 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?

Yes, I think so. cov2cor(vcov(ml.obj)) is easy enough.

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

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

Which there shouldn't be...
 
> 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.

Sounds nice, but surely you don't need a global assignment there? A
superassign ("<<-") perhaps, but that doesn't need to go to
.GlobalEnv. 


> 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.

All fine. The last one could be important as I had a case where the
spline method went rather badly wrong (the data it happened with are
still rather heavily embargoed I'm afraid). 

> 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

Hmm.. I see the effect with the current version too. Depending on
temperament, it is the labels rather than the order that is wrong...  

> 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)

The last one should be trivial, no?

mll2 <- function(a,b,c,d) mll1(log(a),atan(b),c,d)
 
Also: code for combination of likelihoods (i.e. summing likelihoods
for independent subexperiments involving the same parameters;
integrating out nuisance variables.)

MUCH, MUCH HARDER: Figure out what it would take to include higher
order asymptotic inference (as per Brazzale/Bellio's hoa bundle) in a
generic likelihood setting...



-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907