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.

Group Analysis in Functional MRI - a brief overview

Previously, I mentioned subject-level or first-level GLM analysis that is performed on individual subject data. This step gives us within-subject contrasts across parameter estimates, e.g. estimates of brain activity due to visual stimulation vs. rest. In this blog entry, a second-level analysis will be discussed which provides group-level inference on whether they are significantly different from zero, e.g., group status, behavioral performance, etc. Statistically speaking, the goal of the standard group analyses is to do inference of a set of acquired sample to a larger population. How can we perform group analysis in neuroimaging data? There are two ways we can model the "effects" of a behavioral study on the brain, that is, fixed effects and random effects. A statistic that takes both variations into account is called mixed-effects modeling in psychology or behavioral sciences. 

Consider the following scenario. You want to know the distribution of hair length between adult men and women, but doing this for all world population is just impossible. You start to select four men and women. For every single person, there is bound to be differences in hair length which represents a within-subject variability. In addition, different people in our study will have different hair cut or style, a form of between-subject variability. You can sample a strain of hair of each person or multiple hairs.

1. Fixed effects model: we assume that any systematic variability introduced in our measurement is due to within-subject variability only. In other words, the four men (or women) are more or less equal in terms of hairstyle.
2. Random effects model: this model assumes a more realistic scenario where a man (or woman) is randomly sampled from a population. The model takes into account another factor called "subject", and any unpredictable variation in the data is contributed by an individual or between-subject differences. 

In functional MRI, we employ mixed-effects modeling as accurate statistics when analyzing group data. This is done by taking sample voxel's time series per subject. Theoretically speaking, a fixed-effects model is used when a single subject has multiple runs. In order to assess multiple subjects that have multiple runs, we have to use mixed-effects modeling because every person is understood to be unique and we want to make a statement that can generalize to all people.

Fixed Effects Analysis
Suppose Y_i = X_ib_i + ε_i  is the fixed effects model of our neuroimaging data of N-subjects, each is acquired with n time points. The subject-level GLM consists of p parameter estimates. The within-subject residual ε_i has a multivariate normal distribution with a mean vector 0 and variance-covariance matrix of Σ_i = σ_w².I. Under the null hypothesis, H₀: c^Tb_G = 0, we can then compute:

which obeys a t-distribution with # dof = N(n – p). Note how we incorporate the group-level PE from the average of individual subject's PEs, and within-subject error variance from the pooled individual error variance.

Random Effects Analysis
Our group modeling will now have both fixed and random effects. Suppose Y_i = X_ib_i + ε_i  is the model of our neuroimaging data of N-subjects, each is acquired with n time points. The within-subject parameter estimates and residual are b_i and ε_i respectively. In the random-effects model, we further treat the individual subject as coming from a set of the population with its own between-subject variability. Thus, b_i can be further decomposed into b_G + εGi. the group parameter estimates and group error of which the subject belongs. Under the null hypothesis, H₀: c^Tb_G = 0, we can then compute:

which obeys a t-distribution with # dof = N – 1. Note that the numerator is the same as the fixed effects model, but the denominator contains more terms. The denominator in random effects is comparatively larger than that in fixed-effects model. So, if you use fixed effects for group analysis, you tend to get unnecessary activation, that is, a false positive.

Notes on Linear Mixed Model
The statistics mentioned in this post are also known by linear mixed model (LMM). With mixed modeling, the model will fit the average intercept and slope as a fixed effect, and each subject is assigned a different intercept/slope. Hence, LMM is an extension of GLM, but this is unlike a standard linear regression model which has only fixed effects. LMM approach is powerful to model data where there is a hierarchical or nested structure, or we deal with multilevel modeling. LMM can also be used in the case of repeated-measures design, or cases when there are some missing data.

Let's take another classic example of LMM: modeling distribution of performance of college students in New York City. Variation in the performance can primarily be due to a random variation at the individual (each student) level. But suppose the city has many colleges. Each college can also be seen as a contributing factor to the performance of each of the individuals at that school, say, due to the good teachers and attitude of the school. Hence those observations cannot be treated as fully independent of each other, but dependent on a higher level grouping called (which college?) - breaking the major assumptions of more traditional linear models. This example reflects different sources of randomness which are in a hierarchy, i.e. individuals-classes-college. This resembles our fMRI dataset, as we move up from an individual subject data to a group-based inference.

