Looking at PCA in Action

There are different approaches to perform a PCA operation. Software tools (such as Matlab and R) offer two common approaches that we are going to discuss here.

1. Using Eigen-decomposition

The first approach is based on what was written in the previous blog post on PCA. The step-by-step example of PCA below is in Matlab as it is easier to work with matrices. The first step is to mean-center the original data. Suppose you have 2D matrix called data, with [x₁, x₂, ..., x₁₀₀] and [y₁, y₂, ..., y₁₀₀] with 100 observations.

% Step-1: Mean-center the original data.
data2(:,1) = data(:,1) - mean(data(:,1));
data2(:,2) = data(:,2) - mean(data(:,2));

Next, we find the covariance matrix to capture all possible relationships of the variables in the dataset.

% Step-2: Calculate covariance matrix of the mean-centered data.
C = cov(data2);

Once the covariance matrix C is computed we then perform eigen-decomposition on it to find the eigenvectors and eigenvalues according to CV = VD, leading to C = VDV^T. (I think we have to flip the matrix V left and right).

% Step-3: Produce eigenvalues and eigenvectors: a diagonal matrix D of eigenvalues and matrix V
% whose columns are the corresponding eigenvectors so that C*V = V*D.
[V,D] = eig(C);

Remember now the eigenvalues measure the "amount of variation" retained by each principal component. So, we rank the eigenvalues according to the magnitude or "importance". For example, suppose the columns of matrix V and the diagonal elements of D contain the following (see right):

The diagonal elements contain two eigenvalues λ₁ and λ₂. After ranking these two, it is seen that  λ₁ > λ₂, which means that the eigenvector v₁ corresponds to the first component, and the one that corresponds to the second is v₂.
Finally, the last step is to rotate or re-orient the data by [new data] = [original data].[eigenvector]^T. The new dataset is called the component scores.

% Step-4: The new axes is formed using V, to rotate the mean-centered data.
newdata = data2 * V';
newdata = fliplr(newdata) 

What is the percentage of information explained or accounted for by each component? It turns out that this can be found by dividing the corresponding eigenvalue by the total sum of all eigenvalues.

% Step-5: Get the variance accounted for, use matrix D containing the eigenvalues.
explained = D/sum(D(:))

Look at Fig-1. The original dataset in the left panel shows correlated variables. The middle panel shows how the dataset can be represented in terms of two principal components. PC1 axis is the first principal direction along which the data show the biggest variation. PC2 axis is the second most important direction and it is orthogonal to the PC1 axis. The new dataset in the right panel is a "projection" onto the principal components, where projection onto PC1 (X-axis) is uncorrelated with projection onto PC2 (Y-axis). 

Fig-1: The process of rotating the original data into a new dimension of principal components.

Where is the dimensionality reduction here? Take a look at the right panel above. The transformed data have just as many dimensions as the original data. By visual inspection, most variability is along the PC1 (X-axis) so you can just throw away the data plotted on the Y-axis. In other words, our 2D data can be reduced to 1D by projecting each sample onto the first principal component. Do we lose anything? Well, a little bit, depending on the variation that the PC1 and PC2 explain. In data science, data variation means information. Based on the above example, we find that PC1 and PC2 carry 1.284/(0.049+1.284) = 96%; and 0.049/(0.049+1.284) = 4% of the variance of the data respectively.

2. Using Singular Value Decomposition
PCA is intimately related to another statistical analysis called the Singular Value Decomposition (SVD). Although notations may differ in different textbooks, SVD essentially involves the factorization of matrix X = UDV^T, where U is an n × n matrix, with columns as orthogonal unit vectors of length n (left singular vector of X); D is an n × p rectangular diagonal matrix, the singular values of X; V is a p × p, with columns as orthogonal unit vectors of length p (right singular vectors of X).

The new data can be obtained through Y = XV = UD.  The transpose of V is sometimes called the whitening and can be used as a preparation in ICA. Columns of V multiplied by the square root of corresponding eigenvalues, i.e. eigenvectors scaled up by the variances, are the coefficients or loadings.

Now PCA can also be done in Matlab directly using a function:
[coeff,newdata,latent,tsquared,explained,mu] =  pca(data)

This Matlab function automatically mean-centers the data. The coefficients are also known as loadings. The loadings can be understood as the weights for each original variable when calculating the new scores. Each column of coeff contains coefficients for each component, and the columns are in descending order of the variance. If you type coeff and V in the Matlab terminal, you may find similarity. Recall that the eigenvectors are like the coefficients to predict the new scores from the original data (see Step-4 above!). But eigenvectors scaled up by the variances are called loadings. The terms can be confusing because the Matlab command already does normalization by default.

The newdata refers to the scores, which is a 2 × 100 matrix that represents the transformed data in the principal component space. The maximum number of principal components equals the number of variables in X. The component variances or latents are actually the eigenvalues of the covariance matrix of X. The latent values are arranged in descending order according to the magnitude of variance, i.e. ranking the eigenvalues. The function also returns explained, i.e. the proportion of the total variance explained by each component; and mu, the estimated mean of each variable in X.

To recover the original data, simply multiply the component scores with the transpose of matrix V. Remember to put back the mean value to each of the variables.

% Recover the original data
origdata = fliplr(newdata) * V' + mu;

In R, PCA can be performed automatically using a function princomp (using eigen-decomposition) or prcomp (using the SVD operation, for better accuracy in the results) from the built-in stats package. We will discuss one example in R below.

3. Case-study Using PCA
One nice package to easily extract and visualize the results of PCA is by using functions from the "factoextra" package that is based on "ggplot2". As the demo dataset, I am using the decathlon2 dataset from the "factoextra". Briefly, it contains a list of 27 athletes who have performed several different types of sports such as javelin throw, long jump, and so on. We call the athletes as individuals and the sports performance as the numeric variables. The dataset can be divided into:
- Active individuals (row. 1 to 23).
- Active variables (col. 1 to 10), to be used in PCA.
- Supplementary individuals (row. 24 to 27).
- Supplementary variables (col. 11 to 13), used in later prediction.

We first perform the PCA on the Active individuals/variables with the scaling option. This is an important step when you have values of different scales, e.g. metres and feet. After that, we can use the spree plot to visualize the different principal components ranked according to the proportion of variance explained. Recall that the eigenvalues reflect the amount of variance retained by the components. Eigenvalues are large for the first PC and small for the subsequent PCs. The first PC corresponds to a direction with the maximum amount of variation in the data.

We can ask other questions arising from the PCA results:

