Bayesian Workflow for Generative Modeling in Computational Psychiatry

Model space

The set of models was based on hierarchical Gaussian filtering (HGF; Mathys et al., 2011, 2014), a Bayesian learning model comprising a hierarchy of computation inspired by predictive coding theories of brain function (Friston, 2005; Rao & Ballard, 1999). In simple terms, HGF assumes that the brain implements hierarchical Bayesian belief updating to predict the current state of the world. The layers of the HGF represent this hierarchy and track the probability of events (also called tendency) and how quickly this probability changes (called volatility). HGF is a generative modelling framework that implements learning, i.e. belief updating, using precision-weighted prediction errors at multiple levels. Crucially, this hierarchy allows for dynamic adaptation of the learning rate. It describes the evolution of (hidden) states and how these give rise to the sensory inputs (u) an agent receives. We used a 3-level enhanced HGF (eHGF) for binary inputs for the current application, where the states evolve as Gaussian random walks (GRW) at all but the first level and the step-size of the GRW at any given level depends on the next higher level (variance-coupling). The probability of the binary states (stimulus categories) at the lowest level originates from a sigmoid transformation of the quantity at the second level. The state at the second level describes the tendency towards the two stimulus categories and the state at the third level captures the log-volatility of this tendency. The generative model of the eHGF is visualized as graphical model in the upper part of Figure 3 and described by the following set of equations

Figure 3

Further details and equations of the generative model can be found in Mathys et al. (2011).

The 3-level eHGF for binary inputs represents a concrete implementation of the meta-Bayesian ‘observing the observer’ framework, where an observation or response model serves to specify a mapping from inferred beliefs of an agent to observed responses as recorded during our experiment (Daunizeau et al., 2010a; Daunizeau et al., 2010b). The response model uses the perceptual model indirectly via its inversion (Mathys et al., 2014). In this study, we augmented the perceptual model (eHGF) with a novel set of response models combining binary responses and continuous-valued response times.

In the eHGF, variational inversion of the agent’s model of the world (perceptual model) using Variational Bayes (VB) under a mean-field approximation and a quadratic approximation to the variational energies give rise to a set of analytical trial-by-trial update equations that resemble the general structure of Reinforcement Learning (RL) models but explicitly represent inferential uncertainty. For example, the update of the belief at the second level of the eHGF about the cue-outcome contingency takes the form of

where superscripts refer to trial indices and subscripts to the level in the hierarchy of the model. In other words, the prediction (${\mu }_2^{(k)}$) at trial k is equal to the sum of its value at the previous trial ${\mu }_2^{(k-1)}$ and a learning rate ${\sigma }_2^{(k)}$ multiplied with a prediction error term ${\delta }_1^{(k)}$. Here the prediction error ${\delta }_1^{(k)}$ describes the mismatch between the actual ($u^{(k)}$) and the predicted sensory input (${\widehat{\mu }}_1^{(k)}$)

where the predicted sensory input is simply a sigmoid transformation of the belief at the second level

Notably, in the HGF, the learning rate ${\sigma }_2^{(k)}$ is dynamic and corresponds to a function of the temporally evolving uncertainty of the agent’s belief about the cue-outcome contingency ${\mu }_2^{(k)}$. A full account of the update equations describing the evolution of beliefs at all levels of the 3-level eHGF for binary inputs can be found in the analysis plan (Appendix A3).

We defined a set of combined response models that are all paired with a 3-level eHGF for binary inputs (perceptual model). The combination of binary responses and continuous response times in these response models is accomplished by summing the log-likelihood of two individual response models, assuming independence of the two response data modalities conditional on the parameters of the perceptual model. All seven models in our model space are listed in Table 1.

As binary part of the response model, we used a unit-square sigmoid function mapping from the belief ${\widehat{\mu }}_1^{(k)}$ at the lowest level of the eHGF that the next outcome will be 1 onto the probabilities $p\hspace{0.17em}(y_{\mathrm{bin}}^{(k)}=1)$ and $p\hspace{0.17em}(y_{\mathrm{bin}}^{(k)}=0)$ that the agent will choose response 1 or 0 (for simplicity, in the following equation we use $m^{(k)}\stackrel{\text{def}}{=}{\widehat{\mu }}_1^{(k)}$ and we omit time indices on y and m):

Here, the parameter ζ determines the steepness of the sigmoid function and is referred to as inverse decision temperature or inverse decision noise. For reasons of comparability, the binary part was kept constant across all response models in our model space.

For the continuous part of the response models, we used seven different variants of a general linear model (GLM) to predict the log-transformed RTs on a trial-by-trial basis. The regressors of models 1–4 (M1–M4) are belief trajectories from the perceptual model (eHGF), e.g., an individual’s estimated uncertainty about the outcome (${\widehat{\sigma }}_1^{(k)}$) at trial k. This was inspired by the approaches used in previous work where reaction time data had been modelled using the HGF framework (Lawson et al., 2017, 2021; Marshall et al., 2016). A set of alternative combined response models (M5–M7) predicted RTs independently of the perceptual model (eHGF); including these models allowed us to test whether informing RTs by inferred states from the eHGF would improve the models at all.

Our model space can thus be divided into 2 families based on the continuous part of the response models (since the perceptual model and the binary part of the response models were held constant). M1–M4 are ‘informed’ RT models and M5–M7 are ‘uninformed’ RT models. Note that Bayesian model comparison is based on the (log) evidence and therefore requires the data to be the same across all models. In our case, it is thus not possible to include a model that does not model RT at all in the model comparison.

The GLM equations of M1–M7 are listed in Table 1, and a graphical representation of M1 is presented in Figure 3. Concerning the first (informed) set of models, the GLM of M1 was adapted from a model used by Lawson et al. (2017) and tailored towards the specific properties of the SPIRL task. Models M2–M4 are custom built response models with the following underlying motivation: Firstly, the goal was to reduce the number of regressors entering the GLM and hence the complexity of the models compared to M1. Secondly, the idea was to include different types of uncertainty from different levels of the hierarchy described by the perceptual model (eHGF). The eHGF accommodates various forms of uncertainty. Two sources of uncertainty, informational and environmental uncertainty, respectively, are represented in the update equations of the eHGF. M2 includes estimates of informational uncertainty (${\sigma }_2^{(k)}$) and phasic volatility (${\widehat{\mu }}_3^{(k)}$). M3 includes estimates of informational uncertainty at the outcome level (${\widehat{\sigma }}_1^{(k)}$) and at the second level (${\sigma }_2^{(k)}$) as well as environmental ($e^{({\kappa }_2{\mu }_3^{(k-1)}+{\omega }_2)}$) uncertainty whereas M4 only includes estimates about the two latter quantities. For a more detailed discussion of the different forms of uncertainty captured by the eHGF, please see Mathys et al. (2014).

By contrast, the second (uninformed) set of models, M5–M7, have response time GLMs where the regressors are independent of the perceptual model. M5 is a null model, representing the assumption that log RTs are simply noise around a constant intercept. M6 allows for a linear decay over time in the response times, and M7 corresponds to a null model with two different intercepts for correct and incorrect response trials (referring to the binary response data).

Prior specification

Prior elicitation. Since the central questions of the present study are addressed by Bayesian model comparison (whose central quantity, the model evidence, takes into account the models’ priors), the specification of priors is a critical step. We aimed to specify priors that lay in the range of actual human behaviour for the present task, while at the same time minimising the subjectivity involved in the choice of priors. In order to avoid problems of double-dipping, we made use of an independent data set, i.e. the 20 participants from our pilot study, for the elicitation of prior densities. Recently, similar approaches have been used for comparable modelling endeavours (Harrison et al., 2021; Schöbi et al., 2021). We used a two-step procedure for elicitation of prior densities that were subsequently used in the analysis of the main data set (for the results, see Supplementary Material S2). We estimated the sufficient statistics (mean and variance) of the prior densities as follows (all priors take the form of a normal distribution):

1. Inversion of M1–M7 given the pilot data set (N_pilot = 20) using initial prior means and variances. These initial prior means and variances represent the default prior configurations of the eHGF in the HGF Toolbox.

