{ "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" }, "name": "" }, "name": "fixed_and_random_FX_exercises", "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "GLM, fixed versus random effects in fMRI" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Plan" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Remind what is the issue with random and fixed effects\n", "2. Construct some data set from one bold series\n", " 1. Extract a time series \n", " 2. Add some signal (3 conditions)\n", "3. Solve the GLM for one subject\n", "4. Solve the GLM for 7 subjects\n", " 1. Fixed\n", " 2. Random\n", "5. Testing. What is the distribution of our contrast of beta values ?" ] }, { "cell_type": "code", "collapsed": true, "input": [ "# - \n", "from __future__ import print_function, division\n", "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# - set gray colormap and nearest neighbor interpolation by default\n", "plt.rcParams['image.cmap'] = 'gray'\n", "plt.rcParams['image.interpolation'] = 'nearest'\n", "\n", "# - import numpy.linalg with a shorter name\n", "import numpy.linalg as npl" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": true, "input": [ "# - \n", "np.set_printoptions(precision=3)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "1. The issue with random and fixed effect" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "See http://www.fil.ion.ucl.ac.uk/spm/doc/books/hbf2/pdfs/Ch12.pdf" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "2. Construct some data from a BOLD time series" ] }, { "cell_type": "code", "collapsed": true, "input": [ "# - nibabel\n", "import nibabel as nib" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": true, "input": [ "# - change this to the path of this image on your system\n", "my_path_to_file = \"/home/jb/code/pna2015/pna2015/day10/ds114_sub009_t2r1.nii\"" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# load this image, get the 32,32,15 voxel time series without the 5 first points.\n", "# plot the time course" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - normalize the time series to have mean 100. call result 'bold'\n", "\n", "bold = np.asarray(img.get_data()[x,y,z,to_rm:])\n", "bold_std = bold.std()\n", "bold_mean = bold.mean()\n", "print(bold_std, bold_mean)\n" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - Normalized to have mean 100\n", "bold = bold * 100. / bold_mean\n", "bold_std = bold.std()\n", "bold_mean = bold.mean()\n", "sig_mag = 2*bold_std\n", "print(bold_std, bold_mean, sig_mag)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - divide the total duration time in 12, into 12 intervals of \"cond_duration\" TR\n", "# we will fill the time with 3 conditions, (eg A,B,C)\n", "# A, B will be our experimental conditions, C is for control condition (baseline)\n", "\n", "cond_duration = bold.shape[0] // 12 # integer division\n", "\n", "# make a signal for two conditions: A, B. Each block of A or B lasts for cond_duration, \n", "# and these blocks occur regularely, eg if the time is divided in 12 intervals,\n", "# 'A' blocks will occur at intervals 0, 3, 6, 9, 'B' at intervals 1, 4, 7, 10, and the \n", "# remaining intervals will be baseline \n", "\n", "ones = np.ones((cond_duration,1))\n", "a = np.kron(np.eye(3), ones) \n", "\n", "print(a.shape)\n", "signal_A = np.kron(np.ones(4), a[:,0]).flatten()\n", "signal_B = np.kron(np.ones(4), a[:,1]).flatten()\n", "baseline = np.ones(shape=(bold.shape[0]))\n", "plt.plot(signal_A)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - construct the design matrix with condition A, B and baseline using np.vstack\n", "X_ = np.vstack((signal_A, signal_B, baseline)).T\n", "plt.imshow(X_,aspect='auto',interpolation='nearest')\n", "print(X_.shape)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - add some signal with a signal to noise ratio of 2, \n", "# sig_mag = 1.0*bold_std\n", "\n", "sig_mag = 1.0\n", "c = [sig_mag*bold_std, 0, 0]\n", "Y = X_.dot(c) + bold\n", "time = np.arange(X_.shape[0])\n", "plt.plot(time, Y, time, X_.dot(c)+100.)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# Add a trend to X_, call this new matrix X" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Now, recall the GLM and how to test a contrast " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Looking again at our time series, we should add a trend..." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# - estimate the parameters of the model X_ and X using the OLS (see below)\n", "\n", "pinvX_ = npl.inv(X_.T.dot(X_)).dot(X_.T) # this is to match the formula we have seen, \n", "pinvX = npl.inv(X.T.dot(X)).dot(X.T) # in practice do not use this, use npl.pinv(X2)\n", "\n", "beta_ = pinvX_.dot(Y)\n", "beta = pinvX.dot(Y)\n", "print(beta_, beta)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We note our estimates 'beta' $\\hat{\\beta}_1$, $\\hat{\\beta}_2$, $\\hat{\\beta}_3$, etc" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Now, we need to _test_ if our signal of interest $\\hat{\\beta}_1$ is significant" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will need 2 things: signal, and standard error of the signal." ] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "First, what is our observed signal of interest ?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# express our signal of interest as c.dot(beta) - print the signal magnitude. " ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Second, what is our observed standard error ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall the linear model for one observation i :\n", "\n", "$$ Y_i = X \\beta + \\epsilon_i $$\n", "\n", "and in matrix form:\n", "\n", "$$ Y = X \\beta + \\epsilon $$\n", "\n", "with $ \\epsilon_i $ independant gaussian variable with variance $\\sigma^2$, such that the vector $\\epsilon$ has variance $\\sigma^2 I$, $I$ the identity matrix.\n", "\n", "Our OLS estimates are:\n", "\n", "$$ \\hat\\beta = (X^tX)^{-1}X^t Y $$\n", "\n", "Our fitted data : \n", "\n", "$$ \\hat{Y} = X\\hat\\beta = X(X^tX)^{-1}X^t Y = M_X Y$$\n", "\n", "And our residuals:\n", "\n", "$$ r = \\hat{\\epsilon} = Y - \\hat{Y} = Y - M_X Y = (I - M_X) Y = R_X Y $$\n", "\n", "$R_X$ or $R$ when there is no ambiguity, is the __residual forming__ matrix. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Let's now compute $ \\mathrm{var}(c^T \\hat\\beta) $. We have\n", "\n", "$$ \\mathrm{var}(c^T \\hat\\beta) = \\mathrm{var}(c^T (X^tX)^{-1}X^t Y) \\qquad (A) $$ \n", "\n", "Let's call $ X^+ = (X^tX)^{-1}X^t $.\n", "\n", "We have assumed that $\\mathrm{var}(Y) = \\sigma^2 I$\n", "with $I$ the identity matrix.\n", "\n", "$$ (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 $$ \n", "\n" ] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "How do we know $\\sigma$ ? We don't, we have to estimate it. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we will need the residual sum of square:\n", "\n", "Our residuals, the estimated noise, are $ r = \\hat{\\epsilon} = Y - \\hat{Y} $, with $ \\hat{Y} = X \\hat{\\beta}$\n", "\n", "The residual sum of squre is :\n", "\n", "$$ \\textrm{RSS} = \\sum_{i=1}^{n} (Y_i - \\widehat{Y_i})^2 $$\n", "\n", "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 \n", "\n", "$$\\hat\\sigma^2 = \\textrm{RSS} / (n - p)$$\n", "\n", "Where p here is 4, and n = 196.\n", "\n", "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.\n", "\n", "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} $\n", "\n", "The t test is then : $t = {c^T \\hat\\beta}/{\\textrm{SE}} $\n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# - Compute the estimate of sigma: compute the residuals, take the sum of square, and \n", "# divide by n-p, \n", "# n is the number of observations, p is the number of (free) parameters that we estimate\n", "\n", "n,p = X.shape # ok because X has no linear dependencies for its columns !\n", "r = Y - X.dot(beta)\n", "RSS = (r**2).sum()\n", "sigma2_hat = RSS/(n - p)\n", "\n", "print(RSS, sigma2_hat)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - what can you say of this compared to bold_std ? \n", "# Compute the standard error of the signal\n", "\n", "V_1 = npl.inv(X.T.dot(X))\n", "\n", "SE = np.sqrt( sigma2_hat * c.T.dot(V_1).dot(c))[0,0]\n", "print(V_1, \"\\n\", SE)\n" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - compute the t statistics, that follows a student distribution with n-p df\n", "\n", "t_observed = signal/SE\n", "print(t_observed)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - \n", "signal, sigma2_hat, SE" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "How far from the t distribution centered on 0 ?" ] }, { "cell_type": "code", "collapsed": true, "input": [ "# - \n", "import scipy.stats as sst" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# compute the p of the t statistics, use sst.t(df) to define a t_variable\n", "# then find how to compute the p-value (see cdf or sf functions of t)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "F statistics and model selection" ] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Review the F test" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The F-test is best understood in terms of model selection.\n", "The question we are asking is : should all the regressor be in the model ? \n", "\n", "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). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{eqnarray} \n", "F_{\\nu_1, \\nu_2} & = & \\frac{(\\hat\\epsilon_0^t \\hat\\epsilon_0 - \\hat\\epsilon^T\\hat\\epsilon)/ \\nu_{1} }{\\hat\\epsilon^T\\hat\\epsilon/\\nu_{2}} \\\\ \n", " & = & \\frac{(Y^TR_{X_0}Y - Y^TR_{X}Y)/ \\nu_{1} }{Y^TR_{X}Y/\\nu_{2}} \\\\ \n", " & = & \\frac{(Y^T(I-P_{X_0})Y - Y^T(I-P_{X})Y)/ \\nu_{1} }{Y^T(I- P_{X})Y/\\nu_{2}} \\\\ \n", "& = & \\frac{(\\textrm{SSR}(X_0) - \\textrm{SSR}(X))/\\nu_1}{\\textrm{SSR}(X)/\\nu_2}\n", "\\end{eqnarray}\n", "$$\n", "\n", "What are $\\nu_{1}$ and $\\nu_{2}$ ?\n", "\n", "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. " ] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Let's apply it in our simple case :" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# compute the F test for regressor 0 in our previous model\n", "# what is the difference with the t-test" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - define a glm function that returns beta, Yfitted, residuals\n", "\n", "def glm(X, Y):\n", " \"\"\" a simple GLM function returning the estimated parameters and residuals \n", " \n", " parameters:\n", " -----------\n", " X: np array\n", " The desing matrix\n", " Y: np array\n", " The data\n", " \n", " returns:\n", " --------\n", " betah : np array\n", " estimated beta\n", " Yfitted: np array\n", " fitted data \n", " resid: np array\n", " residuals\n", " \"\"\"\n", " \n", " betah = npl.pinv(X).dot(Y)\n", " Yfitted = X.dot(betah)\n", " resid = Y - Yfitted\n", " return betah, Yfitted, resid\n", "\n", "\n", "# define an F test function:\n", "\n", "def F_test(X, X0, Y):\n", " \"\"\"\n", " parameters:\n", " -----------\n", " X: np array\n", " The design matrix\n", " X0: np array\n", " The reduced design matrix\n", " Y: np array\n", " The data\n", " \n", " returns\n", " --------\n", " F: float\n", " the F statistics\n", " nu1, nu2: floats\n", " the degrees of freedom\n", " \"\"\"\n", " betah, Yfitted, resid = glm(X,Y)\n", " betah0, Yfitted0, resid0 = glm(X0,Y)\n", " \n", " nu1 = X.shape[1] - X0.shape[1] # assume X full rank (no linear dependencies)\n", " nu2 = X.shape[0] - X.shape[1]\n", "\n", " numerator = ((resid0**2).sum() - (resid**2).sum())/nu1 \n", " denominator = (resid**2).sum()/nu2\n", " F = numerator / denominator\n", " return F, nu1, nu2\n", "\n", "# define a t test function:\n", "\n", "def t_test(X, c, Y):\n", " \"\"\"\n", " \"\"\"\n", " assert c.ndim == 2\n", " assert c.shape[1] == 1\n", " betah = npl.pinv(X).dot(Y)\n", " signal = c.T.dot(betah)[0,0]\n", " RSS = ((Y - X.dot(betah))**2).sum()\n", " p = npl.matrix_rank(X)\n", " n = X.shape[0]\n", " df = n-p\n", " V_1 = npl.pinv(X.T.dot(X))\n", " SE = np.sqrt( (RSS/df) * c.T.dot(V_1).dot(c))[0,0]\n", " return signal/SE, df \n", " " ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# compute the F test for regressor 0 in our previous model \n", "# use the F_test function\n", "# use sst.f to find the p value, compare it to the p_ttest\n", "# why is p_ttest half of the p_Ftest ?" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": true, "input": [ "# - F test can be computed by stacking several one dimensional contrasts \n", "# into a \"contrasts matrix\" - no demonstration provided here but see reference\n", "# Poline and Brett, Neuroimage 2012." ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - now, test A and B :\n", "X0 = X[:,2:]\n", "f_obs, nu1, nu2 = F_test(X, X0, Y)\n", "\n", "F_variable = sst.f(nu1, nu2)\n", "p_FtestAB = F_variable.sf(f_obs)\n", "\n", "print(p_FtestAB, p_Ftest)\n", "# is this expected ?" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Let's create 5 subjects with the same design" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# - \n", "N = 5\n", "XN = np.kron(np.eye(N), X)\n", "plt.imshow(XN, aspect='auto',interpolation='nearest')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - create N subject data, add noise to the previous data\n", "YN = np.kron(np.ones(N), Y)[:,np.newaxis]\n", "YN += np.random.normal(0, 0.7, size=YN.shape)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - create contrast\n", "c = np.asarray([1,0,0,0]).reshape(4,1)\n", "cN = np.kron(np.ones(shape=(1,N)), c.T).T\n", "cN.shape" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - compute t stat for fixed effect analysis\n", "t_Nfix, df = t_test(XN, cN, YN)\n", "print(t_Nfix, df, sst.t.sf(t_Nfix, df)*2)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Create a random effect model:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# compute signal (c.T.dot(beta)) for each subject" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# What is the random effect t-test ? Create a design matrix and a contrast, use t_test()" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - check that we get same between t_test and scipy stat ttest_1samp\n", "sst.ttest_1samp(signals, 0)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 4, "metadata": {}, "source": [ "Compare the two statistics: which one is \"better\" ?" ] }, { "cell_type": "code", "collapsed": true, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }