[R] About object of class mle returned by user defined functions
Prof Brian Ripley
ripley at stats.ox.ac.uk
Fri Jul 22 15:47:39 CEST 2005
confint has to be able to re-fit the function to form the profile
likelihood. The fit you return refers to values inside the function you
used, and those are not available in the environment you call confint
from. You need to ensure that those values are substituted and not
referred to.
Compare
> fitA at call
mle(minuslogl = llgamma, start = list(shape = (mean(isi1)/sd(isi1))^2,
scale = sd(isi1)^2/mean(isi1)), method = "L-BFGS-B", lower = c(1e-05,
1e-05))
> fit1 at call
mle(minuslogl = minusLogLikelihood, start = initial.para, method =
"L-BFGS-B",
lower = optim.lower, upper = optim.upper)
and the difference should be clear.
You can fix this up by constructing a call containing the values and not
the names (one way is to use substitute()) and then eval() it. Another
way is something like
Call <- quote(mle())
Call$minuslogl <- minusLogLikelihood
Call$start <- initial.para
...
On Fri, 22 Jul 2005, Christophe Pouzat wrote:
> Guys,
>
> I apologize for being slightly misleading in my previous e-mail.
>
> First, I generated some confusion between the scale and rate parameters in
> the gamma distribution. My direct call to mle use a minuslogl function
> "working" with a scale parameter while my call to mle from my function used a
> minuslogl function "working" with a rate parameter!... To add to the
> confusion I had simulated data with a scale ( = 1/rate) value of 1... I
> really hope that none of you lost time with that.
>
> Second, some "^2" in my original function definition got converted into
> exponents on the e-mail, meaning that if some of you tried to copy and paste
> it you must have gotten some insults from R while sourcing it. In order to
> avoid that I attach an ".R" file. In principle if you source it and then type
> the following commands you should get (exactly):
>
>
>> coef(fitA)
> shape scale
> 2.2230421 0.8312374
>> coef(fit1)
> shape scale
> 2.2230421 0.8312374
>> vcov(fitA)
> shape scale
> shape 0.08635158 -0.03228829
> scale -0.03228829 0.01518126
>> vcov(fit1)
> shape scale
> shape 0.08635158 -0.03228829
> scale -0.03228829 0.01518126
>> logLik(fitA)
> 'log Lik.' -146.6104 (df=2)
>> logLik(fit1)
> 'log Lik.' -146.6104 (df=2)
>> confint(fitA)
> Profiling...
> 2.5 % 97.5 %
> shape 1.6985621 2.853007
> scale 0.6307824 1.129889
>> 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)
>
> Here fitA is obtained by a direct call to mle (I mean from the command line)
> while fit1 is obtained by the same call but within a function: newFit.
>
> The fundamental problem remains, I don't understand why confint does work
> with fitA and not with fit1.
>
> Christophe.
>
> PS: my version info
>
> platform i686-pc-linux-gnu
> arch i686 os linux-gnu system i686, linux-gnu
> status major 2 minor 1.1
> year 2005 month 06 day 20
> language R
>
>
>
> Christophe Pouzat wrote:
>
>> 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
>
>
--
Brian D. Ripley, ripley at stats.ox.ac.uk
Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel: +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UK Fax: +44 1865 272595
More information about the R-help
mailing list