[R] About object of class mle returned by user defined functions

[Christophe Pouzat]

Hi,

There is something I don't get with object of class "mle" returned by a 
function I wrote. More precisely it's about the behaviour of method 
"confint" and "profile" applied to these object.

I've written a short function (see below) whose arguments are:
1) A univariate sample (arising from a gamma, log-normal or whatever).
2) A character string standing for one of the R densities, eg, "gamma", 
"lnorm", etc. That's the density the user wants to fit to the data.
3) A named list with initial values for the density parameters; that 
will be passed to optim via mle.
4) The method to be used by optim via mle. That can be change by the 
code if parameter boundaries are also supplied.
5) The lowest allowed values for the parameters.
6) The largest allowed values.

The "big" thing this short function does is writing on-fly the 
corresponding log-likelihood function before calling "mle". The object 
of class "mle" returned by the call to "mle" is itself returned by the 
function.

Here is the code:

newFit <- function(isi, ## The data set
                   isi.density = "gamma", ## The name of the density 
used as model
                   initial.para = list( shape = (mean(isi)/sd(isi))^2,
                     scale = sd(isi)^2 / mean(isi) ), ## Inital 
parameters passed to optim
                   optim.method = "BFGS", ## optim method
                   optim.lower = numeric(length(initial.para)) + 0.00001,
                   optim.upper = numeric(length(initial.para)) + Inf,
                   ...) {

  require(stats4)
 
  ## Create a string with the log likelihood definition
  minusLogLikelihood.txt <- paste("function( ",
                                  paste(names(initial.para), collapse = 
", "),
                                  " ) {",
                                  "isi <- eval(",
                                  deparse(substitute(isi)),
                                  ", envir = .GlobalEnv);",
                                  "-sum(",
                                  paste("d", isi.density, sep = ""),
                                  "(isi, ",
                                  paste(names(initial.para), collapse = 
", "),
                                  ", log = TRUE) ) }"
                                  )

  ## Define logLikelihood function
  minusLogLikelihood <- eval( parse(text = minusLogLikelihood.txt) )
  environment(minusLogLikelihood) <- .GlobalEnv

 
  if ( all( is.infinite( c(optim.lower,optim.upper) ) ) ) {
      getFit <- mle(minusLogLikelihood,
                    start = initial.para,
                    method = optim.method,
                    ...
                    )
  } else {
    getFit <- mle(minusLogLikelihood,
                  start = initial.para,
                  method = "L-BFGS-B",
                  lower = optim.lower,
                  upper = optim.upper,
                  ...
                  )
  }  ## End of conditional on all(is.infinite(c(optim.lower,optim.upper)))
 
  getFit
 
}


It seems to work fine on examples like:

 > isi1 <- rgamma(100, shape = 2, scale = 1)
 > fit1 <- newFit(isi1) ## fitting here with the "correct" density 
(initial parameters are obtained by the method of moments)
 > coef(fit1)
    shape     scale
1.8210477 0.9514774
 > vcov(fit1)
           shape      scale
shape 0.05650600 0.02952371
scale 0.02952371 0.02039714
 > logLik(fit1)
'log Lik.' -155.9232 (df=2)

If we compare with a "direct" call to "mle":

 > llgamma <- function(sh, sc) -sum(dgamma(isi1, shape = sh, scale = sc, 
log = TRUE))
 > fitA <- mle(llgamma, start = list( sh = (mean(isi1)/sd(isi1))^2, sc = 
sd(isi1)^2 / mean(isi1) ),lower = c(0.0001,0.0001), method = "L-BFGS-B")
 > coef(fitA)
      sh       sc
1.821042 1.051001
 > vcov(fitA)
            sh          sc
sh  0.05650526 -0.03261146
sc -0.03261146  0.02488714
 > logLik(fitA)
'log Lik.' -155.9232 (df=2)

I get almost the same estimated parameter values, same log-likelihood 
but not the same vcov matrix.

A call to "profile" or "confint" on fit1 does not work, eg:
 > confint(fit1)
Profiling...
Erreur dans approx(sp$y, sp$x, xout = cutoff) :
    need at least two non-NA values to interpolate
De plus : Message d'avis :
collapsing to unique 'x' values in: approx(sp$y, sp$x, xout = cutoff)

Although calling the log-likelihood function defined in fit1 
(fit1 at minuslogl) with argument values different from the MLE does return 
something sensible:

 > fit1 at minuslogl(coef(fit1)[1],coef(fit1)[2])
[1] 155.9232
 > fit1 at minuslogl(coef(fit1)[1]+0.01,coef(fit1)[2]+0.01)
[1] 155.9263

There is obviously something I'm missing here since I thought for a 
while that the problem was with the environment "attached" to the 
function "minusLogLikelihood" when calling "eval"; but the lines above 
make me think it is not the case...

Any help and/or ideas warmly welcomed.

Thanks,

Christophe.

-- 
A Master Carpenter has many tools and is expert with most of them.If you
only know how to use a hammer, every problem starts to look like a nail.
Stay away from that trap.
Richard B Johnson.
--

Christophe Pouzat
Laboratoire de Physiologie Cerebrale
CNRS UMR 8118
UFR biomedicale de l'Universite Paris V
45, rue des Saints Peres
75006 PARIS
France

tel: +33 (0)1 42 86 38 28
fax: +33 (0)1 42 86 38 30
web: www.biomedicale.univ-paris5.fr/physcerv/C_Pouzat.html