# -
from __future__ import print_function, division
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
# - import numpy.linalg with a shorter name
import numpy.linalg as npl
# -
np.set_printoptions(precision=3)
# - nibabel
import nibabel as nib
# - change this to the path of this image on your system
my_path_to_file = "/home/jb/code/pna2015/pna2015/day10/ds114_sub009_t2r1.nii"
# load this image, get the 32,32,15 voxel time series without the 5 first points.
# plot the time course
# - normalize the time series to have mean 100. call result 'bold'
bold = np.asarray(img.get_data()[x,y,z,to_rm:])
bold_std = bold.std()
bold_mean = bold.mean()
print(bold_std, bold_mean)
# - Normalized to have mean 100
bold = bold * 100. / bold_mean
bold_std = bold.std()
bold_mean = bold.mean()
sig_mag = 2*bold_std
print(bold_std, bold_mean, sig_mag)
# - divide the total duration time in 12, into 12 intervals of "cond_duration" TR
# we will fill the time with 3 conditions, (eg A,B,C)
# A, B will be our experimental conditions, C is for control condition (baseline)
cond_duration = bold.shape[0] // 12 # integer division
# make a signal for two conditions: A, B. Each block of A or B lasts for cond_duration,
# and these blocks occur regularely, eg if the time is divided in 12 intervals,
# 'A' blocks will occur at intervals 0, 3, 6, 9, 'B' at intervals 1, 4, 7, 10, and the
# remaining intervals will be baseline
ones = np.ones((cond_duration,1))
a = np.kron(np.eye(3), ones)
print(a.shape)
signal_A = np.kron(np.ones(4), a[:,0]).flatten()
signal_B = np.kron(np.ones(4), a[:,1]).flatten()
baseline = np.ones(shape=(bold.shape[0]))
plt.plot(signal_A)
# - construct the design matrix with condition A, B and baseline using np.vstack
X_ = np.vstack((signal_A, signal_B, baseline)).T
plt.imshow(X_,aspect='auto',interpolation='nearest')
print(X_.shape)
# - add some signal with a signal to noise ratio of 2,
# sig_mag = 1.0*bold_std
sig_mag = 1.0
c = [sig_mag*bold_std, 0, 0]
Y = X_.dot(c) + bold
time = np.arange(X_.shape[0])
plt.plot(time, Y, time, X_.dot(c)+100.)
# Add a trend to X_, call this new matrix X
Looking again at our time series, we should add a trend...
# - estimate the parameters of the model X_ and X using the OLS (see below)
pinvX_ = npl.inv(X_.T.dot(X_)).dot(X_.T) # this is to match the formula we have seen,
pinvX = npl.inv(X.T.dot(X)).dot(X.T) # in practice do not use this, use npl.pinv(X2)
beta_ = pinvX_.dot(Y)
beta = pinvX.dot(Y)
print(beta_, beta)
We note our estimates 'beta' $\hat{\beta}_1$, $\hat{\beta}_2$, $\hat{\beta}_3$, etc
We will need 2 things: signal, and standard error of the signal.
# express our signal of interest as c.dot(beta) - print the signal magnitude.
Recall the linear model for one observation i :
$$ Y_i = X \beta + \epsilon_i $$and in matrix form:
$$ Y = X \beta + \epsilon $$with $ \epsilon_i $ independant gaussian variable with variance $\sigma^2$, such that the vector $\epsilon$ has variance $\sigma^2 I$, $I$ the identity matrix.
Our OLS estimates are:
$$ \hat\beta = (X^tX)^{-1}X^t Y $$Our fitted data :
$$ \hat{Y} = X\hat\beta = X(X^tX)^{-1}X^t Y = M_X Y$$And our residuals:
$$ r = \hat{\epsilon} = Y - \hat{Y} = Y - M_X Y = (I - M_X) Y = R_X Y $$$R_X$ or $R$ when there is no ambiguity, is the residual forming matrix.
Let's now compute $ \mathrm{var}(c^T \hat\beta) $. We have
$$ \mathrm{var}(c^T \hat\beta) = \mathrm{var}(c^T (X^tX)^{-1}X^t Y) \qquad (A) $$Let's call $ X^+ = (X^tX)^{-1}X^t $.
We have assumed that $\mathrm{var}(Y) = \sigma^2 I$ with $I$ the identity matrix.
$$ (A) \;\; = c^T X^+ \mathrm{var}(Y) X^{+T} c = \sigma^2 c^T X^+ X^{+T} c = \sigma^2 c^T (X^T X)^- c $$First, we will need the residual sum of square:
Our residuals, the estimated noise, are $ r = \hat{\epsilon} = Y - \hat{Y} $, with $ \hat{Y} = X \hat{\beta}$
The residual sum of squre is :
$$ \textrm{RSS} = \sum_{i=1}^{n} (Y_i - \widehat{Y_i})^2 $$But the mean sum of square is not the simple division by the number of time points, we have to correct for the number of estimation we have done when computing $\widehat{Y}$, so that will be
$$\hat\sigma^2 = \textrm{RSS} / (n - p)$$Where p here is 4, and n = 196.
This estimation is also called the Mean Residual Sum of Square (MRSS), and $\hat\sigma = \textrm{MRSS} = \textrm{RSS} / \it{df} $ where $\it{df}$ is the degrees of freedom $(n-p) = 196 - 4$, the number of observations (here, 196) minus the rank of the design matrix (ie the number of independant columns, here 4). The MRSS is the estimation of our noise level.
Coming back to our t-statistics, we can now calculate the Standard Error (SE) of $c^T \hat\beta$ with $\textrm{SE} = \sqrt{\hat{\sigma}^2 c^T (X^t X)^+ c} $
The t test is then : $t = {c^T \hat\beta}/{\textrm{SE}} $
# - Compute the estimate of sigma: compute the residuals, take the sum of square, and
# divide by n-p,
# n is the number of observations, p is the number of (free) parameters that we estimate
n,p = X.shape # ok because X has no linear dependencies for its columns !
r = Y - X.dot(beta)
RSS = (r**2).sum()
sigma2_hat = RSS/(n - p)
print(RSS, sigma2_hat)
# - what can you say of this compared to bold_std ?
# Compute the standard error of the signal
V_1 = npl.inv(X.T.dot(X))
SE = np.sqrt( sigma2_hat * c.T.dot(V_1).dot(c))[0,0]
print(V_1, "\n", SE)
# - compute the t statistics, that follows a student distribution with n-p df
t_observed = signal/SE
print(t_observed)
# -
signal, sigma2_hat, SE
# -
import scipy.stats as sst
# compute the p of the t statistics, use sst.t(df) to define a t_variable
# then find how to compute the p-value (see cdf or sf functions of t)
The F-test is best understood in terms of model selection. The question we are asking is : should all the regressor be in the model ?
To answer this, we construct a reduced model X0 from X. From this model we look at the RSS of model X (full model) and model X0 (reduced). The residuals are always going to be larger for the reduced model (how do we prove this ? Let's see in class).
What are $\nu_{1}$ and $\nu_{2}$ ?
They are the same things that we had for the t-test. The difference $ \hat\epsilon_0^t \hat\epsilon_0 - \hat\epsilon^T\hat\epsilon $ will have the difference of each term df.
# compute the F test for regressor 0 in our previous model
# what is the difference with the t-test
# - define a glm function that returns beta, Yfitted, residuals
def glm(X, Y):
""" a simple GLM function returning the estimated parameters and residuals
parameters:
-----------
X: np array
The desing matrix
Y: np array
The data
returns:
--------
betah : np array
estimated beta
Yfitted: np array
fitted data
resid: np array
residuals
"""
betah = npl.pinv(X).dot(Y)
Yfitted = X.dot(betah)
resid = Y - Yfitted
return betah, Yfitted, resid
# define an F test function:
def F_test(X, X0, Y):
"""
parameters:
-----------
X: np array
The design matrix
X0: np array
The reduced design matrix
Y: np array
The data
returns
--------
F: float
the F statistics
nu1, nu2: floats
the degrees of freedom
"""
betah, Yfitted, resid = glm(X,Y)
betah0, Yfitted0, resid0 = glm(X0,Y)
nu1 = X.shape[1] - X0.shape[1] # assume X full rank (no linear dependencies)
nu2 = X.shape[0] - X.shape[1]
numerator = ((resid0**2).sum() - (resid**2).sum())/nu1
denominator = (resid**2).sum()/nu2
F = numerator / denominator
return F, nu1, nu2
# define a t test function:
def t_test(X, c, Y):
"""
"""
assert c.ndim == 2
assert c.shape[1] == 1
betah = npl.pinv(X).dot(Y)
signal = c.T.dot(betah)[0,0]
RSS = ((Y - X.dot(betah))**2).sum()
p = npl.matrix_rank(X)
n = X.shape[0]
df = n-p
V_1 = npl.pinv(X.T.dot(X))
SE = np.sqrt( (RSS/df) * c.T.dot(V_1).dot(c))[0,0]
return signal/SE, df
# compute the F test for regressor 0 in our previous model
# use the F_test function
# use sst.f to find the p value, compare it to the p_ttest
# why is p_ttest half of the p_Ftest ?
# - F test can be computed by stacking several one dimensional contrasts
# into a "contrasts matrix" - no demonstration provided here but see reference
# Poline and Brett, Neuroimage 2012.
# - now, test A and B :
X0 = X[:,2:]
f_obs, nu1, nu2 = F_test(X, X0, Y)
F_variable = sst.f(nu1, nu2)
p_FtestAB = F_variable.sf(f_obs)
print(p_FtestAB, p_Ftest)
# is this expected ?
# -
N = 5
XN = np.kron(np.eye(N), X)
plt.imshow(XN, aspect='auto',interpolation='nearest')
# - create N subject data, add noise to the previous data
YN = np.kron(np.ones(N), Y)[:,np.newaxis]
YN += np.random.normal(0, 0.7, size=YN.shape)
# - create contrast
c = np.asarray([1,0,0,0]).reshape(4,1)
cN = np.kron(np.ones(shape=(1,N)), c.T).T
cN.shape
# - compute t stat for fixed effect analysis
t_Nfix, df = t_test(XN, cN, YN)
print(t_Nfix, df, sst.t.sf(t_Nfix, df)*2)
# compute signal (c.T.dot(beta)) for each subject
# What is the random effect t-test ? Create a design matrix and a contrast, use t_test()
# - check that we get same between t_test and scipy stat ttest_1samp
sst.ttest_1samp(signals, 0)