LMM is also common for two explanatory variables that are not independent. For example: suppose we want to predict the pitch level as a function of two independent variables: age and gender. Following GLM notation, we can write down:

pitch ~ nationality + gender + error

Notice how age and gender are dependent as coming from the same person or subject. Taking into consideration fixed and random effects with LMM, we rewrite the equation to be:

pitch ~ nationality + gender + (1 | subject) + error

Analyzing fMRI data with GLM

What is GLM?
GLM (general linear model in statistics) forms the core of the statistical pipeline in fMRI data analyses. Suppose we want to localize a brain activity associated with moving the left arm. And suppose we utilize a simple blocked design consisting of two alternate blocks of rest and move. We can predict the parts of the brain associated with moving the left arm through a simple prediction. The area in which BOLD signal changes in the same manner as our task-block can be the focus. In essence, what GLM does is multiple regression, i.e. predicting an outcome based on several predictors (regressors, a.k.a independent variables, or explanatory variables). In fMRI, the outcome Y is the voxel time series and this operation is done for each voxel in the brain.

So, what is GLM? There are three basic meanings:

    a)  It is a model to predict or estimate the actual BOLD response based on a known task design.

    b)  It is 'linear' because our prediction is formed through a linear combination of many factors. 

    c)  It is 'general' because we use assumptions to proceed with the statistical tests (t-test, F-test).

Basic statistics remain the same. The variable that we want to predict is the dependent variable. In GLM, the BOLD response of each voxel in the brain is our dependent variable, or to be more precise, it is the voxel time series. The predictors or regressors are also called EVs or explanatory variables. GLM tries to model each voxel's time series as a linear combination of many EVs. 

In other words:

    a)  GLM is a model-fitting strategy done one by one for every voxel, a mass-univariate statistics.

    b)  The contribution of each regressor is called parameter estimate (PE) or β coefficient.
    c)  As in normal regression, the difference between the fitted and actual data is called the residual.


Fundamentally, GLM is a simple linear regression, that is, a model that consists of a predictor (explanatory variable) and a coefficient (slope or β-value). The following diagram sums up what happens when you construct a GLM when you analyze your fMRI data and how it relates to the subsequent calculation of the t-value. A good source for this is the FSL Course slides.

GLM in a Diagram

Fig-1: An illustration of how General Linear Modeling of functional MRI data analysis works.

Key Steps in GLM:

• Our goal is to predict the voxel time series (y) during the task phase using our predicted BOLD response model. Voxel that best predicts the observed time series should be responsible for task performance.
• We assume we roughly know the shape of the BOLD response based on our experimental paradigm. The approach is to employ a canonical HRF or hemodynamic response function such as the double-gamma function. Convolve this with the boxcar function to obtain the model of the BOLD response that corresponds to our blocked design. Note, this function is also called the basis function because you want to express multivariate data using a certain function.
• Determine the whole set of regressors, EVs, or predictors. These include the canonical model mentioned above, motion-related variables, all other noise signals as noise regressors, etc. A collection of EVs is also known as the design matrix. We include noise regressors to take into account artefacts corrupting the actual BOLD response (hence the term: to regress them out).
• A good design matrix should not contain an EV that is highly correlated or a linear combination of another EVs. If this happens, the model is rank-deficient. Our estimates may then fail to represent the actual BOLD response. A well-designed study/task, hence, should ensure that the regressors in the GLM are orthogonal. This topic needs a separate discussion.
• In GLM, the best prediction occurs when the residual is minimized. Thus, to find the best fit of our data is to estimate β values for each voxel such that the residual is minimum.

The solution for unbiased estimates of betas can be found by using the Gauss-Markov Theorem but with some assumptions involved: 

• The residual errors ~ N(0, σ².I),  are normally distributed and independent of each other, 
• The homogeneity of variance is held, 
• Each regressor cannot be a linear combination of other regressors. 

Unfortunately, fMRI observations contain some sort of temporal autocorrelation denoted by a matrix V, making the dirty error variance Σ_e = σ².V. To ensure the data satisfy the assumption, such autocorrelation shall be handled properly. Data preparation such as prewhitening using FILM in FSL or coloring (temporal smoothing in SPM) are meant to achieve this.

Sources of autocorrelation in fMRI dataset include neural source, scanner drifts, physiological signals (respiration and cardiac pulsation), and the head movement. Not accounting for this temporal autocorrelation results in spuriously high fMRI signal at one time point that bleed over to the subsequent time points, which increases the chance of false positive results in task-based studies. The word 'pre' means the operation is added or applied to the original data. The word 'whitening' is a term that is motivated by the fact that white light contains all visible frequencies in equal amounts, i.e. making the spectrum flatter.

