Pitfalls in Post Hoc Analyses of Population Receptive Field Data

Data binning involves grouping observations into bins and calculating bin-wise summary statistics. It can cope with overplotting and noise, making it a versatile tool for comparing many observations. However, data binning goes awry if the same observations are used for binning (selection) and contrasting (selective analysis). This creates circularity, biasing noise components and resulting in artifactual changes in the form of regression towards the mean. Importantly, these artifactual changes are a statistical necessity. Here, we use (null) simulations and empirical repeat data to expose this flaw in the scope of post hoc analyses of population receptive field data. In doing so, we reveal that the type of data analysis, data properties, and circular data cleaning are factors shaping the appearance of such artifactual changes. We furthermore highlight that circular data cleaning and circular sorting of change scores are selection practices that result in artifactual changes even without circular data binning. These pitfalls might have led to erroneous claims about changes in population receptive fields in previous work and can be mitigated by using independent data for selection purposes. Our evaluations highlight the urgency for us researchers to make the validation of analysis pipelines standard practice. Graphical abstract Highlights Circular data binning produces artifactual changes in the form of regression towards the mean Analysis type, data properties, and circular data cleaning shape these artifactual changes Circular data cleaning and sorting produce artifactual changes even without circular data binning

Bin-wise fictive body weight data and means for Today and Tomorrow in the same group of adults and different data binning scenarios. Body weight data for Today and Tomorrow were either binned according to body weight data for Today (1 st column) or an Independent test occasion (2 nd column). Fictive body weight data were simulated by sampling the body weight of 1000 adults from a Gaussian distribution (M = 70 kg; SD = 10 kg) and disturbing each adult's body weight with random Gaussian noise (SD = 10 kg), separately for each test occasion (Today, Tomorrow, and Independent). The red horizontal lines indicate the location of the overall mean for Today and Tomorrow. Dark brown colors correspond to lower and dark blue-green colors to higher decile bins. The endpoints of the colorful lines represent individual data points and the colorful dots with the black outline bin-wise means. Note that the graphs displayed here are referred to as Galton squeeze diagrams (Campbell and Kenny, 1999;Galton, 1886;Shanks, 2017).
x-axis condition. Subsequently, we quantify for each voxel the position shift from the Baseline to 90 the Interest condition (see Figure 2 for a single pRF). Finally, we calculate the bin-wise mean 91 shift. This is equivalent to calculating the bin-wise simple means for each condition and  For each vertex (gray matter node on the cortical surface), we obtained 6 estimates: pRF 125 position (x 0 and y 0 coordinates), pRF size (σ), pRF baseline (β 0 ), pRF amplitude (β 1 ), 126 and goodness-of-fit (R 2 ). We first smoothed the resulting parameter maps and delineated 127 V1 hemifield maps manually (for more details, see Supplementary methods, 1. Retinotopic 128 mapping experiment). We then pooled the x 0 and y 0 coordinates across V1 hemifield maps 129 and removed empty data points. and Interest condition, a null effect with eccentricity-dependent noise, and a true effect. We 147 use the term 'condition cross-thresholding' to refer to the pair-wise or list-wise deletion of 148 3 Note that when evaluating data distributions with unequal means, variances, or non-linearity, z-standardization might be necessary to detect regression towards or away from the mean (Campbell and Kenny, 1999;Shanks, 2017). In particular, z-standardization makes data distributions directly comparable. As such, bin-wise means should regress to wherever they intersect the identity line. Here, we always display data in native space, as this is typically done in the pRF literature. However, we use crosshairs to indicate the location of the mean and thus provide a visual guideline.

216
We then performed a 1D binning analysis on eccentricity and a 2D binning analysis 217 on x 0 and y 0 as we did for the simulated data. However, here, the eccentricity bins for 218 the 2D analysis ranged from 0 to 18 dva with a constant bin width of 2 dva. All binning 219 analyses and visualizations (including those on simulated data) were implemented in Matlab 220 2016b (9.1; https://uk.mathworks.com/) using custom code (Data and code availability).

221
The color scheme used for color-coding was an adapted version of the BrBG palette from  To expose the regression artifact, we repeatedly perturbed the x 0 and y 0 values of an empirical 227 visual field map with random Gaussian noise to generate a Baseline and Interest condition. 228 We then converted the x 0 and y 0 values to eccentricity. Subsequently, we binned the eccen- for the condition that was used for contrasting and binning (henceforth circular condition).

241
On the contrary, for the other condition (henceforth non-circular condition), this is not the 242 case.

243
The skew in average noise renders the bin-wise eccentricity means of the circular condition

302
This is evidently true even without circular data binning. As such, the reason why the noise 303 components in our cross-thresholding scenarios are sometimes biased even for the non-circular 304 condition 6 ( Figure 5, B., 2 nd column as well as Figure 5, C., 1 st and 2 nd columns) is because 305 we introduced another layer of circularity.

306
The fact that circular cross-thresholding and circular data binning are somewhat distinct 307 but also highly similar issues can, for instance, be appreciated when comparing the overall 308 instead of the bin-wise means. Without circular cross-thresholding, the overall mean in both 309 6 For reasons of clarity and simplicity, we use the term 'circular condition' or 'non-circular condition' exclusively when referring to circular data binning. However, other circular selection procedures, such as circular data sorting or cleaning, might of course render a condition circular above and beyond circular data binning.   Figure 4, although here, original observations having smaller eccentricities (≥ 0 and < 3 dva) were disturbed by random Gaussian noise with a smaller standard deviation (SD = 0.25 dva) and those having larger eccentricities (≥ 3 dva) by random Gaussian noise with a larger standard deviation (SD = 2 dva). B. The same as in Figure 4, although here, we simulated a true effect, that is, a radial increase in eccentricity of 2 dva in the Interest as compared to the Baseline condition.  2014, 2020) was based on unbinned data, and thus the overall means. As such, the apparent 379 significant effect was likely driven by or inflated due to circular cross-thresholding.  also possible that only one but not the other produces artifactual changes. This potential 383 divergence adds another layer of complexity to the issues we discussed here.

