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.