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

From: Christophe Pouzat <christophe.pouzat_at_univ-paris5.fr>
Date: Fri 22 Jul 2005 - 02:01:36 EST


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@minuslogl) with argument values different from the MLE does return something sensible:

 > fit1@minuslogl(coef(fit1)[1],coef(fit1)[2]) [1] 155.9232
 > fit1@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

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Fri Jul 22 02:28:02 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:33:55 EST