(1) How do we know the variance accounted after PCA?
Get this by extracting using the summary() function, we can see that the PC1 accounts for 41% of the total variance in the data, PC2 accounts for 18% variance, and so on. This is also what fviz_eig() provides. There are up to 10 original variables so we get up to 10 principal components.

(2) How many principal components are useful?
PCA does not tell us # components are useful, so there is no magic answer for this. From the cumulative proportion as shown above (or the spree plot), we can safely say that the first six components (PC1 - PC6) are indeed sufficient as they account for more than 90% of the total variance in the data.

(3) How are the original axes look like in terms of the components?
Before PCA, we have in the dataset the original axes which are formed by the 10 variables. After PCA, we get rotated axes which are formed by the principal components. To extract the rotated variables, we write get_pca_var(res.pca). Note that Dim.xx and PCxx are the same thing!

We can then plot how the original variables look like in terms of the first two components, PC1 and PC2. To visualize can be achieved using fviz_pca_var(), which will output the graph as shown below, Fig-2a on the left panel. The plot shows a correlation plot or relationships between all variables based on the first two PCs. Negatively correlated variables are positioned on opposite sides of the quadrants. The distance between variables and the origin measures the quality of the variables on the new map. Variables that are away from the origin are well represented on the map.

Fig-2: (a) Left: A correlation circle with a radius of 1.0 shows how the original variables relate to one another in the new dimensions of the principal components. (b) Right: a graphical representation of how each variable contributes to the principal component. It can also be seen that the representations of the original variables to PC7 to PC10 are minimal.

(4) What is the representation of each variable to the principal components?

We can use the function from the "corrplot" package to highlight the most contributing variables for each component. The function takes in var$cos2 as the argument. This cos2 shows the quality of representation of the variables on the component map (square cosine, squared coordinates).

Refer to Fig-2b. There are 2 ways to interpret this figure. First, it indicates a degree of representation of the original variables (the rows) on each principal component (the columns). Both panels in Fig-2 are somehow related. If a variable has a high cos2 value, it will be seen as a large and dark blue circle in Fig-2b, the same variable is positioned closer to the circumference of the plot in Fig-2a. A small circle indicates that the variable is not perfectly represented by the PCs. In this case, the variable is close to the center of the circle. Another way is to say how each principal component contributes to the variable. So, for a given variable, the sum of the cos2 on all the principal components is equal to one.


For further reading, see below:
(1) PCA in 6 steps
(2) Tutorial on PCA
(3) Making sense of eigenvector/eigenvalue

Additional Notes on fMRI

TR and TR in MRI

When deciding a suitable MRI sequence, you have to decide the choice between short or long repetition time (TR) and echo time (TE). Different combinations of these two values produce different types of MR images. See below. To recap, the T2* image has transverse magnetization which decays much faster than would be predicted by natural atomic and molecular mechanisms. In practice, T2* can be considered as an effective/observed T2 image, whereas the T2 image is considered the natural or true T2 (i.e. T2* is always less than or equal to T2). The T2* images result principally from inhomogeneities in the main magnetic field.

Statistical power in MRI

Statistical power is defined as the probability of rejecting the null hypothesis when it is false. In most cases, it is dependent on three main factors: (1) the desired effect size; (2) alpha value or Type I error; (3) the sample size tested. Of these 3 values, the sample size is the one that is within the control of the experimenter (Desmond & Glover, 2002). Unlike in usual behavioural tasks, fMRI data of each voxel is typically normalized as percent signal % change between the Task and Rest conditions. Variability in the fMRI data consists of the within-scan and between-scan variability. Increasing scan time and the number of subjects will potentially reduce these two kinds of variation. 

The concept of statistical power in fMRI data is complex because it involves different expected activation patterns, analysis pipelines, and also scanning parameters. Loosely, some rule of thumb to obtain a decent signal, alpha = 0.05 and 80% power. For example:

a)  Total N = 12 - 20 subjects per condition.

b)  Scan time per subject = 10 - 15 min.

c)  Scan volume (time points) = 400 - 500.

Non-parametric Approach to MRI

Like the usual inferential statistics, we want to make a conclusion in our neuroimaging data that apply to the population. As the name implies, non-parametric tests do not depend on the nature of the distribution of the fMRI data (e.g. normality assumption of the residual component). Through resampling, we construct an empirical distribution from the data itself that can accurately estimate the probability of a particular observance. This is the fundamental idea of non-parametric statistics. Some of these approaches have been clearly discussed in Nichols & Holmes, 2001. The use of non-parametric statistics if the data is not normal has been discussed briefly in this blog previously. A more thorough review article that includes the history and logic behind non-parametric tests is by Winkler A, et al. (2014) in the NeuroImage journal.

Basic Concepts
To see why this makes sense, the task-based block design paradigm will be used, i.e. a subject performs a task in one block and rest in the other. If a specific location of the brain is used during tasks, we should expect a task-related activation during the task, but not during the rest period. Here, we reshuffle or scramble the data across TRs. If a particular voxel is involved in the task, we should not expect any strong correlation between the predicted BOLD response and the actual activation because the 'reshuffled data' is just random noise. Reshuffling does not change the mean and variance of the dataset, that is, we don't tamper with the data.

In analyzing task-based fMRI data, for example, you form a design matrix according to your behavioural task and conduct a GLM (using FEAT, let say) to obtain parametric estimates or PEs of the corresponding regressors. To correct for the familywise error, the multiple comparison problem across voxels, we then use GRF theory by defining cluster-forming threshold (Z = 2.30 in FEAT). According to Echlund et al, GRF performance in controlling Type-I Error is poor. In fact, it is recommended to use non-parametric tests such as permutation testing to control for Type-I Error.

Using randomise on my dataset
For my iMac machine purchased in 2011 with Intel i5 Quad-core, 3.1 GHz, it takes me 8-9 minutes to analyze a decent task-based, 100 x 100 x 63 x 250, preprocessed 4D NIFTI file using 1000x permutation. With FEAT-glm using the same GLM design matrix, it takes me 2-3 minutes.

Registration with ANTs

When I first embarked in the fMRI project, I didn't learn about other software package. Recently, I managed to compare the performance of T1 -> MNI normalization done with 3 different popular software package: ANTs, FSL, and AFNI. The performance was measured using "similarity metric" that I took from Nipype, a python-based neuroimaging package. My data came from MRI scans of stroke patients, as part of my project at the Jewish General Hospital. As a reference image, I took the MNI-152 with 1mm resolution. The left panel below shows what happened when I used MNI-152 2mm template instead. In both FSL_1 and FSL_2, I used the 2 mm template to obtain the non-linear transformation matrix. However, in FSL_1, I explicitly specified the 1mm template when I applied the non-linear warping (with applywarp) process.

