[Rd] BIC() in "stats" {was [R-sig-ME] how to extract the BIC value}
Martin Maechler
maechler at stat.math.ethz.ch
Tue May 18 18:01:50 CEST 2010
Adding to my own statements (below) :
>>>>> "MM" == Martin Maechler <maechler at stat.math.ethz.ch>
>>>>> on Tue, 18 May 2010 13:05:27 +0200 writes:
>>>>> "MM" == Martin Maechler <maechler at stat.math.ethz.ch>
>>>>> on Tue, 18 May 2010 12:37:21 +0200 writes:
>>>>> "GaGr" == Gabor Grothendieck <ggrothendieck at gmail.com>
>>>>> on Mon, 17 May 2010 09:45:00 -0400 writes:
GaGr> BIC seems like something that would logically go into stats in the
GaGr> core of R, as AIC is already, and then various packages could define
GaGr> methods for it.
MM> Well, if you look at help(AIC):
>>> Usage:
>>> AIC(object, ..., k = 2)
>>> Arguments:
>>> object: a fitted model object, for which there exists a ‘logLik’
>>> method to extract the corresponding log-likelihood, or an
>>> object inheriting from class ‘logLik’.
>>> ...: optionally more fitted model objects.
>>> k: numeric, the _penalty_ per parameter to be used; the default
>>> ‘k = 2’ is the classical AIC.
MM> you may note that the original authors of AIC where always
MM> allowing the AIC() function (and its methods) to compute the BIC,
MM> simply by using 'k = log(n)' where of course n must be correct.
MM> I do like the concept that BIC is just a variation of AIC and
MM> AFAIK, AIC was really first (historically).
MM> Typically (and with lme4), the 'n' needed is already part of the logLik()
MM> attributes :
>>> AIC((ll <- logLik(fm1)), k = log(attr(ll,"nobs")))
MM> REML
MM> 1774.786
MM> indeed gives the BIC (where the "REML" name may or may not be a
MM> bit overkill)
MM> A stats-package based BIC function could then simply be defined as
>> BIC <- function (object, ...) UseMethod("BIC")
>> BIC.default <- function (object, ...) BIC(logLik(object), ...)
>> BIC.logLik <- function (object, ...)
>> AIC(object, ..., k = log(attr(object,"nobs")))
MM> {well, modulo the fact that "..." should really allow to do
MM> this for *several* models simultaneously}
MM> In addition to that (and more replying to Doug Bates):
MM> Given nlme's tradition of explicitly providing BIC(), and in
MM> analogue to the S3 semantics of the AIC() methods,
MM> - I think lme4 (and "lme4a" on R-forge) should end up having
MM> working AIC() and BIC() directly for fitted models, instead of
MM> having to use
MM> AIC(logLik(.)) and BIC(logLik(.))
MM> The reason that even the first of this currently does *not*
MM> work is that lme4 imports AIC from "stats" but should do so
MM> from "stats4".
MM> --> I'm about to change that for 'lme4' (and 'lme4a').
MM> However, for the BIC case, ... see below
MM> - I tend to agree with Gabor (for once! :-) that
MM> basic BIC methods (S3, alas) should move from nlme to stats.
MM> For this reason, I'm breaking the rule of "do not cross-post"
MM> for once, and am hereby diverting this thread to R-devel
What I *did* find is that the stats4 package has already had
all necessary BIC methods -- S4, not S3.
So for lme4 (and R-forge's "lme4a"),
I've only needed to change the NAMESPACE file to have both
importFrom("stats4", AIC, BIC, logLik)# so S4 methods are used!
and later
export(AIC, BIC,
.....)
and also add 'stats4' to the 'Imports: ' line in DESCRIPTION.
So both (development versions of) lme4 and lme4a now have
working AIC() and BIC(),
and I guess Doug could release a new version of lme4 (not .."a")
pretty soon.
I got private e-mail suggestions for extensive S3 methods for
AIC, BIC and logLik.
I think these should happen more in public (i.e. here, on
R-devel), and while I still advocate that a BIC S3 generic +
simple default methods should be added (as above),
I'd be happy if others joined into the discussion,
(and possibly provided simple patches).
Martin Maechler, ETH Zurich
GaGr> On Mon, May 17, 2010 at 9:29 AM, Douglas Bates <bates at stat.wisc.edu> wrote:
>>>> On Mon, May 17, 2010 at 5:54 AM, Andy Fugard (Work)
>>>> <andy.fugard at sbg.ac.at> wrote:
>>>>> Greetings,
>>>>>
>>>>> Assuming you're using lmer, here's an example which does what you need:
>>>>>
>>>>>> (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
>>>>> Linear mixed model fit by REML
>>>>> Formula: Reaction ~ Days + (Days | Subject)
>>>>> Data: sleepstudy
>>>>> AIC BIC logLik deviance REMLdev
>>>>> 1756 1775 -871.8 1752 1744
>>>>> Random effects:
>>>>> Groups Name Variance Std.Dev. Corr
>>>>> Subject (Intercept) 612.092 24.7405
>>>>> Days 35.072 5.9221 0.066
>>>>> Residual 654.941 25.5918
>>>>> Number of obs: 180, groups: Subject, 18
>>>>>
>>>>> Fixed effects:
>>>>> Estimate Std. Error t value
>>>>> (Intercept) 251.405 6.825 36.84
>>>>> Days 10.467 1.546 6.77
>>>>>
>>>>> Correlation of Fixed Effects:
>>>>> (Intr)
>>>>> Days -0.138
>>>>>
>>>>>> (fm1fit <- summary(fm1)@AICtab)
>>>>> AIC BIC logLik deviance REMLdev
>>>>> 1755.628 1774.786 -871.8141 1751.986 1743.628
>>>>>
>>>>>> fm1fit$BIC
>>>>> [1] 1774.786
>>>>
>>>> That's one way of doing it but it relies on a particular
>>>> representation of the object returned by summary, and that is subject
>>>> to change.
>>>>
>>>> I had thought that it would work to use
>>>>
>>>> BIC(logLik(fm1))
>>>>
>>>> but that doesn't because the BIC function is imported from the nlme
>>>> package but not later exported. The situation is rather tricky - at
>>>> one point I defined a generic for BIC in the lme4 package but that led
>>>> to conflicts when multiple packages defined different versions. The
>>>> order in which the packages were loaded became important in
>>>> determining which version was used.
>>>>
>>>> We agreed to use the generic from the nlme package, which is what is
>>>> now done. However, I don't want to make the entire nlme package
>>>> visible when you have loaded lme4 because of resulting conflicts.
>>>>
>>>> I can get the result as
>>>>
>>>>> (fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
>>>> Linear mixed model fit by REML
>>>> Formula: Reaction ~ Days + (Days | Subject)
>>>> Data: sleepstudy
>>>> AIC BIC logLik deviance REMLdev
>>>> 1756 1775 -871.8 1752 1744
>>>> Random effects:
>>>> Groups Name Variance Std.Dev. Corr
>>>> Subject (Intercept) 612.090 24.7405
>>>> Days 35.072 5.9221 0.066
>>>> Residual 654.941 25.5918
>>>> Number of obs: 180, groups: Subject, 18
>>>>
>>>> Fixed effects:
>>>> Estimate Std. Error t value
>>>> (Intercept) 251.405 6.825 36.84
>>>> Days 10.467 1.546 6.77
>>>>
>>>> Correlation of Fixed Effects:
>>>> (Intr)
>>>> Days -0.138
>>>>> nlme:::BIC(logLik(fm1))
>>>> REML
>>>> 1774.786
>>>>
>>>> but that is unintuitive. I am not sure what the best approach is.
>>>> Perhaps Martin (or anyone else who knows namespace intricacies) can
>>>> suggest something.
>>>>
>>>>
>>>>> Tahira Jamil wrote:
>>>>>> Hi
>>>>>> I can extract the AIC value of a model like this
>>>>>>
>>>>>> AIC(logLik(fm0)
>>>>>>
>>>>>> How can I extract the BIC value if I need!
>>>>>>
>>>>>> Cheers
>>>>>> Tahira
>>>>>> Biometris
>>>>>> Wageningen University
>>>>>>
>>>>>> _______________________________________________
>>>>>> R-sig-mixed-models at r-project.org mailing list
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>>
>>>>>
>>>>> --
>>>>> Andy Fugard, Postdoctoral researcher, ESF LogICCC project
>>>>> "Modeling human inference within the framework of probability logic"
>>>>> Department of Psychology, University of Salzburg, Austria
>>>>> http://www.andyfugard.info
>>>>>
MM> ______________________________________________
MM> R-devel at r-project.org mailing list
MM> https://stat.ethz.ch/mailman/listinfo/r-devel
More information about the R-devel
mailing list