3. Design of experiments

date

Aug 5, 2023

slug

fmri-glm-3

author

status

Public

Types of designs

Factorial designs

지난 시간에 배웠던 t-test에서 활용되는 contrast vector를 세워보는 연습을 해보자

(1) red squares, (2) red circles, (3) green squares, (4) green circles로 자극이 주어졌다고 가정해보자

exp_length = 160
TR = 1 # TR이 1인 경우에는 다운샘플링을 무시할 수 있음
xmax = exp_length // TR

y, X = simulate_signal(
    onsets=np.array(
        [0, 40, 80, 120,  # red squares
         10, 50, 90, 130,  # red circles
         20, 60, 100, 140,  # green squares
         30, 70, 110, 150]  # green circles
    ),
    conditions=['rq'] * 4 + ['rc'] * 4 + ['gs'] * 4 + ['gc'] * 4,
    TR=TR,
    duration=exp_length,
    icept=0,
    params_canon=[1, 2, 3, 4],
    rnd_seed=42,
    plot=False
)

X = X[:, :5]

plt.figure(figsize=(15, 5))
plt.plot(X[:, 1], c='r', ls='--')
plt.plot(X[:, 2], c='r')
plt.plot(X[:, 3], c='g', ls='--')
plt.plot(X[:, 4], c='g')
plt.legend(['red squares', 'red circles', 'green squares', 'green circles'], frameon=False)

plt.plot(y)
plt.xlim(0, xmax)

simulate_signal code

def simulate_signal(onsets, conditions, TR=2, duration=None, icept=0, params_canon=None, params_tderiv=None,
                    params_deriv1=None, phi=None, std_noise=1, osf=100, rnd_seed=None, plot=True):
    """ Simulates a somewhat realistic voxel signal with an associated design matrix. """
    if rnd_seed is not None:
        np.random.seed(rnd_seed)
    
    conds = sorted(np.unique(conditions))
    P = len(conds)
    
    if duration is None:
        duration = np.max(onsets) + 30
    
    if params_canon is None:
        params_canon = np.zeros(P)

    # TMP FIX
    params_tderiv = params_deriv1
    
    if params_tderiv is None:
        params_tderiv = np.zeros(P)

    X = np.zeros((duration * osf, P))
    for ons, con in zip(onsets, conditions):
        X[int(ons * osf), conds.index(con)] = 1
        
    hrf = glover_hrf(tr=1, oversampling=osf)
    hrf_d = glover_time_derivative(tr=1, oversampling=osf)

    t_orig = np.arange(0, duration, 1 / osf)
    t_new = np.arange(0, duration, TR)
    
    Xconv = np.zeros((t_orig.size, P * 2))
    idx = 0
    for i, this_hrf in enumerate([hrf, hrf_d]):
        this_hrf /= this_hrf.max()
        for ii in range(P):
            Xconv[:, idx] = np.convolve(X[:, ii], this_hrf)[:t_orig.size]  
            idx += 1
    
    Xconv = np.c_[np.ones(Xconv.shape[0]), Xconv]
    params = np.r_[icept, params_canon, params_tderiv]
    y = Xconv @ params
    
    resampler = interp1d(t_orig, y)
    y = resampler(t_new)       
    
    if phi is None:
        noise_cov = std_noise ** 2 * np.eye(y.size)
    else:
        noise_cov = std_noise ** 2 * phi ** toeplitz(np.arange(y.size))
    
    y = y + np.random.multivariate_normal(np.zeros(y.size), noise_cov)

    Xconv_ds = np.zeros((t_new.size, Xconv.shape[1]))
    for i in range(Xconv.shape[1]):
        resampler = interp1d(t_orig, Xconv[:, i])
        Xconv_ds[:, i] = resampler(t_new)
    
    est_betas = np.linalg.lstsq(Xconv_ds, y, rcond=None)[0]
    if plot:
        plt.figure(figsize=(15, 5))
        plt.plot(y)
        plt.plot(Xconv_ds @ est_betas)
        plt.xlim(0, y.size)
        plt.legend(['y', 'y-hat'])
        plt.show()

    return y, Xconv_ds

Questions

Question 1.

GLM 모델을 만들어라.

print(y.shape, X.shape) # X: intercept 포함되어 있음

from numpy.linalg import inv