There are two main types of cost function: intra-modal (least squares and normalised correlation) and inter-modal (correlation ratio and mutual information-based options).


The left panel shows the registration accuracy of the subject's T1 image to the standard MNI 1mm template of 3 different software packages. The right panel also shows what happened when the MNI 2mm template was used instead.

As a side note, the MNI template has been widely used in various neuroimaging software packages. It was historically produced in a report for the International Consortium for Brain Mapping (ICBM) back in 2001. The atlas was based on linear mapping of 152 healthy adults ranging from 18-44 years old, mostly white people. It has gone through 2 major revisions lately to increase the accuracy: the first one was in 2006 that made use of a non-linear mapping, they call this the MNI-152 6th generation. FSL uses this template. Another revision was the most recent one in 2009.

Recent development in functional MRI

More recent contributions in functional MRI have been based on graph theory (Bullmore & Sporns, 2009), which underlines the importance of seeing the brain as complex networks. This analysis of MRI data has been widely used in both healthy and dysfunctional brain networks. 

It is important for understanding normal brain development and functions, the networks involved in affect regulation, motor control and execution, and learning and memory (e.g. motor learning study, see Sami and Miall, 2013). In addition, the graph theory has been beneficial to understand the implication resulting from e.g. Alzheimer's disease, schizophrenia, and stroke. For example, Buckner et al. found a correlation between the site of the targeted regions and the location of major hubs in Alzheimer’s disease (2009). He et al. (2007) demonstrated how the frontoparietal network is implicated in the spatial neglect of stroke patients.

In order to establish a brain network using modern graph theory, there are a number of steps to be taken: 

- define the network nodes (usually from resting-state fMRI), 

- estimate association/correlation between nodes, 

- compile pairwise associations between nodes and generate an association matrix, 

- and, calculate the network characteristics. 

Yet, despite this emerging literature, many researchers using graph theoretical approaches for fMRI data have not closely examined whether their data violate the assumptions of graph-theory. One important assumption that has to be met is the assumption of stationarity. Time series or time-varying data is said to fulfill the stationarity principle if the statistical properties of a time series are time-invariant, meaning, the mean and variance do not change across time. Generally speaking, stationarity implies that the statistic parameter of interest does not change over time. Such assumption is also essential for analyzing fMRI time series in the frequency domain, as the Fourier transform is suitable for stationarity. For more than a decade, most scientists have always assumed that fMRI signals are stationary, however, one may refer to some recent works that warrant doing fMRI analyses with caution (e.g. Muhei-aldin et al., 2014; Ou et al, 2014).

Logistic Regression: with MLE and more

The following post outlines the concept of logistic regression with the well-known Generalized Linear Model or GLM. Detailed maths is out of scope and we'll just make use of R to aid our learning. To begin, logistic regression is used to model and predict data with the following properties:
     • A categorical dependent variable, e.g. cure vs sick, old vs young (binary values)
     • A continuous, quantitative independent variable

In the discussion herewith, the response variable has 2 possible values, e.g. male/female, success/failure, dead/alive, etc. Our main goal is to predict the probability of either value occurring. A series of such binary responses (Bernoulli process) follows a binomial distribution. Probability density function (pdf) tells us the probability of a particular value for a fixed or known model or distribution, e.g. binomial distribution. Suppose the probability of having Y is denoted as P(Y) or p, where 0 ≤ p ≤ 1. The ratio of having Y and not having Y, p/(1 – p) is called odds. Thus, instead of simply predicting Y from X as in the linear regression, we predict the odds of Y for a given X. 

Logistic Regression as a Special Case of GLM

A simple logistic regression related to P(Y) can be expressed in the form of Equation (1). The word 'logistic' comes from the logistic function of probability shown on the left-hand side. This function can be transformed to the natural logarithmic function of odds called logit as shown on the right-hand side. The transformed model now takes logit as the dependent variable, which reminds us of the linear regression equation. Although the assumption of linearity (and others e.g. normality, homogeneity of variance) does not hold for logistic regression, the linearity between independent variable and logit has to be maintained! An extended version with multiple predictors can be seen in Equation (2).

Fig-1: (Left) An example of a logistic function σ(t) that takes any value between 0, 1. The function has a unique sigmoidal shape. (Right) An example of logistic regression with different values of betas [figure taken from ips8e, supp-Chapter 14, Macmillan].

When we talk about the simple linear regression, we assume X and Y to be linearly related more or less, the random error or residual follows a normal distribution, and the assumption of independence holds true. What happens when some of the assumptions collapse? "Generalized Linear Model" attempts to provide a more general form of regression when such assumptions are not fulfilled. It's sometimes confused with the general linear model found in the ANOVA and neuroimaging analyses. Now, there are three components of GLM:
     • An exponential family distribution for residuals.
     • A set of linear predictors that represents a systematic component.
     • A function that links or connects the expected value or mean of (1) and (2).

In GLM, the outcome variable takes up an exponential family distribution, e.g. normal or Gaussian, beta, gamma, exponential, chi-squared, Poisson, Bernoulli, and Dirichlet distribution. In the case of the ordinary linear regression, it is the Gaussian distribution. Suppose we predict the expected value of the outcome variable as shown by its mean μ. GLM allows a function of μ rather than the mean itself. This function, called the link function, is denoted as g(μ) because it links the mean of the outcome variable to the set of k predictors or explanatory variables.

      In summary, GLM formula states that:    g(μ)  =  a +  b₁X₁  + b₂X₂  + b₃X₃ + .... +  b_kX_k

In the case of linear regression, the link function is the identity link where g(μ) = μ. For binary data such as success/failure, a logit link is used instead. This is the case of the logistic regression. In R, logistic regression can be done by using glm command, instead of lm( ). Now, of so many exponential families, is there a unified way we estimate the parameters in GLM? Yeah, using a technique described in the next section, MLE.

Maximum Likelihood Estimation (MLE)
How do we estimate b? Whereas linear regression uses the least-squares method to minimize the errors, logistic regression uses the maximum likelihood method to arrive at the solution. Key to MLE is really just the "likelihood". Here, we first assume there is a model and this model is valid/correct. If it's correct, we can find the most likely estimate that explains the observed data.