2. For every model m, we defined:

   1. a new prior mean pE_m as the robust mean estimate over Maximum a posteriori (MAP) estimates obtained in step 1.

   2. a new prior variance pC_m as the robust variance estimate over MAP estimates obtained in step 1.

Due to the limited size of the pilot data set, robust mean and variance estimates were computed. We used a variant of the minimum covariance determinant (MCD) method, namely the FAST-MCD algorithm developed by (Rousseeuw & Van Driessen, 1999), as implemented in the ‘robustcov’ function in MATLAB R2019b. We refer to this new set of priors, that we estimated via the above described procedure, as empirical priors (as informed by the pilot data set). Please note that our procedure is distinct from commonly employed ‘empirical Bayesian’ procedures that estimate priors and parameter values from the same data set, using a hierarchical model. The advantage of our procedure is that the data which inform the choice of priors is fully independent from the data that inform parameter estimation.

Prior predictive checking. In order to get an intuition for the range of behaviours which can be generated by our models under the initial and empirical priors, we looked at their prior predictive distributions. More specifically, we randomly sampled 100 parameter values from the respective prior densities and simulated belief trajectories and responses. Adequacy of our chosen prior distributions was determined in a qualitative manner (visualisations) as well as in a quantitative fashion in the form of different recovery analyses (see paragraph ‘Validation of model inversion’). More detailed information describing our choice of priors and prior predictive checking can be found in the analysis plan.

Model inversion

We inverted the generative models using approximate Bayesian inference as implemented in the HGF Toolbox. MAP estimates are computed as the minimum of the negative log joint using gradient-based optimisation techniques. By default, the HGF Toolbox uses an optimisation algorithm from the quasi-Newton methods family (Broyden-Fletcher-Goldfarb-Shanno algorithm; BFGS) (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970). The covariance of the posterior was obtained under a Laplace approximation to the negative log joint at the MAP. This allowed for the specification of credible intervals on the obtained parameter estimates as well as the calculation of an approximate log model evidence (LME) as a measure of ‘model goodness’. For the present analyses, a multi-start optimisation approach was used to alleviate issues with the optimisation algorithm getting stuck in local extrema of the objective function. Specifically, 400 different starting points were used for the inversion of each subject and model. Of these 400 starting points, one always corresponded to the prior mean values of the different parameters (representing the default setting in the HGF Toolbox), whereas the other 399 starting points were randomly sampled values from the respective prior density of the parameters.

It should be noted that this choice of inference technique means that we are dealing with a simplified case of BDA. That is, we only obtain a point estimate of the posterior and approximate the posterior uncertainty. An advantage of our model inversion procedure is the computational efficiency and the avoidance of concerns about convergence, as in the case of Markov Chain Monte Carlo (MCMC) procedures. Having said this, for the specific purpose of our paper, the choice of model inversion approach is not critical.

Validation of model inversion

To examine the identifiability of parameters from the models in our model space, we performed a set of recovery analyses using synthetic data (for a tutorial introduction to recovery analyses, albeit in the context of maximum likelihood estimation, see Wilson and Collins (2019)). Examining the identifiability (or recoverability) of parameters and models can also be seen as a test of face validity, i.e. asking whether the models actually do what they are supposed to do: allowing for veridical parameter estimates and representing a distinct explanation of observed data that can be distinguished from other explanations. This step is an important part of Bayesian workflow because it establishes a boundary between the type of questions that can be addressed using the model space at hand, and those questions for which a meaningful answer cannot be expected. In other words, we assessed whether in principle (i.e., knowing the ‘true’ data-generating models and parameter values), we would be able to identify the data-generating model as the model that explained the data best (model identifiability) and the parameter values of the data generating model (parameter recovery), respectively. In a second step, we assessed whether data-generating model families could be recovered using family-level comparisons.

For each model, we generated a synthetic data set (N_sim = 100) by randomly drawing 100 samples from the empirical prior densities and plugging these values into the likelihood function. We then fit the synthetic data sets by each of the models in our model space under the respective empirical priors. Parameter recovery was assessed by visually comparing simulated parameter values to MAP estimates obtained using the data-generating model and by calculating Pearson correlation coefficients r. Model identifiability was quantified both as the proportion of correctly identified models according to approximate LME scores in a classification analysis, and calculating protected exceedance probabilities (PXP) as part of random-effects (RFX) Bayesian model selection (BMS) (Rigoux et al., 2014; Stephan et al., 2009). Moreover, we computed a balanced accuracy score for the LME winner classification and compared it to the upper bound of a 90%-CI when assuming each model to be selected with equal probability. Family level recovery was quantified using approximate LME scores obtained during the model identifiability analysis to create 5’000 synthetic data sets of N_famsim = 60 subjects each, with different ratios of data generating families. We compared true family frequencies with expected posterior family frequency (Ef) and family exceedance probabilities (XP) resulting from family-level RFX BMS (Penny et al., 2010). More details on the background of BMS methods can be found in the next section (‘Model comparison’) and in the referenced literature.

Model comparison

After careful investigation of our model space and validation of the chosen inference algorithm, we inverted the models given the main data set in order to address the following aims:

Aim 1. Perform a quantitative comparison of response models that utilise binary and continuous-valued response data in different ways. More specifically, we tested whether the family of informed RT models explained the collected data better than the uninformed RT model family, or vice-versa. To this and, we performed family-level RFX BMS as implemented in the VBA Toolbox. Family-level inference serves to reduce uncertainty about aspects of model structure other than the characteristic of interest (Penny et al., 2010). In the present study, we had no specific hypothesis about which of our candidate models explained the data best. Instead, we aimed to show that informing our response time models using quantities from the perceptual model (informed RT family) resulted in a better explanation of the data than using different variants of response time models that explain response times independent of the perceptual model (uninformed RT family). RFX BMS represents a hierarchical approach to model selection, treating the model as a random variable among subjects, allowing for inference on posterior family and model probabilities. RFX BMS accounts for group heterogeneity (different participants may be using different winning models/families) and provides robustness against outliers as opposed to fixed-effects (FFX) procedures (Stephan et al., 2009).

For the family-level RFX BMS, we specified a uniform prior at the family level to avoid biasing our inference (Penny et al., 2010). We computed the posterior family probabilities which correspond to the posterior belief that family k generated the data. Additionally, we computed family XPs which correspond to the belief that family k is more likely than any other of the K families, given the data from all participants.

Aim 2. A quantitative assessment of whether and how parameters characterising subject-specific learning behaviour were associated with an individual’s measured response times. To this end, we first needed to identify a single winning model – i.e. the model that best describes the measured data – by conducting a second RFX BMS analysis on the entire model space, this time at the level of individual models.

Once this model had been identified and in case it belonged to the family of informed RT models, we analysed the influence of the individual regressors of its response time GLM in a second step by examining the posterior parameter estimates of the winning model. To quantify the importance of each free parameter, we tested its significance against the initial prior mean (i.e. 0) using one-sample t-tests at a significance threshold of p < 0.05, Bonferroni-corrected for the number of performed tests (i.e. free parameters).

Model evaluation

After we had fit our models to empirical data and performed a set of statistical tests to answer our research questions, it still remained to be evaluated whether our model(s) actually provided a good explanation of the data, i.e. an assessment of model quality in absolute terms, using posterior predictive checks (van de Schoot et al., 2021). This approach is distinct from the purpose of (Bayesian) model comparison which assesses the quality of candidate models in relative terms.

Specifically, we performed qualitative posterior predictive checks (van de Schoot et al., 2021) by simulating new data conditional on the obtained subject-specific posteriors. For this, we used the Hessian of the negative log joint at the MAP estimate obtained during the optimisation to determine the posterior covariance matrix. This allowed us to construct an approximate multi-dimensional posterior for each model and participant. We then sampled parameter values (N_ppc = 100) from the calculated approximate posterior densities of each participant and model and used the sampled parameter values to generate synthetic response data (binary and continuous). Subsequently, the generated data were compared to the empirical data at an individual subject level. In order to determine the quality of the posterior predictions of our model, we pre-specified the following criteria in our analysis plan which we deemed critical characteristics of the collected data set: Calculated adjusted correctness of binary responses, and the plausibility of predicted log RT trajectories assessed in a qualitative fashion.