Skip to content

suggested modification to the 'mle' documentation?

17 messages · Spencer Graves, Duncan Murdoch, Gabor Grothendieck +4 more

Spencer Graves
# 
Hello: 

      I wish to again express my appreciation to all who have 
contributed to making R what it is today. 

      At this moment, I'm particularly grateful for whoever modified the 
'mle' code so data no longer need be passed via global variables.  I 
remember struggling with this a couple of years ago, and I only today 
discovered that it is no longer the case. 

      I'd like to suggest that the 'mle' help file be modified to 
advertise this fact, e.g., by adding one of the two examples appearing 
below. 

      Best Wishes,
      Spencer Graves
################################
x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
#  Pass data via function arguments rather than global variables
ll.5 <- function(ymax=15, xhalf=6, x., y.)
         -sum(stats::dpois(y., lambda=ymax/(1+x./xhalf), log=TRUE))
(fit.5 <- mle(ll.5, start=list(ymax=15, xhalf=6),
              fixed=list(x.=x, y.=y)))

ll3 <- function(lymax=log(15), lxhalf=log(6), x., y.)
  -sum(stats::dpois(y.,
         lambda=exp(lymax)/(1+x./exp(lxhalf)), log=TRUE))
(fit3 <- mle(ll3, start=list(lymax=0, lxhalf=0),
             fixed=list(x.=x, y.=y)))
Peter Dalgaard
#
Spencer Graves wrote:
In a word: No!!! That is not the design. A likelihood function is a 
function of its parameters, and the "fixed" argument is for holding some 
parameters fixed (e.g. during profiling).

To include data, just make a closure, e.g.

poissonLike <- function(x., y.){
    function(ymax=15, xhalf=6)
      -sum(stats::dpois(y., lambda=ymax/(1+x./xhalf), log=TRUE))}