384
Taken together, the heterogeneity in manifestation we exposed here makes it hard to spot  in the Independent condition, cross-thresholding did not bias the noise components in the 419 Independent condition. As such, although the simple bin-wise means in the Baseline and independent criteria as this is where selection starts. Importantly, the intricacies we just discussed do not only apply to circular data binning, but also circular sorting of change 460 scores and circular condition cross-thresholding.

461
The application of circular data binning, circular sorting of change scores, and/or circular

631
The canonical HRF we adopted is implemented in SamSrf7.

Supplementary figures
x-axis nterest Baseline Figure S1. Interpretation of changes in eccentricity. The same as Figure 2, although here, the pRF shifts from one visual field quadrant to another in the Interest compared to the Baseline condition. This can happen due to noise or when visual field maps partially cover the ipsilateral hemifield. In such cases, an increase or decrease in eccentricity does not necessarily correspond to an outwards or inwards shift in the traditional sense. For instance, imagine that a pRF sits at x0 = -0.18 dva and y0 = -0.18 dva in the Baseline condition (ecc = 0.25 dva) but at x0 = 0.36 dva and y0 = 0.36 dva in the Interest condition (ecc = 0.51 dva). This would result in an increase in eccentricity, which might be interpreted as an outwards shift, although the pRF shifts effectively radially inwards until it reaches the origin and then outwards. We can likewise imagine that the pRF shifts horizontally to x0 = 0.36 dva and y0 = -0.36 dva in the Interest condition. Importantly, removing such shifts would bias noise components (see condition cross-thresholding in the main text and Figure 5, B. and C. as well as Figure S2-S3) and therefore, does not seem a valid option. Dva = Degrees of visual angle.
Simulated null effect -Cross-thresholding (Baseline)  Figure 7, although here, we simulated a true effect, that is, a radial increase in eccentricity of 2 dva in the Interest as compared to the Baseline condition. Note that the eccentricity bins ranged from 0 to 22 dva here (instead of 0 to 20 dva).
A. Empirical repeat data | 25 th %ile | Dorsal    Figure S7. Empirical 1D post hoc binning analysis on eccentricity | Repeat data | 25 th %ile participant | Posterior -Without and with cross-thresholding. The same as in Figure S6, although here, we used data from the posterior complex (V1-V3).  Figure S8. Empirical 1D post hoc binning analysis on eccentricity | Repeat data | 75 th %ile participant | Dorsal -Without and with cross-thresholding. The same as in Figure S6, although here, we used the 75 th %ile participant of the median R 2 distribution.  Figure S9. Empirical 1D post hoc binning analysis on eccentricity | Repeat data | 75 th %ile participant | Posterior -Without and with cross-thresholding. The same as in Figure S7, although here, we used the 75 th %ile participant of the median R 2 distribution. Figure S10. Empirical 2D post hoc binning analysis on x0 and y0 | Repeat data | 25 th %ile participant | Dorsal -Without and with cross-thresholding. Bin-wise x0 and y0 means in the Interest and Baseline condition for repeat data from the HCP belonging to the 25 th percentile participant of the median R 2 distribution and different data binning scenarios. A. Data from the dorsal complex (V3A/B and IPS0-5) without condition crossthresholding. B. Same as A., but with condition cross-thresholding. To this end, eccentricity values falling outside a certain eccentricity range (≥ 0 and ≤ 8 dva) and below a certain R 2 cut-off (≤ 2.2%) in the Baseline condition were removed from both conditions. C. Same as B., although here, condition cross-thresholding was based on both the Baseline and Interest condition. The x0 and y0 values in the Baseline and Interest condition were either binned according to eccentricity and polar angle values in the Baseline (1 st column in A.-C.) or Interest (2 nd column in A.-C.) condition. The marginal histograms (bin width = 0.5 dva; y-axis: relative frequency) show the x0 and y0 distributions for each condition. Magenta histograms correspond to the Interest condition and gray histograms to the Baseline condition. Note that the range of the marginal y-axis is the same for all histograms. The large magenta dots (arrow tip) correspond to the means in the Interest condition and the endpoint of the gray line (arrow knock) to the mean in the Baseline condition. The gray line itself (arrow shaft) depicts the shift from the Baseline to the Interest condition. The magenta crosshair indicates the location of the overall x0 and y0 means for the Interest condition and the gray crosshair the location of the overall means for the Baseline condition. Note that for subtle differences between the Baseline and Interest condition, the histograms and crosshairs almost coincide (see Figure S11 and Figure S13). The light gray polar grid demarks the bin segments. Polar angle bins ranged from 0°to 360°with a constant bin width of 45°and eccentricity bins from 0 to 18 dva with a constant bin width of 2 dva. The maximal eccentricity of the stimulated visual field area subtended 8 dva. HCP = Human Connectome Project. Dva = Degrees of visual angle. %ile = Percentile. (Condition) cross-thresholding = The pair-wise or list-wise deletion of observations across conditions.