Probabiliy means given a model of an event, what's the distribution of the outcome? Likelihood means given a set of observed outcomes, what are model parameters? 

Undergraduate stats course teaches us binomial probability. For example, suppose we toss an unfair coin 20 times and trials or tosses are independent of each other. In each trial, the probability of a head is p = 0.6, so the probability of not getting a head is 0.4. We try to describe the probability of getting heads r times out of n tosses. In fact, we may obtain a total of r = 0, 1, 2, ..., up to 20 heads from tossing that coin 20 times. We compute,  ₀C²⁰p⁰(1 – p)²⁰ ,  ₁C²⁰p¹(1 – p)¹⁹  , ... up to ₂₀C²⁰p²⁰(1 – p)⁰. Refer to Fig-2. The plot on the left is the binomial probability density function. The shaded area under the curve equals 1.

Fig-2: (Left) Discrete probability density function of tossing a coin 20 times, with the probability of obtaining head in each toss p = 0.6. The X-axis is the number of possible heads from 0 to 20. Obviously, if observing a head in each toss is 0.6 then I would expect 0.6 x 20 = 12 heads in 20 tosses. (Right) Suppose we were to estimate p here based on the observation that there are 13 heads out of 20 tosses. The best estimate of p is the one that has highest likelihood, in this case, the peak is at p = 0.65.





The step shown above is actually generating a possible data set from known parameters (n, p). What if the process is reversed? Twist your mind! We want to estimate p for a given series of n observations or outcomes. A likelihood function is exactly this, i.e. the opposite of probability density function. So, if we observed 13 heads in tossing the coin 20 times, what is the probability of a head in each coin toss? We try to estimate the value of p such that our observation is most likely to happen. Theoretically, if our model is correct, the most likely value of p = 13/20 = 0.65. The right panel Fig-2 shows this is the case. The highest likelihood is when the X-axis or p equals 0.65.

L(data | model) is a function that describes the probability of a model given our observed data set, p(model | data). If independence among trials is assumed, the likelihood function equals the joint-probability density which is none other than the product of each and every p(y_i | x_i), where y_i takes up a binary value.

As the end result of multiplication may be too small to manage, people take the natural logarithm of the likelihood which is then called log-likelihood (LL). Twist your mind again, taking the log of multiplication yields a summation of logs. Look at Equation (3). Now, since the numerical optimizer is designed to find the minimum, not the maximum, we multiply this by –1. Equation (4) is the final equation we want to minimize, just like minimizing sum of squares residuals in linear regression! Numerical optimization is used as there is no closed-form analytical solution. In R, this can be achieved by using optim() function.

Fitting a model to our data in logistic regression means to estimate the coefficient b₀ and b. The first is the intercept, the value of Y when X = 0. But how do we interpret b? The logit part is our dependent variable in the logistic regression. The logit is also easily converted back into the odds, where odds = exp(b₀ + bX). Suppose we have X₁ and X₂ = X₁ + 1; we compare both odds by finding the ratio of odds(X₂) = exp(b₀ + bX₁ + b) and odds(X₁) = exp(b₀ + bX₁), which is just equal to exp(b). Hence, the coefficient b can be interpreted as the log of the relative increase in the odds for a unit increase in score X. 

Assessing the model fit

There are a few common measures of model fit in the logistic regression. The first is the log-likelihood statistics as mentioned previously, i.e. summing the probabilities of the predicted and actual outcomes. This is analogous to the residual sum of squares (SS_R) in multiple regression in the sense that it is an indicator of how much unexplained information there is after the model has been fitted. The larger the value, the poorer the fitting is. A similar measure called deviance statistics is more popular, where deviance = –2 LL. Again, a higher number means a bad fit. We'll talk more about deviance later.

The second measure is Pearson's goodness of fit test. This is a statistical test to investigate how close values predicted by the model are to the observed values. The null hypothesis here is that our data follow the logistic regression model. We first find the chi-squared (c²) value and then the degree of freedom to find the p-value. The third performance assessment is by R², which is related to both the deviance and z² (Wald statistics). This measure is similar to the one we saw in the linear regression but with a different formula. R².can be unnecessarily inflated when we add more predictors. Like in multiple linear regression, we can use BIC or AIC to compute the parsimony measure of R², penalizing for useless extra predictors in the model. AIC/BIC is useful for choosing a model, with the preferred model is the one with a lower AIC/BIC index.

Software packages such as SPSS and R make use of the deviance statistics in the report. More importantly, there are two types of deviance. The null deviance shows how well the response is predicted by a reduced or baseline model. In the case of the simple logistic regression, the model assumes that the outcome can be predicted by nothing but a constant value, that is, the intercept b₀. In other words, it agrees with the null hypothesis that the estimated b = 0. The residual deviance shows how well the response is predicted by the full model when all predictors are included, implying that this is against the null hypothesis. Each deviance carries its own # dof.

The difference between the null and residual deviance is called likelihood ratio, another goodness of fit that also has a chi-squared distribution. The degree of freedom is defined as the difference between two deviances. We then easily get the p-value that suggests whether the model is better than chance at predicting the outcome.

How to assess the importance of the individual predictor in the logistic model? Similar to linear regression, it is assessed by carrying out statistical tests of the significance of the coefficient. Whereas t-test is used in linear regression, Wald test (z²) is used to evaluate the statistical significance of each predictor. It is calculated by taking the ratio of the square of the regression coefficient to the square of the standard error of the coefficient. The null hypothesis is similar, i.e. whether the estimated b is different from zero. As a side note: Wald test has been shown to be less reliable for small sample sizes!

Fig-3: Two R outputs from a logisitc model operation where you have one predictor (intervention) and two (intervention + duration) to predict the same dependent variable. Click to enlarge!


Let's take a look at Fig-3. Suppose there is a data frame containing three fields: whether the disease is cured (dependent variable), whether any intervention is performed (first categorical predictor), and the sickness duration (another, but quantitative, predictor). The model on the left fits the outcome variable (cured vs. not-cured) with whether or not treatment has been performed. What does R show us? First, the coefficients portion suggests there is a significant intervention predictor with estimated b = 1.23, z = 3.07, p < 0.005. This is Wald statistics. A z-value that is sufficiently far from 0 means that the estimate is both precise enough to be statistically different from 0 and large to have an effect on the response.

