Neurocomputational Substrates Underlying the Effect of Identifiability on Third-Party Punishment

Linear mixed model

Participants' TPP decisions were analyzed using a linear mixed model implemented in the R statistical package (version 3.6.1) with the lmerTest package (Kuznetsova et al., 2017). Specifically, participants' TPP decisions were modeled with within-subjects variables of Inequity (computed as differences between dictators and recipients according to the dictators' offers) and Identifiability (identified vs unidentified) as fixed variables and subjects as a random variable with varying intercepts and slopes. p values were determined by likelihood ratio tests between the full model with the effect in question against the model without the effect in question (Brown, 2021).

Computational modeling

To uncover the cognitive underpinnings of participants' TPP decisions, we modeled their choice behaviors within the framework of the inequity aversion model (Fehr and Schmidt, 1999). This widely used model assumes that participants consider both monetary payoffs and inequity when making decisions. In the current study, we used different variants of the inequity aversion model to formally test hypotheses regarding participants' perspectives (egocentric, impartial, or vicarious) in the TPP task (Zhong et al., 2016), as well as potential mechanisms (reference points vs sensitivity to inequity) by which identifiability influenced TPP decisions (Wright et al., 2011; Zhong et al., 2016).

First, an egocentric model (M1) of inequity aversion assumes that people are averse to outcome inequity between self and others only. In this model, consider 3 players indicated by i ∈ [1, 2, 3] for P1 the dictator, P2 the recipient, and P3 the third party (i.e., participants); let d indicate the corresponding TPP decision made, M = [m₁, m₂, m₃] denotes the vector of monetary payoffs of all 3 players, and θ represents a set of parameters in each model. The utility (U) of choice could be represented as follows (Fehr and Schmidt, 1999): U(d;M,θ)=(m3−d)−α·(|(m1−2d)−(m3−d)|+ |m2−(m3−d)|)(1) Where α (range, 0-10) denotes the extent to which one dislikes that others have different monetary payoffs than self. It should be noted that the current design did not allow for the distinctions between advantaged and disadvantaged inequality, because disadvantaged inequality between self and dictators and advantaged inequaltiy between self and recipients were intermixed in the current design. As such, the parameter α reflects general aversion to inequity between self and others.

Second, an impartial model of inequity aversion assumes that third party takes an impartial perspective to reduce general inequity between P1 and P2 with the consideration of his own payoff (Zhong et al., 2016). In other words, a participant cares about his own payoff as well as is averse to general inequity between P1 and P2 but is not directly involved in inequity comparison with P1 or P2 as assumed in the egocentric model. According to this model, the utility of choice could be represented as follows: U(d;M,θ)=(m3−d)−δ·|(m1−2d)−m2|(2) where δ (range, 0-10) represents the extent to which participants dislike the distributional inequity between P1 and P2. Higher δ values indicate a greater sensitivity to distributional inequity between others.

Two variants of the impartial model were used to investigate two potential mechanisms by which identifiability might influence TPP decisions. First, the scaling impartial model (M2) assumes that participants assign different weights to inequity based on identifiability. In other words, identifiability is hypothesized to impact TPP decisions by affecting individuals' sensitivity to inequity (e.g., Wright et al., 2011). In the scaling model, the utility of choice takes the following forms for identified and unidentified contexts, respectively: Ui(d;M,θ)=(m3−d)−δi·|(m1−2d)−m2|(3) Uu(d;M,θ)=(m3−d)−δu·|(m1−2d)−m2|(4) where δ_i represents the sensitivity to inequity in the identified context, whereas δ_u represents the sensitivity to inequity in the unidentified context.

Alternatively, extensive research has documented the influential role of reference points in human decision-making, with reference-dependent preferences being a fundamental concept in prospect theory (Tversky and Kahneman, 1991). In light of this, the additive impartial model (M3) posits that identifiability modifies TPP decisions by modulating reference points of inequity. In other words, participants may hold different prior expectations regarding typical allocations by dictators in identified and unidentified contexts, which serve as reference points for assessing inequity in those respective contexts (e.g., Sanfey, 2009; Chang and Sanfey, 2013). A notable distinction between scaling and additive models lies in the relationship between the effects of identifiability and the level of inequity. In the scaling model, the effects of identifiability increase proportionally with higher levels of inequity. However, in the additive model, the impact of identifiability on punishment remains constant across different levels of inequity as long as they exceed the reference points (Jenkins et al., 2018). Within the additive model, context-specific reference points of inequity are modeled as additive shifting parameters (Zhong et al., 2016; Jenkins et al., 2018), as opposed to being fixed at an equal split between dictators and victims (i.e., 0). Accordingly, the utility of choice takes the following forms: Ui(d;M,θ)=(m3−d)−δ·|(m1−2d)−m2−ηi|(5) Uu(d;M,θ)=(m3−d)−δ·|(m1−2d)−m2−ηu|(6) where η_i and η_u (range, 0-12) represent the reference points of inequity in the identified and unidentified contexts, respectively. In other words, higher η values indicate that participants adhere to a less stringent equity norm. For instance, a η value of 0 indicates that participants consider an equal split as the typical offer, while a η value of 12 (i.e., total amounts of jointly earned money) indicates that participants consider dictators taking all as typical behavior. Negative values of η are considered implausible, as they would assume that participants expected dictators to allocate more to the recipients than to dictators themselves. The reference point parameter has been well validated in previous studies, such that this parameter captures the experimentally manipulated changes in expectations of allocation (Herz and Taubinsky, 2013; Xiang et al., 2013).

Third, a vicarious model of inequity aversion assumes that participants take a vicarious perspective with P2 and only disliked the disadvantageous inequity of P2 with P1. Accordingly, third-party's TPP decisions are driven by reducing the vicarious disadvantageous inequity of P2 with the consideration of his own payoff. Similar to impartial models, the vicarious model assumes that participants are not directly involved in inequity comparison with P1 or P2. A crucial distinction between impartial and vicarious models lies in participants' attitude toward the condition where P2 has more than P1. In impartial models, participants dislike both P1 having more than P2 and P2 having more than P1; hence, any disparities between P1 and P2 result in disutility. In vicarious models, however, participants only object to P1 having more than P2 and remain indifferent to P2 having more than P1. This distinction is behaviorally significant since participants' TPP decisions could lead to a scenario where P2 has more than P1.

Similar to impartial models, two variants of vicarious model were constructed to examine whether identifiability modulated sensitivity to inequity or reference points of inequity. In the scaling vicarious model (M4), it is assumed that identifiability affects TPP decisions through impacting sensitivity to inequity, and the utility of choice takes the following forms: Ui(d;M,θ)=(m3−d)−δi·max((m1−2d)−m2,0)(7) Uu(d;M,θ)=(m3−d)−δu·max((m1−2d)−m2,0)(8)

In the additive vicarious model (M5), it is assumed that identifiability affects TPP decisions through impacting reference points of inequity, and the utility of choice takes the following forms: Ui(d;M,θ)=(m3−d)−δ·max((m1−2d)−m2−ηi,0)(9) Uu(d;M,θ)=(m3−d)−δ·max((m1−2d)−m2−ηu,0)(10)

Finally, based on the utility of each choice, we used the softmax function to evaluate the probability (P) of each choice as follows: P(d;M,θ)=eλ·U(d;M,θ)∑j∈Jeλ·U(dj;M,θ)(11) where λ (range, 0-10) is the softmax inverse temperature capturing the stochasticity in decision-making. The smaller λ is, the more randomly participants behaved, and d^j denotes all the possible TPP decisions (i.e., 0, 2, 4, 6 MUs) that could be implemented in each trial.

Parameter estimation

Parameter estimation was conducted using hierarchical Bayesian analysis, which enables more stable and precise estimation compared with maximum likelihood estimation (Ahn et al., 2013). The hierarchical Bayesian analysis was implemented with the RStan package in R. The RStan package depends on an interface to R provided by Stan (Carpenter et al., 2017), which adopts Markov Chain Monte Carlo sampling methods to execute fully Bayesian inference (see also Ahn et al., 2017). Following the approach in the hBayesDM package (Ahn et al., 2017), it was assumed that individual-level parameters were drawn from group-level normal distribution Normal (μ₀, σ₀) where μ₀ and σ₀ refer to group-level mean and SD, in other words, the hyper parameters of individual-level parameters. Weakly informative priors were adopted for the priors of the group-level normal means and SDs: μ₀, Normal (0, 10); and σ₀, half-Cauchy (0, 5). This was to minimize the impact of priors on the posterior distributions. It should be noted that closely similar results of parameter estimates were revealed with an even more weakly informative hyper-priors (μ₀, Normal (0, 100) and σ₀, half-Cauchy (0, 10)). The ranges of parameters were set based on theoretical/empirical plausibility as well as previous studies using similar models. For the parameters of inequity aversion (α or δ) and the softmax inverse temperature (λ), the value 0, as the lower bound, was set based on theoretical plausibility. Negative values of inequity aversion parameters are implausible, as they would assume that participants preferred inequity between self and others or between others in the framework of inequity aversion model. Likewise, negative values of softmax inverse temperature are implausible, as they would assume that participants preferred options with lower utility than those with higher utility. The upper bound (i.e., 10) of these parameters was set based on previous studies using similar models, which have either explicitly set similar ranges or empirically demonstrated that parameter estimates are within the range of 0-10 (Wright et al., 2011; Xiang et al., 2013; Zhong et al., 2016). In accordance, our data also showed that all α/δ/λ estimates are <10 (see also Results). For the reference points of inequity (η), we have mentioned above that the range between 0 (participants expected the dictator to equally split the money units) and 12 (participants expected the dictator(s) take all money units) includes all plausible values. The constraint was implemented with inverse probit transform following the approach of the hBayesDM package (Ahn et al., 2017). Models were defined in a Stan file, which was then compiled in the R environment. We fit each model with three Markov Chain Monte Carlo chains, each of which contains 2000 iterations after 2000 iterations for warm-up. Parameter estimates of the winning model (i.e., M5, the additive vicarious model, see also Results) were used for the comparisons between contexts, as well as correlational analyses with self-reported ratings, personality traits, and neuroimaging measures. The details of these analyses are illustrated in Materials and Methods and Results.

Model comparison

The leave-one-out information criterion (LOOIC) and widely applicable information criterion (WAIC) were used as the criteria for model comparison, as recommended in recent studies (Vehtari et al., 2017). Both LOOIC and WAIC represent the estimation of out-of-sample pointwise predictive accuracy using the posterior simulations. LOOIC uses leave-one-out cross-validation, and WAIC is based on the series expansion of leave-one-out cross-validation (Vehtari et al., 2017). By convention, a lower LOOIC or WAIC indicates better prediction accuracy of candidate models; thus, the model with the lowest WAIC and LOOIC is the winning one. In general, a 10 point difference in LOOIC or WAIC between two models can demonstrate superiority (Burnham and Anderson, 2004). The LOOIC and WAIC of all candidate models were computed using the “loo” package in R.

Based on the differences in LOOIC/WAIC value between each candidate model and the winning model, it is also possible to obtain LOOIC/WAIC weights. These weights indicate the probability of being the best model for each model, given the data and the set of candidate models. Previous studies have provided detailed information on the calculation of LOOIC/WAIC weights (Wagenmakers and Farrell, 2004; Park et al., 2021).

Model validation

To further validate our model's performance, model validation was implemented. First, we correlated the average amounts of TPP between actual behavioral results and model predictions across participants. Second, we examined the correlation between the average amounts of TPP in actual behavioral results and model predictions across trials.

Parameter recovery

We ran parameter recovery analyses to assess how accurate a model estimates true parameter values from the simulated choice data generated from the true parameter values (Palminteri et al., 2017; Wilson and Collins, 2019). For this purpose, the estimated posterior means from actual data were used as the true parameter values, which provided plausible combinations of parameter values for parameter recovery analysis (see also Park et al., 2002). Afterward, simulated choice data were generated by using the true parameter values (for 27 subjects and 80 trials per subject). Finally, we used the simulated choice data to estimate the parameters, and then evaluated the model performance by calculating correlations between the true and predicted parameter values. The simulation was implemented 50 times for each participant. The means and SDs of the Pearson correlation coefficients between actual and simulated parameter estimates were reported.

Univariate analysis

Neuroimaging data analyses were performed with SPM 12 (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/). Preprocessing of functional data included slice timing, realignment through rigid-body registration to correct for head motion, normalization to MNI space, interpolation of voxel sizes to 2 × 2 × 2 mm and smoothing (8 mm FWHM kernel).

A two-level GLM was used to analyze the functional data. For the first level, boxcar regressors were defined for each participant and for each epoch of the time course. The regressors modeled the BOLD response to the epoch of both split phase (4 s) and response phase (3 s). Two regressors were created for the split phase: one for the identified context and the other for the unidentified context, with the levels of inequity (i.e., 0, 2, 8, 10, 12) as a parametric modulator. Another regressor was created for the response phase across contexts. These regressors were convolutions between respective boxcar stimulus function with the canonical hemodynamic impulse-response function (Büchel et al., 1998). The GLM also modeled first-order temporal autocorrelations in the residual. For the second level, parametric modulations of inequity were compared between identified and unidentified contexts with a paired t test to examine the interaction between Inequity and Identifiability. To investigate the neural correlates of individual differences in changes of reference points between the identified and unidentified contexts, as estimated from the winning model (additive vicarious model, M5, as described in Results), participant-specific parameter values were incorporated as a covariate at the second level (Wright et al., 2011). These parameter values were included in the analysis of the interaction contrast, which compared the effects of inequity between the identified and unidentified contexts. The main effects of Inequity and Identifiability were also assessed.

For false positive control, we used whole-brain cluster correction with a cluster-defining threshold of p < 0.001 and a family-wise error corrected threshold of p < 0.05 unless otherwise stated (Woo et al., 2014; Eklund et al., 2016). A small volume correction was also conducted within ROIs, focusing on those brain areas consistently implicated in costly punishment in light of recent meta-analyses (Gabay et al., 2014; Feng et al., 2015, 2021a). A small volume of interest was created to include these regions in both hemispheres according to the automated anatomic labeling (aal) atlas (Tzourio-Mazoyer et al., 2002): aal insula ROIs were used for insula; aal Cingulum_Ant ROIs were used for ACC; mPFC was constructed as a combination of the separate aal ROIs Frontal_Sup_Medial, Frontal_Med_Orb, and rectus; dlPFC was generated with aal ROIs Frontal_Mid and Frontal_Sup; aal caudate ROIs were used for caudate; and aal putamen ROIs were used for putamen (for similar ROI constructions, see also Strobel et al., 2011; Corradi-Dell'Acqua et al., 2013; Feng et al., 2017). Statistical significance was established through either whole-brain cluster-level inference or small volume correction and was stated explicitly in Results.

Multivoxel pattern analysis (MVPA)

To examine brain regions that were differentially involved in the dissociation of fair and unfair events in the identified and unidentified contexts, we performed the MVPA implemented in The Decoding Toolbox (Hebart et al., 2014). For each subject, we estimated a GLM in which trials of four conditions (i.e., fair offers in the identified context, unfair offers in the identified context, fair offers in the unidentified context, unfair offers in the unidentified context) were individually modeled. The estimated β maps were then fed to a support vector machine classifier with a fixed cost parameter (c = 1) and searchlight decoding scheme (Kriegeskorte et al., 2006; Corradi-Dell'Acqua et al., 2011). For each voxel in the individual participants' brain image in the native space, a sphere with a radius of four voxels was defined. For images associated with each condition in each run, the parameter estimates for each of the N voxels in a given sphere were then extracted to represent an N-dimensional pattern vector. Pattern vectors in one run (training set) were used to train the support vector machine to discriminate between fair and unfair offers separately for identified and unidentified contexts, and the performance of the classifier was measured on the other independent run (test set). The d′ was computed as a metric of the sensitivity of the classifier to discriminate fair and unfair offers in the test set (Corradi-Dell'Acqua et al., 2011, 2016). The procedure was repeated twice, with each run being the test set once. The mean d′ across two folds was calculated and assigned to the central voxel of the sphere. The classification was conducted for each voxel, resulting in a 3D d′ map for each context and each subject. At the group level, the classification performance was compared between identified and unidentified contexts, and the statistical significance was determined using a permutation test (N = 5000) implemented with the SnPM toolbox of SPM (Corradi-Dell'Acqua et al., 2011, 2014). The MVPA focused on classifying unfair and fair offers rather than decoding the magnitude of inequity, because of the relatively small and unbalanced number of trials available for each level of inequity.

The MVPA and univariate analyses share the use of mean activation, but MVPA also leverages voxel-level variability, even if a single voxel does not significantly change across conditions. In this study, both univariate and MVPAs were reported for several reasons. First, univariate analysis has been the standard and conventional approach in the literature. Reporting the results of univariate analyses allows for direct comparisons with previous findings and facilitates future meta-analytic studies. Second, MVPA and univariate approaches are sensitive to different variances, noise sources, and different types of effects (Davis et al., 2014). As a result, the two approaches can be used in a complementary manner to provide a more comprehensive understanding of the data (see also Weaverdyck et al., 2020).

Connectivity analysis

The activation and MVPAs identified the involvement of left vAI, dAI, and dACC in the identifiability effect (see also Results). We examined whether these regions worked together with other brain regions to underlie the identifiability effect, with an analysis of psychophysiological interaction (PPI) (Friston et al., 1997) using the revealed areas as ROIs. Specifically, we used the generalized PPI toolbox (http://www.nitrc.org/projects/gppi, version 13.1) (McLaren et al., 2012) with fMRI signal time courses individually extracted from a given ROI as the seeding signals. These seeding signals were then deconvolved with the canonical HRF, resulting in estimates of underlying neuronal activity. Subsequently, the interactions of these estimated neuronal time-series and vectors representing each of the onsets for each regressor and corresponding parametric modulators as in the univariate activation analysis were computed. Last, these interaction terms were reconvolved with the HRF and entered into a new GLM along with the vectors for the onsets of each event (i.e., the psychological terms). Group-level analysis of the PPI data was similar to that of the activation data, except that the β values used were derived from the PPI regressors.