betas_todo = inv(X.T @ X) @ X.T @ y

Question 2.

빨간색 자극이 녹색 자극보다 더 많은 활동을 일으킨다는 가설을 세우려고 한다. 이에 맞는 contrast vector를 만들어라.

cvec_red_green = np.array([0, 1, 1, -1, -1])

Question 3.

원 자극이 사각형 자극보다 더 많은 활동을 일으킨다는 가설을 세우려고 한다. 이에 맞는 contrast vector를 만들어라.

cvec_circle_square = np.array([0, -1, 1, -1, 1])

Question 4.

빨간색 자극이 녹색 자극보다 더 많은 활동을 일으키면서, 원 자극이 사각형 자극보다 더 많은 활동을 일으킨다는 가설을 세우려고 한다. 이에 맞는 contrast vector를 만들어라.

cvec_interaction = np.array([0, -1, 1, 1, -1])

Parametric designs

변수는 categorical뿐만 아니라 continuous, ordinal할 수 있음

이렇게 지속적으로 영향을 미치는 변수를 가지고 가설로 설정하는 것을 parametric design이라고 부름

parametric design은 두 가지 방식으로 복셀 반응에 영향을 미친다고 가정함

An “unmodulated” response
A parametric modulation of the response

위의 예시에서 가설은 보상 규모에 의한 반응을 보이는 복셀 영역이 있는지 찾는 것임

np.random.seed(42)
exp_time = 60 * 10  # i.e., 14 minutes in seconds
n_trials = 30
reward_onsets = np.arange(4, exp_time, exp_time / n_trials).astype(int) # size: 30

reward_magnitudes = np.random.randint(1, 11, size=reward_onsets.size) / 10

plt.figure(figsize=(15, 5))
plt.plot(reward_magnitudes, marker='o')

위의 가설에 따라 자극이 없을 때 0과, 자극이 있을 때 1을 가진 회귀식을 만들 것임

이 때, “unmodulated” response, parametric modulation of the response을 만족하는 조건을 가지는 두 개의 회귀식을 만듦

regressor of “unmodulated” response

hrf = glover_hrf(tr=1, oversampling=1)

x_unmod = np.zeros(exp_time)
x_unmod[reward_onsets] = 1
x_unmod = np.convolve(x_unmod, hrf)[:exp_time] # “unmodulated” response

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

regressor of parametric modulation of the response

보상의 크기 - 보상의 평균 크기를 활용해서 회귀식을 만듦
unmodulated regressor에서 modulated regressor만큼 decorrelate 하기 때문임

def parametric_modulation(reward_onsets, reward_magnitudes, exp_time, hrf, ans):
    reward_magnitudes_ms = reward_magnitudes - reward_magnitudes.mean() # parametric modulation of the response
    x_mod = np.zeros(exp_time)
    x_mod[reward_onsets] = reward_magnitudes_ms
    x_mod = np.convolve(x_mod, hrf)[:exp_time]

modulated regressor: 자극에 대한 일반적인 반응이 아니라 modulation of response에 대한 분산만 고려함

Design variance/efficiency

Optimizing effects

Psychological factors

“Are subjects really doing/thinking the way I want them to?”

“Am I really measuring what I think I’m measuring?”

Design factors

design variance

Optimizing design variance

design variance가 높으면 효율성이 낮고, design variance가 낮으면 효율성이 높음 (design variance는 각 변수의 variance와 covariance로 구성됨)

효율적인 설계는 다음과 같은 요인을 고려하면 됨

high variance (i.e. they should vary a lot relative to their mean)
low covariance (i.e. they should not correlate between each other a lot)

여기서 는 분산-공분산 행렬을, 는 contrast vector를 의미함

큰 효과(t)를 얻기 위해서는 세가지가 필요함

A large response/effect (β)
An efficient design (high variance, low covariance)
Low noise/unexplained variance (low SSE / DF)

def conceptual_tvalue_calculator(signal, noise, design_variance):
    """ Calculate the effect (t-value) from the signal, noise, and design variance components.
    
    Parameters
    ----------
    signal : int/float
    noise : int/float
    design_variance : int/float
    
    Returns
    -------
    effect : float
    efficiency : float
    """
    efficiency = 1 / design_variance
    effect = signal / (noise * design_variance)
    return effect, efficiency

