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?
Ben Bolker <bolker@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@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]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@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 argumentsWhich 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 orderHmm.. 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@biostat.ku.dk) FAX: (+45) 35327907