2. Using the GLM to fMRI: estimation

2.1 Recap of linear regression

2.1.1 Notation

• (선형) 회귀 모델을 생성할 때, 우리는 종속 변수와 독립 변수를 지정함

• 위의 식에서 y는 종속 변수, x는 독립 변수임 (N은 observations의 수, P는 predictors의 수임)

np.random.seed(1)

# 가상 데이터셋 만들기
prespecified_covariance = np.array([[1, .8], [.8, 1]]) # y와 x는 0.8의 상관관계 지님

# x 평균은 3, y 평균은 7이면서, 서로의 상관관계는 0.8인 데이터 30개 만들기
data = np.random.multivariate_normal(mean=[3, 7], cov=prespecified_covariance, size=30)
X = data[:, 0, np.newaxis] # (30, 1)
y = data[:, 1, np.newaxis] # (30, 1)

# intercept modeling
n_obs = y.size
intercept = np.ones((n_obs, 1)) # (30, 1)

tuple_with_arrays = (intercept, X)
X_with_icept = np.hstack(tuple_with_arrays) # (30, 2)

2.1.2 Parameters in linear regression

• linear regression은 다음과 같이 표현함
• 여기서, 각 (독립) 변수가 종속변수에 영향을 미치는 정도를 beta-parameter라고 함

• beta-parameter를 추정하는 방법으로 Ordinary Least Squares(OLS, 최소제곱법) 사용함
• y_hat은 실제 y(종속변수)가 아니라 선형회귀식으로 얻은 종속변수의 추정치임
• 실제 y(종속변수)와 y_hat의 차이의 제곱을 줄이는 것을 목표로 함

from numpy.linalg import inv

est_betas = inv(X_with_icept.T @ X_with_icept) @ X_with_icept.T @ y

print(est_betas)
>>>
[[4.25897963] # parameter of intercept
 [0.88186203]] # parameter of x

2.1.3 Residuals and model fit

• residual은 다음과 같이 정의함
• 실제 종속변수에서 예측치를 뺀 값임

• 아래의 MSE, R^2는 회귀 예측 정도를 나타내는 척도임

mse = np.mean((y - y_hat) ** 2) # Mean Square Error

numerator = np.sum((y - y_hat) ** 2) # 관측값 - 추정값(SSR, Residual Sum of Squares)
denominator = np.sum((y - np.mean(y)) ** 2) # 관측값 - 관측값의 평균(SST, Total Sum of Squares)
r_squared = 1 - numerator / denominator

print('The R² value is: %.3f' % r_squared)
>>> The R² value is: 0.549 # 1에 가까울수록 좋은 성능

2.2 GLM in fMRI analyses

2.2.1 Variables

• 우리는 regression하기 위해 독립변수와 종속변수를 지정할 필요가 있음

• 독립변수는 (시각) 자극 시점이며, 독립변수로부터 예측하고자 하는 변수(종속변수)는 voxel signal임

• fMRI 데이터로는 discrete한 변수이므로 continuous하게 plot하면 다음과 같이 그릴 수 있음

voxel_signal = np.load('example_voxel_signal.npy')
plt.figure(figsize=(25, 5))
plt.plot(voxel_signal) # discrete value to continuous line
plt.show()

• 50초 간격으로 1초간 사각형 화면과 원 화면을 번갈아가며 자극을 줄 때 다음과 같이 표시할 수 있음

onsets_squares = np.array([10, 110, 210, 310, 410, 510, 610, 710], dtype=int) # 사각형 화면 자극
onsets_circles = np.array([60, 160, 260, 360, 460, 560, 660, 760], dtype=int) # 원 화면 자극

onsets_all = np.concatenate((onsets_squares, onsets_circles))
print(onsets_all)
>>> [ 10 110 210 310 410 510 610 710  60 160 260 360 460 560 660 760]

predictor_all = np.zeros(800)
predictor_all[onsets_all] = 1
print(predictor_all)
>>>
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 1., ...]

• 우리가 만든 독립변수(stimulus_predictor, 800)는 실제 신호 시간(400)과 다른 시간 단위를 가짐
• fMRI는 대략 2초동안 아래부터 위까지 층으로 올라가면서 스캔하기 때문임

• 여기에서는 2초마다 하나의 데이터 포인트로 down-sampling 통해 해결 가능

• Python에서 scipy.interpolate 패키지 안의 interp1d 함수를 통해 리샘플링 가능
• interp1d 함수 사용 시에 resampler 변수 안에 원하는 스케일을 저장하면 그 스케일에 맞춰 리샘플링 가능
• (+ 실제로는 정밀하게 계산하기 위해 밀리초 규모로 사용함, [initial upsampling of predictors] 부분 확인할 것)

from scipy.interpolate import interp1d

original_scale = np.arange(0, 800, 1)  # from 0 to 800 seconds
print(original_scale)
>>> 800

resampler = interp1d(original_scale, predictor_all)

desired_scale = np.arange(0, 800, 2)
print(desired_scale)
>>> 400 # (0, 2, 4, ... 800, in volumes)

predictor_all_ds = resampler(desired_scale)