mll <-  poissonLike(x, y)
mle(ll, ....

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Gabor Grothendieck
    

    

      
      #
    

  


  

      
The closure only works if you are defining the inner function yourself.
If you are not then its yet more work to redefine the environment of
the inner function or other workaround.

        
On Dec 6, 2007 6:01 PM, Peter Dalgaard <p.dalgaard at biostat.ku.dk> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
Gabor Grothendieck wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        
At this point I'd just like to advertise the "bbmle" package
(on CRAN) for those who respectfully disagree, as I do, with Peter over
this issue.  I have added a data= argument to my version
of the function that allows other variables to be passed
to the objective function.  It seems to me that this is perfectly
in line with the way that other modeling functions in R
behave.

  (My version also has a cool formula interface and other
bells and whistles, and I would love to get feedback from other
useRs about it.)

   cheers
    Ben Bolker

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

      
Ben Bolker wrote:

            

      
      
      
    

        
This is at least cleaner than abusing the "fixed" argument. As you know,
I have reservations, one of which is that it is not a given that I want
it to behave just like other modeling functions, e.g. a likelihood
function might refer to more than one data set, and/or data that are not
structured in the traditional data frame format. The design needs more
thought than just adding arguments.

I still prefer a design based a plain likelihood function. Then we can
discuss how to construct such a function so that  the data are
incorporated in a flexible way.  There are many ways to do this, I've
shown one, here's another:

            

      
      
      
    

        
Call:
mle(minuslogl = f, start = list(lambda = 10))

Coefficients:
 lambda
12.3402

It is not at all an unlikely design to have mle() as a generic function
which works on many kinds of objects, the default method being
function(object,...) mle(minuslogl(obj)) and minuslogl is an extractor
function returning (tada!) the negative log likelihood function.

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Gabor Grothendieck
    

    

      
      #
    

  


  

      
On Dec 7, 2007 8:10 AM, Peter Dalgaard <P.Dalgaard at biostat.ku.dk> wrote:

            

      
      
      
    

        
The explicit environment manipulation is what I was referring to but
we can simplify it using proto.  Create a proto object to hold
f and x then pass the f in the proto object (rather than the
original f) to mle.  That works because proto automatically resets
the environment of f when its added to avoiding the evalq.

            

      
      
      
    

        
Call:
mle(minuslogl = p[["f"]], start = list(lambda = 10))

Coefficients:
  lambda
12.46000

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Duncan Murdoch
    

    

      
      #
    

  


  

      
On 12/7/2007 8:10 AM, Peter Dalgaard wrote:

            

      
      
      
    

        
We should allow more general things to be passed as data arguments in 
cases where it makes sense.  For example a list with names or an 
environment would be a reasonable way to pass data that doesn't fit into 
a data frame.

            

      
      
      
    

        
We really need to expand as.environment, so that it can convert data 
frames into environments.  You should be able to say:

environment(f) <- as.environment(d)

and get the same result as

environment(f)<-evalq(environment(),d)

But I'd prefer to avoid the necessity for users to manipulate the 
environment of a function.  I think the pattern

model( f, data=d )

being implemented internally as

environment(f) <- as.environment(d, parent = environment(f))

is very nice and general.  It makes things like cross-validation, 
bootstrapping, etc. conceptually cleaner:  keep the same 
formula/function f, but manipulate the data and see what happens.
It does have problems when d is an environment that already has a 
parent, but I think a reasonable meaning in that case would be to copy 
its contents into a new environment with the new parent set.

Duncan Murdoch

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Gabor Grothendieck
    

    

      
      #
    

  


  

      
On Dec 7, 2007 8:43 AM, Duncan Murdoch <murdoch at stats.uwo.ca> wrote:

            

      
      
      
    

        
Something close to that is already possible in proto and its cleaner in proto
since the explicit environment manipulation is unnecessary as it occurs
implicitly:

1. In terms of data frame d from Peter Dalgaard's post the code
below is similar to my last post but it replaces the explicit
manipulation of f's environemnt with the creation of proto object
p on line ###.  That line converts d to an anonymous proto object
containing the components of d, in this case just x, and then
creates a child object p which can access x via delegation/inheritance.

library(proto)
set.seed(1)
f <- function(lambda) -sum(dpois(x, lambda, log=T))
d <- data.frame(x=rpois(100, 12.34))
p <- proto(as.proto(as.list(d)), f = f) ###
mle(p[["f"]], start=list(lambda=10))

2. Or the ### line could be replaced with the following line
which places f and the components of d, in this case just x,
directly into p:

p <- proto(f = f, envir = as.proto(as.list(d)))

again avoiding the explicit reset of environment(f) and the evalq.

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

        
Gabor Grothendieck wrote:

            

      
      
      
    

        
Agreed.

            

      
      
      
    

        
Fair enough.

            

      
      
      
    

        
Well, my current design specifies a named list: I *think* (but am not
sure) it works gracefully with a data frame as well.  Hadn't thought of
environments -- I'm aiming this more at a lower-level user to whom that
wouldn't occur.  (But I hope it would be possible to design a system
that would be usable by intermediate users and still useful for experts.)

            

      
      
      
    

        
My version still allows a plain likelihood function (I agree that
there will always be situations that are too complicated to encapsulate
as a formula).

            

      
      
      
    

        
HEAR, HEAR.

I think the pattern

            

      
      
      
    

        
OK.

            

      
      
      
    

        
*** I still feel very strongly that end users shouldn't have
to deal with closures, environments, protos, etc. --  I want
mle to LOOK LIKE a standard modeling function if at all possible,
even if it can be used more creatively and flexibly by
those who know how. ***

            

      
      
      
    

        
Agreed.  This would work for formulas, too.

  Have any of you guys looked at bbmle?  The evaluation stuff is
quite ugly, since I was groping around in the dark.  I would love
to clean it up in a way that made everyone happy (?) with it and
possibly allowed it to be merged back into mle.

   Ben

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.6 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iD8DBQFHWWbpc5UpGjwzenMRApxZAJwLYuW+9beykCO1fJvBO4ICZxbEJwCfXgYR
F0nNR+/+/xy11xav9uDZSBE=
=bgiY
-----END PGP SIGNATURE-----

  
  


  


      

    

    

      
      

        


  

        

      Luke Tierney
    

    

      
      #
    

  


  

      
On Fri, 7 Dec 2007, Gabor Grothendieck wrote:

        

            

      
      
      
    

        
I make extensive use of lexical scoping in my programming and I NEVER
use explicit environment manipulaiton--for me that is unreadable, and it
is not amenable to checking using things like codetools.  In the
example above all that is needed is to define x directly, e.g.

     > f <- function(lambda) -sum(dpois(x, lambda, log=T))
     > x <- rpois(10000, 12.34)
     > mle(f, start=list(lambda=10))

     Call:
     mle(minuslogl = f, start = list(lambda = 10))

     Coefficients:
     lambda
     12.337

