AIC/BIC?

[DunovanK]

Hi all,

Is there an easy way to calculate AIC and BIC for model comparison in HDDM? A reviewer has requested that I include these in my manuscript. 

Thank you!

Kyle

[Imri Sofer]

AIC and BIC are not appropriate for hierarchical models.

You can use the DIC (which is related to the AIC) by  using model.print_stats()

--
You received this message because you are subscribed to the Google Groups "hddm-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hddm-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[DunovanK]

Hi Imri,

Thanks for the quick response. Yes, for the hierarchical models, I have reported the DIC.

What about individual models fit using the Bayesian MCMC approach? Or models fit using quantile optimization (hddm.optimize('ML'))?

Thanks!

Kyle

[DunovanK]

Also, I would be very grateful for any recommendations for papers to cite stating that AIC and BIC are not appropriate for comparing hierarchical models.

Thanks!

Kyle

[Imri Sofer]

regarding the citation: not sure. maybe the original DIC paper. also check out the DIC page of the BUGS software.

one problem with them is that it is not clear what is the effective number of parameters you have, when using a hierarchical model.

For the individual models:

optimize('ML') does not use quantiles, it use ML on the original data. I didn't implement AIC/BIC but it should be fairly easy. To get the log likelihood you can just take the observed nodes and sum their logp:

maybe something like this:

sum([x.logp for x in model.get_observeds()['node']])

Then use the equations from wikipedia.

also, a general comment:

all information criterions come with their own assumptions. A much better approach is cross validation. You cannot wrong with that. But sometimes its better not to argue with the reviewer...

[DunovanK]

Imri,

Thank you very much, that helps me out a lot. I'll post again if I run into any issues, but I think that should take care of it.

Thanks again,

Kyle

[Thomas Wiecki]

Hi Kyle,

I implemented AIC and BIC in separate code. You could inherit from hddm.models.HDDM and add these properties or just rewrite them to take as input the model instead of self.

    @property
    def aic(self):
        k = len(self.get_stochastics())
        return 2 * k - self.mc.logp

    @property
    def bic(self):
        k = len(self.get_stochastics())
        n = len(self.data)
        return -2 * self.mc.logp + k * np.log(n)

Thomas

Thomas Wiecki

PhD candidate, Brown University

Quantitative Researcher, Quantopian Inc, Boston

[Imri Sofer]

I think that self.mc.logp returns that logp of the whole model and not only of the observed nodes, so it would not give you what's you need

[Thomas Wiecki]

Thanks, Imri.

Here is the updated code:

    @property
    def aic(self):
        k = len(self.get_stochastics())

        logp = sum([x.logp for x in self.get_observeds()['node']])
        return 2 * k - logp

    @property
    def bic(self):
        k = len(self.get_stochastics())
        n = len(self.data)

        logp = sum([x.logp for x in self.get_observeds()['node']])
        return -2 * logp + k * np.log(n)

[Thomas Wiecki]

Joachim alerted me to the fact that AIC is missing a factor. The updated AIC and BIC code should thus be:

    @property
    def aic(self):
        k = len(self.get_stochastics())
        logp = sum([x.logp for x in self.get_observeds()['node']])

        return 2 * k - 2 * logp

    @property
    def bic(self):
        k = len(self.get_stochastics())
        n = len(self.data)

        logp = sum([x.logp for x in self.get_observeds()['node']])
        return -2 * logp + k * np.log(n)

[DunovanK]

Sorry it's taken me so long to respond. Thank you so much for the information everyone. T'is very much appreciated. 

Best,

Kyle