• resampling을 거친 독립변수와 종속변수를 다음과 같이 시각화할 수 있음

plt.figure(figsize=(25, 10))
plt.plot(voxel_signal)
plt.plot(predictor_all_ds + voxel_signal.mean(), lw=3)
plt.show()

2.2.2 Regression on fMRI data & interpretation parameters

• 우리의 예시 데이터를 활용해 회귀식을 만들었을 경우 다음과 같은 결론을 얻을 수 있음

def test_regression_signal_simple(mse_ans, r2_ans, X, y):
    X = np.hstack((np.ones((X.shape[0], 1)), X))
    b = np.linalg.lstsq(X, y, rcond=None)[0]
    y_hat = X.dot(b)
    c_ans_mse = np.mean((y - y_hat) ** 2)
    c_ans_r2 = pearsonr(y, y_hat)[0] ** 2
    try:
        np.testing.assert_almost_equal(mse_ans, c_ans_mse)
    except AssertionError as e:
        try:
            print('Your MSE is %.3f; I expected %.3f' % (mse_ans, c_ans_mse))
            raise e
        finally:
            e = None
            del e

    else:
        print('MSE is correct')
    try:
        np.testing.assert_almost_equal(r2_ans, c_ans_r2)
    except AssertionError as e:
        try:
            print('Your R2 is %.3f; I expected %.3f' % (r2_ans, c_ans_r2))
            raise e
        finally:
            e = None
            del e

    else:
        print('R2 is correct')
    print('Well done!')

test_regression_signal_simple(mse_simple, r2_simple, predictor_all_ds, voxel_signal)

def plot_signal_and_predicted_signal(y, X, x_lim=None, y_lim=None):
    """ Plots a signal and its GLM prediction. """
    des = np.hstack((np.ones((y.size, 1)), X))
    betas_simple = np.linalg.lstsq(des, y, rcond=None)[0]
    plt.figure(figsize=(15, 5))
    plt.plot(y)
    plt.plot(des @ betas_simple, lw=2)
    plt.xlabel('Time (in volumes)', fontsize=15)
    plt.ylabel('Activity (A.U.)', fontsize=15)
    
    if x_lim is not None:
        plt.xlim(x_lim)
    
    if y_lim is not None:
        plt.ylim(y_lim)
    
    plt.legend(['True signal', 'Predicted signal'], loc='upper right', fontsize=15)
    plt.title("Signal and predicted signal", fontsize=25)
    plt.grid()
    plt.show()

plot_signal_and_predicted_signal(voxel_signal, predictor_all_ds, x_lim, y_lim)

2.2.3 Using the BOLD response in GLM models

• 우리가 현실에서 BOLD 신호의 즉각적 변화를 목격하는 것은 비현실적임
• 진정한 신경 활동은 약 5-7초 후에 최고조에 달하기 때문임

• 따라서 우리는 canonical HRF와 convolution 계산을 통해서 BOLD 신호와 유사한 그래프를 만들 것임

• fMRI 데이터의 BOLD 신호 단변량 분석(univariate analyses): 복셀은 항상 같은 방식으로 자극에 반응한다고 가정 → using “canonical haemodynamic response function (HRF)”

• canonical HRF: formalization of how the a voxel is going to respond to a stimulus

• 가장 많이 사용하는 HRF function은 double-gamma HRF임(하나는 overshoot을 모델링하고 다른 하나는 undershoot을 모델링함)

TR = 2 # 스캔 지속 시간
osf = 2 # 오버샘플링
length_hrf = 32  # sec

canonical_hrf = glover_hrf(tr=TR, oversampling=osf, time_length=length_hrf, onset=0) # glover_hrf는 double-gamma HRF 함수임
canonical_hrf /= canonical_hrf.max()
print(canonical_hrf.size)
>>> 32

t = np.arange(0, canonical_hrf.size)
plt.plot(t, canonical_hrf)

• HRF와 자극 사건의 콘볼루션 연산은 자극과 BOLD response의 지연되는 모양을 예측할 수 있음

predictor_conv = np.convolve(predictor_all.squeeze(), canonical_hrf)
print(predictor_conv)
>>> (831,)

predictor_conv = predictor_conv[:predictor_all.size]
print(predictor_conv.shape)
>>> (800,)

predictor_conv = predictor_conv[:, np.newaxis]
print(predictor_conv.shape)
>>> (800, 1)

plt.plot(predictor_all)
plt.plot(predictor_conv)

• BOLD response를 (주어진 형태의) 선형으로 확장하는 것

plt.figure(figsize=(15, 5))

N = 40
one_stim = np.zeros(N)
one_stim[5] = 1
one_stim_conv = np.convolve(one_stim, canonical_hrf)[:N]
two_stim = np.zeros(N)
two_stim[[5, 6]] = 1
two_stim_conv = np.convolve(two_stim, canonical_hrf)[:N]
three_stim = np.zeros(N)
three_stim[[5, 6, 7]] = 1
three_stim_conv = np.convolve(three_stim, canonical_hrf)[:N]

plt.plot(one_stim_conv)
plt.plot(two_stim_conv)
plt.plot(three_stim_conv)
plt.show()

