how to extract the BIC value

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

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.

you may note that the original authors of AIC where always
allowing the AIC() function (and its methods) to compute the BIC,
simply by using 'k = log(n)' where of course n  must be correct.

I do like the concept that BIC is just a variation of AIC and
AFAIK, AIC was really first (historically).

Typically (and with lme4), the 'n' needed is already part of the logLik()
attributes :

AIC((ll <- logLik(fm1)), k = log(attr(ll,"nobs")))

REML 
1774.786 

indeed gives the BIC (where the "REML" name may or may not be a
bit overkill)

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


--
Martin Maechler, ETH Zurich

>> On Mon, May 17, 2010 at 5:54 AM, Andy Fugard (Work)

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

>>>> 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
    >>> 
    >>> _______________________________________________
    >>> R-sig-mixed-models at r-project.org mailing list
    >>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
    >>> 
    >> 
    >> _______________________________________________
    >> R-sig-mixed-models at r-project.org mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
    >> 

    GaGr> _______________________________________________
    GaGr> R-sig-mixed-models at r-project.org mailing list
    GaGr> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models