Next, the model fit. The null deviance = 154.08 and residual deviance = 144.16. The fact that residual deviance is lower suggests that one more predictor is better than just a constant in reducing the residual error. So adding "intervention" improves the fit by 9.93 unit. The difference in # dof of 1 so we can find the p-value of the chi-squared. Accordingly, c²(1) has p-value = 0.00163. The model on the right introduces another predictor, i.e. total duration. Based on the deviance values, there is no benefit of adding one more predictor. The book says that anova() can be used to compare both model fits. The test yields no significant difference between the left and right models.

Special example: Psychometric curve
A fundamental concept in psychophysics, the psychometric function relates a parameter of a sensory stimulus to a subjective response of a participant. We learnt previously how the psychometric curve has a sigmoidal shape and it's getting much clearer why. Take an example of a two-alternative-forced-choice (2AFC) task. What we're interested in is how the probability of responding with one of the two choices varies with the stimulus. Here, we will consider a situation where the varying sensory stimulus is the movement direction θ of a participant's arm to the left or right of the body midline. This becomes the continuous independent variable X. The desired response is binary, either left or right. This response is our dependent variable Y.

The model suitable for this is the simple logistic regression shown in Equation (1). Let p = P(y = right | X) is the probability of responding "right" given a certain direction X. So, 1 – p is equivalent to P(y = left | X). We typically represent the response in numerical format, 1 = right and 0 = left. The value X in which p = 0.5 corresponds to the perceptual threshold that in ideal case is the body midline, but one may have a perceptual bias. From our experiment, we obtain an array of direction X (in degree) and an array of binary responses R (in number 0, 1). First construct our cost function, name it NLL. This is then fed into a numerical optimizer to estimate b₀ and b using MLE. An initial guess for each estimate has to be given, e.g. (0.1, 0.1).

NLL <- function(B,X,R) {
          y = B[1] + B[2]*X
          p = 1/(1+exp(-y))
          NLL = -sum(log(p[R==1])) - sum(log(1-p[R==0]))
  }

out = optim(par=c(-.1,.1),NLL,X=X,R=R)    #Numerical optimization! 
Bfit = out$par    #Retrieve parameter estimates of our model

#Let's construct our model with Bfits above and plot the data...
Xp = seq(-15,15,1)
myModel = logistic(Bfit[1] + Bfit[2]*Xp)

PS: I have recently found another nice R package called "quickpsy" for computing with psychometric functions.


Main references
(1) Andy Field, et al, "Discovering Statistics using R" (2012).
(2) Gribble's note on MLE, Psychology_9041B course.

Linear Regression - The Basics and More

Sigh, I should have studied this even before learning about the fMRI statistical pipeline. In itself, the scope of regression and GLM is super wide. In this blog post, the emphasis would be on concepts, not calculation, since statistical software such as R/SPSS/Matlab is able to help us.

Simple Linear Regression
Linear regression is one of the most basic tools to model data. In its simplest form, the model can be represented by Equation (1). This is because the equation can be used to find or “predict” the observed variable Y based on X, and so independent variable X is also known as a predictor. The point where the line crosses the vertical Y-axis is called intercept β₀ of the model. The parameter beta or β is also called the gradient or slope of the equation. The slope represents how much the Y changes resulting from a unit change in X. Both coefficient β and β₀ are called regression coefficients. If there is more than one predictor, the term multiple regression is used instead, see Equation (2). Here, we attempt to predict Y based on several predictors such as in the case of functional MRI statistical analysis. Using multiple regression performed at every brain voxel altogether, I have discussed how a task-activation map can be produced.

(1) How can we estimate the coefficients?
Performing a linear regression will also mean finding the best fit model or Ŷ (the ^ hat sign indicates that the variable is an estimated or predicted value). We position or "fit" the line where the deviation from each data point is the least. The regression line used to model the actual data set does not obviously pass every data point. Whereas variance is like the deviation from the data center (mean), deviations from the regression line are called residuals. A method. called the method of least squares, which aims to get a numerical estimate of β that minimizes these deviations or residuals. The Gauss-Markov Theorem has been presented when talking about fMRI statistics. The estimated beta is sometimes named parameter estimate, PE. Briefly, the formula is shown in Equation (3). Assumptions are made for this formula to be valid. In fact, for an accurate model, several assumptions have to be held:

• The outcome variable has to be quantitative, continuous, and unbounded.
• All predictors must be quantitative or categorical (at least two categories). They should not have zero variance. 
• The relationship between predictors and the outcome variable is linear.
• Residuals Σe = σ².I are uncorrelated and independent of each other. They are normally distributed with zero means.
• Residuals maintain homogeneity of variance.

(2) How well does the model fit my data?
The model performance can be quantified by using three different sums of squares defined as follow:

1. Residual sum of squares, SS_R = deviation of the observed data from the fitted line. This is actually an error or a measure for the inaccuracy of the model. Hence, it represents a portion of the data that cannot be well predicted by the line. This portion is minimized through the method of least squares. Refer to the red solid line in the figure above.
2. Model deviation sum of squares, SS_M = deviation of the predicted Ŷ based on the model from the mean of the data set. Refer to the green solid line in the figure above.
3. Total deviation sum of squares, SS_T = deviation of the observed data Y from their mean. This is simply the variability of the data set around its mean. 

Variance partitioning says that SS_T = SS_M + SS_R . The total variance in our observed data can be decomposed into two parts: the portion explained by the model and by the residuals. Hence, SS_M signifies the merits of using the regression line. The ratio of SS_M  and SS_T is called the R², or r-squared. It is a statistical measure of how close the data are to the fitted regression line. Generally speaking, between 0.16 ‒ 0.25 is considered a good model fit. There is also another measure called the F-ratio which is defined as follow:

F-ratio is a ratio between how good the model is with how bad it is (residual). Mathematically, it is the ratio between mean-squared error of the model MS_M, and mean-squared error from the residual, MS_R. F-ratio is used a lot in analysis of variance (ANOVA) and the significance of the value can be assessed using a F-distribution table.

Actually, a mean-squared error is the sum of square divided by a degree of freedom (dof). For a simple linear regression with just one predictor, the model has # dof equals 1 and the whole data set has N ‒ 1. The # dof of the residual is, therefore, N ‒ 2.

Residuals are important to identify poor model fit. Mathematical software such as R is able to give us a summary of the linear regression, lm. Residual can be quickly computed using resid. The sum of all residuals and the expected value E[res] equals zero. Lastly, residuals can be thought of as outcome variable y having removed the linear relationship with the predictor x.

