In some areas, it is common to fit a model separately to multiple clusters in a data set, for instance fitting a cognitive model separately to data from each participant in an experiment.
Model comparison is a bit more complicated in this scenario, since rather than having a single deviance, AIC, BIC, or Bayes Factor per model, we have one score per model per participant.
A way that works for Bayes Factors
Stephan et al (2009; NeuroImage) discuss this issue with an emphasis on Bayesian analysis of MRI data, and identify two approaches.
In the Fixed Effects approach, we assume that all subjects’ data are generated by the same model, and so the Bayes Factor for the group is just the product of individual participants' Bayes Factors:
$$ BF_{\text{Group}} = \prod_i BF_i $$
In the Random Effects approach, introduced in that paper, we assume that there is a true distribution of models in the population that is, some participants' data is generated by model 1, some by model 2, etc., the distribution of these models is described by a multinomial distribution with model probability parameters $r$, and the posterior distribution over these model probability parameters is described by a Dirichlet distribution with concentration parameters $\alpha$.
$$ \begin{align} \text{Data}_i &\sim \text{Model}_i\\ \text{Model}_i &\sim \text{Multinomial}(r)\\ r &\sim \text{Dirichlet}(\alpha) \end{align} $$
Estimating this model allows us to infer useful quantities such as how likely it is that a specific model generated the data of a randomly chosen subject, and what is the probability that model M is the most prevalent in the population. In practice, these parameters are estimated using either model Bayes Factors, marginal likelihood from variational methods, or BIC scores.
Doing the same for AIC
My question is whether similar methods exist for evaluating collections of AIC scores obtained by fitting various models to various participants?
It seems pretty reasonable to calculate
$$ \begin{align} AIC_{\text{Group}} &= \sum_i^n AIC_i\\ &= \sum_i^n 2k_i + (\sum_i^n -2\text{ln}(\hat L_i)) \end{align} $$
since this is the same as calculating AIC for a single model with $nk$ parameters and log-likelihood of $\sum_i^n \text{ln}(\hat L_i)$.
I've also seen some papers just running t-tests on the individual AIC scores:
t.test(aic.score ~ model, paired=T, data=aic.scores)
Is there any more principled solution to this problem?
Notes
- An R implementation of the Stephan et al (2009) random effects procedure can be found here.