SIGNAL = 0.5
NOISE = 2.8
DESIGN_VARIANCE = 0.02

design_var = c @ inv(X.T @ X) @ c.T
efficiency = 1.0 / c @ inv(X.T @ X) @ c.T

The effect of predictor variance on design variance/efficiency

변수의 값이 평균에서 더 많이 벗어날수록, 분산은 커짐

큰 분산을 가진 데이터는 큰 MSE를 가지며, 높은 t-value를 가짐

큰 분산을 가진 데이터는 작은 분산을 가진 데이터보다 design variance가 적다는 것을 의미함

The effect of predictor covariance on design variance/efficiency

Multicollinearity

변수들 간의 공분산을 multicollinearity라고 부름

변수들 간의 연관성이 높은 경우, 어떤 변수에 의한 변화인지 정확하게 측정할 수 없어서 efficiency가 낮음

(if correlation is high, then the ‘ambiguity, about the source of an effect is also high)

from scipy.stats import pearsonr

corr = pearsonr(X[:, 1], X[:, 2])[0]
eff = 1.0 / cvec.dot(np.linalg.inv(X.T.dot(X))).dot(cvec.T)

Evaluating multiple contrasts

하나의 기준이 아닌 여러 기준(A, B, A-B)에 대한 변화를 측정하고 싶을 때 사용함

아래의 예시는 contrast vector는 A, B, A-B에 대한 변화를 측정함

def overall_eff(X):
    cvecs = [np.array([0, 1, 0]), np.array([0, 0, 1]), np.array([0, 1, -1])]
    overall_eff_ans = len(cvecs) / np.sum([cvec @ inv(X.T @ X) @ cvec.T for cvec in cvecs])

How should we design our experiment to maximize efficiency?

Estimation vs. detection

1) [Majority research] You want to know whether different conditions activate voxel activity differently.

→ HRF-amplitude 에만 집중하면 됨 (detection)

2) [Small portion] You want investigate how different conditions influence voxel activity not only by investigating the “amplitude” parameter of the HRF, but also parameters relating to other properties of the shape of the HRF (like width, lag, strenght of undershoot, etc.).

→ HRF 모양이 나이 등에 따라 달라질 수 있음

→ HRF 함수에서 다양한 파라미터를 사용해야 함

Event-related vs. blocked designs

Event-related designs: 다른 자극이 무작위로 시행

Blocked designs: 같은 자극이 블록으로 그룹화되어 시행

np.random.seed(2)

N = 180
dg_hrf = glover_hrf(tr=1, oversampling=1)

# Blocked designs
blocked_pred1_onsets = list(range(10, 30)) + list(range(90, 110))
blocked_pred2_onsets = list(range(50, 70)) + list(range(130, 150))
N_stim = len(blocked_pred1_onsets)
blocked_pred1, blocked_pred2 = np.zeros(N), np.zeros(N)
blocked_pred1[blocked_pred1_onsets] = 1
blocked_pred2[blocked_pred2_onsets] = 1

icept = np.ones((N, 1))
X_blocked = np.hstack((
    icept,
    np.convolve(blocked_pred1, dg_hrf)[:N, np.newaxis],
    np.convolve(blocked_pred2, dg_hrf)[:N, np.newaxis]
))

# Event-related designs
er_stims = np.arange(N)
er_pred1_onsets = np.random.choice(er_stims, N_stim, replace=False)
er_stims_new = np.array([o for o in er_stims if o not in er_pred1_onsets])
er_pred2_onsets = np.random.choice(er_stims_new, N_stim, replace=False)
er_pred1, er_pred2 = np.zeros(N), np.zeros(N)
er_pred1[er_pred1_onsets] = 1
er_pred2[er_pred2_onsets] = 1

X_er = np.hstack((
    icept,
    np.convolve(er_pred1, dg_hrf)[:N, np.newaxis],
    np.convolve(er_pred2, dg_hrf)[:N, np.newaxis]
))

우리가 detection에 포커싱을 하고 싶다면, covariance가 낮고(correlation), variance가 높은(spread) blocked design을 사용해야 한다.

corr_blocked = pearsonr(X_blocked[:, 1], X_blocked[:, 2])
corr_er = pearsonr(X_er[:, 1], X_er[:, 2])
print(corr_blocked, corr_er)
>>> -0.238, -0.268

