{'fraction_non_alphanumeric': 0.024678837052062204, 'fraction_numerical': 0.009519234190954058, 'mean_word_length': 4.7367561498417885, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Getting Blood from a Stone: Improving Neural Inferences without Additional Neural Data

In recent years, the cognitive neuroscience literature has come under criticism for containing many low-powered studies, limiting the ability to make reliable statistical inferences. Typically, the suggestion for increasing power is to collect more data with neural signals. However, many studies in cognitive neuroscience use parameters estimated from behavioral data in order to make inferences about neural signals (such as fMRI BOLD signal). In this paper, we explore how cognitive neuroscientists can learn more about their neuroimaging signal by collecting data on behavior alone. We demonstrate through simulation that knowing more about the marginal distribution of behavioral parameters can improve inferences about the mapping between cognitive processes and neural data. In realistic settings of the correlation between cognitive and neural parameters, additional behavioral data can lead to the same improvement in the precision of inferences more cheaply and easily than collecting additional data from subjects in a neuroimaging study. This means that when conducting an neuroimaging study, researchers now have two knobs to turn in a design analysis: the number of subjects collected in the scanner and the number of behavioral subjects collected outside the scanner (in the lab or online).


Introduction
One factor which limits progress in human cognitive neuroscience is statistical precision and power (Munafò et al., 2019;Cremers et al., 2017;Yarkoni, 2009). Many reported effects are small relative to the noise in neuroimaging signals, a problem compounded by the small samples sizes used in most studies. Within a hypothesis testing framework, a lack of precise estimates of some measurements can make it difficult to replicate reported effects. For example, even quite strong, known effects such as the relationship between motor movements and fMRI BOLD activation in sensorimotor regions can be unlikely to reach significance in current standard sample sizes of around 25-30 subjects . Having lower precision also means that any estimated effect that is significant will also have a high probability of the estimate being much greater in magnitude or even having the wrong sign compared to the true effect (Gelman and Tuerlinckx, 2000;Cremers et al., 2017;Yarkoni, 2009). However, due to a combination of factors (including lack of sufficient funding for large sample studies), research in human cognitive neuroscience tends to be under-powered in many cases .
Several remedies have been proposed for this situation, the most obvious being to simply collect more neuroimaging data Munafò et al., 2019). Because the funding available to most labs is relatively small (compared to the cost of a large scale neuroimaging study), this might require a move towards working in larger consortia around critical topics. Another approach might be to leverage the power of publicly available open datasets, which can allow for meta-analyses across several smaller experiments . However, increasing sample sizes is not the only option for dealing with a lack of precision. For example, increasing the signal-to-noise ratio of neuroimaging measures by creating more detailed statistical models of the neuroimaging signal (Lindquist et al., 2009), improving experimental design  or even improving the measurement process itself (Feinberg and Yacoub, 2012;Lombardo et al., 2016) are all measures that are likely to help. In addition, it has been suggested that a focus on computational cognitive modeling and behavior can also allow for extracting more signal (Krakauer et al., 2017;Palmeri et al., 2017) by creating better models of the underlying cognitive and neural processes. In that theme, this paper explores a potentially under-studied option for improving the statistical precision of human cognitive neuroscience studies: collecting additional behavioral data without neural recordings.
We begin by noting that many common statistical analyses in cognitive neuroscience involve behavioral data, typically collected from the subjects who also provided neural data. In one common scenario, the interest is in the relationship between individual differences in behavioral traits and more static measures of neural activity, such as performance in a task and functional connectivity (Rosenberg et al., 2016). In another typical analysis, researchers are interested in how within-subject changes in a cognitive state variable over time relate to changes in the neural signal from a particular brain region, such as the prediction error in a reinforcement learning task (O'Doherty et al., 2003). One interesting, but often overlooked, point is that these behavioral measurements (i.e. survey questionnaires, behavioral responses to stimuli, summaries of task performance or estimates of cognitive model parameters) are themselves noisy estimates of an individual's underlying cognitive processes and the distributions of these patterns that exist in the population. Thus our ability to infer relationships between brain and behavior is limited by sampling error and estimation error. Sampling error is a result of the fact that a small sample may not be representative of the population, so estimates of the relationship based on the sample may be far from the true relationship. A small sample therefore prevents us from being confident about the value of a parameter in a statistical model. Estimation error is a form of measurement error that arises because researchers cannot measure latent cognitive variables directly and instead use models to estimate their values. The relationship between neural recordings and cognitive variables measured with error will necessarily not be as strong as the relationship involving the true values, a notion sometimes known as regression dilution or attenuation. Ignoring these two sources of error can prevent us from learning as much as possible from neural data.
In the following, we demonstrate that one way to reduce both sources of error and improve statistical estimates of the relationship between the brain and behavior is to collect behavioral data from subjects without collecting corresponding neuroimaging data. This idea is somewhat counterintuitive in light of the fact that most common approaches in neuroimaging analyses assume that the only relevant data is from individuals who's neural activity was recorded. However, behavioral data (especially data collected over the Internet) is often significantly less resource-intensive to collect than neuroimaging data, in terms of money and the time and effort of both researchers and human subjects. We show how using modern statistical methods that can leverage information from the additional behavioral data can help mitigate the two types of error. To anticipate one consequence of our analysis, we are able to provide a cost-benefit analysis (in terms of statistical power and precision) of collecting data from subjects with only behavior against collecting more expensive additional data from subjects with both behavior and neural recordings. This analysis shows that in some cases, due to its relatively low cost, a more effective use of resources is to collect more behavioral data instead of additional neuroimaging subjects. This result has implications for designing studies in cognitive neuroscience both prospectively and retrospectively (i.e., the re-analysis of archived open-science data sets.)

Common paradigms for relating brain and behavior in neuroimaging studies
For the purposes of articulating our thesis, it will be useful to describe two common types of analyses used in cognitive neuroscience that use behavioral data to interpret neural signals and lay out some notation for the rest of the paper. While this is by no means an exhaustive list, we claim that many analyses can be understood as falling into two broad categories.
Regression (Conditional Distribution) In many cognitive neuroscience studies, researchers are interested in how the conditional distribution of neural measures changes as a function of variation in cognitive states induced by a task. For instance, researchers might be interested in regions of the brain that track whether or not it is viewing a face (Kanwisher et al., 1997), the valence of the currently presented stimulus (Chikazoe et al., 2014) or the prediction error in a reinforcement learning model of decision making (O'Doherty et al., 2003). More formally, a subject i presented with a set of stimuli s i that occur on trials (or times) t ∈ {1...T i }. Each stimulus presentation s it is assumed to induce some latent cognitive state c it that is associated with some neural activity n it . Researchers can identify neural signals related to the cognitive states by fitting a regression model where where it is uncorrelated noise (for the purposes of this discussion, we ignore features of particular signals like the fMRI BOLD hemodynamic response). However, this is not quite so straightforward because researcher need to derive estimates of the cognitive state for these analyses. The most common assumption is that the relevant states are closely tied with objective features of the stimulus, such as the color of an image or experimental conditions like the reward offered for a choice, so c it can reasonably be replaced by the stimulus feature or experimental condition itself.
In other cases, the cognitive states of interest are tied to subjective features of the stimulus, such as the valence associated with s it or the typicality of s it for a particular category. In some cases, researchers collect behavioral ratings b it which depend on these latent states. Researchers then assume that the ratings come from a cognitive process C such that b it ∼ C(c it , s it ) and obtain the c it from that model. In the simplest case, this can mean simply using the valence ratings given by subjects (Chikazoe et al., 2014) or using the average typicality from a group of subjects tested separately (Wilson-Mendenhall et al., 2015).
More recently, "model-based" neuroimaging analyses have investigated cognitive states that are involved in a particular cognitive process C such that in addition to the current stimulus s it , c it may also depend on θ i are latent cognitive "trait" variables that do not vary over the relevant time as well as the past states c i 1:(t−1) . Researchers can then use the behavioral data to fit the model b it ∼ C(θ i , c it , s it ) and obtain the relevant c it . These states may be, for instance, the reward prediction error or predicted value in a reinforcement learning model of sequential decision making (O'Doherty et al., 2003).
Because the latter two types of analyses involve behavioral data, those will be our focus in this paper. In particular, we will focus on the model-based case but nearly all of our analyses will also apply to the second case.
Modern multivariate pattern analysis techniques (Polyn, 2005;Haxby, 2001) invert this regression and try to predict the cognitive state from neural signals. Typically, this means predicting stimulus properties or experimental conditions but in some cases the cognitive state will be inferred from behavioral data (e.g. Mack et al., 2013;Kragel and Polyn, 2016). While the example analyses in the rest of the paper do not directly apply to this setting, similar principles of collecting additional behavioral data may hold there as well.
Correlation (Joint Distribution) In other cases, researchers are interested in the joint distribution of behavioral measures and neural measures (often across subjects). For instance, is there a relationship between individual differences in behavior and individual differences in static measures of the brain? Researchers often use a direct measure of a subject's behavior b i such as answers to a questionnaire (Treadway et al., 2013) or performance on a behavioral task (such as a sustained attention task Rosenberg et al., 2016). More recently (e.g. Homan et al., 2019), researchers have also used cognitive "trait" variables which can be related to behavioral data through cognitive models (i.e. θ i above). Static measures can include measures like functional connectivity or average resting state activity that do not depend on time. Typically we only get one n i per subject and the measure is not always collected simultaneously to the task where b i is recorded. In this setting, researchers are typically interested in characterizing the joint distribution of neural measures n i and behavioral (b i ) or latent cognitive (θ i ) trait variables. In particular, if the neural measures and behavioral variables can be approximated by a bivariate normal, researchers want to find neural measures where the correlation ρ is non-zero.
Another common analysis that involves correlation is Representational Similarity Analysis (Kriegeskorte, 2008). In these analyses, researchers correlate the similarities of the neural response to pairs of stimuli with the similarities of the two stimuli according to a model. Most often, these models are created based on stimulus properties but they are sometimes derived from behavioral ratings (Chikazoe et al., 2014;Bruffaerts et al., 2013). In these cases, our analyses about the benefits of additional behavioral data will apply as well.
We argue that these two categories of research design and inferential techniques cover a vast array of studies in cognitive neuroscience and neuroimaiging in particular. Each of these analyses can be thought of as attempting to learn about the conditional distribution of neural signals given behavioral data (regression) or the joint distribution of neural and behavioral data (correlation). The focus of this paper asks whether we can improve inferences for these two types of analyses by learning more about the marginal distribution of behavioral data, that is, collecting more behavioral data without neural data.

Two Sources of Error in Neuroimaging Analyses
As the Introduction laid out there are two key sources of error in neuroimaging analyses. In this section, we describe each of the two sources of error in more detail as well as how they can affect the precision of inferences about the relationship between neural signals and behavior.

Sampling Error
In a correlation analysis, we want to know what ρ is in the general population but estimate it by using the sample correlation r. With a small number of subjects, r can be quite far from the true ρ. An approximate confidence interval for ρ can be constructed as arctanh(r) ± 1 √ N −3 , where N is the number of subjects or, more generally, the number of data points used to estimate r (Fisher, 1915). As this makes clear, the precision of the estimate (and thus finding a ρ that we can say is non-zero with confidence) is constrained by the size of our sample. In many neuroimaging studies this is often relatively small because neuroimaging data is expensive and time consuming to collect. In Figure 1, we show an example of where the smaller sample size can lead to noisy estimates and larger error bars. 1 Sampling error can also be an issue in regression analyses as well. However, because in that case we are only interested in the conditional distribution of the neural signals, we cannot reduce sampling error by collecting more behavioral data so do not discuss it here.

Estimation Error
In the model-based regression analysis described above, we want to find a set of neural signals that have a linear relationship with a subject's latent cognitive states c i . 2 The true model we are interested in is where in particular we are interested in which neural signals have a significant β ci parameter (typically, at the group level). In practice, because we don't have direct access to c it we typically replace it with estimatesĉ it by fitting the model C(θ, c it , s it ) to the behavioral data b it . It is typically assumed that the model doesn't perfectly account for behavior and that there is some variability in responses beyond what can be predicted by the model. For instance, when fitting a reinforcement learning model, a decision noise parameter is often used to deal with variability beyond what can be accounted for by the latent representations in the model. Because we don't have infinite data, this typically means that several sets of parameters could account for the data even if the model is well identified. In other words, the parametersĉ it are estimated with some uncertainty. This may seem obvious to many readers but standard techniques that are typically used in cognitive neuroscience do not easily allow this uncertainty to be taken into account in the above regression. Instead of the model above, researchers often simply plug in the estimatedĉ i into the regression, i.e.
In this model, the value of the slope estimate using the estimated latent variables, βĉ i , will not necessarily be the same as when using the true parameter β ci . In fact, it can be shown that on average, the least squares estimate ofβĉ i will be where ρ(c i ,ĉ i ) is the correlation between the true c i and the estimatesĉ i (Wilson and Niv, 2015). This is perhaps not always relevant since the neuroimaging signal can often only be known up to a multiplicative constant. However, Wilson and Niv (2015) also showed that, as long as the estimates are uncorrelated with the regression error term it , the expected t-statistic for theβĉ i iŝ where T i is the number of trials from subject i used to estimate the regression and σ 2 is the variance of the residuals it , due to noise in the neuroimaging signal. We can see that this t-statistic gets closer to 0 as ρ(c i ,ĉ i ) gets further from 1. A smaller t-statistic means that the standard error (or any multiple of it such as a confidence interval) is more likely to overlap with 0. Thus having a lower correlation between the estimates and the true values, or, in other words, greater estimation error, can make it harder to find relevant signals in the typical regression setting. Estimation error can also be relevant for correlation analyses looking at correlations between latent cognitive trait variables and neural signals. If we estimate the trait variable θ i from behavioral data, Katahira (2016) showed that the correlation between estimatesθ and structural neural signals n, can be decomposed as ρ(n,θ) = ρ(n, θ)ρ(θ,θ) Therefore, the correlation is upper bounded by the true correlation between the neural signals and latent trait variables ρ(n, θ) and will decrease as the correlation between the estimates and the true values ρ(θ,θ) decreases (or as estimation error increases). Because the size of confidence intervals for a correlation coefficient only depend on the total number of subjects, the intervals for a correlation with a noisier estimate will be more likely to overlap with zero. This means that again, greater estimation error will make it harder to find relevant signals.

The Benefits of Additional Behavioral Data
Having laid out why these two errors can deteriorate our inferences, we now demonstrate how collecting additional behavioral data can help. In addition, for each case, we lay out a tradeoff between collecting cheaper behavioral data and more expensive imaging data in terms of how much each improves our inferences. These analyses show in which cases it may actually be more valuable to collect behavioral data than neuroimaging data, even if the main object of interest is the relationship between brain and behavior. To preface, section 3.1 demonstrates how little known estimators are able to make use of additional behavioral data to reduce sampling error for correlation analyses. We then follow with simulations to understand where this estimator improves over standard ones. Section 3.2 then documents how the use of hierarchical modeling can allow for combining datasets and improving estimation for common regression analyses as well as correlation analyses. Further analytical work shows the benefits of decreasing estimation error for making group-level inferences about the relationship between cognitive processes and neural signals.

How to Decrease Sampling Error: Behavioral Data as Neuroimaging Data with a Missing Variable
Typically we are interested in generalizing beyond our sample and characterizing the joint distribution of neural signals and behavioral variables in the population. However, as mentioned, most neuroimaging studies are limited by getting a relatively small sample. One somewhat counter-intuitive notion from the field of statistics focused on missing data is that we can often get better estimates of a parameter in the population (in terms of a lower mean squared error) by including data points that do not have every variable recorded in our analysis (Little and Rubin, 2019). This is because we are often assuming that variables have some joint distribution. If the correlation between two variables is not zero then we can know something about the value of what a missing variable must be from knowing the value of the other variable. If we have collected neural and behavioral data from N nb subjects and only behavioral data from an additional N b subjects, Anderson (1957) showed that the maximum likelihood estimates of the mean and standard deviation of n, as well as the correlation between n and b, ρ, are not the same when we include all of the data points as they would be if we only include the data with all variables recorded (the complete data). We can see this by rewriting the likelihood for a bivariate normal. If where µ = (µ n , µ b ) and Σ = σ 2 then, we can factor the likelihood as The advantage of writing the joint distribution as we did above is that the first term is just the regression of n on b which is only defined for the data points that have both variables recorded. Therefore, we can get the maximum likelihood likelihood estimate for α nb , β nb and σ nb using standard least squares regression estimates with the N nb data points. Maximizing the likelihood of the second term only depends on the b data so we can use standard estimates of the mean and variance using all of the data. We can write the correlation as Then plugging in the maximum likelihood estimatesβ nb ,σ nb andσ b , we get the maximum likelihood estimate for ρ using all of the data: Because the maximum likelihood estimate of σ b using all N nb + N b data points will not in general be the same as only using the N nb data points, the maximum likelihood estimate of ρ using all the data will not be the same as using only the complete data.
We can now try to illuminate how this new estimate is different from the correlation estimate based on only the N nb complete cases. Little and Rubin (2019) derive the missing data estimatorρ in a different way, showing that where r,σ b ,σ n are the estimates of the correlation (i.e. Pearson's r), the variance of the behavioral variable and the variance of the neural variable among the N nb complete cases.σ n is the variance of n after adjusting for theσ b estimate based on all of the data. Little and Rubin (2019) derive the formula forσ n aŝ We can instead write this in terms of r and rearrange terms, i.e.
Plugging this back in, we get To our knowledge, we are the first to show that the Anderson (1957) estimate,ρ, only depends on the correlation estimate in the complete cases, r and the ratio of the variance of the behavioral variable in the complete cases and the full data set,σ  (1957) estimate (ρ) of the correlation using all of the data as a function of Pearson's r using only N nb complete cases and the difference in the log of the variance of the behavioral data in the N nb complete cases and all N b cases. In generalρ is modified up or down by the relationship between the variance of the behavioral data on the full data set and that for the complete cases.
marginal distribution of the behavioral data contributes to the estimated measure of association. In figure 3, we plot howρ changes as a function of these two variables. In general, if the variance increases when including the additional data, the correlation estimate increases and visa versa.
3.1.1 When does it work? Garren (1998) showed analytically that asymptotically, the maximum likelihood estimateρ using all the data will have lower error than estimates for r using only the complete data. This result suggests that there can be cases where, if one variable is more expensive to collect than the other one (for instance neuroimaging data and behavioral data), the same quality of estimates of the correlation can can be achieved for lower cost by collecting some data with the more expensive variable missing (Hocking and Smith, 1972). However, the result is an asymptotic result, meaning that it is only guaranteed to apply in the limit of infinitely large data sets. It is therefore still unclear when collecting more behavioral data will be better than collecting more neuroimaging data in practice. In the following, we conduct a simple simulation to explore this tradeoff. Following the description above, suppose we have collected a behavioral variable (b) and a neuroimaging variable (n) from N nb subjects. In addition, we collected just behavioral data from N b subjects. For simplicity, we assume that we have a single behavioral measure and a single neuroimaging measure from each subject.
In order to specify the generative model with missing data, we can create a matrix of data BN from our neuroimaging subjects, an N nb -by-2 matrix with n as the first column and b 1:N nb as the second. The data matrix is generated for i ∈ [1, N nb ] according to where Σ nb1,1 = σ 2 n , Σ nb2,2 = σ 2 b and Σ nb1,2 = Σ nb2,1 = ρσ 2 n σ 2 b . The additional behavioral subjects then have distribution For the purposes of our simulation, we set µ nb = 0 and σ 2 n = σ 2 b = 1. Having generated data according to the above model, we now compare two ways of estimating the correlation, standard Pearson's r using only data with both variables and the Anderson (1957) estimate of the correlation using all of the data.
We compare estimates with for four values of N nb (25, 50, 100 and 200) and ten values of ρ from 0 to .9. To investigate the benefit of the additional behavioral data for methods that can make use of the information, for each value of N nb , we test five values of N b (0,25,50,100,200). Finally, to compare the benefit of the additional behavioral data to collecting an additional neuroimaging subject, we compare all estimates to Pearson's r with one additional nueroimaging subject. For each set of (N nb , N b , ρ), we generate 10,000 datasets.
Of note here is that collecting neuroimaging data can be quite expensive. For instance in 2022, collecting fMRI data at New York University's Center for Brain Imaging cost researchers around $450 per hour of subject time. In addition to the price, collecting neuroimaging data can have many other nonmonetary costs such as requiring the subject to come into the lab, requiring the time of trained researchers or technicians to run the machine and only being able to run one subject at a time. In contrast, many behavioral studies can be easily run in parallel online without significant monitoring by researchers for a cost of around $11 per hour of subject time (Hara et al., 2018;Crump et al., 2013). Thus, if 40 behavioral subjects can improve inferences more than a single fMRI subject, this might be a worthwhile tradeoff in many experimental designs. Figure 4 plots how additional behavioral data can influence the mean squared error of the correlation estimate using Anderson's estimator. Comparing the two methods of correlation estimation, we find that the effect of the additional behavioral data can be nonlinear in the size of the underlying neuroimaging dataset and the true value of ρ. Having too few neuroimaging data points to estimate a weak correlation (e.g. only 25 subjects for a correlation of .1) makes the behavioral data less useful and it can be better to collect more neuroimaging data. This finding resembles similar simulation results from Garren (1998). At low amounts of neuroimaging data and lower true values of ρ, the Anderson estimate of the correlation using all of the behavioral data can actually perform We also compare to estimates from a neuroimaging dataset with N nb + 1 subjects. The orange line is generally below the blue line showing the improved estimate that comes from one additional neuroimaging subject. The curved green line shows the Anderson (1957)ρ changes in the estimate as the number of addition behavior only subject are added to the analysis.
worse than the standard Pearson estimate. However, as the true ρ and number of neuroimaging subjects increase, the value of the additional data increases such that collecting additional behavioral data can be better than collecting a smaller amount of additional neuroimaging data. In some cases, collecting just 25 behavioral subjects can improve precision more than another neuroimaging subject. At the cost levels described above, this can mean that behavioral data can be a better investment than equivalently priced neuroimaging data. Therefore, a design analysis that takes costs into account has the potential to recommend collecting behavioral data rather than neuroimaging data. However, as shown in Figure 4, that is not true across the whole parameter range and, perhaps in particular, it isn't true for low sample sizes and low true correlation values, which is exactly where we might hope to gain in terms power. Why is is that the Anderson estimator performs worse here? To answer this question, we can decompose the mean squared error into the bias and the variance, i.e.

M SE(ρ, ρ)
where Bias is defined as As with mean squared error, we can also look at how various methods perform in terms of bias, i.e. is the estimate of ρ equal to ρ on average? It is well known that Pearson's r is biased downward such that, on average, correlation estimates will be lower than their true value (Olkin and Pratt, 1958). In Figure 5, we can see that the Anderson estimator of ρ always decreases the bias in estimating ρ, relative to an N b of 0 where it is equivalent to the Pearson r. This means that the increase in mean squared error in the Anderson estimate for lower values of ρ and N nb is due to the behavioral data adding variance to the estimate. It may be that the adjustment shown in Figure 3 is only likely to be in the right direction when there is enough data to estimate ρ and σ b precisely. We hope that future statistical work will address this question and potentially develop new estimators that can control this variance better in the low sample size, low correlation regime. For now, the choice of whether to use the Anderson (1957) estimator will depend on the research question. One somewhat counter-intuitive consequence of the above discussion is that Anderson's estimator can help not only in studies of individual differences across subjects but also in studies using representational similarity analyses where the matrix of item-item similarities based on a model is correlated with the similarities based on the multivariate neural signal (Kriegeskorte, 2008;Edelman et al., 1998). If researchers simply increase the number of similarities computed for the model, the above analyses suggest that there may be cases where this increases the power to detect a correlation with neural measures. In practice, since representational similarity matrices are often created from behavioral ratings (e.g Chikazoe et al., 2014;Bruffaerts et al., 2013), this is yet another way in which collecting additional behavioral data can increase power. While it is common in this literature to use Spearman's rank correlation to assess similarities, we can put this within the linear correlation framework by first ranking the two variables and then computing the correlation. For Pearson's r, this is exactly equivalent to computing the Spearman correlation.

How to Decrease Estimation Error: Improving Within-Subect Model Estimates by Collecting Data from Additional Subjects
We now investigate how collecting additional behavioral data can also help reduce estimation error, making it applicable to the regression analyses described above. In our previous description of estimation error, we were agnostic as to how exactly parameters were estimated from behavioral data. The most common method of estimating cognitive model parameters for use in a modelbased regression is to use maximum likelihood estimation (Myung, 2003). The maximum likelihood estimate of parameters θ i for subject i is the set of parameters in a cognitive model C that maximizes the probability of that subject's behavioral data b i i.e.θ One problem is that many neuroimaging experiments are relatively short so not many trials can be collected per subject. In complex cognitive models, when the number of data points is small, many sets of parameters can be consistent with the data and, therefore, two datasets generated from the cognitive model with the same θ i could result in very different maximum likelihood estimatesθ i , i.e. ρ(θ i , θ i ) can be low on average. As demonstrated in the above estimation error section, this will result in lower power for estimating correlations and regression slopes when attempting to relate such variables to neural recordings.
Two methods have been used in the literature to increase ρ(θ i , θ i ) and ρ(ĉ it , c it ): sharing information across subjects or constraining the possible pa-rameter values. One way to use information from other subjects data is to simply assume that all subjects have the same θ values which means that they can be estimated using all of the behavioral data collected in the experiment (Daw, 2011). This increases the data used to fit the parameters by a factor of n. The c it will still vary by subject because of the dependence on the stimuli in the experiment. If the variability in true parameters θ i across subjects is small relative to the estimation error, this can be an effective, heuristic way to increase ρ(ĉ it , c it ). However, this is not an option in a correlation analysis when we are often interested the variation across subjects itself. In addition, many cognitive processes do have large individual differences so this may actually decrease ρ(ĉ it , c it ).
Another option is to constrain parameter values by adding a prior, P (θ i |φ). If φ is chosen correctly, using a prior can decrease the likelihood of implausible parameter values such as a learning rate of 0 (Daw, 2011). We can then find the parameters that maximize the posterior, i.e.
This is known as the maximum a posteriori (or MAP) estimation. By using the prior, we can force θ i to be closer to more plausible values. However, it can be difficult to know in general if your priors are constraining the parameters to be in the correct space. If your priors are wrong, this could also decrease ρ(θ i , θ i ).
One principled way to incorporate information from other subjects and use priors is hierarchical modeling (Gelman, 2006). This is an increasingly popular method in cognitive science and cognitive neuroscience (Ahn et al., 2011;Rouder and Lu, 2005) that assumes that there is a potentially nonuniform distribution P (θ i |φ) of parameters in the population. Functionally, this constrains the possible values of θ i as in the prior above but because this is describing the population distribution, we can now estimate φ from all of the subjects data rather than setting φ by hand. The resulting estimate of θ i will be a weighted average of the estimate from the population with the estimate that just uses an individual's data. The weights are determined by the strength of the individual data. For a subject whose data provides less constraint on the parameters of the cognitive model (i.e. if they had fewer trials or made more inconsistent choices), their estimates will be closer to the prediction from the prior. If a subject made very consistent choices, we might already have a lot of information about the parameters of their cognitive processes and we do not need to use as much information from the population. Thus, hierarchical modeling provides a way to decide from the data how much we should pool information from the population in creating each individual estimate. Katahira (2016) demonstrates that estimates from a hierarchical model strictly dominate both individual maximum likelihood estimates and population maximum likelihood estimates in terms of bias and variance.
While there exist frequentist methods for estimating hierarchical models (such as restricted maximum likelihood as implemented in popular packages for fitting hierarchical linear models like lmer (Bates et al., 2015)), for arbitrary cognitive models, fitting in this way often requires a complicated derivation or numerical integration which can be challenging in high dimensions. With modern Bayesian inference tools like Stan (Carpenter et al., 2017), it is usually much more straightforward to place a prior on φ and use Bayesian inference to compute a joint posterior over φ and θ. The posterior for the hierarchical model is then In addition to Stan, many R or python packages have made hierarchical Bayesian versions of specific cognitive models particularly easy to fit (e.g. (hddm (Wiecki et al., 2013) for drift diffusion models or hBayesDM (Ahn et al., 2017)) for many types of decision making models). One feature of hierarchical modeling is that individual hierarchical estimates can often be improved simply by sampling more data from the population. That is, the cognitive model parameter estimates for subjects whose data you already have, i.e. subjects who contributed neuroimaging data, can be improved by adding data from new subjects. This is because the error of the estimates of φ (or the standard deviation of the prior) will, for most models, converge towards the correct estimates with increasing the number of subjects.
To demonstrate this, we conduct a simulation with a very simple hierarchical model to show that additional behavioral data from new subjects can improve inference for latent cognitive parameters θ i for subjects i ∈ [1, (N nb + N b )], N nb of whom were collected in an original set neuroimaging experiment and N b who were collected in a separate behavioral experiment. We assume each subject's latent cognitive parameter θ i is drawn from a Gaussian population distribution with mean µ θ = 0 and variance σ 2 θ = 1. Each subject then provides behavioral data b it and to keep things simple, we assume that the cognitive state variables c it are essentially equal to b it such that where σ bi is drawn from a uniform distribution between .5 and 1 for each subject. Further simplifying, we summarise b it by its meanb i and assume that there are enough trials t that we can treat σ bi as known. For the hierarchical model, we need to put a prior on σ θi . While in practice it is often better to use a weakly or strongly informative prior, for generality, in these simulations, we assume an improper uniform prior, i.e. P (σ 1 ) ∝ 1. Using the derivation in (Gelman et al., 2013), we construct a grid approximation to the posterior for θ given b. We now investigate how the correlation between the true subject parameters θ i and their posterior mean estimateθ i for the N nb neuroimaging subjects changes as we add a number of new behavior-only subjects, N b , to the experiment. Note that we are only computing the correlation within the original N nb subjects. As we demonstrated in the section on estimation error, this correlation is what limits the power to detect true correlations or regression slopes between model estimates and neural signals. Figure 6 shows that across a range of N nb values, the addition of N b new behavioral subjects can increase the correlation between model estimates for the N nb neuroimaging subjectsθ [1:N nb ] and true parameters θ [1:N nb ] . This implies that one way to increase ρ(θ i , θ i ) (and therefore also increase ρ(ĉ it , c it )) for subjects in a neuroimaging experiment is to simply collect behavioral data from more subjects. However, the exact rate at which more subjects improve the model estimates will in general depend on the distribution of uncertainty in the individual θ i as well as the particular form of the model you use.
Beyond it's use in improving hierarchical estimates of cognitive model parameters, collecting additional behavioral data can improve neuroimaging regression inferences in several other ways. While not computationally feasible yet for most neuroimaging analyses, joint modeling (Turner et al., 2013) is a recently proposed technique in which the cognitive model and the neuroimaging model are jointly fit, hierarchically, to both the neuroimgaging and behavioral data. This can allow for even more efficient use of information by allowing both sources of data to inform the cognitive model parameters. In addition, it allows the uncertainty in the cognitive model parameters to be propagated into the estimates of the regression or correlation, resulting in more powerful estimates without the biases mentioned above. Turner et al. (2016) demonstrated how joint modeling can be used to improve multimodal inferences even when we don't have data for every mode for every row in the dataset (i.e. some subjects' data were collected with EEG and others with fMRI). This can be generalized to the case we discuss where we only have rows in the dataset with only one type of data (e.g. behavior).
In addition, to improving estimates within a particular model, collecting additional behavioral data can be useful in many other ways for improving model estimates. For instance, collecting additional behavioral data in experimental conditions or on novel stimuli not in the neuroimaging experiment could assist in model selection. Choosing the correct model could potentially have a much larger impact on the correlation between model estimates and true parameters.
Finally, we can view reducing estimation error as relevant not only for modelbased regression but also for regression approaches where the stimulus features are subjective or latent, such as the valence of a stimulus. The covariates in this case are often created from ratings based either on subjects in the experiment or outside behavioral subjects. Similar to their use with cognitive models, hierarchical models with additional behavioral data can reduce error in these stimulus ratings as well. Even just aggregating the average rating however, will be improved through greater amounts of data per item. Indeed, going back to at least Kent and Rosanoff (1910), it is already quite common in cases like this to use large "norming" studies to understand stimuli, particularly for word stimuli.

When does it work?
Having demonstrated that it is possible to reduce estimation error in a neuroimaging dataset by collecting more behavioral data, we might now wonder when this is likely to help. As we mentioned above, the amount that hierarchical modeling can help will in general depend on the model and the experimental design. In addition, collecting additional behavioral data can improve estimates in other ways such as through model selection. In order to say something general about how improvements in estimation improve statistical power to find neuroimaging effects, we will simply assume that it is possible to increase the correlation between model estimates and true parameters using behavioral data by a certain amount and investigate how that changes estimates about the relationship with neural data. To begin, we will build on the derivation from Wilson and Niv (2015) where they showed how the least squares estimate of the regression slope is impacted by misestimation (i.e., equations 4 and 5). Wilson and Niv (2015) only solved this for the case where the analyst is trying to infer a single regression slope from the data (i.e. a non-hierarchical regression model). But in model-based neuroimaging analyses, it is rarely assumed that all subjects have the same effect size. Therefore, the most common analysis method in model-based analysis is to use a two stage approach where a first-stage regression is fit to each subject (or each run) and then the subject regressions are aggregated at the group-level by approximating a hierarchical model, testing whether the average effect is significantly different from zero (Beckmann et al., 2003;Woolrich et al., 2004). This means that the benefit of collecting an additional neuroimaging subject cannot be captured by the Wilson and Niv (2015) equations alone. Following Friston et al. (2002) and Beckmann et al. (2003), we assume a Gaussian hierarchical model of population effects: where the noise distributions are both Guassian: Using the same notation as above, for subject i on trial t, c it is that subject's cognitive state at that time and n it is the associated neural signal. β ci represents the relationship between the cognitive state and the neural signal for a particular subject and β represents that average relationship in the population. σ 2 G represents the variance of the magnitude of the relationship in the population and σ 2 represents the additional variance in the neural signal that is not due to the cognitive state (either due to noise in the measurement or the signal being related to other unrelated cognitive processes as well.) From Wilson and Niv (2015), we know the distribution of estimates in the first level regression, i.e.
Given the population distributions above and using standard formulas for linear combinations of Gaussian random variables, we can get We can now write the distribution of first level estimates aŝ The parameter we are interested in is β G , the average effect in the population. While there are more efficient but computationally expensive estimators (Beckmann et al., 2003;Woolrich et al., 2004), the most common way to estimate β G is the OLS method, i.e.β β G is then just the mean of independent Gaussian random variables and will therefore have a Gaussian distribution with mean ρ(c i ,ĉ i )β G and variance The first term is a constant. In the second term, because β c is a Gaussian random variable, the normalized squared sum has a non-central chi-squared distribution, and a non-central chi-squared multiplied by a constant will have a generalized chi-squared distribution. Plugging in the mean of that distribution, the expectation of the variance is now: Combining this with the expected estimate, we can get an t-statistic at the expected variancê Rather than parameterizing this in terms of the raw effect size, β G , it can be useful for understanding to standardize the variables by the total standard deviation, σ 2 + σ 2 G , similar to the derivation in Wilson and Niv (2015). For clarity, we define two new variables. Following Wilson and Niv (2015), we call the standardized effect size, β G √ σ 2 +σ 2 G , the contrast-to-noise ratio (CNR). In multilevel modeling, an important statistic is the intra-class correlation coefficient or ICC which is the proportion of the total variance that is explained by the subject level variance, i.e. σ 2 G σ 2 G +σ 2 (Chen et al., 2018;Shrout and Fleiss, 1979). If we multiply the right hand side of the above equation by √ σ 2 +σ 2 G √ σ 2 +σ 2 G , we can now write the equation in terms of these two variables, i.e., Having derived the expected t-statistic, we can now compute the power under that t-statistic using standard formulas. This allows us to ask the question "how much would additional behavioral data have to improve our latent variable ICC is set to 2 3 corresponding to a between subjects variance that is twice as large as the neuroimaging noise. In general, power increases as the number of subjects increases, as the number of trials increases and as accuracy of the model estimates (the inverse of the estimation error) increases. estimates in order to improve power as much as another neuroimaging subject?" That is, what is the increase in ρ(c i ,ĉ i ) that is equivalent to increasing N nb by one? More formally, if P (t, N ) is the power of the one-sample t-test with a t-stat of t and N subjects, we can write and ask what value of δ ρ makes this equation true?
Since δ ρ must be between 0 and 1 − ρ(c i ,ĉ i ), this is a straightforward constrained root finding problem that we can solve exactly using numerical methods. We can solve this for several values of the other parameters (CNR, ICC, T and N nb ) using Brent's method (Brent, 1972) as implemented in the scipy package (SciPy 1.0 Contributors et al., 2020) in python 3.7.7. We plot the power as a function of CNR, T and N nb values in Figure 7.
Which settings of the parameters in equation 39 are reasonable for neuroimaging experiments? Wilson and Niv (2015) derive estimates for CNR from two different studies with a very wide range of .4 to 11. Their CNR definition was only normalized by σ and not the total standard deviation including the group variance as we have done here. Therefore, to translate those values to values in equation, we need to take into account the ICC. In for this simulation, we set the ICC to 2 3 , corresponding to a between subjects variance that is twice as large as the neuroimaging noise, which seems to be quite large. With that value of ICC, Wilson and Niv (2015)'s CNR of .4 corresponds to a CNR of .4 √ 3 which is approximately .231. We can see in Figure 7 that a value of .2 for the CNR leads to somewhat implausible power estimates for most cognitive neuroscience task-based studies (Button et al., 2013;Turner et al., 2018). Therefore, we will assume that the CNR is closer to .1 in the following simulation. However, with the above derivation, it is easy to compute this for other values of the generative model and task parameters. Also of note in these plots is that for a particular value of N nb the power function will asymptote as a function of ρ(c i ,ĉ i ). This means that given a good enough model estimate, it is impossible to improve power more by improving the estimate than by collecting more subjects. Figure 8 shows the solution to equation 39 in terms of δ ρ as a function of the correlation between estimate and true parameters with one more neuroimaging subject. Overall, as noted above, there is a range of parameter settings where improving model fit (potentially by collecting more behavioral subjects) cannot help improve the power more than collecting another neuroimaging subject. Wilson and Niv (2015) showed that, for a reinforcement learning model, it is fairly easy to obtain a high ρ(c i ,ĉ i ), even if the latent subject-level parameters are severely misestimated. This suggests that in that setting, model-fitting "isn't necessary", or at least improving the model fit may be unlikely to significantly affect power. However, this may not be the case for other cognitive models, although it seems unlikely that a model that is fit to the data will have a correlation as low as .2 or .3.
While in this paper we have assumed that data collected with and without a neuroimaging recording device (such as in an fMRI scanner or with an EEG cap on) are exchangeable, that is unlikely to be exactly true in practice. This does not make the described methods impossible to use but does potentially require more complex models that allow for differences between the data collected under different modalities. That being said, there is relatively little work documenting the magnitude of these purported differences. If these differences were found to be large, this would likely have implications for learning about cognition through neuroimaging in general, not just with the statistical methods described here. We strongly encourage both future empirical work documenting these differences and future statistical work on how to aggregate inferences across data collection modalities.
In future work, we hope to explore this in the context of other cognitive models. Assuming that correlations between are reasonably high and that a smaller δ ρ is easier to achieve, Figure 8 shows that improving the model fit with additional behavioral data is most useful when there is already a reasonably sized N nb and T is smaller. While even 50 subjects is already significantly above average for a neuroimaging study, most use numbers of trials that are much higher than 50. One thing that this suggests is that using extra behavioral data may be particularly useful for neuroimaging studies with small numbers of trials. These are not common because of the low power they achieve but methods like those mentioned in the previous section may make them more feasible. This opens the door to more neuroimaging studies that are difficult to run with high numbers of trials such as long-term memory studies where subjects are unlikely to be able to remember a large number of items over a long delay. Finally, in this section, we focus on inferring the conditional relationship between cognitive model parameters and neural signals, assuming that the linear model is correct. However, there are many situations where it might be statistically reasonable or even preferable to use a linear model even if it is known that the neural signals are not exactly a linear function of the cognitive parameters (White, 1980;Buja et al., 2019). Recent work by Azriel et al. (2021) has shown that in this case, inference for regression parameters can be improved by additional unlabeled data as we show in the above work. Future work should also explore the value of behavioral data for model-based fMRI in this approximation regime.

Discussion
In this paper, we have presented two ways in which collecting additional behavioral data outside of a neuroimaging experiment has the potential to improve the precision of an estimate of the relationship between brain measures and cognition. We also attempted to quantify a tradeoff between collecting each type of data, showing how this could potentially affect the design of future studies or reanalyses of past work. In both cases, behavioral data had the potential to help but not in all cases. Unfortunately, it was particularly unhelpful when both true effect sizes and the number of subjects were small, which is exactly  Figure 8: Solution to equation 39 as a function of parameters ρ(c i ,ĉ i ) for a range of T and N nb . CNR is set to .1 and ICC is set to 2 3 . Lines end when it is no longer possible to achieve better performance by increasing the correlation. the case when more power is needed. However, our analyses pointed to possible future statistical work, in particular the potential to improve the Anderson (1957) correlation estimate by appropriately reducing the variance. In addition, our analyses suggest new regimes where behavioral data could be helpful; for instance, the potential for improving model fits with behavioral data could allow for feasibly running neuroimaging studies with small numbers of trials. In many cases, such as in long term memory studies, it is not feasible to run large numbers of trials because participants cannot reasonably remember hundreds of items. Both of these directions will require more work to be put into practice but we hope this paper will inspire collaborations between cognitive neuroscientists and statisticians to work more on ways of fusing data sources in order to learn about brain and behavior correlations.
Finally, while we have mostly discussed collecting additional behavioral data in the context of designing a new study, it could also be a consideration in a study that resuses an old dataset. It could be that a model is not identifiable with a particular dataset but it could become identifiable by collecting behavioral data in more conditions. With the advent of more open neuroimaging data sets, and particularly with the current difficulty in collecting new data, these methods point to a way to learn new things from old datasets simply by collecting more cheap behavioral data.