article: Computational modeling plays an important role in modern neuroscience research. Much previous research has relied on statistical methods, separately, to address two problems that are actually interdependent. First, given a particular computational model, Bayesian hierarchical techniques have been used to estimate individual variation in parameters over a population of subjects, leveraging their population-level distributions. Second, candidate models are themselves compared, and individual variation in the expressed model estimated, according to the fits of the models to each subject. The interdependence between these two problems arises because the relevant population for estimating parameters of a model depends on which other subjects express the model. Here, we propose a hierarchical Bayesian inference (HBI) framework for concurrent model comparison, parameter estimation and inference at the population level, combining previous approaches. We show that this framework has important advantages for both parameter estimation and model comparison theoretically and experimentally. The parameters estimated by the HBI show smaller errors compared to other methods. Model comparison by HBI is robust against outliers and is not biased towards overly simplistic models. Furthermore, the fully Bayesian approach of our theory enables researchers to make inference on group-level parameters by performing HBI t-test.
  Across different areas of neuroscience, researchers increasingly employ computational models for experimental data analysis. For example, decision neuroscientists use reinforcement learning (RL) and economic models of choice to analyze behavioral and brain imaging data in reward learning and decision-making tasks [1,2]. The field of computational psychiatry uses these models to characterize patients and people at the risk of brain disorders [3–6]. Neuroimaging studies use models of neural interaction, such as dynamic causal modeling [7,8], as well as abstract models to analyze brain signals [2,9]. The success of these efforts heavily depends on statistical methods making inference about validity and robustness of estimated parameters across individuals, as well as making inference on validity and generalizability of computational models. A key theoretical and practical issue has been capturing individual variation both in a model’s parameters and additionally in which of several candidate models a subject expresses, which may also vary from subject to subject. Computational models usually rely on free parameters, such as learning rate in RL models, which often capture quantities of scientific interest but typically vary across individuals and must be estimated from data. A dataset includes a number of subjects, and often the question of interest is to characterize parameters in a population: Is choice consistency altered in patients with attention-deficit hyperactive disorders? Do cognitive enhancers, such as Ritalin, enhance the learning rate at the population level? These questions are most naturally framed in terms of hierarchical models, which characterize both the population distributions over a model’s parameters and also each individual subject’s parameters given the population distribution. Since these two levels are mutually interrelated, they are often estimated simultaneously, using methods like expectation maximization or sampling (MCMC). For example, the hierarchical parameter estimation (HPE) procedure [10,11] regularizes individual estimates according to group statistics, producing better individual estimates and permitting reliable group-level tests. Because subjects typically share underlying structure, hierarchical Bayesian approaches can leverage this structure to yield better individual estimates and to provide better predictions for unseen data, compared to approaches that fit each subject separately [12]. A second, and seemingly logically prior, question is which of several candidate models provides the best explanation for the data. This is important both for providing the setting within which to do parameter estimation, and also for investigating questions of scientific interest. Are rodents’ reaction times best explained by independent or competing accumulators? Do compulsive gamblers rely more on model-free RL compared to controls? Importantly, in principle (and apparently in practice) the model expressed might also vary from subject to subject; thus modern model comparison techniques rely on estimating which of several models obtains for each subject [13]. Estimating such variation is important since the prior assumption that the same model obtains across all individuals (treating model identity as a fixed effect) is a very strong (and in most cases potentially unwarranted) assumption, which makes model comparison very sensitive to outliers [13]. To estimate this variation, in turn, depends on the likelihood of each subject’s data given each model (and, thus, on each subject’s parameters for each model). Intuitively, evaluating whether a model is a good model for a subject’s data precedes estimation of its specific parameter values; and indeed, previous research has used separate tools to solve these two problems. But statistically, the two questions are actually interconnected, because individual parameters and hence individual fit depend on which subjects belong to the population that expresses the model. Here, we address this challenge from a fully Bayesian viewpoint. This work addresses issues of statistical inference over both parameters and models, which have remained elusive with the previous hierarchical methods. Notably, although it is accepted (for the reasons discussed above) that the best-fitting model may vary from subject to subject, hierarchical parameter estimation (conducted separately) has typically assumed that the given model is expressed over all subjects, i. e. that it is a fixed effect (and if multiple models are compared, these are each fit to the entire population). This assumption biases parameter estimation, at both individual and group levels, because it entails that the estimated parameters for each individual subject equally affect group-level estimates, even though some members of the population may be better understood as expressing altogether different models. This same bias, in turn, affects the estimation of which subjects are best fit by each model. In this work, we introduce a hierarchical and Bayesian inference method, which solves these problems by addressing both model fitting and model comparison within the same framework using variational techniques. Furthermore, our fully Bayesian approach enables us to assess uncertainty and provide a rigorous statistical test, HBI t-test, for making inference about parameters of a model at the population level, an issue that has not been addressed in some previous hierarchical models. This paper is structured as follows. First, we highlight the main theoretical advances in our approach. A full formal treatment is given in Materials and methods and S1 Appendix. We then apply the proposed method to synthetic choice datasets as well as empirical datasets to demonstrate its advantages over previous methods.
  Consider a typical computational modeling study in which data of a group of subjects have been measured and a set of candidate models are considered as possible underlying computational mechanisms generating those data. Such studies have generally two main goals: 1) to compare model evidence across competing models; 2) to estimate free parameters of models for each individual and their group-level distributions. All this is typically characterized in terms of inference in a hierarchically structured model of the data, which captures how each subject’s observations depend on their parameters and the individual parameters on their group distribution. The HPE procedure [10,11] employs a hierarchical approach to define the priors based on statistics of the group. This method typically assumes that for a particular model k, all individual parameters are normally distributed, p (h k n) = N (h k n | μ k, V k), where hkn is a vector of the free parameters of the kth model for subject n, μk and Vk are the mean and variance parameters, respectively, indicating the prior distribution over hkn. It is important to distinguish the statistical model itself from the algorithms or approximations used to estimate it. HPE uses the expectation-maximization algorithm [14], a well-known iterative procedure, for obtaining estimating group parameters μk and Vk and individual parameters hkn. Every iteration of this algorithm alternates two steps: 1) an expectation step in which the individual parameters are estimated in light of the group-level distribution; and 2) a maximization step in which the group parameters, μk and Vk, are updated given the current estimates of the individual parameters. Importantly, reflecting the assumption that all subjects express model k, this update weights the individual subjects’ estimates equally; for instance, the update for μk is given by the average of subject level mean estimates (denoted θkn) across all subjects: μ k = 1 N ∑ n θ k n, where N is the number of subjects. Although HPE characterizes variation across subjects in the model parameters hkn (that is, it treats those parameters as random effects), a critical assumption of the procedure is that the parameters for model k are estimated assuming that the same model is responsible for generating data in all subjects. That is, the model identity is taken as a fixed effect, in contrast to the random effects approach that assumes different models might be responsible for generating data in different subjects. The fixed effects assumption has two important implications: 1) for parameter estimation, group parameters, the group mean μk and variance Vk, are influenced equally by all subjects, even those who would be better fit by some other candidate model j ≠ k; 2) for model comparison, the straightforward procedure (e. g. iBIC from [10,11]) is to compare models according to the sum of individual model evidences over all subjects, i. e. again treating the model identity as a fixed effect. Note that while it is possible to submit individual model evidence values (per subject and model) derived from HPE to a separate model comparison procedure that treats model identity as a random effect (such as random effects model selection [13]), these will be biased both from having been fit under the fixed effects assumption and also due to the optimization of the free group-level parameters. For this reason, HPE has typically been accompanied by fixed-effects model comparison [10,11,15], whereas attempts to study subject-subject variation in model identity [13] have typically been conducted using a different, non-hierarchical parameter estimation procedure. Altogether, violations of the fixed effects assumption can adversely influence both parameter estimation and model comparison. Here, we extend HPE’s generative model with another level of the hierarchy, specifying for each subject which model generated their data. This is governed by a subject-specific multinomial random variable, itself drawn from a distribution controlling the proportion of each model in the population. This, in effect, merges the Bayesian model selection model from Stephan et al. [13] with HPE. To accomplish inference in this model, we then lay out a procedure for joint inference over model identities and parameters, including quantifying the probability that each model is responsible for generating data for each subject. To achieve this goal, we adopt a fully Bayesian framework in which the group parameters for each model, μk and Vk, are also random variables. This also gives us a straightforward way to quantify the level of certainty in group-level estimations. We use mean-field variational Bayes [16,17], an extension of expectation-maximization [18], which is able to deal with multiple latent variables in a probabilistic model. Since HBI is a mean-field variational framework, the resulting algorithm (see Materials and methods) is an iterative algorithm. On every iteration, the HBI performs 4 steps: calculates the summary statistics, updates its estimates of the posterior over group parameters, updates its estimate of the posterior over each individual parameter and finally updates its estimates of responsibility of each model in generating each individual data. The algorithm and other important mathematical issues are given in Materials and methods. Here, we highlight three main results. The mathematical proofs are given in S1 Appendix. As noted above, the HBI method estimates the probability of each subject’s dataset being generated by each model, or the responsibility of model k for generating data for subject n, rkn, which is expressed as (expected) probability. Larger values of rkn (i. e. close to 1) indicate that model k is likely to be the true underlying model of the nth subject. In contrast, smaller values of rkn (close to 0) indicate that model k is unlikely to be the underlying model for the nth subject. Based on the responsibilities, it is then possible to estimate the number of subjects explained by each model, N¯ k: N ¯ k = ∑ n = 1 N r k n. Thus N¯ k is always less than the number of subjects and indexes the predominance of model k in the population. Furthermore, the fraction N¯ k / N is called model frequency, which always lies between 0 and 1 and is a useful and intuitive metric for model comparison. In practice, in many situations, researchers are interested in selecting a single best model (rather than relative comparisons among several) even in the face of variation in model identity across subjects. One way to accomplish this goal is to compute the exceedance probability of each candidate model, a metric commonly used for model selection [13]. Exceedance probability is the probability that model k is more commonly expressed than any other model in the model space. Furthermore, the random effects approach enables us to quantify how likely the observed differences in model evidence is simply due to chance [19]. In this case, model selection is not statistically supported, as there is no meaningful difference between models. A metric called protected exceedance probability [19], which typically is more conservative than the exceedance probability, takes into account this possibility (see Materials and methods). Altogether, the random effects approach results in a more robust model comparison and model selection, one less driven by outliers than fixed-effects methods. Note that previous attempts to do model selection at group level using exceedance probability assumed no hierarchy for parameter estimation, thus did not deal with the issue that parameter estimation was not properly conditionalized by group distributions based on model identity. We noted above that an issue with the HPE is that the influence of subjects on the group parameters is equal, due to the assumption that the model is a fixed effect. However, by virtue of its random effects structure, the comparable parameter in our approach, the mean of posterior distribution over μk, denoted by ak, shows an important property: Algorithmically, a subject’s effect on this parameter depends on the degree to which the model is estimated to be the underlying model for that subject. Specifically, this parameter, ak, is updated at each iteration as: a k = 1 1 + N ¯ k (a 0 + ∑ n r k n θ k n), where θkn is the mean of the individual posterior and a0 is the prior mean over μk. The important point in this equation is that ak is a weighted average of individual parameters, in which the weights are the corresponding responsibilities, rkn. This is not specific to the group mean, but it is rather a general feature of our approach: contribution of model k to group parameters is weighted according to the responsibility of model k in generating data in the nth subject, rkn. As mentioned above, another issue that has been incompletely treated in HPE is related to inference on parameters of a fitted model at the population level. Statistically, one needs the uncertainty of the estimated group mean, μk, to be able to make inference on the corresponding parameter at the group level. Since parameters fitted by the HPE are not independent but instead regularized according to the variance given by data, one cannot employ regular statistical tests, such as t-test, to test whether a specific model parameter is “significantly” different from zero. Using those tests on such parameters is biased in favor of generating a significant p-value (more false positives). The HBI framework solves this problem by quantifying the uncertainty of the posterior over the group parameter, resulting in a statistical test similar to the t-test, which we call it HBI t-test. Specifically, it is possible to show that the posterior over the ith group parameter in model k, μki, takes the form of standard Student’s t-distribution centered at the corresponding group mean, aki, with n k = 1 + N¯ k as degrees of freedom. The resulting t-value takes an intuitive form: t = μ k i − a k i s k i / n k, where ski is the empirical deviance statistics for the ith parameter of model k. Therefore, s k i / n k plays the role of standard error, which we call it hierarchical error. Note that the degrees of freedom of the test depend on the number of subjects (i. e. evidence) in favor of model k given by N¯ k, not the total number of subjects. Other group statistics, aki and ski, are also weighted according to the responsibilities of model k in generating data of each subject (as formally obtained in Materials and methods). Using this marginal distribution for population-level group parameters, the HBI t-test enables researchers to determine whether a parameter is significantly different from an arbitrary value, notably 0. For example, the parameter is significantly different from 0 at P < 0. 05 if 0 does not fall within the 95% credible interval. In this section, we apply the proposed HBI method to synthetic datasets and compare its performance with that of HPE, as well as with a non-hierarchical inference (NHI) method estimating parameters for each subject independently according to some fixed, a priori Gaussian priors [20–23]. Importantly, these methods differ in their statistical assumptions about the generative process of data. The NHI assumes no hierarchy in parameter estimation. We then used the individual-level evidence approximated by the NHI (S1 Text) to subsequently perform random effects model comparison using the procedure introduced by Stephan et al. [13,19]. This means that whereas the NHI procedure assumes no hierarchy across parameters, it does (via the Stephan procedure [13]) allow for a hierarchical structure over model identity. In contrast, the HPE procedure, as introduced by Huys et al. [10,11], assumes a hierarchy over parameters, but no hierarchy over model identity: we accordingly, use it with a fixed-effects model comparison procedure. The HBI assumes that both parameters and model identities are generated hierarchically in turn. Note that related approximations, as similar as possible, have been used for making inference in these methods, which allows for a fair comparison (S1 Text) since our main points concern the statistical structure of the methods, not the estimation techniques. In particular, HPE builds upon NHI’s Bayesian inference of per-subject parameters to condition these on additional group level parameters, by using expectation-maximization [14]; and HBI extends that algorithm to condition these on an additional level of model identity variables, by using variational Bayes [16,17]. We also use the same (Laplace) approximation to marginalize the subject-level variables in all three methods. The HBI algorithm has been given in Materials and methods and details of implementing the NHI and HPE have been given in S1 Text. The details of simulation analyses and parameters used in simulations have also been given in S1 Text. The HBI is general and could be applied to any type of data, such as choice data, reaction times, physiological signals and neural data. Since we are primarily interested in models of choice data, we focus on decision-making experiments. We then applied the HBI method to an empirical choice dataset from 31 subjects performing the two-step Markov decision task. The data used for this analysis have been reported elsewhere [28]. In Fig 13, we have plotted protected exceedance probabilities of each model, model frequencies and estimated group means and corresponding hierarchical errors. According to this analysis, the hybrid model is the most likely model across the group. We also performed further analysis testing whether individual differences found by the HBI generalize to individual differences in conceptually related, yet independent, data. We reasoned that subjects showing a hybrid strategy might be slower in their choice, as the hybrid model requires combining of model-based and model-free values (which in some trials might be in conflict). Therefore, we looked at the median of response time across all first-level choices for each subject and tested whether there is a difference in response times between those subjects who (according to the separate analysis of choices) employed a hybrid strategy vs. those who employed a model-based strategy as estimated by the HBI. The subgroup attributed to the hybrid model by the HBI showed slower response time compared to those subjects attributed to the model-based account (P = 0. 03, Wilcoxon test). These results suggest that HBI reveals meaningful individual differences generalizing to unseen data. We applied the HBI to another choice dataset of Parkinson’s disease (PD) patients (N = 31), who performed a probabilistic reward and punishment learning task with binary choices (160 trials), which has been used previously for studying maladaptive learning in PD patients. All patients tested on medication. The dataset used here has been reported elsewhere [15]. Previous studies proposed that positive and negative prediction errors might be communicated through different dopaminergic receptors or striatal pathways [25,26,29], and thus the PD patients might have different learning rate parameters for learning from positive and negative prediction errors [29]. Therefore, we considered a model space including the RL model, the dual- RL model and a simple strategy that selects actions based on the most recent outcome. In both RL models, we also included a perseveration parameter, which models the tendency to repeat or avoid the same choice regardless of the value [15,30]. This analysis showed that the dual-α RL model was more likely across the group. Protected exceedance probabilities, model frequencies and estimated group means and corresponding hierarchical errors are plotted in Fig 14A. We then considered data from matched control participants (N = 20), who performed the same task. The analysis with the HBI showed that the RL model is more likely for the control group (Fig 14B), suggesting that PD (or dopaminergic medication in PD) increases the discrepancy between the learning rates for positive and negative prediction errors. We finally performed a permutation test to formally test the significance of this difference (1000 permutations). For each permutation, all participants were randomly divided into control and PD groups with the same size as the real control and PD groups. The HBI was then used to fit the same model space to each random group. The relative model frequency statistics (RL vs. dual-α RL) was quantified for each permutation. This permutation test confirmed that the dual-α RL was significantly more likely than the RL model in PD patients compared with controls (P <0. 001).
  In this work, we have introduced a novel method, a hierarchical and Bayesian inference framework, for parameter estimation and model comparison. The HBI framework is hierarchical in the sense that parameters at the individual level are regularized by statistics across all individuals in the group. The HBI framework is Bayesian in the sense that all uncertainties at both individual and group levels are represented by probability distributions. The HBI framework has major theoretical advantages over current state-of-the-art methods, mainly because it combines, in a single hierarchical structure, two sorts of inference (about model identity and model parameters), which are interdependent but have previously been treated separately. Our simulation results demonstrated these advantages experimentally. In this work, we took an empirical Bayes approach [31,32], in which priors are constructed based on data. In other words, parameters at the individual level are regularized by statistics across all individuals in the group. Furthermore, we took a so-called random effects approach to model identity [13], which indicates that different models might underlie data in different subjects. This is in contrast to previous hierarchical methods for model fitting, which assume the same model underlie data in all subjects (fixed effects assumption [10,11]). The random effects approach to hierarchical inference has important consequences for both parameter estimation and model comparison. Moreover, we took a fully Bayesian approach by quantifying uncertainty at the group level, which enabled us to develop statistical tests about group parameters and to quantify corresponding statistical errors. Empirical Bayes methods play an increasing role in modern statistics. These methods essentially take a hierarchical approach, by assuming that individual data are generated based on the probabilistic properties of the population. This hierarchical approach has important consequences. The most important consequence is that they provide a promising solution to the classical problem of priors in Bayesian statistics by providing informative, yet objective, priors at the individual level. Furthermore, by partly sharing parameters across subjects, they reduce overfitting relative to non-hierarchical models, which in turn allows them to confidently fit more complex models with a smaller penalty for overfitting. This is because non-hierarchical methods assume that the extra parameters of a complex model are independent. For example, consider a model space in which the more complex model has one extra free parameter and there are 40 subjects in the dataset. Fitting the dataset with the complex model using non-hierarchical methods introduces 40 additional independent free parameters, driving the danger of overfitting, and accordingly an excessive penalty to account for this possibility in assessing the evidence for the model. The hierarchical approach, however, assumes that the individual parameters are dependent, as they are all generated according to the same distribution, sharing a single mean parameter and smaller deviations from it. Modeling this hierarchical dependency enables those methods to avoid penalizing complex models as excessively. Our simulation results demonstrate this point experimentally (Fig 3D). While the non-hierarchical method failed to select the correct model with one additional parameter, evidently because the overfitting penalty was too extreme, the HBI was successful in selecting the correct model (Fig 3D). The HBI method introduced in this paper is built based on the random effects view that different models might underlie data in different subjects. Taking this view enabled us to address problems caused by taking the model identity as a fixed effect in some hierarchical parameter estimation procedures. For parameter estimation, the fixed effects assumption biases the group parameters because it assumes that all subjects contribute equally to the group parameters. The proposed HBI framework solves this problem by weighting contribution of each subject to group statistics by the degree to which that model is likely to be the true underlying model for that subject (Figs 1 and 3). For model comparison, the fixed effects assumption leads to oversensitivity to outliers as the evidence across the group is driven by the sum of individual evidence. Our simulation results (Fig 2) showed that only a few outliers lead to incorrect model selection inference made by the fixed effects assumption. The proposed HBI method solves this problem by normalizing individual evidence across all candidate models. Specifically, the HBI framework quantifies the responsibility of each model k in generating each subject data, a metric lying between 0 and 1. For every subject, the responsibility sums up to 1 across all candidate models as it partitions probability space among those models (see [13,19] for a similar non-hierarchical approach). It is then easy to compare models by enumerating responsibilities across the group in favor of each model or by estimating the most likely model. Another major contribution of this paper is to provide a statistical test, HBI t-test, to the inference problem at the group level using hierarchically fitted parameters. For models fitted by a non-hierarchical method, such as maximum likelihood or Laplace approximation, it is statistically valid to use classical statistical tests on fitted parameters to make inference at the group level. However, for datasets fitted by a hierarchical method in which the individual fits are regularized according to statistics of the group data, conventional statistical tests are not valid, because the parameter estimates are non-independent from subject to subject. Our fully Bayesian approach enabled us to address this issue. Our method provides an intuitive solution to this problem in the form of a t-statistic, in which all the group statistics are computed according to the estimated responsibilities of the corresponding model in generating each individual data. Thus, the HBI quantifies the uncertainty of the group parameters and thereby the corresponding hierarchical errors. Our analysis showed that the HBI performed better than both the NHI and HPE in detecting true effects and also that it was well calibrated, displaying the appropriate number of false positives when effects were absent. Therefore, the HBI framework enables researchers to make statistical claims about parameters at the group level. It is important, however, to note that the foundation of the HBI t-test is completely different from the classical t-test, as it is a Bayesian (in contrast to frequentist) test using posterior probabilities. In particular, this test is based on the posterior distribution of the statistics of interest (i. e. group mean) marginalized over all other parameters (e. g. group variance), which is given by a Student’s t-distribution (Eq 24). Statistically, the precise claim of the HBI t-test is that whether a specific point is outside of a credible interval, which is the interval that the group parameter value falls with a particular subjective probability. For example, if the HBI t-test indicates that a parameter is significantly different from 0 at P <0. 05, it means that 0 does not fall within the 95% credible interval. One important difference between Bayesian credible intervals and classical (frequentist) confidence intervals (used in classical Student’s t-test) is that Bayesian credible intervals depend on priors. However, since we used minimally-informative priors (statistically proper priors with very little effects on posteriors, see Materials and methods), the HBI t-test almost entirely depends on data. In fact, that is the reason that under the null, the HBI t-test generates p-values uniformly as shown by simulations (Fig 11). Notably, the same Student distribution can also be used to accept the null hypothesis for example using a “region of practical equivalence” procedure described by Kruschke [33]. It is also possible to employ the more common way and make inference in favor of the null hypothesis using model selection. In this case, one needs to perform a model selection between a model in which the group-level mean of the parameter of interest is fixed at the null value (null model) and compare that with a full HBI with no restriction (alternative model) using Bayes factor (i. e. difference log model evidence). In addition to model comparison, the HBI framework can also be used for model selection in situations where the goal is to select one of the models as the best model across the group. Exceedance probability is a metric proposed [13] to perform model selection using a random effects approach. An important revision of this metric called protected exceedance probability [19] also takes into account the null possibility that none of the models in model space is supported sufficiently by data, i. e. the differences in model evidence are due to chance. As the HBI framework treats model identity as a random effect, it is possible to compute exceedance and protected exceedance probabilities (Eqs 26–28). Note that if this procedure indicates that models’ ability to explain data are not different (i. e. their difference is likely to be due to chance), one cannot rely on estimated parameters, as those are also dependent on estimated model frequencies. In this situation, we recommend to obtain parameters by fitting models separately to data using the HBI, which makes sense as there is no evidence that models are differently expressed across subjects. In our analysis with simulated and empirical data, however, we never encountered this situation as probability of the null (P0 in Eq 28) was always very small. In this study, we compared the performance of the HBI with two alternative methods with different statistical assumptions about the generative process of data. The NHI assumes a hierarchy in model identity for generating individual data. The HPE assumes that parameters are generated in a hierarchical fashion, but assumes no hierarchy regarding model identities. The HBI assumes that both model identity and parameters are generated hierarchically. Importantly, the inference procedure for all these methods is very similar, which allows a fair comparison of them largely based on their statistical assumptions. In particular, the three methods all employ Laplace approximation for making a quadratic approximation of individual-level posteriors. Furthermore, the HBI is based on variational Bayes, which is an extension to the case of multiple latent variables of the expectation-maximization procedure used previously for implementing the HPE [10,11] (see also [34] for a variational implementation), which itself extends the one-level Bayesian inference of NHI. There are other ways to make an inference, for example using Markov chain Monte Carlo methods. Future studies should investigate the pros and cons of those methods, compared with the variational Bayes used here, for making inference in HBI. There are increasing efforts to exploit advances in computational modeling for understanding mental disorders [3–6]. Recent works, however, have started to tackle challenges related to quantifying uncertainty in diagnosis and also in the evaluation of treatment effects. For example, hierarchical unsupervised generative modeling, have used Monte-Carlo and variational methods to identify a cluster of subjects showing similar patterns of neural connectivity [35,36]. HBI also offers a promising solution by quantifying uncertainty in model attribution to individuals. Our simulation analyses showed that the HBI performs better than other alternatives in model attribution. This can help us to move towards better diagnosis and precise evaluation of different treatments [37]. In summary, the HBI framework proposed in this work rests on a hierarchical view of both hypothesis testing (i. e. model comparison) and parameter estimation for multi-subject studies and thus provides a generic framework for statistical inference. Moreover, the HBI framework runs fully automatically and it does not rely on hand tuning of parameters. Therefore, we expect this method to be useful for a wide range of studies testing different hypotheses in a multi-subject setting. This includes not only computational models of learning and decision making but also any statistical models of brain or behavior.
  We begin by describing the probabilistic model of the HBI. Consider an observed dataset X = {x1, …, xN} where xn is the dataset (e. g. choices) of nth subject and N indicates the number of subjects and a model-space including K candidate models, M1…MK. Moreover, suppose that the prior probability of each model in the population is given by m = {m1, …, mK}. For each dataset, xn, we assume that there is a latent variable zn comprising a 1-of-K binary random vector, in which zkn is one if xn generated is by the kth model. Thus, the probability of the latent variable across all subjects, Z = {z1, …, zN}, is assumed to have a multinomial distribution, p (Z | m) = ∏ n ∏ k m k z k n. (1) Each model Mk in the model-space is supposed to compute the probability of a given dataset (e. g. a set of choices) given a set of parameters, hkn. For example, the reinforcement learning model computes the probability of choices using two parameters: a learning rate and a decision noise parameter. The number of models and their structures depend on specific scientific questions. Here, we take a general approach by making no specific assumption about the number of models, K. Thus, the kth model in the model-space, Mk, computes the probability of dataset xn given the parameter vector hkn, which is denoted by p (xn|hkn, Mk). Note that the number of parameters in model k, denoted by Dk, might be different across models. Since data for each subject is generated by one of the models, which is denoted in the binary vector zn, the probability of the observed dataset given the model-space is p (X | H, Z) = ∏ k ∏ n p (x n | h k n, M k) z k n, (2) where H denotes all the parameters across all participants and models. The parameters of kth model are assumed to have a multivariate normal distribution with mean μk and precision matrix Tk, p (H | Z, μ, T) = ∏ k ∏ n N (h k n | μ k, T k − 1) z k n, (3) where Tk is a diagonal matrix with positive elements. We also introduce a distribution over model frequencies, m. We use the Dirichlet distribution, which forms the conjugate prior for the multinomial distribution, as the prior: p (m) = Dir (m | α 0) = C (α 0) ∏ k = 1 K m k α 0 − 1, (4) where C (α0) is the normalizing constant for the Dirichlet distribution. We also take group parameters μ and T as random variables, which allows us to evaluate their posterior distribution given data. We introduce conjugate priors for these variables, a Gaussian-Gamma prior in which the distribution over μk depends on Tk: p (μ | T) = ∏ k = 1 K N (μ k | a 0, (b T k) − 1) p (T) = ∏ k = 1 K ∏ i = 1 D k G (τ k i | v, s), where G (.) denotes Gamma distribution. Here, τki is the ith diagonal element of Tk. Assuming that τk is a vector containing τki, by defining Tk = diag (τk), in which diag (.) is an operator outputting a diagonal matrix with elements given by τk, we can write these two equations in a compact form: p (μ, τ) = ∏ k = 1 K N (μ k | a 0, diag (b τ k) − 1) G (τ k | v, s), (5) where we have defined: G (τ k | v, s) = ∏ i = 1 D k G (τ k i | v, s), in which v is a scalar and s is a vector with Dk elements all equal to s. The full probabilistic model is given by, p (X, H, Z, μ, τ, m) = p (X | H, Z) p (H | Z, μ, τ) p (Z | m) p (μ | τ) p (τ) p (m). (6) The task of Bayesian inference is to compute the posterior probabilities of latent variables given data, p (H, Z, μ, τ, m|X). Since the inference is intractable for the probabilistic model outlined in the previous section, we employ variational inference to compute approximate posteriors. We take a so-called mean-field approach [16,17] by assuming that the posterior is partially factorized as follows: q (H, Z, μ, τ, m) = q (H, Z) q (μ, τ, m). (7) Note that we force no factorization in the posterior between latent variables, Z and H. Using a quadratic approximation of the conditional posterior, q (H|Z), we prove in S1 Appendix that these posteriors are given by, q (H, Z) = ∏ k ∏ n r k n z k n N (h k n | θ k n, A k n − 1) z k n (8) q (μ, τ, m) = Dir (m | α) ∏ k q (μ k, τ k) (9) q (μ k, τ k) = N (μ k | a k, diag (β k τ k) − 1) G (τ k | ν k, σ k), (10) where 0 ≤ rkn ≤ 1 is the responsibility of model k for nth subject, θkn and Akn are the subject-level mean and precision, νk and βk are scalars and σk is a vector with the same size as τk. In the next section, we provide the HBI algorithm, which iteratively updates the parameters of these distributions, rkn, θkn, Akn, α, ak, νk, βk, and σk. After initializing the individual parameter estimates, θkn and Akn and responsibilities rkn for all subjects and models, as well as setting prior parameters a0, b, s, v and α0 (which will be defined later), the HBI algorithm performs these steps: Calculate the summary statistics: N ¯ k = ∑ n r k n (11) θ ¯ k = 1 N¯ k ∑ n r k n θ k n (12) V ¯ k = 1 N¯ k ∑ n r k n (θ k n θ k n ⊤ − θ ¯ k θ ¯ k ⊤ + A k n − 1). (13) Update parameters of q (μ, τ, m) for all models: a k = 1 N ¯ k + b (N¯ k θ ¯ k + b a 0) (14) β k = b + N¯ k (15) σ k = s + 1 2 diag (N ¯ k V ¯ k + b N ¯ k b + N ¯ k (θ ¯ k − a 0) (θ ¯ k − a 0) ⊤) (16) ν k = v + 1 2 N¯ k (17) α k = α 0 + N¯ k. (18) Update the individual posterior parameters θkn, Akn and fkn, by obtaining a quadratic approximation of the function, ℓkn (h), with respect to h: ℓ k n (h) = p (x n | h, M k) N (h | E [ μ k ], E [ T k ] − 1), (19) where E [ μ k ] = a k and E [ T k ] − 1 = 1 ν k diag (σ k). This approximation can be written as ℓ k n (h) ≃ f k n exp (− 1 2 (h k n − θ k n) ⊤ A k n (h k n − θ k n) ). (20) Note that any quadratic approximation can be used here. For example, using a Laplace quadratic approximation (which is a very common approximation for analyzing behavioral and neural data [20,22,23]), θkn, Akn and fkn are given by the mode, Hessian of log ℓkn and the maximum value of ℓkn, respectively: θ k n = arg max hlogℓ k n (h) A k n = − ∇ ∇ logℓ k n (h) | θ k n f k n = ℓ k n (θ k n). Update responsibilities, r k n = ρ k n ∑ j = 1 K ρ j n, (21) where logρ k n = logf k n + 1 2 D klog2 π − 1 2 log| A k n | + λ k + E [ logm k ] (22) λ k = D k 2 (ψ (ν k) − logν k − 1 β k) (23) E [ logm k ] = ψ (α k) − ψ (∑ k = 1 K α k), in which ψ (.) is the digamma function. Terminate if stopping criteria are met, otherwise go to 1. An important goal of computational modeling studies is to compute the distribution of parameters given data across the whole population. From a Bayesian viewpoint, this is given by the marginal posterior over the mean of group parameters, μk, which reads p (μ k | X) ≃ ∫ q (μ k, τ k) d τ k = ∫ N (μ k | a k, (β k τ k) − 1) G (τ k | ν k, σ k) d τ k = S t (μ k | a k, η k, n k), where n k = 2 ν k = 2 v + N¯ k is the number of degrees of freedom of the Student distribution and η k = ν k β k σ k − 1 is the inverse-scale parameter. Therefore, the random variable t = η k 1 2 (μ k − a k) takes a form of standard Student distribution with nk degrees of freedom. By defining s k i 2 = 2 β k σ k i, in which s k i 2 corresponds to empirical variance (c. f. Eq (16) ), we can write this result in an intuitive form, p (μ k i | X) = S t (μ k i − a k i s k i / n k | n k). (24) Noting the similarity between s k i / n k and the standard error of the mean, we called s k i / n k the hierarchical error. Note that if we assume v = 0. 5 (which is reasonable as explained later), we obtain n k = 1 + N¯ k. In many situations, researchers are interested to fit a new dataset to a particular model and find corresponding parameters. In Bayesian statistics, this is called the predictive distribution and it is given by marginalizing over group parameters. Suppose that x* and h k * denote the new dataset and its corresponding parameters for model k. The marginal distribution p (x *, h k * | z k * = 1, X) is the predictive distribution given the observed dataset X assuming that the new data is generated by the kth model. This distribution is given by: p (x *, h k * | z k * = 1, X) = ∫ p (x * | h k *, M k) p (h k * | μ k, τ k, z k * = 1) p (μ k, τ k | X) d μ k d τ k = p (x * | h k *, M k) S t (h k * | a k, (1 + β k) − 1 η k, n k), where ηk and nk have been defined in the previous section. This distribution can also be written in terms of standard Student distribution with nk degrees of freedom. Furthermore, if we assume that b = 2v, which is a reasonable assumption (see the next section), this distribution is given by p (x *, h k * | z k * = 1, X) = p (x * | h k *, M k) S t (diag (s k) − 1 (h k * − a k) | n k), where sk is a vector of corresponding empirical deviance parameters, defined in the previous section. Using this joint distribution, one can use sampling methods to obtain the posterior over parameters, p (h k * | z k n = 1, X, x *), or to obtain the maximum-a-posteriori parameters, θ k *, given by θ k * = arg max hlogp (x * | h, M k) S t (diag (s k) − 1 (h − a k) | n k). (25) Note that for many degrees of freedom due to large values of N¯ k, the Student distribution tends to a Gaussian with mean ak and deviance matrix diag (sk). However, small values of N¯ k lead to a small number of degrees of freedom and heavier tailed distributions than Gaussians, which are more robust against outliers. As the mean-field variational inference is an iterative framework, it also depends on the initialization of the parameters. In this section, we provide priors that do not bias the final solution and also provide some intuitive criteria for the initialization. We initialize the parameters θkn and Akn by fitting all models separately to all participants (with some initial Gaussian prior), i. e. , assuming as if zkn = 1. These values are then used to calculate summary statistics according to Eqs (11) – (13). Furthermore, we need to define prior parameters. The free parameter α0 indicates prior frequency of each model. We take uninformative priors on frequency of models, which is given by α0 = 1 for all models. The prior mean, a0k, is assumed to be zero. Given Eq (15), we see that b can be interpreted as the effective number of prior samples associated with models. Also, given Eq (17), v could be interpreted as the half of the effective number of prior samples associated with models. Assuming that the priors account for one sample, which is a common assumption in Bayesian statistics, we take b = 1 and v = 1 2. Finally, since s has always an additive effect on σk according to Eq (16), we assume a small positive value for s, allowing that σk to be driven dominantly by data. In all our analyses, we assumed s = 0. 01. It is also important to note that by choosing a small value for s, we ensure that if a model loses entirely (takes no responsibility), its corresponding parameters at the individual level converge to the prior mean, a0k, with a very small variance. Finally, the HBI algorithm presented above requires stopping criteria. In our analyses, we terminated the algorithm if the change in normalized value of parameters between two consecutive iterations, j − 1 and j, defined as d ^ = 1 K ∑ k 1 D k ∑ i (θ ^ k i j − θ ^ k i j − 1) 2, was smaller than 0. 01. Here, θ ^ k i j is defined according to summary statistics of parameters on the jth iteration: θ ^ k i j = θ ¯ k i / V ¯ k i 1 2, where θ ¯ k i and V ¯ k i are the ith element of θ ¯ k and V ¯ k defined in (12 and 13), respectively. In our analyses, we also set 50 as the maximum number of iterations, although almost always the algorithm stopped before hitting this number. Using the posterior over m, one can also derive the so-called exceedance probability and protected exceedance probability, as defined in previous works [13,19]. We reproduce the equations here for completeness. The exceedance probability of kth model, ϕk, is defined as the probability that model Mk is more likely than any other model in the model-space and it is given by ϕ k = Prob (m k > m j | α), ∀ j ≠ k. (26) Computing protected exceedance probabilities, as defined in [19], also requires to run the HBI under the (prior) null hypothesis, H0, that there is no difference between models (i. e. α0 → ∞). The alternative hypothesis, H1, is the original case, in which α0 = 1. If we define L and L0 as the log-likelihood (actually the variational lower bound as its approximation) of all data given the model-space under H1 and H0, respectively, then the protected exceedance probability of kth model, ϕ ˜ k, is defined as: ϕ ˜ k = ϕ k (1 − P 0) + 1 K P 0, (27) where P 0 = 1 1 + exp (L − L 0). (28) Note that if P0 is close to 1, then model frequencies should be ignored, as the difference between models in the model space is due to chance. Furthermore, if data does not support any model, i. e. P0 is close to 1, then parameters should be estimated by fitting each model separately using the HBI.
