[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