[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 13:05:27 CEST 2010


>>>>> "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")))

 {well, modulo the fact that "..." should really allow to do
  this for *several* models simultaneously}

In addition to that (and more replying to Doug Bates):

Given nlme's tradition of explicitly providing BIC(), and in
analogue to the S3 semantics of the AIC() methods,

- I think lme4 (and "lme4a" on R-forge) should end up having
  working  AIC() and BIC() directly for fitted models, instead of
  having to use
	 AIC(logLik(.))   and    BIC(logLik(.))

  The reason that even the first of this currently does *not*
  work is that lme4 imports AIC from "stats" but should do so
  from "stats4".
  --> I'm about to change that for 'lme4' (and 'lme4a').
    
  However, for the BIC case, ... see below


- I tend to agree with Gabor (for once! :-)  that
  basic BIC methods (S3, alas) should move from  nlme to stats.
  
  For this reason, I'm breaking the rule of "do not cross-post"
  for once, and am hereby diverting this thread to R-devel

Martin


    MM> --
    MM> 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
    >>>>



More information about the R-devel mailing list