It isn't necessary to go through the data frame or environment
munging.  If you want to be able to work with likelihoods for several
data sets at once then you can either use diferent names for the
variables, like

     x1 <- rpois(10000, 12.34)
     f1 <- function(lambda) -sum(dpois(x1, lambda, log=T))
     x2 <- rpois(10000, 12.34)
     f2 <- function(lambda) -sum(dpois(x2, lambda, log=T))

If you are concerned that x, x1, x2 might have been redefined if you
come back to f1, f2 later (not an issue with typical useage inside a
function but can be an issue at top level) then you can create a
closure that cpatures the particular data set you are using.  The
clean way to do this is with a function that creates the negative log
likelihood, e.g.

     makePoisonNegLogLikelihood <- function(x)
 	function(lambda) -sum(dpois(x, lambda, log=T))

Then you can do

     f <- makePoisonNegLogLikelihood(rpois(10000, 12.34))
     mle(f, start=list(lambda=10))

which I find much cleaner and easier to understand than environment
munging.  Once you are defining a likelihood constructor you can think
about things like making it a bit more efficient by calculating
sufficient statistics once, for example

     makePoisonNegLogLikelihood <- function(x) {
         sumX <- sum(x)
         n <- length(x)
 	function(lambda) -dpois(sumX, n * lambda, log=T)
     }

Best,

luke

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Luke Tierney
    

    

      
      #
    

  


  

      
On Fri, 7 Dec 2007, Duncan Murdoch wrote:

        

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Luke Tierney
    

    

      
      #
    

  


  

      
On Fri, 7 Dec 2007, Duncan Murdoch wrote:

        

            

      
      
      
    

        
For working at the general likelihood I think is is better to
encourage the approach of definign likelihood constructor functions.
The problem with using f, data is that you need to mathc the names
used in f and in data, so either you have to explicitly write out f
with the names you have in data or you have to modify data to use the
names f likes -- in the running example think

     f <- function(lambda) -sum(dpois(x, lambda, log=T))
     d <- data.frame(y=rpois(10000, 12.34))

somebody has to connext up the x in f with the y in d. With a negative
log likelihood constructor defines, for example, as

     makePoisonNegLogLikelihood <- function(x)
         function(lambda) -sum(dpois(x, lambda, log=T))

this happens naturally with

     makePoisonNegLogLikelihood(d$y)

            

      
      
      
    

        
Both (simple) bootstrapping and (simple leave-one-out) crossvalidation
require a data structure with a notion of cases, which is much more
restrictive than the conext in which mle can be used.  A more ngeneric
aproach to bootstrapping that might fit closer to the level of
generality of mle might be parameterized in terms of a negative log
likelihood constructor, a starting value constructor, and a resampling
function, with a single iteration implemented soemthing like

     mleboot1 <- function(nllmaker, start, esample)  {
 	newdata <- resample()
 	newstart <- do.call(start, newdata)
 	nllfun <- do.call(nllmaker, newdata)
 	mle(fnllfun, start = newstart)
     }

This would leave decisions on the resampling method and data structure
up to the user. Somehing similar could be done with K-fold CV.

luke

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

        
Luke Tierney wrote:

            

      
      
      
    

        
I hate to sound like a jerk, but I do hope that in the end we come
up with a solution that will still be accessible to people who don't
quite have the concept of writing functions to produce functions.  I
feel it is "natural" for people who have multiple data sets to have the
variables named similarly in different data sets.  All of the
constructor stuff is still accessible to anyone who wants to use the
function that way ... is there any way to do a cheesy default
constructor that is just equivalent to taking the likelihood function
and arranging for it to be evaluated in an environment containing
the data?  That way if "nllmaker" below were just a formula
or a log-likelihood function it could still work ...

  [snip]

            

      
      
      
    

        
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.6 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iD8DBQFHWcS1c5UpGjwzenMRAig2AJ9iTzhI1p8tBb7Q15jgT4nA+Zds+gCgggc2
sI2que28Hl1M5cVGa+anEL0=
=hCiS
-----END PGP SIGNATURE-----

  
  


  


      

    

    

      
      

        


  

        

      Antonio, Fabio Di Narzo
    

    

      
      #
    

  


  

      