What are the steps? 

• First, we fit the GLM ignoring the autocorrelation, estimate the residual errors which contain the correlated voxels and have variance Σ_e. 
• Use the obtained residual to estimate the autocorrelation matrix V. FSL performs a local voxelwise estimation (instead of the global in SPM), which can be modelled quite well with an AR(1) function. Then smoothen this estimate V to correct for any biases using a Tukey Taper, hoping to downweight noisy estimates at higher lags. 
• Once V is obtained, we find a new matrix K such that KVK' = I, where I is the identity matrix. Subsequently, prewhitening is achieved by multiplying both sides of the original GLM equation with this matrix K. 
• We then repeat the GLM again to obtain the correct beta weights. The resulting beta estimates are now unbiased! 

Hypothesis Test: Contrast?

Okay, suppose the whitened GLM is completed. What's next? We want to know which voxels are highly activated due to task performance. How? If the resulting β from the GLM is big and the residual is small, highly likely the voxel is significantly activated during the task. As a quick note: in simple linear regression, β can be treated as a slope. It represents the change in the dependent variable Y resulting from a unit change in the predictor. In statistics, we assess how good the predictors are by comparing whether the PEs are significantly different from zero. Why? If the β = 0, it means the predictor does not contribute anything to the outcome Y.

We are now in the position to perform a t-test on every voxel to know which ones carry significant βs. The formula is given below, the denominator denotes the var[c^T β]. This variance depends on the residual and the critical assumption here is that the residual follows a normal distribution ~ N(0, σ²I).

What is the # dof? On each voxel, we have n time points and n is usually large. A few minutes scan can give you 500 volumes, In our t-test, the degree of freedom = n ─ p. With such a huge # dof, our t-statistics is approximately a z-distribution. Our statistical test gives a map called "Statistical Parametric Map". Because of the number of voxels in the brain, and knowing in reality that each voxel isn't that independent, the problem of multiple comparisons is challenging. Another post will discuss how this can be managed, e.g. through Random Field Theory.

The null hypothesis that a particular voxel isn't significantly related to our task is H₀: c^T β  = 0, where c represents a set of contrast. This contrast is to test the β values and serves as a linear combination of regressors. For p-number of regressors (EVs), the contrast can be written as a p × 1 matrix, each is associated with a beta. In our blocked design, the first regressor is associated with the predicted BOLD response. This allows us to put the first contrast as 1, leaving the rest of unwanted regressors as zeros, c = [1 0 0 ... 0]. This does not look like a normal contrast but as long as there is no singularity in the GLM, the hypothesis test is deemed valid.

Another example, say, if you use two different and independent tasks, e.g. visual and auditory tasks, you want to localize voxels related to your visual, auditory, or both on average. You may then set the contrast to be c₁ = [1 0 0 ... 0], then c₂ = [0 1 0 ... 0], c₃ = [1 1 0 ... 0]. If you want to localize voxels that are significantly more active in the visual than in the auditory task, then c₄ = [1 -1 0 ... 0]. The importance of assigning correct contrasts should not be ignored. Note that we can also perform an F-test, just like ANOVA, to test whether any of the contrast is significant.

Fig-2: An example of a rank-deficient GLM. The design matrix contains two task-related regressors from a blocked experiment (no noise regressors for simplicity). We can equally well use either EV1 or EV2 to explain what we see in the data. In other words, Y can be fit with any linear combination of the EVs, as long as β1 + β2 = 0.9. This yields a problem in hypothesis tests with [1 0] or [0 1] contrasts, but we can still get away with that when testing [1 1] contrast. The computation may still be successful, but your inference hereafter won't be accurate. 

For more detailed explanation on GLM, I find the following sources useful:

      [1]   A nice review article by Martin Monti (Front. Hum. Neurosci., 2011).      
      [2]   Chapter 5, in "Statistical Analysis of fMRI Data", a book by Ashby FG (2011).
      [3]   A more classic Chapter 9, in "Functional MRI: An introduction to methods" by K. Worsley (2001).
      [4]   For this post, I consulted Chapter 7, in "fMRI Techniques and Protocols", Woolrich M, et al. (2009).

