[R-sig-ME] how to extract the BIC value
Martin Maechler
maechler at stat.math.ethz.ch
Tue May 18 12:37:21 CEST 2010
>>>>> "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
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
>>>
>>> _______________________________________________
>>> 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
More information about the R-sig-mixed-models
mailing list