<jupyter_start><jupyter_text>"Uncorrelation" demoFrom <NAME>., <NAME>., & <NAME>. (2020). Inverse Transformed Encoding Models–a solution to the problem of correlated trial-by-trial parameter estimates in fMRI decoding. Neuroimage, 209, 116449.<jupyter_code>import numpy as np
from scipy.linalg import toeplitz, sqrtm
from scipy.stats import pearsonr
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Encoding modelsLet's specify a two-stage data simulation scenario. We assume that a set of $P$ trials across $N$ time points &mdash; summarized in a "design matrix" $X_{1}$ &mdash; "activate" a voxel according to weights $\beta$, plus some additive noise $\epsilon$, resulting in a particular $N \times 1$ time series $y$:\begin{align}y = X_{1}\beta + \epsilon\end{align}Importantly, there might be some inherent trial-by-trial correlation between the $P$ variables, which is common in single-trial fMRI analyses (as we assume here), which we summarize in a $P \times P$ covariance matrix $U$, generated through an AR1 process with a particular base correlation $\rho$.Now, we also assume that the "activations", $\beta$, are in turn generated by a set of $Q$ properties of, e.g., our stimulus set &mdash; summarized in a "design matrix" $X_{2}$ &mdash; weighted by a set of parameters $\phi$:\begin{align}\beta = X_{2}\phi\end{align}<jupyter_code>N = 600
P = 120
Q = 5
ρ = 0.8
U = ρ ** toeplitz(np.arange(P))
σ2 = 2

X1 = np.random.multivariate_normal(np.zeros(P), U, size=N)
X2 = np.random.normal(0, 1, size=(P, Q))
ϕ = np.random.normal(0, 0.1, size=Q)
β = X2 @ ϕ<jupyter_output><empty_output><jupyter_text>We can run a "first-level" model to estimate the "activation" parameters $\hat{\beta}$ (assuming uncorrelated errors, $\epsilon$):\begin{align}\hat{\beta} = (X_{1}^{T}X_{1})^{-1}X_{1}^{T}y\end{align}To estimate the set of parameters $\hat{\phi}$, we can run a "second-level" model in which we, for now, assume uncorrelated errors:\begin{align}\hat{\phi} = (X_{2}^{T}X_{2})^{-1}X_{2}^{T}\hat{\beta}\end{align}However, as is shown in the Soch et al. paper, it's more efficient to incorporate knowledge about the trial-by-trial covariance ($U$) as induced by the first-level design matrix $X_{1}$ using generalized least squares:\begin{align}\hat{\phi} = (X_{2}^{T}U^{-1}X_{2})^{-1}X_{2}^{T}U^{-1}\hat{\beta}\end{align}where, in practice, $U$ needs to be estimated from the design matrix:\begin{align}U = (X_{1}^{T}X_{1})^{-1}\end{align}Alternatively, it is possible to "prewhiten" the second-level design matrix ($X_{2}$) and the second-level data ($\hat{\beta}$):\begin{align}\hat{\beta}^{*} = U^{-\frac{1}{2}}\hat{\beta} = (X_{1}^{T}X_{1})^{\frac{1}{2}}\hat{\beta} \\X_{2}^{*} = U^{-\frac{1}{2}}X_{2} = (X_{1}^{T}X_{1})^{\frac{1}{2}}X_{2}\end{align}which outputs can subsequently be used in an ordinary OLS model (which is equivalent to running GLS as outlined above):\begin{align}\hat{\phi} = (X_{2}^{*T}X_{2}^{*})^{-1}X_{2}^{*}\hat{\beta}\end{align} Below, we simulate our data ($y$) iteratively and estimate the second-level parameters $\hat{\phi}$ with and without consideration of the first-level trial-by-trial covariance ($U$), showing that the former leads to a more efficient estimation of the second-level parameters and to a better fit in general:<jupyter_code>model = LinearRegression()
iters = 100

ϕ_naive = np.zeros((iters, Q))
ϕ_uncor = np.zeros((iters, Q))
r2_naive = np.zeros(iters)
r2_uncor = np.zeros(iters)
for i in range(iters):
    ϵ = np.random.normal(0, σ2, size=N)
    y = X1 @ β + ϵ

    βh = np.linalg.lstsq(X1, y, rcond=None)[0]
    ϕ_naive[i, :] = np.linalg.lstsq(X2, βh, rcond=None)[0]
    r2_naive[i] = cross_val_score(model, X2, βh, cv=6).mean()
    
    U_inv = X1.T @ X1
    U_sqrt = sqrtm(U_inv)
    βh_pw = U_sqrt @ βh
    X2_pw = U_sqrt @ X2
    ϕ_uncor[i, :] = np.linalg.inv(X2_pw.T @ X2_pw) @ X2_pw.T @ βh_pw
    # Line below: GLS formulation of uncorrelation (mathematically the same as line above)
    #ϕ_uncor[i, :] = np.linalg.inv(X2.T @ U_inv @ X2) @ X2.T @ U_inv @ βh
    
    r2_uncor[i] = cross_val_score(model, X2_pw, βh_pw, cv=6).mean()
    
