Martin Maechler
2010-May-18 11:05 UTC
[Rd] BIC() in "stats" {was [R-sig-ME] how to extract the BIC value}
>>>>> "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 >>>>
Martin Maechler
2010-May-18 16:01 UTC
[Rd] BIC() in "stats" {was [R-sig-ME] how to extract the BIC value}
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