summary: Computational modeling of brain and behavior plays an important role in modern neuroscience research. By deconstructing mechanisms of behavior and quantifying parameters of interest, computational modeling helps researchers to study brain-behavior mechanisms. In neuroscience studies, a dataset includes a number of samples, and often the question of interest is to characterize parameters of interest in a population: Do patients with attention-deficit hyperactive disorders exhibit lower learning rate than the general population? Do cognitive enhancers, such as Ritalin, enhance parameters influencing decision making? The success of these efforts heavily depends on statistical methods making inference about validity and robustness of estimated parameters, as well as generalizability of computational models. In this work, we present a novel method, hierarchical Bayesian inference, for concurrent model comparison, parameter estimation and inference at the population level. We show, both theoretically and experimentally, that our approach has important advantages over previous methods. The proposed method has implications for computational modeling research in group studies across many areas of psychology, neuroscience, and psychiatry.
keywords: learning
  medicine and health sciences
  decision making
  neurodegenerative diseases
  applied mathematics
  social sciences
  neuroscience
  learning and memory
  simulation and modeling
  algorithms
  cognitive psychology
  mathematics
  probability distribution
  cognition
  research and analysis methods
  statistical distributions
  kalman filter
  probability theory
  movement disorders
  psychology
  normal distribution
  neurology
  parkinson disease
  biology and life sciences
  physical sciences
  cognitive science
title: Hierarchical Bayesian inference for concurrent model fitting and comparison for group studies