This notebook illustrates a possible mechanims through which the apparent difference in value coding between conditions in figure 5 of Seamans et al. (https://doi.org/10.1101/2024.06.07.597894) may arise in the absence of any actual difference in value coding, due to differences in the distribution of ival across conditions.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | import numpy as np import seaborn as sns from sklearn.linear_model import LinearRegression import pylab as plt # Generate two distributions of ival where one is biased towards high values and the other towards low values. ival1 = np.hstack([np.random.rand(500), 0.7+0.3*np.random.rand(500)]) # Biased to high values ival2 = np.hstack([np.random.rand(500), 0.3*np.random.rand(500)]) # Biased to low values # Generate two sets of neural activity which are linearly dependent on ival but noisy activity1 = ival1 + np.random.randn(1000) activity2 = ival2 + np.random.randn(1000) # Plot the activity against ival with linear fit. plt.figure(1) sns.regplot(x=ival1, y=activity1, label='high_ival') sns.regplot(x=ival2, y=activity2, label='low_ival') plt.ylabel('Activity') plt.xlabel('Ival') plt.legend(); |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | # Perform a regression predicting ival from activity. pred_ival1 = LinearRegression().fit(activity1.reshape(-1, 1), ival1).predict(activity1.reshape(-1, 1)) pred_ival2 = LinearRegression().fit(activity2.reshape(-1, 1), ival2).predict(activity2.reshape(-1, 1)) # Compute reiduals residuals1 = ival1 - pred_ival1 residuals2 = ival2 - pred_ival2 # Plot residuals as a function of ival. plt.scatter(ival1, residuals1, label='high_ival') plt.scatter(ival2, residuals2, label='low_ival') plt.xlabel('Ival') plt.ylabel('Residuals') plt.legend(); |
We can see that although the true statistical relationship ival and activity and ival is the same in both conditions, there is a systematic offet in the residuals when they are plotted as a function of ival, due to the difference in the distribution of ival values in the two datasets, as observed in figure 5C.
To see why this occurs it is useful to plot the linear regression predicting ival from activity, for which the above plot shows the residuals.
1 2 3 4 5 6 7 | # Plot the ival against activity with linear fit. plt.figure(1) sns.regplot(x=activity1, y=ival1, label='high_ival') sns.regplot(x=activity2, y=ival2, label='low_ival') plt.ylabel('Ival') plt.xlabel('Activity') plt.legend(); |