• convolved predictor와 BOLD signal의 스케일 맞추기 위해 resampling을 다시 함

original_scale = np.arange(0, 800)
resampler = interp1d(original_scale, np.squeeze(predictor_conv))

desired_scale = np.arange(0, 800, 2)
predictor_conv_ds = resampler(desired_scale)

• 자극은 딱 떨어지는 초 단위로(10s 등) 나타나지 않음

• 보다 정확한 초 단위로 자극 예측기를 만들기 위해 100으로 오버샘플링을 함

onsets2 = np.array([3.62, 16.26, 34.12, 42.98])  # 정밀한 자극 시간
duration = 50  # 자극 간 사이 시간
osf = 100  # osf = OverSampling Factor

pred2 = np.zeros(duration * osf)
onsets2_in_msec = (onsets2 * osf).astype(int)
pred2[onsets2_in_msec] = 1

osf = 100 * TR # TR = 2
hrf_ms = glover_hrf(tr=TR, oversampling=osf, time_length=32)
hrf_ms /= hrf_ms.max()

t = np.arange(0, hrf_ms.size)

plt.figure(figsize=(12, 5))
plt.plot(t, hrf_ms)

pred2_conv = np.convolve(pred2, hrf_ms)[:pred2.size]
plt.plot(pred2_conv)

• Initial upsampling을 하는 경우 다시 downsampling을 통해 차원의 크기를 맞춤

resampler = interp1d(original_scale, pred2_conv)
desired_scale = np.arange(0, 50, TR)
pred2_conv_ds = resampler(desired_scale)

intercept = np.ones((predictor_conv_ds.size, 1))
X_conv = np.hstack((intercept, predictor_conv_ds))
betas_conv = inv(X_conv.T @ X_conv) @ X_conv.T @ voxel_signal

plt.subplot(2, 1, 1)
plt.plot(voxel_signal)
plt.plot(X_conv @ betas_conv, lw=2)

X_simple = np.hstack((icept, predictor_all_ds))
betas_simple = np.array([1000.64701684, 1.02307437])

plt.subplot(2, 1, 2)
plt.plot(voxel_signal)
plt.plot(X_simple @ betas_simple, lw=2)

print(betas_conv[1], betas_simple[1])
>>> (8.181, 1.023)

plot_signal_and_predicted_signal_zoom(voxel_signal, predictor_conv_ds, x_lim=(25, 65), y_lim=(995, 1010))

from numpy.linalg import lstsq  # numpy implementation of OLS, because we're lazy

y_hat_conv = X_conv @ betas_conv
y_hat_orig = X_simple @ lstsq(X_simple, voxel_signal, rcond=None)[0]

MSE_conv = ((y_hat_conv - voxel_signal) ** 2).mean()
MSE_orig = ((y_hat_orig - voxel_signal) ** 2).mean()

print(MSE_conv, MSE_orig)
>>> (6.022, 10.327)

R2_conv = 1 - (np.sum((voxel_signal - y_hat_conv) ** 2) / np.sum((voxel_signal - voxel_signal.mean()) ** 2))
R2_orig = 1 - (np.sum((voxel_signal - y_hat_orig) ** 2) / np.sum((voxel_signal - voxel_signal.mean()) ** 2))

print(R2_conv, R2_orig)
>>> (0.41919, 0.00388)

• HRF 함수를 이용하면 원 신호보다 자극이 느리게 나타날 수 있음

• 이를 해결하기 위해 Temporal basis function을 이용함(ex. double-gamma basis functions)

• double-gamma basis functions을 1계 미분하면 시간의 지연을 조정할 수 있고(temporal), 2계 미분하면 너비의 차이를 조정할 수 있음(dispersion) (거리-속도-가속도의 관계랑 유사함)

tderiv_hrf = glover_time_derivative(tr=2, oversampling=2)
tderiv_hrf /= tderiv_hrf.max()

plt.plot(canonical_hrf, lw=3)
plt.plot(tderiv_hrf, c='tab:orange', lw=3)

y_tbs, X_tbs = simulate_signal(~)

# regression with HRF only
beta_canon = inv(X_canon_only.T @ X_canon_only) @ X_canon_only.T @ y_tbs
yhat_canon = X_canon_only @ beta_canon

# regression with HRF + temporal deriv.
beta_td = inv(X_tbs.T @ X_tbs) @ X_tbs.T @ y_tbs
yhat_td = X_tbs @ beta_td

Summary

Linear Model

BOLD response in GLM

canonical_hrf = glover_hrf(tr, oversampling, time_length, onset)
predictor_conv = np.convolve(predict_all, canonical_hrf)

original_scale = np.arange(0, 800)
resampler = interp1d(original_scale, np.squeeze(predictor_conv))

desired_scale = np.arange(0, 800, 2)
predictor_conv_ds = resampler(desired_scale)

onsets2 = np.array([3.62, ...])
duration = 50
oversampling = 100
pred2[onsets2_in_msec] = 1

tderiv_hrf = glover_time_derivative(tr=2, oversampling=2)
tderiv_hrf /= tderiv_hrf.max()