2007/12/7, Ben Bolker <bolker at ufl.edu>:

            

      
      
      
    

        
I don't really agree with this.
I found really natural writing functions which builds other functions,
for handling
in a clean way the data-dependency problem, much more than
manipilating function environments.
As a useR, I think that if I'm able to write a likelihood function myself:

data <- whatever
negloglik <- function(theta)
  a + very * complicated / function - of %% theta %*% and %o% data

to be used in mle, I'm also good at abstracting it a bit this way:

nllmaker <- function(data)
  function(theta)
    a + very * complicated / function - of %% theta %*% and %o% data

negloglik <- nllmaker(whatever),

don't you think? I use this kind of tricks routinely for simulations.
In general, I think it should be more emphatized functional style in R coding.
In fact, I like a lot the recent introduction of some higher order
functions in the base package (Reduce, Filter, Map).

Bests,
Antonio, Fabio.

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

      
Luke Tierney wrote:

        
[misc snippage]

            

      
      
      
    

        
[more snippage]

That's not really worse than having to match the names in a model 
formula to the names of the data frame in lm(), is it?

The thing that I'm looking for in these matters is a structure which 
allows us to operate on likelihood functions in a rational way, e.g. 
reparametrize them, join multiple likelihoods with some parameters in 
common, or integrate them. The join operation is illustrative: You can 
easily do 

negljoint <- function(alpha, beta, gamma, delta)
    negl1(alpha, beta, gamma) + negl2(beta, gamma, delta)

and with a bit of diligence, this could be the result of Join(negl1, 
negl2). But if the convention is that likelihods have their their data 
as an argument, you also need to also automatically define a data 
argument fot negljoint, (presumably a list of two) and organize that the 
calls to negl1 and negl2 contains the appropriate subdata. It is the 
sort of thing that might be doable, but you'd rather do without.

-pd

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Luke Tierney
    

    

      
      #
    

  


  

      
On Fri, 7 Dec 2007, Ben Bolker wrote:

        

            

      
      
      
    

        
Any programming language has some idioms and conventions that are
worth lerning if you are going to make most effective use of the
language.  R is no exception.  One can use R as interactive C but the
results aren't likely to be too satisfactory.  Three ideas worth
learning in R (not an exclusive list) are how to use vectorized
arithmetic, how to use the apply family of functions, and how to take
advantage of lexical scope.  None of hese is hard to learn, and basic
lexical scope may be the easiest of the three, and the small
investment will pay off.

Best,

luke

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Luke Tierney
    

    

      
      #
    

  


  

      
On Sat, 8 Dec 2007, Peter Dalgaard wrote:

        

            

      
      
      
    

        
Yes and no.

If the likelihood is simple engough to include in line, is in

     d <- data.frame(y=rpois(100,12.34))
     mle(function(lambda) -sum(dpois(d$y, lambda, log = TRUE)),
         start = list(lambda=10))

or neaarly in line, eg in a with or local construct, like

     with(d, {
         f <- function(lambda) -sum(dpois(y, lambda, log = TRUE))
         mle(f, start = list(lambda=10))
     })

or

     local({
         y <- d$y
         f <- function(lambda) -sum(dpois(y, lambda, log = TRUE))
         mle(f, start = list(lambda=10))
     })

then I think it is essentially the same.  But if the function is
complex enough that you will want to define and debug it separately
then you will probably want to be able to reuse your code directly,
not with copy-paste-edit.  At that point things are different.

In a sense this difference also exists with model formulas as well. We
usually write formular in line, rather than something like

     f <- y ~ x
     lm(f)

With simple formulalas that is reasonable. But it would be nice to be
able to abstract out common patterns of more complex fomulas for
simple reuse. A simple-minded example might be to be able to define a
splitPlot formula operator so one can write

     yield ~ splitPlot(whole = fertilizer, sub = variety)

This sort of thing would become more useful in more complicated
multi-level models.  I could be wrong but I don't think BUGS has the
ability to abstract out submodel patterns in this way.  Don't know if
any of the other multi-level modeling systems provide this.  Might be
worth looking into; it's not unrelated to the issues you raise below.

luke