2. Using the GLM to fMRI: inference

2.3 How to compute t-values and p-values

2.3.1 t-values

• t-values: 베타 변수의 표준 오류로 베타 변수를 나눈 값(uncertainty-normalized)

• 표준 오류는 모델의 설명할 수 없는 분산(noise)와 설계된 분산(design)으로 나뉨

• noise는 제곱 잔차(SSE)를 구하고 모델의 자유도(DF)로 나눔

• design variance는 아래의 코드를 참고하길 바람

def design_variance(X, which_predictor=1):
    ''' Returns the design variance of a predictor (or contrast) in X.
    
    Parameters
    ----------
    X : numpy array
        Array of shape (N, P)
    which_predictor : int or list/array
        The index of the predictor you want the design var from.
        Note that 0 refers to the intercept!
        Alternatively, "which_predictor" can be a contrast-vector
        (which will be discussed later this lab).
        
    Returns
    -------
    des_var : float
        Design variance of the specified predictor/contrast from X.
    '''
    
    is_single = isinstance(which_predictor, int) # isinstance(): 객체가 특정 클래스의 인스턴스인지 확인
    if is_single:
        idx = which_predictor
    else:
        idx = np.array(which_predictor) != 0
    
    c = np.zeros(X.shape[1])
    c[idx] = 1 if is_single == 1 else which_predictor[idx]
    des_var = c.dot(np.linalg.inv(X.T.dot(X))).dot(c.T)
    return des_var

# 선형 회귀 식
beta_meter = inv(Xn.T @ Xn) @ Xn.T @ y
y_hat_meter = Xn @ beta_meter

N = y.size
P = Xn.shape[1]
df = (N - P)
print(df)
>>> 98

sigma_hat = np.sum((y - y_hat_meter) ** 2) / df # noise
print(sigma_hat)
>>> 0.005

design_variance_weight = design_variance(Xn, 1) # design variance
print(design_variance_weight)
>>> 0.000334

t_meter = beta_meter[1] / np.sqrt(sigma_hat * design_variance_weight) # t-test of beta
print(beta_meter, t_meter)
>>> (0.013, 10.431)

2.3.2 p-values

• scipy 패키지 안의 stats.t.sf(t_value, df) 함수를 통해서 p-value를 쉽게 구할 수 있음

• (주의) 이 함수는 항상 right-tailed p-value를 구함, t-value가 음수라면 left-tailed p-value를 구해야함(np.abs(t_value) 를 통해 구할 수 있음)

from scipy import stats

p_value = stats.t.sf(np.abs(t_meter), df) * 2 # two-tailed p-value는 2를 곱하면 구할 수 있음
print(df, t_meter, p_value)
>>> 98, 10.431, 0.00000000

2.4 Contrasts

2.4.1 t-tests

• 가설 세우기
1. Using a contrast of a parameters “against baseline”
• (null-hypothesis, i.e. there is no activation compared to baseline)
• (alternative hypothesis, i.e. blue activates relative to baseline)
1. Using a contrast between parameters
• (null-hypothesis, i.e. there is no difference)
• (alternative hypothesis, i.e. there is some difference)

• Simple example of hypothesis in this section
• condition 1: images of male faces with a happy expression
• condition 2: images of male faces with a sad expression
• condition 3: images of male faces with a neutral expression
• condition 4: images of female faces with a happy expression
• condition 5: images of female faces with a sad expression
• condition 6: images of female faces with a neutral expression

data = np.load('data_contrast_example.npz')
X, y = data['X'], data['y']

betas = inv(X.T @ X) @ X.T @ y
betas = betas.squeeze()
print(betas)
>>> array([ 0.08208567, -0.21982422, -0.16284892,  0.53208935,  0.26214462,
        0.38945094,  0.21565532])

• using contrast-vectors: 복잡한 가설 검증을 할 때 유용하게 활용 가능

# 가설 3을 검정하고 싶을 때 
contrast_vector = np.array([0, 0, 1, 0, 0, 0, 0])
contrast = (betas * contrast_vector).sum()
print(contrast)
>>> -0.16284892

# (가설 3 - 가설 2)을 검정하고 싶을 때
contrast_vector = np.array([0, -1, 1, 0, 0, 0, 0])
contrast = (betas * contrast_vector).sum()
print(contrast)
>>> 0.0569753

# (감정에 관계 없이) 남성 얼굴과 여성 얼굴과의 차이를 알고 싶을 때
contrast_vector = [0, 1, 1, 1, -1, -1, -1]
male_female_contrast = (contrast_vector * betas).sum()
print(male_female_contrast)
>>> -0.71783467

• “generalized” version of the t-test of contrast/hypothesis

# (가설 3 - 가설 2)을 검정하고 싶을 때의 t-score
c_vec = np.array([0, 1, -1, 0, 0, 0, 0])  # contrast vector
dvar = design_variance(X, which_predictor=c_vec # design variance
sigma_hat = np.sum((y - y_hat_meter) ** 2) / df # noise

def test_contrast(X_new, y_new, betas, tval_ans, pval_ans):
    cvec = np.array([0, -1, 1, 0, -1, 1, 0]) # contrast vector
    this_dvar = design_variance(X_new, cvec) # design variance
    y_hat = X_new.dot(betas) # y_hat
    this_sse = ((y_new - y_hat) ** 2).sum() / (X_new.shape[0] - X_new.shape[1]) # sigma_hat
    tval_todo = cvec.dot(betas) / np.sqrt(this_sse * this_dvar) # t-value
    pval_todo = stats.t.sf(np.abs(tval_todo), X_new.shape[0] - X_new.shape[1]) * 2
    np.testing.assert_almost_equal(tval_todo, tval_ans) # score of t-value
    np.testing.assert_almost_equal(pval_todo, pval_ans) # score of p-value

2.4.2 F-tests