# "Uncorrelation" demo
From Soch, J., Allefeld, C., & Haynes, J. D. (2020). Inverse Transformed Encoding Models–a solution to the problem of correlated trial-by-trial parameter estimates in fMRI decoding. Neuroimage, 209, 116449.


```python
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
```

## Encoding models
Let'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}


```python
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 @ ϕ
```

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:


```python
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")
```

And plot second-level model fit ($R^{2}$):


```python
plt.hist(r2_naive)
plt.hist(r2_uncor)
plt.legend(['naive', 'uncorrelation'])
plt.xlabel('R2')
plt.ylabel('Freq')
```

## Decoding
We 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).


```python
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 @ ϕ
```

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}


```python
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()
```

6ae7d484dbb134e8935191c30e7432ea5cd6443d
ipynb
Jupyter Notebook
uncorrelation_simulation_demo.ipynb
lukassnoek/random_notebooks
d7df507ce2b6949726c29de0022aae2d0dc583ac
MIT
2018-05-28T13:45:11.000Z
2021-08-31T11:41:34.000Z
uncorrelation_simulation_demo.ipynb
lukassnoek/random_notebooks
d7df507ce2b6949726c29de0022aae2d0dc583ac
MIT
uncorrelation_simulation_demo.ipynb
lukassnoek/random_notebooks
d7df507ce2b6949726c29de0022aae2d0dc583ac
MIT
2018-05-28T13:46:05.000Z
2018-06-11T15:25:59.000Z
Qwen/Qwen-72B
1. YES
2. YES