plt.figure(figsize=(15, 10))
plt.plot(ϕ, c='k', ls='--', lw=3)
plt.plot(ϕ_naive.mean(axis=0))
plt.fill_between(np.arange(Q),
                 ϕ_naive.mean(axis=0) + ϕ_naive.std(axis=0),
                 ϕ_naive.mean(axis=0) - ϕ_naive.std(axis=0))
plt.fill_between(np.arange(Q),
                 ϕ_uncor.mean(axis=0) + ϕ_uncor.std(axis=0),
                 ϕ_uncor.mean(axis=0) - ϕ_uncor.std(axis=0))
plt.legend(['true', 'naive', 'uncorrelation'])
plt.ylabel("Parameter value")
plt.xlabel("Parameter index")<jupyter_output><empty_output><jupyter_text>And plot second-level model fit ($R^{2}$):<jupyter_code>plt.hist(r2_naive)
plt.hist(r2_uncor)
plt.legend(['naive', 'uncorrelation'])
plt.xlabel('R2')
plt.ylabel('Freq')<jupyter_output><empty_output><jupyter_text>DecodingWe can also apply uncorrelation in a decoding analysis, where we try to predict $X_{2}$ from the estimated activity ($\hat{\beta}$) of a pattern of $K$ voxels. Here, we simulate some data in which $X_{2}$ is a matrix with two categorical variables (i.e., class 0 and class 1, e.g., faces vs. houses).<jupyter_code>K = 50
Q = 2
X2 = np.zeros((P, 2))
#X2[:(P // 2), 0] = 1
#X2[(P // 2):, 1] = 1
X2[::2, 0] = 1
X2[1::2, 1] = 1
ϕ = np.random.normal(0, 0.1, size=(Q, K))
β = X2 @ ϕ<jupyter_output><empty_output><jupyter_text>Below, we iteratively simulate data for 50 voxels:\begin{align}Y = X_{2}\beta + \epsilon\end{align}where $Y$ is an $N \times K$ matrix. Then, in a 10-fold cross-validation analysis, we fit either a "naive" (OLS-based) decoding model, which estimates weights $\hat{W}$:\begin{align}\hat{W} = (X_{2}^{T} X_{2})^{-1}X_{2}^{T}\hat{\beta}\end{align}or a decoding model with "uncorrelation" (effectively a GLS model):\begin{align}\hat{W} = (X_{2}^{T}U^{-1}X_{2})^{-1}X_{2}^{T}U^{-1}\hat{\beta}\end{align}where $U^{-1}$ refers to the covariance matrix of the first-level design matrix:\begin{align}U^{-1} = X_{1}^{T}X_{1}\end{align}<jupyter_code>from sklearn.model_selection import StratifiedKFold
cv = StratifiedKFold(n_splits=10)
target = X2.argmax(axis=1)
iters = 200

Uinv = X1.T @ X1

acc_naive = np.zeros(iters)
acc_uncorr = np.zeros(iters)
for i in range(iters):
    ϵ = np.random.normal(0, 2, size=(N, K))
    Y = X1 @ β + ϵ
    βh = np.linalg.inv(X1.T @ X1) @ X1.T @ Y
    for train_idx, test_idx in cv.split(βh, target):
        X_train, X_test = βh[train_idx, :], βh[test_idx, :]
        y_train, y_test = target[train_idx], target[test_idx]
        Wh = np.linalg.inv(X_train.T @ X_train) @ X_train.T @ y_train
        preds = np.round(X_test @ Wh)
        acc_naive[i] += np.mean(preds == y_test)
 
        Uinv_c = Uinv[train_idx, :][:, train_idx]
        Wh = np.linalg.inv(X_train.T @ Uinv_c @X_train) @ X_train.T @ Uinv_c @ y_train
        preds = np.round(X_test @ Wh)
        acc_uncorr[i] += np.mean(preds == y_test)
        
    acc_naive[i] /= cv.get_n_splits()
    acc_uncorr[i] /= cv.get_n_splits()
    
plt.hist(acc_naive)
plt.hist(acc_uncorr)
plt.legend(["naive", "uncorr"])
plt.xlabel("Accuracy")
plt.ylabel("Freq")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17178502879078694, 'fraction_numerical': 0.0255004112969564, 'mean_word_length': 4.493222891566265, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '396950', 'license': 'MIT', 'n_tokens_mistral': 2961, 'n_tokens_neox': 2762, 'n_words': 710, 'path': 'uncorrelation_simulation_demo.ipynb', 'repo_name': 'lukassnoek/random_notebooks', 'size': '7306'}