Authors: Brendan Williams, Lily FitzGibbon, Daniel Brady, Anastasia Christakou
Categories: Original Manuscript, Reliability, Test retest, Sample size, Reinforcement learning, Computational modelling, Reversal learning, Cognitive flexibility
Source: Behavior Research Methods
Authors: Brendan Williams, Lily FitzGibbon, Daniel Brady, Anastasia Christakou
Intraclass correlation coefficients (ICCs) are a commonly used metric in test–retest reliability research to assess a measure’s ability to quantify systematic between-subject differences. However, estimates of between-subject differences are also influenced by factors including within-subject variability, random errors, and measurement bias. Here, we use data collected from a large online sample (N = 150) to (1) quantify test–retest reliability of behavioural and computational measures of reversal learning using ICCs, and (2) use our dataset as the basis for a simulation study investigating the effects of sample size on variance component estimation and the association between estimates of variance components and ICC measures. In line with previously published work, we find reliable behavioural and computational measures of reversal learning, a commonly used assay of behavioural flexibility. Reliable estimates of between-subject, within-subject (across-session), and error variance components for behavioural and computational measures (with ± .05 precision and 80% confidence) required sample sizes ranging from 10 to over 300 (behavioural median N: between-subject = 167, within-subject = 34, error = 103; computational median N: between-subject = 68, within-subject = 20, error = 45). These sample sizes exceed those often used in reliability studies, suggesting that sample sizes larger than are commonly used for reliability studies (circa 30) are required to robustly estimate reliability of task performance measures. Additionally, we found that ICC estimates showed highly positive and highly negative correlations with between-subject and error variance components, respectively, as might be expected, which remained relatively stable across sample sizes. However, ICC estimates were weakly or not correlated with within-subject variance, providing evidence for the importance of variance decomposition for reliability studies.
The online version contains supplementary material available at 10.3758/s13428-025-02599-1.
The study of learning and decision-making processes in psychology and neuroscience research typically relies on the use of conditioning tasks to assay behaviour. These include instrumental conditioning tasks where subjects learn associations between actions and outcomes, such as pulling a slot machine lever and winning money, through experience. These associations either increase or decrease the likelihood of a given action being made in the future, depending on whether the outcome was rewarding or punishing. However, these associations also need to be updated flexibly for an individual to respond adaptively to dynamic and changing environments. Supporting this flexible updating of action–outcome associations relies on cognitive flexibility. Cognitive flexibility, broadly defined, is a complex competence that enables the maintenance of a given goal (such as winning money) while appropriately updating the required actions to achieve that goal (Dajani & Uddin, 2015), by identifying, selecting, and executing the optimal response strategy (Yu et al., 2019).
A commonly used measure of cognitive flexibility is the reversal learning task (Izquierdo et al., 2017). In this task, actions are associated with differing probabilities of rewards and/or punishments. These varying probabilities mean that different actions will, on average, have outcomes that make them more or less favourable than other actions. During the task, the probabilities associated with the available actions change such that previously favourable actions become less favourable, and vice versa. When this occurs, a cognitively flexible agent shifts action selection towards the newly favourable and away from the previously favourable outcomes to maintain goal-directed behaviour (assuming the agent’s goal is to maximise gains). Performance in the reversal learning task can be indexed using summary measures of behaviour, such as choice accuracy, reaction time, and the pattern of win–stay/lose–shift behaviour, or by deriving latent descriptors of performance using computational modelling. For these latter indices, reinforcement learning models can be fitted to choice behaviour to estimate parameters that describe features of behaviour, such as the rate of learning, or to what degree behaviour is driven by estimates of expected value.
The reliability of reversal learning behaviour has been previously reported in several papers using behavioural and computational modelling (Schaaf et al., 2023; Waltmann et al., 2022), in schizophrenia (Reddy et al., 2016), and using functional magnetic resonance imaging (Freyer et al., 2009). At the group level, reversal learning performance has been shown to produce reliable activation in regions of prefrontal and parietal cortices and the cingulate gyrus (Freyer et al., 2009), and reliable behavioural effects in both individuals with schizophrenia (Reddy et al., 2016) and healthy non-clinical populations (Schaaf et al., 2023; Waltmann et al., 2022). However, there appears to be less agreement in the current literature about the reliability of parameters derived using computational modelling.
Two recent studies are particularly pertinent to this issue, both exploring the test–retest reliability of computational model parameters, but drawing somewhat different conclusions. Waltmann et al. (2022) assessed the reliability of behavioural and computationally derived measures of performance during reversal learning. They used several approaches, including the calculation of intraclass correlation coefficients (ICCs), split-half reliability, variance decomposition, and simulation work. Behavioural measures of performance (e.g., reaction time, accuracy, stay behaviour, and perseveration) had good to excellent reliability (based on ICC and correlation coefficients). Waltmann et al. (2022) also used variance decomposition to demonstrate that some behavioural measures (such as accuracy and staying after losses) had high between-subject variability while others (such as reaction time measures) had high within-subject variability. Similarly, parameter estimates from the best-fitting computational model were found to have good to excellent reliability. Schaaf et al. (2023) assessed the reliability of behaviour and computational measures of performance during both a reversal learning and a two-armed bandit task. Using the same ICC coefficient interpretations as Waltmann et al. (2022), Schaaf et al. (2023) found fair reliability for accuracy and lose–shift behaviour and good reliability for win–stay behaviour in the two-armed bandit task. Additionally, Schaaf et al. (2023) found good reliability for accuracy and lose–shift behaviour and excellent reliability for win–stay behaviour in the reversal learning task. Schaaf et al. (2023) also found that parameter estimates showed good identifiability (meaning equivalent likelihoods cannot be produced by different sets of parameter estimates during fitting (Gershman, 2016)) for simulated behaviour but not for human subject data. Thus, this suggests that although computational models can recover parameters reliably when behaviour is truly stable (i.e., generated by an artificial agent), they struggle with behaviour from real subjects, which is variable and influenced by context. Indeed, Schaaf et al. (2023) also demonstrate that momentary mood can influence model parameters estimated from behaviour, with happiness and stress being associated with decreased and increased sensitivity to negative feedback, respectively, in the two-armed bandit task.
ICCs are widely used in reliability studies as an indicator of a method’s ability to measure systematic differences between subjects. However, a method’s ability to capture systematic differences between subjects is influenced by factors including within-subject variability, random errors, and measurement bias (Liljequist et al., 2019). Generally, an ICC is calculated by taking the ratio of between-subject variance and the total amount of variance for a given measure (McGraw & Wong, 1996). However, because an ICC is a ratio, the individual contributions of variance components to the overall coefficient cannot be accounted for. A second limitation for only using ICCs to determine reliability is that ICC calculations are more heavily influenced by increases in between-subject variance than session variance effects (Barnhart et al., 2007, 2016; Gorgolewski et al., 2013). Therefore, stable between-subject effects over time are essential to minimise disproportionate biasing of ICC estimates. Variance decomposition, by contrast, keeps the variance of a given measure in its composite parts; in the case of ICC(A,1), this includes measures of within-subject session variance, between-subject variance, and error variance.
Waltmann et al. (2022) used variance decomposition to demonstrate that some behavioural measures (such as accuracy and staying after losses) and some computational model parameters (learning rate) had high between-subject variability, while others (such as reaction time measures and reinforcement sensitivities for wins and losses) had high within-subject variability. Schaaf et al. (2023) did not use variance decomposition to assess the reliability of their computational modelling parameters, but they measured subjects’ mood at each time point and found some evidence that changes in mood can explain within-subject variability in some model parameters.
Important factors to consider when assessing the reliability of a variable are the precision and accuracy of the metric used to quantify its variability. One factor that influences precision and accuracy is sample size, with smaller sample sizes generally producing less precise estimates of reliability (Chen et al., 2023; Clarke & Wheaton, 2007; Maas & Hox, 2004, 2005; Neubauer et al., 2020). Previous work has also suggested that the precision of variance component estimates can be improved by increasing sample size (Paccagnella, 2011), and increasing sample sizes can reduce bias in estimates of reliability measures within subjects (Neubauer et al., 2020). Here, we present results on the reliability of reversal learning using a sample size larger than in previous reports (Freyer et al., 2009; Reddy et al., 2016; Schaaf et al., 2023; Waltmann et al., 2022).
To complement previous research, we firstly replicate the analytical pipeline of Waltmann et al. (2022) using data from a large sample of subjects (N = 150) on a similar reversal learning task. We then build upon previous research by investigating the effect of sample size on estimates of reliability using synthetic datasets based on the statistical properties of the collected “ground truth” data. To pre-empt our results, our replication analysis showed that behavioural variables and computational modelling parameters had similar reliability patterns to Waltmann et al. (2022) as assessed using ICCs. We then used simulation to investigate the effects of sample size on estimates of variance components, and the association between these variance components and ICC values. We generated synthetic datasets based on the underlying distributions and associations between sessions from our collected data. We generated 1000 synthetic datasets with a sample size of 300 for each task performance measure. For each synthetic dataset we calculated estimates of variance components and ICC measures at each sample size from 10 to 300. We then determined at what sample sizes the proportions of variance components stabilised in our synthetic dataset, based on our ground truth data. Our results show that the critical sample size for stable estimates of variance components for our task performance measures change as a function of the level of precision (half-width) and confidence (point of stability percentile). Stable estimates of variance components required sample sizes ranging from 10 to over 300 participants (behavioural median N: between subjects = 167, within subjects = 34, error = 103; computational median N: between subjects = 68, within subjects = 20, error = 45) with ± 0.05 precision and 80% confidence, which exceeds sample sizes typically observed in test–retest research. Thus, our results suggest that a larger than usual sample size (and/or data density) is required to accurately estimate variance components and, in turn, infer reliability.
The Prolific online recruitment platform (https://www.prolific.co/) was used to enrol eligible subjects (filters described in supplementary materials) in this study over two waves (wave August–September 2021; wave March–April 2022). A total of 251 subjects completed the experiment in the first part of the study. After the first phase, two subjects were excluded because they failed instructional attention check questions, and 14 were excluded because they failed nonsensical attention check questions or careless/insufficient effort (C/IE) responding checks (further described below). A total of 222 subjects completed the second part of the experiment. After the second phase, one subject was excluded because they failed instructional attention check questions, and eight were excluded because they failed nonsensical attention check questions or the C/IE responding checks. The mean interval between the two sessions was 12.53 days (SD = 1.89, range = 11.45–22.25 days). Lastly, we used a binomial test to identify 58 subjects performing at or below chance level in the reversal learning task in at least one session, and excluded them from further analyses (Zorowitz et al., 2023). Our final sample used for statistical analysis included 150 subjects (mean age = 35.450, SD = 13.30, range = 19–73, female = 97), and demographic information for these participants is summarised in Table 1. Table 1Demographic information for the participants included in this studyMeanSDMinMaxAge35.4513.301973Prolific approvals511.27442.50322672RangeFrequencyAge19–254326–395940–593960–758NA1SexMale53Female97StudentNo106Yes32NA12Prolific approvals0–10015101–20020201–30024301–40017401–50018501–75023751–1000161001–1500111501–200042001–250002501–30002Prolific rejections057142313225486253Country of residenceUnited Kingdom143Portugal3Greece3Poland1Employment statusFull-time74Not in paid work13Other10Part-time28Unemployed (and job-seeking)10Due to start a new job within the next month2NA13
Subjects who successfully completed the first phase of the study were reimbursed £1.75 for their time, framed as a basic pay rate of £1.25 plus a 50p bonus based on their performance in the reversal learning task. This was done to maximise the likelihood that subjects would remain focused while completing the task. Subjects who successfully completed the second phase of the study were reimbursed £2.50 for their time, framed as a basic rate of £1.25 plus a performance bonus of £1.25. This payment bonus was larger during the second phase of the study than the first to encourage subjects to complete both parts of the study. This study was approved by the research ethics committee of the University of Reading (2021–50-AC).
Subjects were presented with two visually distinguishable abstract stimuli that would appear randomly on screen, one left of centre, and one right of centre. Subjects selected one of these stimuli by pressing the ipsilateral arrow key on their keyboard. Subjects were given up to two seconds to make a valid choice response. After stimulus selection, subjects were presented with the outcome of their choice. Choice outcomes were either the gain or loss of a single point. Initially, one stimulus was randomly assigned as the “correct” stimulus. Selection of the correct stimulus meant subjects had a 75% chance of gaining a point and 25% chance of losing a point. The “incorrect” stimulus had the inverse probabilities for gains and losses. Outcomes for the correct and incorrect stimuli were pseudorandomised, so the assigned outcome probabilities were true over contiguous blocks of four trials. Subjects experienced nine reversals during the task, with each reversal occurring every 15 ± up to 3 trials (uniform distribution). At the point of reversal, the identity of the correct stimulus was changed, such that the correct stimulus became the incorrect stimulus and vice versa. If subjects did not make a valid choice response within the two-second time limit, then they were told they were too slow, and lost a single point. Subjects completed 150 trials of the reversal learning task. This task was created using the JavaScript library jsPsych (https://www.jspsych.org/) version 6.1.0, and custom JavaScript code (Fig. 1). The JavaScript code for running the task is available in the following GitHub (https://github.com/bwilliams96/SR_Online_task). Testing was hosted on the Gorilla online platform (https://gorilla.sc).Fig. 1Reversal learning task overview. Subjects completed 150 trials. On each trial they were presented with two abstract stimuli, one of which was associated with reward on 75% of trials, and the other on 25% of trials (true over consecutive blocks of four trials). After subjects made their choice, they received either reward (+ 1 point) or punishment (− 1 point). The assigned reward probabilities reversed nine times during the task, and this reversal occurred every 15 ± up to 3 trials, based on a uniform distribution
The reversal learning task included automated careless/insufficient effort (C/IE) responding checks that terminated the task prematurely if met. These conditions were (1) not making a valid choice over five consecutive trials, or (2) not making a valid choice for over 5% of the total number of trials. As part of the study, subjects completed the 12-item version of the Intolerance of Uncertainty Scale (Carleton et al., 2007) (data not reported here). To check for C/IE, we added two instructional attention check questions and one infrequency attention check question (see supplementary materials), following best practice guidelines for online research (Huang et al., 2015; Zorowitz et al., 2023). Subjects were made explicitly aware of the use of attentional checks, and that their responses to the Intolerance of Uncertainty Scale questions would not influence their bonus payment.
Behavioural measures of task performance were derived from choices made in the reversal learning task. Accuracy is the probability of selecting the “best” (most likely to be rewarding) choice on trial n, regardless of whether a reward was obtained or not. Perseveration is a measure of persistence in choosing the previously “best” choice after reversal, and is defined as the probability of selecting the “worst” (least likely to be rewarding) choice on trial n, after receiving two losses when making that choice following a reversal. Stay/switching behaviour is determined as the probability of making the same choice as on the previous trial, and is defined as the probability of staying both overall and after either a win or loss on the previous trial. Reaction time is defined as the amount of time the participant took to make a choice on trial n and, like staying behaviour, is calculated both overall and whether a win or loss was experienced on the previous trial. We also calculated the difference in reaction time following a win and loss (win RT − loss RT). Behavioural measures were calculated using (1) simple means (referred to as “mean”), estimated from mixed-effects logistic/linear (dependent on whether the variable was binary) regression models that either (2) calculated estimates for each session using separate regression models (referred to as “separate”), or (3) calculated estimates for each session using a single model which explicitly modelled the effect of session (referred to as “joint”).
Reinforcement learning models, which are commonly used for modelling reward learning tasks, were fitted to reversal learning behaviour, replicating the methods of Waltmann et al. (2022). We fit two families of models, differentiated by how the expected value determined choice. The first family used an inverse temperature parameter (β) in the softmax choice function to define choice stochasticity by determining the steepness of the softmax function. The second family used a reinforcement sensitivity parameter (ρ) to determine the maximum difference in expected values between choices, placing a lower bound on choice stochasticity. To determine the best-fitting model within each family, we fitted a range of models with combinations of parameters that are commonly used in the reinforcement learning literature. Models had either a single inverse temperature (β)/sensitivity parameter (ρ) or separate inverse temperature (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\beta }_{win/loss}
\usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\rho}_{win/loss} $$\end{document}ρwin/loss) for wins and losses. Expected values for actions were updated using prediction errors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left({\lambda }_{t}-{V}_{t}^{k}\right) $$\end{document}λt-Vtk, the difference in value of the actual (λ) and expected outcome (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ V $$\end{document}V) of action \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ k $$\end{document}k on trial \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ t $$\end{document}t (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {V}_{t}^{k} $$\end{document}Vtk). The rate of expected value updating was captured by the learning rate, with models having either a single learning rate (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \alpha $$\end{document}α) or separate learning rates for wins and losses (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\alpha }^{+/-} $$\end{document}α+/-, dual learning rate models). Models either updated the expected value of only the chosen action (single update models) or of both the chosen and unchosen actions (dual update models) using the inverse outcome for updating the unchosen action. The dual update models were fitted with and without a discount weight (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \kappa $$\end{document}κ) for the unchosen action. An in-depth explanation of model variants can be found in the supplementary materials. #### Model fitting Our model fitting replicated the approach taken by Waltmann et al. (2022). Briefly, each model described in the previous section was fitted to the data with one of three estimation methods, maximum likelihood (ML), maximum a posteriori estimation with uninformative priors (MAP0), and maximum a posteriori estimation with empirical priors inferred from the multivariate distribution of parameter estimates across subjects (expectation–maximisation [EM]). The best-fitting model for each family was determined using the integrated Bayesian information criterion (iBIC) from the EM fitting approach, with a lower iBIC indicating a better model fit (Huys et al., 2011, 2012; Waltmann et al., 2022). #### Synthetic dataset generation and reliability assessment We firstly replicated the ICC, correlational, split-half, and variance decomposition-based reliability assessment of behavioural and computational measures as reported by Waltmann et al. (2022) (for a detailed overview see supplementary materials and the original publication). We assessed the effect of sample size on components of between-subject, within-subject session, and error variance for our behavioural and computational modelling measures of task performance (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ X $$\end{document}X). To do this, we used a regression-based approach to generate statistically related and plausible synthetic two-session data based on the underlying statistical properties of our collected data (Fig. 2). Firstly, we regressed session 2 measures of task performance (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {X}_{2} $$\end{document}X2) onto session 1 measures (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {X}_{1} $$\end{document}X1), then extracted the residuals (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ e $$\end{document}e) from this regression for each participant. We used histograms to estimate the probability density function of our task performance measure in session 1 and the residuals from the regression model, using the Freedman–Diaconis estimator (a robust estimator that accounts for data variability and sample size) to determine the optimal number of bins. We used the Fitter Python package (Cokelaer et al., 2024) to identify the best-fitting probability distribution function and calculated its parameters (from the 80 available in SciPy using sum of squares error as the fitting metric) for our task performance measure in session 1 (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {X}_{1} $$\end{document}X1) and the residuals (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ e $$\end{document}e) from the regression model. We then drew *n* random samples from the probability distribution functions of task performance measure in session 1 (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {X}_{1synth} $$\end{document}X1synth), and the residuals from the regression model (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {e}_{synth} $$\end{document}esynth). Lastly, we took the intercept (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\beta }_{0} $$\end{document}β0) and scalar (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\beta }_{1} $$\end{document}β1) coefficients from regressing our task performance measure from session 2 onto session 1, and generated session 2 values using the following \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {X}_{2synth}={\beta }_{0}+{\beta }_{1}{X}_{1synth}+{e}_{synth} $$\end{document}X2synth=β0+β1X1synth+esynth. We used this regression-based approach to generate 1000 synthetic datasets with sample sizes of *n* = 300.Fig. 2Regression-based method for generating plausible, statistically related two-session synthetic data from our “ground truth” behavioural data For each synthetic dataset we iterated through sample sizes ranging from 10 to 300 and calculated variance components and ICC(A,1) values. We used the following formulae to calculate components of between-subject, within-subject session, and error variance. Between-subject variance was calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ Between=\frac{k}{n-1}\sum_{i=1}^{n}{\left({\overline{x} }_{i}-\overline{x }\right)}^{2} $$\end{document}Between=kn-1∑i=1nx¯i-x¯2, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ k $$\end{document}k is the number of sessions and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ n $$\end{document}n is the number of subjects. Within-subject session variance was calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ Within=\frac{n}{k-1}\sum_{j=1}^{k}{({\overline{x} }_{j}-\overline{x })}^{2} $$\end{document}Within=nk-1∑j=1k(x¯j-x¯)2. Error variance was calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ Error=\frac{SS-\left(n-1\right)MSr-\left(k-1\right)MSc}{(n-1)(k-1)} $$\end{document}Error=SS-n-1MSr-k-1MSc(n-1)(k-1), where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ SS $$\end{document}SS is the total sum of squares, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ MSr $$\end{document}MSr is the between-subject variance, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ MSc $$\end{document}MSc is the within-subject session variance. We then corrected within-subject session variance to control for varying sample sizes by multiplying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ Within $$\end{document}Within by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \frac{k}{n} $$\end{document}kn (as is the case when calculating ICC(A,1)), and normalised variance components so they summed to one. Lastly, we calculated ICC(A,1) for each synthetic dataset. This reliability assessment procedure was carried out using our three methods for estimating behavioural measures of task performance and the computational model parameter estimates from the best-fitting model. ICCs were interpreted following guidelines from Cicchetti (1994; but also see Gell et al., 2023, for an overview of how even a small change in reliability can affect accuracy), with ICC < 0.4; 0.4 ≤ ICC < 0.6; 0.6 ≤ ICC < 0.75; 0.75 ≤ ICC. To test whether the effect of sample size on variance component estimates changed as a function of between-session variance, we also generated synthetic datasets with noise added to the session 2 estimates. To ensure noise was equally scaled across measures, we first generated statistically related two-session data as above (Fig. 2), then *z*-scored the data from the two sessions (mean = 0, *SD* = 1) and added Gaussian noise to values from session 2 sampled from a normal distribution (mean = 0) with varying levels of noise (*SD* = [0.25, 0.5, 1, 1.5, 2]). Finally, values for both sessions were reverse *z*-scored. Using this procedure we generated 1000 independent synthetic datasets with sample sizes of *n* = 300 for each behavioural and computational modelling variable, for each noise level (*SD* = [0.25, 0.5, 1, 1.5, 2]). To determine the critical sample size at which variance component estimates stabilised, we determined the *point of stability* using the method described by Schönbrodt and Perugini (2013). The *point of stability* was determined based on two values. These were the parameters *g*, which was the ground truth of a given measure, and *w,* which was the half-width of a *corridor of stability* around *g*, (CoS = *g* ± *w*). Here, g was defined using variance component proportions calculated from our observed data collected from participants, and we had corridors of stability with half-widths *w* = [0.025, 0.05, 0.75, 0.1, 0.15, 0.2]. The *point of stability* was defined as the smallest sample size in a vector of values (in our case these were variance component estimates for sample sizes ranging from 10 to 300 for one iteration of our simulation) from which all subsequent sample sizes had values that fell within the defined *corridor of stability*. This procedure was performed for each synthetic dataset for each variance component, meaning a distribution of *points of stability* could be generated. Critical sample sizes were then calculated for a given percentile (here, we use 80th, 90th and 95th percentiles) of the distribution of *points of stability*(Fig. 3). If a vector of variance components never stabilised within the *corridor of stability* then the *point of stability* was defined as the maximum sample size (300), while a vector of variance components that never deviated outside the *corridor of stability* had the *point of stability* defined as the minimum sample size (10), as in Schönbrodt and Perugini (2013). We determined critical sample sizes for our synthetic datasets that were generated both with and without additional Gaussian noise (Figs. 4, 5 and 6).Fig. 3Actual (thick black line) and simulated (transparent grey lines, each line represents one iteration of the simulation) variance proportions for a synthetic dataset, where the true variance proportion (*g*) is equal to 0.5. For a given width (0.05), the *n*th percentile of the *point of stability* (PoS) can be calculated within a *corridor of stability* (CoS range = g ± w), with the *corridor of stability* half-width in this example equal to 0.05 (CoS range = 0.45–0.55). More explicitly, this is the sample size at which the *n*th percentile of simulations had variance proportions that no longer exceeded the range of the *corridor of stability*Fig. 4Distributions of between-subject variance proportions for simulated measures of behavioural performance generated using our regression-based approach. These data were generated using behavioural measures estimated using the “joint” regression model, which explicitly modelled the effect of session. The mean proportion of variance for each sample size (purple), 90th inter-percentile range (dark pink), and interquartile range (light pink) are overlaid on individual data points (blue). Distributions of simulated data generated using behavioural measures from the “separate” regression model and using simple means can be found in the supplementary figuresFig. 5Distributions of within-subject variance proportions for simulated measures of behavioural performance generated using our regression-based approach. These data were generated using behavioural measures estimated using the “joint” regression model, which explicitly modelled the effect of session. The mean proportion of variance for each sample size (purple), 90th inter-percentile range (dark pink), and interquartile range (light pink) are overlaid on individual data points (blue). Distributions of simulated data generated using behavioural measures from the “separate” regression model and using simple means can be found in the supplementary figuresFig. 6Distributions of error variance proportions for simulated measures of behavioural performance generated using our regression-based approach. These data were generated using behavioural measures estimated using the “joint” regression model, which explicitly modelled the effect of session. The mean proportion of variance for each sample size (purple), 90th inter-percentile range (dark pink), and interquartile range (light pink) are overlaid on individual data points (blue). Distributions of simulated data generated using behavioural measures from the “separate” regression model and using simple means can be found in the supplementary figures ## Results ### Reversal learning performance reliability First, we assessed reversal learning task performance using data collected from our participants. Behavioural performance in the reversal learning task was measured using accuracy, perseveration, staying behaviour (making the same choice on trial *t* as *t* − 1), and reaction time. Simple means and regression models were used to estimate behavioural measures from trial-wise measures; both a separate regression model for each session (separate) and a single model which explicitly estimated the effect of session (joint) were used. Accuracy, perseveration, and stay behaviour were calculated as proportions, while reaction time measures were calculated in milliseconds. Summary statistics for mean, standard deviation, and range are presented in Table 2 for all behavioural estimates. Overall, behavioural estimates of task performance show small but statistically significant increases in performance between sessions 1 and 2 (e.g., increased accuracy, reduced perseveration, increased staying after wins; paired-samples *t*-tests, all *p* values < [0.05/27]), although staying after losses also significantly increased between sessions 1 and 2. A mixed-effects logistic regression revealed a main effect of previous feedback on staying behaviour, with subjects switching more after losses than wins (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \beta =3.665, z=31.756, p<0.001 $$\end{document}β=3.665,z=31.756,p<0.001), while a mixed-effects linear regression revealed a main effect of previous feedback on reaction time, with subjects responding more quickly after wins than losses (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \beta =-23.748, t\left(149.4\right)=-6.84, p<0.001 $$\end{document}β=-23.748,t149.4=-6.84,p<0.001). Table 2Summary statistics for behavioural measures of reversal learning task performance in sessions 1 and 2Behavioural measureEstimation methodSessionMean*SD*MinMaxAccuracyMean10.690.070.570.8520.730.060.60.87Joint10.70.050.610.8120.720.050.620.82Separate10.70.050.610.820.730.040.640.82PerseverationMean10.130.060.010.2820.110.060.010.27Joint10.130.050.030.2620.110.050.030.25Separate10.130.050.040.2620.10.050.030.25Reaction timeMean1546.2111.4319.9917.12509.699.89349.5915.1Joint1540.5108.8316.6894.82506.798.48349908.5Separate1544.8109.8317.9911.72507.198.56351.3906.4Reaction time winsMean1528.8109.3312.28912499.499.75334.8933.4Joint1528.7105.7316.28752499.497.19338.4910.2Separate1528.8105.8315.1876.42499.496.86340.5913Reaction time lossMean1560.7120.9313.710282514.6106.3354.4892.8Joint1558117.1317.1986.92517.6103.4359.4905.9Separate1560.8116.6320.8994.82514.7102.5360.2899.8Reaction time win − lossMean1 − 31.8554.39 − 230.986.512 − 15.2246.63 − 181123.3Joint1 − 29.337.58 − 180.664.522 − 18.231.94 − 148.169.05Separate1 − 32.0336.53 − 171.444.832 − 15.3830.51 − 124.671.86StayMean10.730.120.450.9320.770.110.510.93Joint10.830.110.460.9720.860.090.560.98Separate10.770.120.420.9520.820.120.480.96Stay winsMean10.930.070.66120.960.050.681Joint10.940.060.690.9920.960.050.691Separate10.930.060.680.9920.960.050.70.99Stay lossMean10.440.230.030.8320.490.230.050.85Joint10.440.210.070.8220.490.220.070.83Separate10.440.210.070.8120.490.210.080.83Measures were estimated either using simple means (mean), or from mixed effects regression models calculated using separate regression models (separate) or a single model which explicitly modelled the effect of session (joint). Accuracy, perseveration, and stay behaviour are presented as proportions, and reaction times are presented in milliseconds. Test–retest reliability of task performance was assessed by replicating the statistical methods of Waltmann et al. (2022). These results are briefly summarised here, but see supplementary materials for a detailed overview. Accuracy reliability estimates ranged from good when jointly estimating the effect of session to poor when using mean values. Staying behaviour showed good reliability across all estimation methods, and showed good to excellent reliability for staying after losses, but only fair (mean and separate estimation) to good (joint estimation) reliability for staying after wins. Perseveration reliability was poor when using mean and separate estimation methods, and fair when using joint estimation. Reaction time showed good reliability when estimated using mean and separate estimation methods, and excellent reliability when using joint estimation, and the same pattern was observed for reaction time after a win. Reaction time after a loss had good reliability across all estimation methods; poor reliability was found for the difference in reaction between wins and losses when estimated using separate and separate estimation measures, but was good when joint estimation was used. Computational measures of task performance were derived by fitting reinforcement learning models to choice behaviour during the task. We fit two families of models that either altered the softmax choice function to define choice stochasticity (by varying the inverse temperature parameter β), or used a reinforcement sensitivity parameter (ρ) to determine the maximum difference in expected values between choices, placing a lower bound on choice stochasticity. Comparisons of model fit were performed using each model’s integrated Bayesian information criterion (iBIC) from the EM fitting approach (supplementary Fig. 3), which provided the best model fit scores and reliability of parameter estimates. The best-fitting model for the softmax family was the dual update model with a discount weight for the unchosen option with separate learning rates plus softmax temperatures for wins and losses (DU-2β2ακ). The best-fitting model for the reinforcement sensitivity family was the dual update model with separate reinforcement sensitivities for wins and losses plus a single learning rate parameter (DU-2ρα), which was also the best-fitting model overall (Table 3). This model had good to excellent estimates of reliability for all parameters (learning rate, reinforcement sensitivity for win, reinforcement sensitivity for loss) when model parameters were fitted using the EM approach. The parameters from the best-fitting model from the softmax family had reliability that ranged from poor to good. Model parameters had poor reliability when estimated using maximum likelihood methods (Supplementary Figs. 4 and 5). Table 3Summary statistics for parameter estimates from the best-fitting computational model in sessions 1 and 2Mean*SD*MinMaxReinforcement sensitivity win (session 1)2.050.920.364.43Reinforcement sensitivity win (session 2)2.511.140.425.5Reinforcement sensitivity loss (session 1) − 0.70.43 − 1.560.56Reinforcement sensitivity loss (session 2) − 0.710.41 − 1.620.58Learning rate (session 1)0.710.180.350.96Learning rate (session 2)0.710.150.340.95Parameters were estimated from a dual update model from the reinforcement sensitivity family, with separate reinforcement sensitivities for wins and losses, and a single learning rate. Parameters were fitted using expectation–maximisation. Overall, the reliability of both our behavioural and computational modelling measures of reversal learning performance were in line with those presented by Waltmann et al. (2022). Moreover, we found that the best-fitting models from both the softmax and reinforcement sensitivity families were the same reported by Waltmann et al. (2022). Next, we tested how robust these reliability effects were across varying sample sizes by investigating the effect of sample size on the individual components of variance used when calculating ICCs. ### Effects of sample size on behavioural measures of task performance We generated synthetic two-session data using our regression-based approach to investigate the effect of sample size on estimates of variance components for behavioural performance measures. These behavioural performance measures were derived from our three estimation methods. Figure 4 show the variance component estimates for data simulated using the distributions of behavioural measures, themselves estimated using a regression model that explicitly modelled the effect of session on behavioural performance (for figures using simple mean estimates and estimates using separate regression models for each session, see supplementary Figs. 8–13). Each point represents a given sample size between 10 and 300 for each iteration from the simulation, and overlaid are mean estimates for each sample size (purple), interquartile range (light pink), and the 90th inter-percentile range (dark pink). We then calculated at what sample size variance component estimates stabilised. Using variance component proportions from our collected data as the ground truth, we determined the *point of stability* for variance component estimates from each iteration of our simulation, using a range of half-widths (0.025, 0.05, 0.75, 0.1, 0.15, 0.2) for the *corridor of stability* (Fig. 3). For each half-width we determined the critical *point of stability* by calculating the 80th, 90th, and 95th percentiles of the calculated *points of stability* (Fig. 7, “joint” estimation method; see supplementary Figs. 14 and 15 for *points of stability* using data from the “separate” and “mean” estimation methods). As the half-width of the *corridor of stability* increased, the sample size required for variance component estimates decreased (Table 4). At narrower *corridor of stability* half-widths (i.e., at more precise estimates of variance component proportions), the sample size where variance component estimates stabilised were, on average, between sample sizes of 126.84 and 193.99 for half-widths of 0.05 (ground truth variance proportion ± 0.05) and were between 250.15 and 290.41 for half-widths of 0.025 (ground truth variance proportion ± 0.025) (4).Fig. 7Critical *point of stability* of variance component estimates for synthetic behavioural measures across a range of *corridor of stability* half-width values, calculated based on the distribution of data from the “joint” regression model. Critical *point of stability* calculations using behavioural measures from the “separate” regression model and using simple means can be found in the supplementary figuresTable 4Summary statistics for critical *points of stability* for between-subject, within-subject, and error variance components for synthetic behavioural measures across a range of *corridor of stability* half-width values and *point of stability* percentilesMeasureVariance component80th Percentile90th Percentile95th Percentilew = 0.025w = 0.05w = 0.075w = 0.1w = 0.15w = 0.2w = 0.025w = 0.05w = 0.075w = 0.1w = 0.15w = 0.2w = 0.025w = 0.05w = 0.075w = 0.1w = 0.15w = 0.2MedianBetween subjects3001676232.21410300245103.1522315300300146.1803722Within subjects189.23415101010252.3492416101029868.138231210Error296.210329.2151010300159.248.3251110300216.0591.05451911MeanBetween subjects294.44187.04117.1851.5117.2211300229.47150.3788.1928.715.34300258.26175.36122.5842.3521.58Within subjects195.3348.0217.8912.2210.0210237.6276.3829.2418.6711.7810.11274.36108.3643.7928.0114.7811.01Error260.67145.44115.6972.1319.4710.89288.91180.27128.51100.4233.4414.67296.89215.37146.58121.1249.6820.46*SD*Between subjects12.4993.72102.2440.127.81.49077.18105.9472.0614.785.54062.9893.8493.64228.37Within subjects76.4938.017.873.150.06056.4758.8913.177.472.860.3135.5679.6420.4412.35.671.7Error51.67111.84130.4896.6518.772.5127.4195.17121.68112.9340.929.68.881.61109.63120.6562.1517.65 ### Effects of sample size on computational modelling parameter measures of task performance To test for the effects sample size on estimates of variance components for computational modelling parameter estimates (for our best-fitting model), we generated synthetic two-session data using our regression-based approach (as above for our behavioural data, Fig. 8). We then calculated the critical *point of stability* for variance components of our simulated computational modelling parameters in the same way as we did for our behavioural measures (Fig. 9). As observed for our behavioural data, as the half-width of the *corridor of stability* increased, the sample size required for variance component estimates decreased. At narrower *corridor of stability* half-widths (i.e., at more precise estimates of variance component proportions), the sample size where variance component estimates stabilised were, on average, between sample sizes of 51.67 and 106.59 for half-widths of 0.05 (ground truth variance proportion ± 0.05), and were between 173.31 and 249.58 for half-widths of 0.025 (ground truth variance proportion ± 0.025) (Table 5).Fig. 8Distributions of within-subject, between-subject, and error variance proportions for parameter estimates from our best-fitting computational model estimated for different sample sizes generated using our regression-based approach. The mean proportion of variance for each sample size (purple), 90th inter-percentile range (dark pink), and interquartile range (light pink) are overlaid on individual data points (blue)Fig. 9Critical *point of stability* of variance component estimates for synthetic computational modelling parameters across a range of *corridor of stability* half-width valuesTable 5Mean critical *points of stability* for between-subject, within-subject, and error variance components for synthetic computational modelling parameters across a range of *corridor of stability* half-width values and *point of stability* percentilesMeasureVariance component80th Percentile90th Percentile95th Percentilew = 0.025w = 0.05w = 0.075w = 0.1w = 0.15w = 0.2w = 0.025w = 0.05w = 0.075w = 0.1w = 0.15w = 0.2w = 0.025w = 0.05w = 0.075w = 0.1w = 0.15w = 0.2MedianBetween subjects264.268291710103009542.1271410300124.0560.05361913Within subjects47201110101081.132.11811.11010126.054623161010Error203451610101029070291510103009342251210MeanBetween subjects232.538435.3319.6710.331027612654.0331.031510.67295.05169.3878.3742.720.713.68Within subjects94.7324.3312.6710.671010135.3738.0719.3313.710.3310160.3554.3827.0218.3312.3310.33Error192.6746.6718.3311.671010263.170.3731.3316.671010293.3396.0243.332512.3310*SD*Between subjects71.6146.3118.127.590.47033.9471.6127.4612.614.550.94793.7941.6817.837.043.32Within subjects85.0413.823.090.9400107.8920.088.224.480.470102.9229.6112.597.933.30.47Error53.9815.976.342.360045.320.868.734.64009.4322.5810.665.722.050 ### Effects of sample size on ICC(A,1) and variance component associations We tested for associations between individual variance components and ICC(A,1) measures of reliability. To do this, we calculated Spearman’s correlation coefficient and *p* values (Bonferroni-corrected) between each variance component measure and ICC(A,1) across the 1000 synthetic datasets generated at each sample size *n*, for each behavioural (Fig. 10, Fig. 11, Fig. 12) and computational modelling parameter (Fig. 13) measure of task performance. Overall, we observe that between-subjects variance had a strong positive correlation with ICC(A,1), error variance had a strong negative correlation with ICC(A,1), and within-subject variance had a weak/no correlation with ICC(A,1). Secondly, these associations are true for most behavioural and computational measures of task performance. Yet, where there are inconsistencies (e.g., staying after wins for between-subject variance and accuracy for within-subject variance), we find that the strength of these associations is consistent even when the behavioural measures were estimated using different methods (see supplementary Figs. 16–21). We also find that the strength of these associations appears relatively stable across different sample sizes, as indexed by the consistency of Spearman correlation coefficients.Fig. 10Effects of sample size on the association between ICC coefficients and variance component estimates for behavioural measures. For each behavioural measure and at each sample size, we took the set of 1000 simulated datasets and calculated correlation coefficients to measure the strength of the association between each dataset’s respective ICC(A,1) and variance component estimates. The point estimate for the correlation coefficient and its statistical significance (coloured green for significant, red for non-significant; Bonferroni-corrected) were then plotted. Overall, between-subject variance was strongly positively correlated with ICC(A,1). These plots were generated from synthetic data generated using the distribution of behavioural measures from the “joint” regression approach, which explicitly modelled the effect of session. See supplementary figures for plots generated using data from the “separate” regression approach and using simple meansFig. 11Effects of sample size on the association between ICC coefficients and variance component estimates for behavioural measures. For each behavioural measure and at each sample size, we took the set of 1000 simulated datasets and calculated correlation coefficients to measure the strength of the association between each dataset’s respective ICC(A,1) and variance component estimates. The point estimate for the correlation coefficient and its statistical significance (coloured green for significant, red for non-significant; Bonferroni-corrected) were then plotted. Overall, within-subject variance was weakly or not correlated with ICC(A,1). These plots were generated from synthetic data generated using the distribution of behavioural measures from the “joint” regression approach, which explicitly modelled the effect of session. See supplementary figures for plots generated using data from the “separate” regression approach and using simple meansFig. 12Effects of sample size on the association between ICC coefficients and variance component estimates for behavioural measures. For each behavioural measure and at each sample size, we took the set of 1000 simulated datasets and calculated correlation coefficients to measure the strength of the association between each dataset’s respective ICC(A,1) and variance component estimates. The point estimate for the correlation coefficient and its statistical significance (coloured green for significant, red for non-significant; Bonferroni-corrected) were then plotted. Overall, error variance was strongly negatively correlated with ICC(A,1). These plots were generated from synthetic data generated using the distribution of behavioural measures from the “joint” regression approach, which explicitly modelled the effect of session. See supplementary figures for plots generated using data from the “separate” regression approach and using simple meansFig. 13Effects of sample size on the association between ICC coefficients and variance component estimates for computational modelling parameters. For each parameter and at each sample size, we took the set of 1000 simulated datasets and calculated correlation coefficients to measure the strength of the association between each dataset’s respective ICC(A,1) and variance component estimates. The point estimate for the correlation coefficient, and its statistical significance (coloured green for significant, red for non-significant; Bonferroni-corrected) were then plotted. Overall, between-subjects variance was strongly positively correlated with ICC(A,1), error variance was strongly negatively correlated with ICC(A,1), and within-subjects variance was weakly or not correlated with ICC(A,1) ### Effects of noise on variance estimates Lastly, we tested how adding noise to our simulated data influenced the calculation of *points of stability*. We generated synthetic two-session data using our regression-based synthesis approach as described above, with an additional step of adding Gaussian noise to the calculated value for session 2. We generated multiple datasets with varying levels of noise; by sampling noise from normal distributions with mean = 0 and *SD* = [0.25, 0.5, 1, 1.5, 2], we were able to investigate the effects of increasing levels of noise on *point of stability* calculations, in comparison to the originally simulated data (*SD* = 0). Adding noise to simulated data caused a monotonic change in the *point of stability* of all variance components with increasing levels of noise (Fig. 14). Across a range of percentile values (80th, 90th, 95th) the *point of stability* is equal to the largest sample size (300), indicating that a *point of stability* was not reached for a majority of variance component estimates, meaning that variance component estimates did not stabilise as simulated noise increased. These results were consistent across *point of stability* percentiles, *corridor of stability* half-widths (w), and behavioural estimation approaches (“joint”, “separate”, and “mean” modelling approaches; see supplementary Figs. 22 and 23; included code can be used to recreate figures with different corridor widths).Fig. 14Effects of noise on *point of stability* calculations for synthetically generated data. As the amount of noise added to synthesised behavioural and computational measures of task performance is varied, a monotonic change in the *point of stability* for variance components is observed. For the majority of variance components, increasing amounts of simulated noise cause variance components to fail to reach a *point of stability* before the largest sample size is reached, meaning that variance component estimates remain unstable. These plots were generated from synthetic data generated using the distribution of behavioural measures from the “joint” regression approach, which explicitly modelled the effect of session ## Discussion In this paper we investigated the test–retest reliability of behavioural and computational measures of performance on a reversal learning task, and the effects of sample size on such estimates of reliability. We calculated the reliability of these measures using several approaches, including ICCs and variance decomposition, replicating a previous study on this topic (Waltman et al., 2022). The retest–reliability of ICCs for our behavioural measures were good to excellent for staying and reaction time behaviour, while accuracy and perseveration were less reliable between sessions. ICCs for parameter estimates from our best-fitting computational model (single learning rate and separate reinforcement sensitivity parameters for wins and losses) showed good reliability when estimated using an expectation–maximisation (EM) model fitting approach. These results were broadly in line with the previous findings of Waltmann et al. (2022). Using our behavioural and computational measures of task performance, we then investigated the effects of sample size on individual components of variance. We used a regression-based approach to generate statistically related and plausible synthetic two-session datasets based on the underlying statistical properties of our collected data. Sample size influenced estimates of within-subject, between-subject, and error variance for all behavioural measures of task performance and all computational modelling parameter estimates. Importantly, we demonstrate that variance component estimates do not stabilise until sample sizes much greater than those often used in test–retest research are achieved. We also tested the effects of sample size on the association between measures of reliability, as assessed using ICC(A,1) and individual variance components. This is important to understand because individual variance component calculation decomposes an ICC coefficient into its constituents and enables the investigation of how sources of variance contribute to the summary reliability statistic represented by ICCs (the ratio of between-subject variance over the total amount of variance). ICCs had significant and large positive correlations with between-subject variance and negative correlations with error variance across a broad range of performance measures, and these correlations remained relatively stable across sample sizes. However, ICCs were either weakly or non-significantly correlated with within-subject variance, a finding that also remained relatively stable across sample sizes. These results suggest that ICCs may be insufficient for discriminating reliability, particularly for within-subject variance. Lastly, we demonstrate that increasing noise between data from sessions 1 and 2 increases the critical *point of stability*, meaning that data from more participants are required before stable estimates of variance components are reached. Therefore, studies of reliability should ensure they use variance decomposition methods alongside larger sample sizes to more informatively measure the reliability of task performance. The results presented here highlight the importance of sample size considerations for test–retest reliability work, and support existing work indicating that greater sample sizes than are often used are required for reliable (individual differences) research (Button et al., 2013; de Winter et al., 2016; Hedge et al., 2018; Hirschfeld et al., 2014; Kretzschmar & Gignac, 2019; Marek et al., 2022; Paccagnella, 2011; Schönbrodt & Perugini, 2013). For instance, Schönbrodt and Perugini (2013) suggest there are few scenarios where personality psychology using correlational methods should justify sample sizes smaller than 150, while Marek et al. (2022) and Gell et al. (2023) suggest that thousands of participants may be required for reliable brain–behaviour associations, although this depends on the reliability and validity of the measures used (see also Spisak et al., 2023, for discussion). A viable alternative to collecting large samples is to instead improve reliability by increasing the density of data collected within a set of individuals (Kraus et al., 2023; Smith & Little, 2018; Tiego et al., 2023). In precision research, each subject acts as their own replication unit, with a large amount of data collected within small/single-subject units. This may be particularly useful in situations where practical impediments, such as time and funding restrictions or specialist populations, would prevent collection of data from hundreds or thousands of individuals. This approach may also enhance the prediction of momentary factors that influence the rank order of a given data point. For instance, intensive longitudinal designs (Lydon-Staley et al., 2019) could be used to enhance estimates of both within- and between-subject effects. This would have the added benefit of providing insight into how momentary changes in cognitive and affective state influence behaviour and model parameter estimates, which are missed in large-*N* studies with a single time point, since temporal dynamics cannot be modelled. One relevant example of the utility of this approach comes from Schaaf et al. (2023), who found that the current state of an individual significantly influences their reward learning (using data from two time points). Yet, this is a nascent field of research, and few insights into temporal aspects and predictors of reward learning behaviour exist. One important similarity between the work presented here and previous work looking at the reliability of reversal learning task performance (Schaaf et al., 2023; Waltmann et al., 2022) is the consistency of the best-fitting models despite differences in task structure. For instance, the best-fitting model of Waltmann et al. (2022) was the same as our best-fitting model here (dual update, single learning rate, and separate reinforcement sensitivity parameters for wins and losses). Similarly, although Schaaf et al. (2023) only fit models from the softmax family, our best-fitting model in the softmax family (dual update, separate learning rates for wins and losses, and an update discount weight for the unchosen option) matched the best-fitting model of both Schaaf et al. (2023) and Waltmann et al. (2022). This consistency suggests these models are useful for approximating latent processes underlying reversal learning. The commonality of counterfactual updates across all these models makes sense in the context of two-choice reversal learning, where only one of the two outcomes can be optimal at any given moment. Therefore, when an agent updates the expected value of an action based on an outcome, simultaneously updating the expected value of the unchosen action using a counterfactual outcome makes behaviour more responsive and able to rapidly adjust in response to a change, in line with Bayesian state inference approaches to reversal learning (Bartolo & Averbeck, 2020; Costa et al., 2015). With experience, this dual updating should reduce perseverative responding as the agent learns with greater fidelity when transitions in the optimal choice assignment occur. This assumption is supported by our finding of lower levels of reliability for both accuracy and perseveration between sessions (which typically increase and decrease, respectively, for subjects) versus other behavioural measures, and suggests that subjects get “better” at the task. It may also explain why we observed lower reliability for our learning rate parameter estimates relative to other model parameters, as the rate at which expected values are updated is refined through experience, while other parameters (e.g., reinforcement sensitivity) may be less experience-dependent. Another similarity is that our work here shows that reliability, as assessed using ICCs, is broadly in line with previous work (Schaaf et al., 2023; Waltmann et al., 2022). For instance, our confidence intervals around ICC measures for our collected behavioural data overlapped with those presented by Waltmann et al. (2022) for all behavioural measures (except for lose–stay behaviour estimated from the joint session regression model) and parameter estimates from our best-fitting model using the EM approach (dual update, single learning rate, and separate reinforcement sensitivity parameters for wins and losses). Given our larger sample size and narrower confidence intervals relative to Waltmann et al. (2022), we suggest that our estimates of reliability presented here may be more representative of the true underlying reliability of reversal learning performance measures. One important sample size-related consideration for test–retest reliability studies is ensuring they are sufficiently powered, and this could be achieved by using expected ICC values in power calculations. This approach is similar to the use of effect sizes for a priori sample size calculations. In the latter case, effect sizes from previously published work are used to identify the sample size required to achieve the given effect size at a particular level of statistical power. For ICCs, Doros and Lew (2010) present a method where the width of confidence intervals is used to estimate appropriate sample sizes for reliability studies. While the work presented here provides valuable insight into the reliability of reversal learning task performance measures, there are several limitations worth mentioning. Firstly, it is important to highlight that these behavioural data were collected from an online sample. Although online data collection enables large amounts of data to be more readily collected, as researchers we are unable to control the environment in which each subject completed the reversal learning task. To mitigate this we included response checks in our task and questionnaire measures to identify and exclude subjects that were clearly inattentive. However, there may be nuanced and non-systematic differences in the behaviour of our subjects that influenced how reliable their performance was over the two testing sessions. To account for some of these challenges, future work could, as previously mentioned, collect data from the same subjects over a greater number of testing sessions or could use tools such as WebGazer (Papoutsaki et al., 2016) to track gaze directions for compliance monitoring. A second limitation of the work presented here is that point estimates of model parameters across subjects were taken using the best-fitting model at the group level. However, the best-fitting model for a given subject will not necessarily be the same as that for the group, and alternative approaches to model fitting can be used to infer the best-fitting model at the level of both the subject and the group (Piray & Daw, 2020; Piray et al., 2019; Williams & Christakou, 2022). Coupling individualised model fits with momentary measures of individual state could, again, improve the explanatory power of model parameters. Finally, it may be the case that task performance or within-subject variance changes as subjects gain further experience on the task. Although we provide subjects with explicit instruction about the task’s general structure, such as one choice is better than the other and that the better choice will change throughout the task, subjects may still be refining their understanding of aspects of the task structure (the statistical relationships between actions and outcomes, their time estimate since last reversal, etc.), even after two sessions. Therefore, measuring reliability while subjects are still learning these representations will place an upper bound on reliability, and will be dependent on how quickly a stable representation of the task environment is generated, which in turn will vary between individuals. In the future, it may be worth considering how reliable performance is after sufficient overtraining on a given task. In summary, we assessed the reliability of reversal learning behaviour using data collected from a large online sample and found good reliability of behavioural and computational model parameters at the group level, in line with findings from previous literature. However, our results also suggest that while behaviour may appear stable at the group level based on ICC values, sample size contributes significantly to variability in the estimates of variance components that underlie ICCs. Moreover, associations between estimates of variance components and calculated ICC values appear to remain relatively stable across sample sizes, with between-subject variance being highly positively and error variance being highly negatively correlated, but within-subject variance being weakly or non-significantly correlated with ICC values. This effect for within-subject variance challenges traditional practices in assessing test–retest reliability, and demonstrates the importance of understanding individual factors that contribute to (un)reliability. For instance, within-subject variance could be due to momentary differences in cognition and/or affect, and future work should aim to address how momentary state influences behaviour. These results also hold for the effects of sample size on estimates of computational modelling parameters, suggesting further characterisation of the stability of model parameters within and across tasks over time is needed before point estimates of parameter values can be considered stable trait-like measures. ## Supplementary Information Below is the link to the electronic supplementary material.Supplementary file1 (DOCX 9440 KB)