The “paradox” of efficiency

우리는 살아있는 인간의 fMRI 데이터를 scan하므로, 효율적인 Blocked designs에만 집중하는 것은 오류를 낳을 수 있음

fMRI 실험 자체의 특징(지겨움, 예측가능성, 습관, 불필요한 움직임 등)들에 의해서 오류에 빠질 수 있음

특히 불필요한 움직임으로 인해 LTI(Linear time-invariant)가 생겨 BOLD 신호를 방해함

심리적 요인 등의 가능성을 줄이기 위해, semi-random event-related designs을 활용할 필요가 있음

그렇다면, event-related designs을 어떻게 디자인해야 최적화할 수 있을지 고민해볼 필요 있음

Improving design efficiency for event-related designs using jittering

자극이 없는 동안의 짧은 기간을 ISI(inter-stimulus interval)이라고 부름

아래 예시는 8초의 ISI(자극의 간격)를 가짐

def simulate_data_fixed_ISI(N=420):

    dg_hrf = glover_hrf(tr=1, oversampling=1)
    
    # Create indices in regularly spaced intervals (9 seconds, i.e. 1 sec stim + 8 ISI)
    stim_onsets = np.arange(10, N - 15, 9)
    stimcodes = np.repeat([1, 2], stim_onsets.size / 2)  # create codes for two conditions
    np.random.shuffle(stimcodes)  # random shuffle
    stim = np.zeros((N, 1))

    c = np.array([[0, 1, 0], [0, 0, 1]])

    # Fill stim array with codes at onsets
    for i, stim_onset in enumerate(stim_onsets):
        stim[stim_onset] = 1 if stimcodes[i] == 1 else 2
    
    stims_A = (stim == 1).astype(int)
    stims_B = (stim == 2).astype(int)

    reg_A = np.convolve(stims_A.squeeze(), dg_hrf)[:N]
    reg_B = np.convolve(stims_B.squeeze(), dg_hrf)[:N]
    X = np.hstack((np.ones((reg_B.size, 1)), reg_A[:, np.newaxis], reg_B[:, np.newaxis]))
    dvars = [(c[i, :].dot(np.linalg.inv(X.T.dot(X))).dot(c[i, :].T))
             for i in range(c.shape[0])]
    eff = c.shape[0] / np.sum(dvars)
    return X, eff

X, eff = simulate_data_fixed_ISI()

위와 달리 ISI가 시도 때마다 다양해질 수 있음 → 이를 ‘jittering’이라고 부름

jittering은 covariance를 줄이고 variance를 높일 수 있음

아래 예시는 2-6초의 (가변) ISI를 가짐

def simulate_data_jittered_ISI(N=420):
    
    dg_hrf = glover_hrf(tr=1, oversampling=1)
    
    stim_onsets = np.arange(10, N - 15, 9)
    stimcodes = np.repeat([1, 2], stim_onsets.size / 2)
    np.random.shuffle(stimcodes)

    # Here, we pick some *deviations* from the standard ISI (i.e., 8),
    # so possible ISIs are (8 - 2, 8 - 1, 8 - 0, 8 + 1, 8 + 2)
    ISIs = np.repeat([-2, -1, 0, 1, 2], repeats=11)
    np.random.shuffle(ISIs)
    
    stim = np.zeros((N, 1))
    c = np.array([[0, 1, 0], [0, 0, 1]])

    for i, stim_onset in enumerate(stim_onsets):
        # We subtract the stim-onset with -2, -1, 0, 1, or 2 (from ISIs)
        # to simulate jittering
        stim[stim_onset - ISIs[i]] = 1 if stimcodes[i] == 1 else 2
        
    stims_A = (stim == 1).astype(int)
    stims_B = (stim == 2).astype(int)
    reg_A = np.convolve(stims_A.squeeze(), dg_hrf)[:N]
    reg_B = np.convolve(stims_B.squeeze(), dg_hrf)[:N]
    X = np.hstack((np.ones((reg_B.size, 1)), reg_A[:, np.newaxis], reg_B[:, np.newaxis]))

    # Loop over the two contrasts
    dvars = [(c[i, :].dot(np.linalg.inv(X.T.dot(X))).dot(c[i, :].T))
             for i in range(c.shape[0])]
    eff = c.shape[0] / np.sum(dvars)
    return X, eff