Voxel vs Cluster-based Thresholding
In general, the main statistical analysis involves performing GLM at every voxel of the brain. That's why it is also called the mass univariate analysis. Once the statistical test is carried out, we get a z-map or t-map with a large #dof. Each voxel is constructed under the null hypothesis that nothing interesting happens because of our behavioral task. In statistics, the most common (frequentist) approach of a hypothesis testing is to obtain a p-value which is then compared against a certain threshold α. This α also represents the chance of a false positive which is capped at 0.05. Simply put, when we have data with 100,000 voxels, we have 5,000 voxels deemed false positive when we use α = 0.05. The error or bias in drawing a conclusion due to multiple comparison problems is also called familywise error because we are essentially repeating the same t-test for all voxels.

In statistics, the most common correction method for familywise error is Sidak-Bonferroni correction. In fMRI, however, the method is highly conservative. Two more popular methods are widely used in the neuroimaging community: the classic Random Field Theory or GRF (Worsley, K. 1995/1996), the False-discovery Rate or FDR (Benjamini & Hochberg, 1995). A more recent development based on non-parametric statistics using the permutation method is also popular. Detailed discussions about these methods are not presented here.

Post-stats thresholding is the last step in the pipeline to draw the conclusion about our neuroimaging data. This step is done with two objectives in mind. First, we want to know whether a certain voxel or a set of voxels is really active due to the task. Second, we draw that conclusion after having a proper correction for multiple comparisons. Technically, there are two different ways of thresholding:

(1) Voxel-wise thresholding
We can correct for multiple comparisons and do thresholding at every voxel by showing which part of the brain is active at a particular significance level. We reject the null hypothesis that there is no activation if t > u_v. This is called voxel-wise thresholding with threshold u_v. The merit of this method is the high specificity but it faces a serious multiple comparison problem. In earlier days when the scientific community was overly excited about fMRI, the results were reported using voxel-wise thresholding using a more stringent Bonferroni correction, say p < 0.0001, focusing on an area of interest, say frontal motor cortex or visual cortex only. Theoretically, Bonferroni correction is not suitable as it assumes adjacent points are not correlated.

(2) Cluster-based thresholding
Cluster-based thresholding is the more accepted way in the community currently. If a brain region is activated, we expect that not only one, but the surrounding set of voxels gets activated too. Instead of dealing with an individual voxel, we look at a set of voxels called a cluster. Cluster-based thresholding is done in two steps. We first create a binary image (pass/fail) of voxels passing the cluster-forming threshold u_k. Following this, we find another thresholding parameter that is based on either cluster size or peak voxel. For example, we set the cluster parameter with a size greater than a certain value p_k. Here, the sensitivity is typically better but worse specificity. Correction for multiple comparisons is more complicated because it involves a series of contagious voxels. This will be discussed in more detailed in my other post (Random Field Theory, RFT).

Case study: My own data
I'd like to write a simple example of how one can perform statistical analysis (GLM) in FSL. A participant performs a very simple motor task inside the scanner with the right-dominant arm. This is a blocked design with two conditions (rest | move | rest | move |.....) lasting for about 7 minutes (TR=1.69 sec). The scan recorded 250 volumes, meaning, the data would have 250 time points.

The following graphs illustrate what we expect from the GLM analysis to identify voxels that are significantly associated with the task. The red line in the top panel is our observed data, the time series of a particular voxel at a certain coordinate. Ideally, the best fit would see a smooth HRF curve that spans over the 'move' block. The purple plot of a full model fit of the actual time series given by the GLM analysis, that is, given all of the regressors, whether any activation in this voxel can be attributed to any of the task conditions. Lastly, the green line represents only the contrast of interest, and is usually only meaningful when having 3-4 task conditions or taking a simple main effect. Another way of looking at full model vs partial mode fit is this. The former shows what happens when all EVs are used in the fit, whereas the latter shows what happens when only some EVs are used (those involved in the contrast).

The resulting z-map shown at the bottom panel is thresholded at Z = 3.0 with a corrected cluster p-value < 0.05. This simply means that among the voxels that satisfy Z > 3.0, we apply the cluster-based thresholding with a correction for multiple comparisons. In FSL, GRF theory is used as the default correction method.

Fig-3: Outputs of the task-based analysis using GLM. The color map overlayed on the brain image is sometimes called the statistical parameter mapping or SPM (after Friston, et al). Voxelwise temporal autocorrelation must be taken care of by default (prewhitening stage) to ensure that the GLM outputs are valid.