(3) How well does the predictor do the job?
To evaluate further, we can also test whether X is a good predictor of Y. Because the slope β represents how much Y changes for every unit of predictor X, a lousy predictor means a flat or very small slope. Using one-sample t-test on the slope β (technically on the parameter estimate of that), we test whether it is significantly different from zero. Refer to Equation (4) in the context of simple linear regression. The denominator is the standard error that tells us the variability β's across different samples and σ_R² is the residual variance. Commonly, a significant p-value of the t-test suggests that the predictor is reliable enough to predict the outcome variable.

Multiple Linear Regression
When there is more than one predictor, multiple regression kicks in. For example, we want to model or predict electricity bills during winter as a function of outdoor temperature and oil price. The basic principle remains, R² represents the proportion of the variance explained by the model with the total variance in the observed data set. Note that more predictors means higher variance explained by the model, artificially inflating the R². Some adjusted measures based on AIC, BIC, etc. are used as parsimony adjusted measures that penalize you for adding more independent variables in the model unnecessarily. As a rule of thumb, higher AIC/BIC means worse fit. Adjusted R² tells you the percentage of variation explained by only the predictors that actually affect the outcome, dependent variable. It is computed using the formula 1 – ((1 – R²)((N – 1) /( N – k – 1)) where k is the number of predictors. The more junk variables, the lesser the value will be.

Fig-1: Two R outputs from a linear model operation where you have one predictor (IV2) and three (IV1, IV2, IV3) to predict the same DV. Click to enlarge! SPSS, like R, is able to produce similar detailed regression results.

In Fig-1, I want to predict DV using one (left) and 3 predictors (right panel). The beta coefficients are shown together as "Estimate" with their standard error, t-statistics, and p-value results. In R, these p-values are two-tailed. On the left panel, a positive slope means that there is ~0.27 increase in IV1 for 1 unit increase in DV. A very small p-value indicates that we cannot retain the null hypothesis that the β equals zero. Also, the summary tells us about the residual component, the multiple and adjusted r-squared. Look at how adding 2 extra predictors reduces our degree of freedom to 19, but does not give many benefits to the model. First, the residual is not dramatically reduced. Second, the left panel has a greater reduction in adjusted r-squared. The last line tells us about the F-statistics together with their p-value. A good model should have a high F-value with a low p-value against the null hypothesis. For a simple linear regression on the left, the p-value of the F-statistics equals the p-value from the t-statistics. This is unlike the regression with 3 predictors on the right.

In the context of multiple regression, we can use F-statistics similar to that used in ANOVA. The null hypothesis is that all β_i equals zero. The alternate hypothesis is that at least one β_i > zero. The F-ratio formula doesn't change, but the degree of freedom of the model and residual would now be k and N ‒ k ‒ 1 respectively, where k is the # predictors.

In multiple regression, the statistical test on the betas can be further analyzed according to different contrast. Such term may be familiar in ANOVA. For example, if there are 3 predictors, the contrast of [~1 0 +1] attempts to test whether the difference in the first and last β values is significant.

(1) Parameter estimates in multiple regression
How to find β-parameters in multiple regression? Exactly what we learnt in the fMRI statistics. We can derive betas by minimizing ∑(Y ‒ X₁β₁ ‒ X₂β₂ ‒ ... ‒ X_k β_k )²  just like before. The formula for R² and F-ratio are similar but DOF of the model and residual are different. Conceptually, the estimate for β₁ is the regression through the origin estimate having [X₂, X₃, ... , X_k β_k ] regressed out of both Y and the X₁. Similarly, the estimate for β₂ is the regression through the origin estimate having [X₁, X₃, ... , X_k β_k ] regressed out of both Y and the X₂. Regression through the origin happens when Y = 0 for X = 0.

The left graph below depicts the first model with IV2 only, the slope is +0.2712. The second model has IV1, IV2, and IV3 to predict DV. Here, I'm plotting IV2 and DV but with the influence of IV1 and IV3 removed. It can be seen that the slope becomes negative, ‒0.382776. Hence, the meaning of coefficient in multiple regression is the expected change in the outcome or dependent variable per unit change in one predictor, holding all other predictors fixed. By holding those predictors constant, it is said that we have controlled for other predictors in the model. Each regression coefficient is also called partial regression coefficient. Note: italicized words in brown carry similar context!

Fig-2: The slope between IV2 and DV is positive (left). When two other predictors are introduced, the same slope changes to negative (right). The result in the right panel is when we partial out IV1 and IV3 from the model. In this way, we say we have "adjusted" the model that only IV2 contributes to DV.

If we are to make a model with several predictors, there are several methods available. In hierarchical regression, the most important predictor is entered first followed by other predictors. But in what sequence? Based on past literature or if we have a strong hypothesis in mind. After the known predictor is included, we can add new predictors either in one go or stepwise manner. Then, the R² is checked each time to see whether adding predictors reliably improves the model. If you want to include all predictors at one go, this is called forced-entry method. The most exploratory method that receives less scientific respect is actually stepwise regression where you start with the whole bunch of predictors, and predictors that contribute the least will be removed step by step. Here, what eventually remains is a set of "good" predictors.

Lastly, there is a predictor that helps to increase the variance explained in the model when placed together with other predictors, although in itself it isn't correlated with the dependent variable. This is called suppressor variable, since it suppresses the residuals.

(2) Adjustment and interaction
Adjustment, is the general idea of putting predictors into a linear model to investigate the role of a third variable on the relationship between another two. Example, suppose you recorded blood pressure (continuous var) based on two types of participants, those who do and don't seek treatment (binary, 0/1). A third, continuous variable x, is bodyweight that is correlated with blood pressure. We may model our data with  Y = β₀ + β₁X + β₂T + e. Refer to Fig-3. Here, the treatment status T and their fitted line are depicted as a different color. The horizontal dotted line represents the mean of each group when we disregard X. Try to regress X out of Y in R, and see whether the distance of two dotted lines is similar to the difference in intercept!

Fig-3: You want to predict blood pressure (Y) from treatment effect (T). A third variable X (body weight) carries a linear relationship with blood pressure as well. Now, each treatment status has the same relationship (slope) with Y. The black line is what you would get if you just fit X and ignored the group.

The predictor X is unrelated to treatment status (color). The groups have different intercept that depends on treatment status. So, β₁ represents the difference in intercept, and β₂ is the slope common to both treatments. Notice that the estimated relationship between the group variable and the outcome doesn’t change much, regardless of whether X is accounted for or not. This is seen by comparing the distance between the horizontal dotted lines and the distance between the intercepts of the fitted lines: blue and red. That the relationship doesn’t change much is ultimately a statement about balance. The nuisance variable (X) is well balanced between levels of the group variable. So, whether you account for X or not, you get about the same answer. One way to try to achieve such balance with high probability is to randomize the group variable. This is especially useful, of course, when one doesn’t get to observe the nuisance covariate.

(3) Issues: outliers, confounds, and multicollinearity
Outliers mislead people in finding parameter estimates and eventually drawing conclusions. That's why, scatter-plot is our good friend in correlation/regression analysis. There are a few diagnostic tools for quantifying an outlier. We can find the distance from an outlier to the fitted line, e.g. Cook's distance, dfbetas() in R, etc. A confound is defined as a variable that is correlated to both your outcome and predictor variables in a way able to explain the relationship. Not accounting for this confound may introduce bias to the estimates.

Multicollinearity arises because there is a strong correlation between two or more predictors in the model. Look at the case of Fig-4 where IV1 and IV2 has a very strong correlation. Each predictor alone is able to reliably predict DV with an almost identical slope. When together, IV1 and IV2 have different slope estimates, but very low p-values > 0.10 because the variance of the unexplained component (residual) is too big. In multiple regression, the problem of multicollinearity makes it difficult to see the isolated effects of each predictor on the dependent variable. In practice, this can be diagnosed by first computing the correlation matrix and see which predictors are correlated one another. It can also be diagnosed systematically using variance inflation factor (VIF), where VIF >> 10 usually means high multicollinearity. Refer to Fig-5 to visualize multicollinearity using a Venn's diagram.

Fig-4: What happens when > 2 predictors have a strong correlation. This concept becomes useful when dealing with fMRI data!

Visualization using Venn's Diagram
I like this concept; learning multiple regression becomes easier using Venn's diagram. Each variable in the regression model is represented by a different circle: a blue circle normally for the dependent variable, while red circles for independent variables/predictors. The radius or size represents the variance. The manner in which the circles overlap illustrates how they behave in the same manner or covary (covariance). E.g., area A in the left panel = cov(X₁, Y). It tells us the variance of Y explained by X₁. Note how r-squared is computed visually and how it becomes tricky when two predictors are present.

Computing parameter estimates is our next topic. The shaded area in the left panel shown as green is relevant when we want to explain Y using X₁ alone. When predicting Y from both predictors altogether, the green shaded area is the only relevant information for estimating the slope of X₁. Why? As mentioned before, β₁ is the regression through the origin estimate having X₂ removed from both Y and the X₁. To visualize this situation, we remove (regress out, partial out) areas A + A" + D" from both Y and X₁ (look at the inset on the right!!).

Fig-5: With Venn's diagram, we can try to explain multiple regression. Example, the blue shaded area on the right panel shows variance of the observed data (y) explained only by x₂, not x₁. Note that this diagram does not provide the sign, meaning there's no way to see negative covariance and parameter estimates, or suppressor effect. The lower left panel shows what happens when two predictors are related. The residual part increases, but the variance of each predictor alone decreases (B >>, C <<). This causes inflation in standard error of each estimates, and a big drop in the t-statistics.

Suppose we start with a predictor X₁ and an outcome variable. The bias in slope estimate only occurs when we fail to take into account another predictor, say X₂, that is correlated with both the dependent variable and one of the predictors (also called the confounding variable). Look at the lower right panel of Fig-5. Technically, bias will not occur if there is no correlation between X₁ and X₂. However, it will give a lower t-statistics. Why? As with multicollinearity, the value of the denominator that's proportional to the variance of the residual σ_R² would go up.


Main references
(1) Andy Field, et al, "Discovering Statistics using R" (2012).
(2) Brian Caffo, notes from online Coursera, Regression.
(3) This link for visualizing multiple regression using Venn diagram.

Notes on Resting-state fMRI Analyses (Part II)

Fieldmap Correction and Coregistration
What is fieldmap correction?
Apart from having a poor resolution, functional images acquired using common EPI sequences suffer distortion due to magnetic field (B₀) inhomogeneity introduced by different tissue types in our heads. Such a thing occurs due to the existence of non-homogeneity in RF receive and transmit of the head coils. The more channels you have, the more inhomogeneity the image may have. The most severe inhomogeneity includes the air-bone or air-brain tissue interfaces in the sinuses in the inferior frontal gyrus and medial temporal lobes. This poses a serious effect on our data, a geometrical distortion and signal loss as depicted in Fig-1.

Magnetic field inhomogeneity can be measured with fieldmap images; which can give us a geometric distortion and signal loss. These values can then be used to compensate for the loss by geometrically unwarping the EPI images, and applying cost-function masking in registrations to ignore areas of signal loss. The correction is most useful during image co-registration as it dramatically improves the registration accuracy. Areas where signal loss has occurred unfortunately cannot be restored with any form of post-processing. In other words, it is impossible to recover time-series data in those locations.

There is no separate sequence for acquiring the fieldmap and different scanners give different images. The sequence can be EPI, Spin-echo, or Gradient-echo sequences, but it isn't recommended to use the EPI-based sequence since it will suffer the same problem. There exist 2 different methods of acquiring fieldmap images for the purpose of correction. 

When you do the fieldmap acquisition, you usually acquire two different images: a pair of magnitude images captured with different echo times, and a phase difference image (Fig-1). The acquisition can also be controlled either in the AP (j+) or PA (j-) direction. These images should be acquired in the same orientation as the target EPIs. The phase difference between the two images is proportional to the difference in echo time (ΔTE) and the B0 inhomogeneity observed. The fieldmap is calculated by taking the difference between the two-phase images, and dividing that by the echo time difference.

Method-2 is called the blip-up blip-down method, which calculates the fieldmap based on the difference in distortion between the two consecutive acquisitions. This method acquires two diffusion-weighted images (DWI) with opposite phase encoding directions, that is, the AP and AP directions. It is assumed that there is no change in the magnetic field and sudden motion during the two acquisitions. You can use TOPUP in FSL to help you do fieldmap processing using this method.

Fig-1: Images obtained from the scanner (left) are converted to get a fieldmap image (right). Red circles show distorted regions that require correction. This is a standard procedure of double gradient-echo performed in Siemens 3T scanner.

How to process this in FSL?
At the MNI, our brain imaging center uses Siemens 3T scanner, which is a good thing as FSL provides a ready-to-use tool, fsl_prepare_fieldmap, to obtain a fieldmap phase image in rad/sec. The magnitude image resembles a lower resolution version of the T1 structural (anatomical) image. FSL FUGUE, which is incorporated in FEAT, helps us to do distortion correction using this method. Both the complete and skull-stripped versions of the magnitude image and the processed phase image (rad/sec) should be defined in FEAT. FSL will then attempt to unwarp the distorted EPI image before mapping it to the structural image. The unwarp direction has to be specified and is typically given by the scanner operator depending on how the fieldmap acquisition is set.

In FSL, fieldmap correction is incorporated as part of the registration (preprocessing) pipeline. The highly accurate functional-to-structural coregistration is also called boundary-based registration or BBR  (Greve and Fischl, 2009). The method is based on changes in the intensity along the white matter boundaries instead of the less reliable grey matter boundaries. This means that an accurate segmentation of the structural image is required and bias-field correction reliably improves the accuracy. Performing BBR registration without a fieldmap correction doesn't give many benefits than the usual 6DOF method with FLIRT (Fig 2-3).

Also, there must be some grey-white intensity contrast in the EPI, though it doesn't have to be good enough for segmentation. The FSL website said since only intensities near the white-matter boundary are used by BBR, it is likely to be more robust to a range of pathologies and artefacts in the EPI or the structural.

Fig-2: Comparison of 3 situations using BBR coregistration with fieldmap correction: when there full magnitude image with the skull wasn't supplied to FSL(left); when the correct magnitude image was used but the unwarping direction was the opposite (middle); the correct BBR-registered image with a superior accuracy (right).

Fig-3: Performing coregistration of a functional image to a structural image using BBR is superior than the usual linear 6DOF registration in FSL. Note that the asterisk ( * ) sign indicates the region with severe signal loss. Without using the fieldmap correction, the corpus callosum mapping becomes inaccurate as denoted by a hex sign (#).

 

Slice Timing Correction?

Scientists more or less agree that the slice timing correction is important.  For the more recent multiband sequence, some experts said that slice timing misalignment may not have a huge impact on the analysis. In the earlier version of the Siemens WIP, I was told that the slice timing information contained in the DICOM files was not correct. This can be retrieved easily with a Matlab function. As a result, I didn't perform this correction in my fMRI paper (MB3, TR=1690 msec).

Until recently, one can deduce the slice timing information based on the CMRR Multiband protocol here. For comparison, I have included the effect of slice timing correction to my resting state data with an MB 3x acceleration measured on a single voxel.

Fig-4: Time series with and without the slice timing correction measured on a single voxel @ MNI coordinate (67,41,49).

Adding additional EVs to GLM

Additional regressors (EVs) can be added to the GLM in FSL FEAT. I find that the GUI is a bit tricky, better write a script for that. First, we have to recreate the design matrix by adding the extra regressors, using either:

       ⁍ Pointing to a file for each regressor by constructing a full model design    

       ⁍ Creating a space-delimited text-file comprising all confound EVs

Fig-5: If you click the "Full model setup", a new GUI will appear as shown on the left. Select an appropriate setting (number 1-3). You don't have to perform another temporal filtering. The temporal derivative is optional too. Another way is to construct a text file and select "Add additional confound EVs" (number 4). 

Refer to the Fig-5 above. If you click "Full model setup", a new GUI will appear. Choose the input file as 1-entry per volume (see number 1), no need to convolve it with the HRF anymore (number 2). Also, you do not have to perform another temporal filtering if the regressors are derived from the prefiltered data (number 3). This method is time-consuming. A better option is to use the additional confound EVs (number 4), just that the text file has to be space-delimited, not a comma-separated file! Using a wrong delimiter will cause FEAT to ignore these additional regressors! 

It is also safer to write a code or script rather than getting restricted with the GUI features. To create a design matrix, use the feat_model command. To manually perform GLM according to the design matrix with prewhitening, use the film_gls command. The command is equipped with sophisticated estimations of autocorrelation and Tukey tapering, which is very important for making statistical inferences in task-based fMRI.

One last note about the design matrix is about Orthogonalization. Most EVs generally are almost orthogonal, so enforcing orthogonality is not gonna help much in the results.

Denoising nuisance components with ICA

As mentioned often, rs-fMRI has one major drawback: the data is recorded at rest so it is prone to noise or artefacts. This is so because we don't have any reference pattern as we do when we perform task-based functional imaging. Hence. proper cleanup is paramount to getting a correct deduction or conclusion. This is even more serious for me who is doing learning-related rs-fMRI. One of the things I'm struggling with is choosing the best cleaning method! Technically for my project, I plan to try out the ICA denoising method since the work of our previous postdoc used a different technique.

As mentioned in the previous post, the subject-level data cleaning steps cover the following:

1. First, you perform ICA, e.g. using FSL MELODIC and identify the nuisance components. From my experience, the tool performs a pretty good job. Let the algorithm choose the best number of ICs.
2. Identify the nuisance components and run fsl_regfilt script in FSL, producing a so-called clean or denoised fMRI dataset. 
3. Now, to assess its performance after cleanup, you can either conduct another round of ICA on the residual image or compute the temporal standard deviation of the GM regions (or a specific ROI in the motor cortex).

The script above effectively removes nuisance components identified from the original dataset. If I conducted another ICA on the denoised data, the resulting components are much cleaner!

How does the script differ from the usual GLM function in FEAT? They are not quite the same. The core function of the GLM in FEAT is film_gls. The fsl_regfilt, on the other hand, uses a simple time-varying GLM function called fsl_glm, which does not perform any sophisticated modelling of temporal autocorrelation and whitening. I finally found this difference after struggling for so long!    

Fig-6: Demeaned time-series of the same voxel obtained from FEAT and ICA denoising tool.

The output of regression is called the residual image, res4d.nii.gz, which is supposed to be as clean as the denoised image, but having the mean removed (just to add back). See Fig-6 for representative time-series of a voxel from the residual outputs of each FEAT and fsl_regfilt. There are some minor differences. I think pre-whitening is not necessary because the nuisance EVs are all spatially independent from the MELODIC. 

On the other hand, I finally found out that the residual output of the fsl_regfilt is in fact the same as the one produced by the fsl_glm. The data has to be demeaned, and the design matrix regressors des_norm has to be normalized into unit variance. The FSL gurus in their forum claimed that both methods use the same GLM methods basically, but I'm not sure which GLM it was. Knowing this similarity is essential because now I can compare different methods of denoising, e.g. using WM/CSF average time-series.

To regress out nuisance components following either ICA or other noise modelling tool (e.g. RETROICOR), you can just use the simple a GLM function. Sophisticated estimation of temporal autocorrelation becomes important when you want to do statistical inference of neural activity or connectivity.