Linear models: Getting the syntax right!

When conducting statistical analysis and mixed-effects modelling, I often either got confused or forgot how to write the syntax properly. I thought this entry will be useful. 

Simple linear models

In simple linear models in R using the lm library, we can write the formula as DV ~ IV (DV is a dependent variable). The command employs a few arithmetic symbols. A '+ sign' indicates more than one main effect or predictor (independent variable, IV) with no interaction. A '* sign' provides a short form of main effects with interaction. Number '1' with a '+ sign' refers to an intercept. There is no need to explicitly write this, as it is always implied and estimated by default. Sometimes, you want to omit this by simply writing '0'.

When we want to include fixed/random factors, we then use, e.g. the lmer library, and the syntax is slightly different.

DV ~ 1 + IV1 * IV2        DV ~ IV | grouping                 

1. The dependent variable (DV) is the response variable to be predicted. The independent variable (IV) is the one whose effects we will assess. By default, IV is treated as a fixed factor.
2. A vertical bar denotes a grouping factor or random factor. It separates expressions for design matrices from the grouping factor. For example: to fit a predictor for each random factor, you can write 1 + A|S, which means the same thing as (A|S).
3. A '/ sign' indicates nesting. So (school/class) means classes are nested within school.
4. Fixed factors can be included without any grouping. You can have additional random factors without any fixed factor (an intercept-only model). In each random/fixed factor, ask yourself whether you allow the intercept to change, the slope, or both!

Hierarchical Linear Modeling

Begin with the first and simplest model which is the intercept-only model. It has no predictor or independent variable (so, no slope!). Suppose we have one random factor, usually the participant factor. The model is also called the unconditional model.

Y ~ (1|S)            ; intercept-only model 

The next model is the random intercept model. We assign a new variable or predictor X as Level 1, which acts as a fixed factor. Here, the intercepts for different subjects will vary but not the slopes. Note that the intercept for X can be omitted as the function understands it to be 1+X. 

Y ~ X + (1|S)        ; random intercept model

Now, we can add complexity by allowing different subjects to have different intercepts and slopes. This is the random intercept + slope model. There are two possibilities: the variations of intercept and slope can be independent or correlated!

Y ~ X + (1|S) + (0 + X|S)    ; independent
Y ~ X + (X|S)                ; correlated

 Source: here

What is Hierarchical Linear Modeling?

This is going to be a very brief conceptual discussion on hierarchical linear modelling (HLM), also known as multi-level modelling in social sciences. People sometimes call this mixed-effects modelling or mixed modelling, but be careful that not all mixed-effects models have a hierarchy. 

Note: Mixed-effects models by definition are models which have both random and fixed effects (see below).

(1)  When do we use it? 

In social sciences, the data types are usually hierarchical in nature, i.e. they have a certain hierarchy that is based on grouping or nesting. We use HLM precisely in cases where observations are nested or clustered in some ways at a different level. It has many applications in cross-sectional and longitudinal research (repeated-measures design). Research studies create data hierarchies through the way the data are sampled. For example:

• Students (S) are nested within a particular Class, which again is nested in a School.
• Time points are nested within participants in a longitudinal study.

Fig-1: Pictorial representation of (a) two-level hierarchy and (b) three-level hierarchy (Source: SAGE Ency. 2018)

Suppose we wanted to examine A-level exam scores this year to predict the university entrance test in the following semester. So entrance exam is the outcome variable [DV] and the A-level score is the predictor [IV]. One way is to randomly sample students from the whole population such that each individual is guaranteed to be independent. However, what if the students recruited only belong to a particular high school? In HLM, the students are treated as the first level or the base, and the school as the second level (2-level hierarchy). In other words, we performed sampling primarily from four different local schools only. To advance further, we have to introduce 'Class' in between both levels (now, 3-level hierarchy), since each school surely has different classes. 

In HLM, independence can no longer be freely assumed. Why? E.g. think of how students of Teacher John of School ABC have a different profile, intelligence, general ability than those attending the class of Teacher Mary from School XYZ. In such a nested condition, assumptions fulfilled for the OLS will be violated: observations are no longer independent of each other, and residual errors are correlated within clusters. The error variance will be different for different clusters. As a result, the standard errors of the regression coefficients will generally be underestimated and the significance level will be wrong, leading to misinterpretation. One way is to bring the students' data of Level 1 up to the next level through some averaging, but this is a poor practice. HM treats levels as something to take into account. It is basically an extension of the ordinary least square (OLS) of the linear regression.   

(2)  How to construct the model? 

For a given A-level exam score x_i in a simple linear regression, we estimate the intercept b₀ and slope b₁. The term e_i is the residual. The equation states that for a unit increase in the A-level exam score, the score on the entrance test y_i increases by the value of b₁. But wait! The outcome dataset (entrance test scores y_i) is clustered by the school. If this is not taken into account, then the data from all students is treated as unique observations.  

Suppose we adopt a 2-level HLM comprising Students and Schools, i.e. the students are grouped or nested within a school. We introduce a group-level variation by changing the suffices, such that y_ij is the score on the entrance test for student i in school j, and u_j is the group-level residual for that j-th school. This u_j is a new random variable, referred to as the Level-2 or group residuals, and is assumed to follow a normal distribution, ~ N(0, σ²). The new HLM equation has the following terms:

• We see b₀ as the grand mean of the university entrance test y.  
• Level 1 contains the basic form of a linear regression, y_ij = b_0j + b₁x_ij + e_ij  ; in turns, Level 2 consists of b_0j = b₀ + u_j.
• The mean of y for the school group j is b₀ + u_j, where u_j is the "school effect", i.e. the difference between the mean of the school group j and the grand mean. With this, each school can have different mean scores on the entrance test. 
• The individual-level residual e_ij is the difference between the value of y for a student i and the individuals group mean b₀ + u_j. This reflects differences in students’ individual test scores from their respective school (or group) mean.

What can this improved model tell us? Some schools will have means that are higher than the grand mean, suggesting that their students perform better on the A-level exam on average, and some schools will have lower cluster mean values. 

(3)  Defining fixed and random terms 

The core of hierarchical models is the assignment of fixed effects and random effects, two terms that already appeared in another blog post on group-level fMRI analysis. You have to define each term of a mixed-effects model to be either a fixed or random effect.

• A fixed effect is a parameter that is fixed across all groups and does not vary in the model. It is the term of interest in experimental manipulation. Estimating a fixed-effect of a term is like estimating a regression slope b. Differences in the slope If all terms in the model above are fixed, it behaves as the usual linear regression.
• In contrast, a random effect allows each group to have a different estimate. In HLM, random effects represent a higher level variable under which data points are grouped. This implies that random effects must be categorical (but cannot be continuous!). For example, the residual errors (u_j and e_ij) are considered random effects. Estimating a random effect is like looking for an effect in our data to come from a large group of normally distributed datasets. 

The most basic model has only an intercept without any predictors x_ij, which is also called the intercept-only model. Another name for it is the Unconditional Model or Null Model. Usually, the model has only participants e_ij as a random error.

The next basic model is the one shown in the equation above. It allows for different groups or schools j to have a varying intercept or mean b_0j = b₀ + u_j. Thus, the model is also known as the random intercept model. 

(4)  The more complete model 

A random intercept model shown above assumes that the relationship between the university test and the A-level exam score (the predictor) is the same for each group. This assumption can be relaxed by allowing for different slopes for the predictor in each group, making it a random slope model as shown in the right panel of Fig-2 below. Each slope will be estimated separately for each group. In this new equation, there is a new term u_1j x_ij , where:

• The intercept for school group j is now b₀ + u_0j. 
• The slope for school j is now b₁ + u_ij , where b₁ is the average slope across groups.
• Subscript “0” differentiates the random effect for the intercept u_0j from the random effect for the slope u_1j. Both random effects are assumed to be normally distributed. 

  Fig-2: Hierarchical models [Source: Univ Bristol]

The variance of the intercept and slope are assumed to be correlated; the covariance between the intercept and slope is estimated as part of the random slopes model. The random slopes model is also commonly known as the random coefficient model or a growth curve model when using repeated measures or longitudinal analysis. 

(5)  Final notes 

How is the implementation in practice? Apart from defining fixed/random terms, we have to know which one nested under which variable.

• When we deal with a dataset we have to begin with the simplest model, which is the intercept-only model.
• Then we add additional terms, the simplest being the random intercept model, i.e. each group has its own group mean. 
• If we have more than one predictor or IV, define the fixed/random effect. A mixture of fixed and random effect slopes is possible. 
• The most complete model is the random slope model, where intercepts and slopes are treated as random effects. Often in a longitudinal dataset, we can add the quadratic term "Time" (Level-1).

So how do we know which model is best? We can use the usual performance metrics such as Bayesian Information Criterion (BIC) or Akaike Information Criterion (AIC). The two most common libraries for HLM analysis in R are lme4 and nlme. In the commands, you have to indicate which ones are the fixed and random effects, and which variable is nested within what. The outputs of the HLM in R is quite similar to what we expect from OLS in linear regression. 

The random slope model often mimics reality, but this requires a large sample size. Parameter estimation in HLM is achieved using either the Maximum Likelihood (ML) or Restricted Maximum Likelihood (ReML). ReML works better if the sample size is small.

Some good references:

[1]  https://rpsychologist.com/r-guide-longitudinal-lme-lmer

[2]  https://www.lcampanelli.org/mixed-effects-modeling-lme4/

[3]  https://m-clark.github.io/mixed-models-with-R/introduction.html

Basic Math for Principal Component Analysis

To begin, Principle Component Analysis (PCA) is a statistical tool for exploratory analysis that is famous for reducing the dimensionality of the data. It is a technique that involves vector space transformation to represent data in another form. It is useful for feature extraction and feature elimination, removing redundant features. The outcome of a PCA operation is a set of orthogonal variables, called the principal components. To learn PCA from a technical angle, I recommend the work by J Shlens (ref.[1]).

1. Intro to Matrix Operations
Before going into the math, there are a few technical terms in matrix algebra we have to be familiar with. First, eigenvector and eigenvalue. An eigenvector v is of a linear transformation (or simply a square matrix A) is a non-zero vector that changes only by a scalar factor (i.e. by the eigenvalue, λ) when that transformation is applied to it. We can say, [matrix].[eigenvector] = [eigenvalue].[eigenvector]; or we can write: Av = λv.  If the eigenvalue is negative, the direction is reversed. Refer to the figure on the right for a geometric description.

An eigenvector, with a real nonzero eigenvalue, points in a direction that is stretched by the transformation and the eigenvalue is the factor by which it is stretched.

Another term is about an orthogonal matrix. Let a square matrix A is invertible (it can be inversed) or non-singular.

A square matrix is orthogonal with real entries if, and only if, its transpose (A^T) is equal to its inverse (A^-1), that is AA^T yields an identity matrix I.  The square matrix is also normal if it satisfies A^TA = AA^T.

Importantly, if the matrix is orthogonal then it is necessarily invertible, unitary, and normal. The determinant of an orthogonal matrix is either +1 or −1. Next, we learn about another special matrix called a diagonal matrix, i.e. a square matrix with non-zero diagonal elements and all-zero remaining (non-diagonal elements). When a square matrix A is diagonalizable, there exists an invertible matrix B (a matrix that can have an inverse) such that B^-1AB produces a diagonal matrix.

Matrix S is symmetric only if it has a characteristic of S^T = S. Such matrix is also orthogonally (orthonormally) diagonalisable.

An example of a square matrix related to PCA is called the covariance matrix. More specifically, in this matrix we have the variances on the diagonal entries of this matrix, and the covariances on the off-diagonal entries. Take a look below. Here, the red box denotes the variance of variable height, whereas the green box shows the covariance between the score and age.

Suppose we have a huge dataset which has too much information. Can we re-express the original dataset optimally in a new form? Now the word 'optimality' means we do not want correlated or redundant features. After a series of matrix operations, the new or transformed dataset will have a new coordinate system. The whole idea of PCA is to force as much of the variation as possible into fewer dimensions, you can throw away the rest without losing much information. Let's relate PCA with the previous ideas on matrices:

In PCA, each eigenvector is a unit vector pointing in the direction of a new coordinate axis. The axis with the highest eigenvalue is said to be the axis that explains the largest variation, meaning an eigenvalue represents variance. The directions with the largest variances are the most “important” or principal. 

2. Technical Definition

Suppose we have a dataset of m × n matrix X, where the n columns are the number of observations and the m rows are the variables. We wish to linearly transform this matrix X into another matrix Y, also of dimension m × n, so that for some m × m transformation matrix P, we will obtain Y = PX.

So, PCA attempts to express correlated m variables into uncorrelated new variables through some sort of transformation. In matrix algebra, all possible correlations of any variable-pairs can be captured using the covariance matrix C_X. How do we decouple the variables? It turns out, we can get this by maximizing the variances or diagonal elements, and minimizing the co-variances or off-diagonal entries of the covariance matrix. So our goal is to obtain a new data Y such that the covariance matrix of the transformed data C_Y is strictly diagonal. To put it in another way, if we can find the transformation matrix P in such a way that C_Y is diagonal, then our objective is met.

The solution is based on an important property of matrix algebra on C_X called the eigen-decomposition. The effort to make C_Y a diagonal matrix is essentially to find P = E^T. How to get this new matrix? It turns out that the orthonormal matrix E can be found by diagonalizing the covariance matrix of the original dataset, C_X = EDE^T , where E is the matrix of eigenvectors which applies the transformation on X; and matrix D has diagonal elements containing the eigenvalues in descending order. In short, the principal components are formed by:
- The “new directions” of our data as indicated by the eigenvectors;
- The “magnitude”, or importance of each direction, as indicated by the eigenvalues.

Often, the normalized eigenvectors by the variance are also called the coefficients or loadings of the principal components. Remember, these coefficients define how the new principal component is formed as a linear combination of the original variables. Each column in the eigenvector matrix contains coefficients for each component in the descending order of the variance.

Dimensionality reduction works as follows. Suppose our original data is an m × n matrix X with m dimension and n data points.

If we want to transform the points to k dimensions (where, k < m), then select the first k eigenvectors of the matrix C_X sorted decreasingly according to the eigenvalues, and form a matrix with them, then use them as k × m matrix P. The resulting matrix is Y with size (k × m)(m × n) = k × n.

3. Closing remarks

PCA does not simply do a rotation. Normalization is always involved prior to PCA, which can be done by mean-centering the data and divide them with the standard deviation. This is important if our data contain variables of different scales, e.g. kilogram and pounds.

PCA and ICA are two slightly different things. While PCA expresses the data into a set of orthogonal variables, ICA simply represents data into an independent set of variables (not necessarily orthogonal). In ICA, the order of importance is irrelevant. Prior to ICA, data whitening is performed and PCA is one method to whiten the data. The figure on the right is common on the internet. PCA is predominantly used as a dimensionality reduction technique, for example, in computer vision, facial recognition, and image compression. It is also used for finding patterns in a high dimensional dataset in the field of finance, data mining, biomechanics, bioinformatics, etc.

References:
(1) "A tutorial on Principal Components Analysis" by J. Shlens (2014). Very nice technical tutorial!
(2) https://towardsdatascience.com/a-step-by-step-explanation-of-principal-component-analysis-b836fb9c97e2
(3) https://stats.stackexchange.com/questions/612/is-pca-followed-by-a-rotation-such-as-varimax-still-pca

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.

Experimental Design - Statistical Tests

Basic inferential statistics.   I decided to separate this from the earlier entry. We shouldn't forget that what we observe is based on a sample, not a population. Inferential statistics has its goal in making a conclusion about a population from just this sample. In doing so, there bound to be some errors as our measure is just an estimate. Suppose we have a population of interest and that we draw all possible samples of size n. And suppose we compute a statistic e.g., a mean, proportion, standard deviation for each sample. The probability distribution of this statistic for all possible samples is called a sampling distribution; the standard deviation of which is called standard error which roughly tells us the spread around the mean of the sampling distribution. Mathematically, statisticians have to make some assumptions in order to model the problem accurately. This will be discussed later on.

How can we infer about a population parameter by using our sample? We use the concept of the sampling distribution (of means) to see where the true population mean lies w.r.t our sample mean. Having obtained the sample mean and standard deviation, we can find a confidence interval within which the true population parameter lies. This interval can be found with a certain confidence level, say 95% or 99%, with the help of a transformation. This data transformation converts our sample into a known sampling distribution, e.g. using a Student t-distribution, by first converting our scores into a t-value. The conclusion of our experiment is obtained by making use of a known sampling distribution and putting that distribution in a certain hypothesis test. The type of statistical test used (e.g. one-sample t-test, ANOVA, or chi-square test) depends on the nature of our variables.

In constructing a hypothesis, there are by definition the null and alternative hypothesis, indicated as H₀ and H₁. In the language of experimental design, the H₀ states that the IV or treatment (marijuana) has no effect on DV (exam score). Conversely, the H₁ states that the treatment has an effect on DV. To do, we first create a sampling distribution assuming the null, H₀ is true. We decide an alpha level α or significant threshold to represent a concept of "very unlikely-ness" or chance that we observe a treatment effect when H₀ is true. We then compare the likelihood of getting our sample statistics p to the alpha level that is usually fixed at 0.05 or 0.01. The final outcome of the test is whether we disapprove or retain the null hypothesis. For example,

• The probability of our sample statistics of 0.01 means there is only 1 out of 100 chance you find no treatment effect in your data when the null H₀ is true  ..... We reject H₀.
• The probability of our sample statistics of 0.80 means there is 8 out of 10 chance you find no treatment effect in your data when the null H₀ is true  ..... We retain H₀.
• Rejecting the null while in reality it's true is a Type I error with a chance α. Accepting the null while in reality it's false is a Type II error with a chance β. Note that these errors are not measurement errors made by mistakes but rather a measure of chance or probability.

Fig 1: Summary of hypothesis test with its associated Type I/II errors. 

Behavioral science is unique in the sense that the test must first assume the null hypothesis to be correct. Also, the test result simply makes a case without even saying which hypothesis, H₀ or H₁ is indeed true or factual. Even if H₀ is indeed true, we don't know how much the difference is with the untreated sample group. For such reason, the effect size d comes into play. The effect size is defined as the difference in means between the experimental group and control group, divided by the standard deviation of the control group.

Lastly, the ability to correctly reject the H₀ or find something significant in your data while H₀ is indeed wrong is called power, and it is simply 1 – β. A powerful design is able to detect a small effect size < 0.20 unit. Power increases or β is reduced when the alpha level increases, the mean difference between two conditions under the null hypothesis increases, and sample size n increases. The concept of power is crucial because it determines the replicability of our study. According to Cohen (1966), most of the behavioral studies have low power, i.e. the chance of replicating the same experiment yielding the same result is 50-50.

A technique called power analysis is used to determine the appropriate sample size for the study. Generally, power = 0.80 is an acceptable value. Tools such as "pwr" in R is able to help us in power analysis, e.g. find the minimum sample size to achieve certain p-value and power of 0.80.

Okay, we move on to talk about assumptions in statistics, a crucial component to make our problems mathematically solvable. Most statistical procedures such as t-test, linear regression, and ANOVA are based on assumptions. For example, a linear regression is valid when the DV and IVs are linearly related to begin with (for categorical variables such as "Yes/No", one uses logistic regression instead). Regression has other important assumptions in general, e.g. the residuals are independent, normally distributed, and their variance is homogenous. The type of test that requires certain assumptions to be met is called parametric statistics. Violating these assumptions causes incorrect inference or decision. The assumptions for parametric tests are presented below.

[Reference #1]: Among some relevant textbooks, I chose the materials from the classic JL Myers' "Fundamentals of Experimental Design" (1979).

The normality assumption.   As mentioned, statisticians make assumptions on their methodology to simplify the maths. Both t-test and ANOVA have to satisfy the following main assumptions:
(a) Samples are independent of each other; the score of Tom is independent of Ni;
(b) The normality assumptions (for ANOVA, it means the distribution within groups);
(c) The homogeneity of variance.
The normality assumption means that the sampling distribution comes from a population that is normally distributed or Gaussian. This is hard to know as we don’t have the luxury of knowing all possible samples related to our population of interest. Out of faith, we assume that if the data we take is normally distributed then the sampling distribution will also be normally distributed. Every stats guru also tells us about the central limit theorem that is, the distribution of the sample means approaches a normal distribution as the sample size N increases, regardless of the shape of the population of interest. As a rule of thumb, N > 30 is generally accepted for the assumption to be valid. The average of our sample means is itself the population mean, and the standard deviation is the standard error.

Fig 2: Comparison between normal data and its QQ-plot (left) and "not-really" normal data (right)

This can be tricky if N is as few as 15 or 20. Now, how do we test for normality?
(1) Visual comparison. Normal data has a bell-shaped distribution. It has a specific kurtosis and skewness characteristic. It is recommended to transform the kurtosis and skew values to z-scores. Another useful visual test is through using Q-Q plot that compares theoretical values against the actual or observed values. An ideal case will be a straight line with 45-deg w.r.t the horizontal. The following figure shows an example how the data that is less normal has the tendency to have Q-Q plot that is not straight, but rather curved away from the perfect straight line. Most of the times, the presence of outliers distort the normality, the skewness or kurtosis of your data. It is quite common to treat outliers by removing the particular data point or by transforming the data.

(2) The Shapiro-Wilk test. In R, you can simply write shapiro.test(variable name); or in Matlab as swtest(x, alpha). This test is super cool as it directly assesses whether a certain distribution is significantly different from the normal distribution with the same mean and standard deviation. If the p is significant (p < 0.01, say), then the data is different from the normal distribution.

Homogeneity of variance.   Suppose you have three groups in your experiment. Homogeneity of variance says that the variability of the data in all groups should be equal, because strictly speaking, they come from the same population. The same notion applies, for example, in the case of repeated measures where you test the same group in two different experimental manipulations. If the homogeneity principle holds true, the variance of the data in either case remains rightfully the same. In most cases, we usually ignore this principle and immediately perform t-test or ANOVA, etc. Refer to the following illustration.

There are two tests that you can use to assess this: Levine’s Test and Hartley’s F-test (good for large sample). Both methods test the null hypothesis that the variance of the two datasets or groups are equal. The Hartley's F-test is also called the equality of variance test, or sometimes variance ratio test since you're actually taking the ratio of. The computed F-statistics will be compared against the standard F-table according to their respective degree of freedom. In general, people have shown that ANOVA is robust against a mild non-normality and heterogenous variance (only for equal sample size among groups!).

Fig 3: Comparison between dataset with homogenous (left) and non-homogeneous variance (right)

Non-parametric statistics.    Non-parametric statistics means the methods don't rely on any assumption discussed above! In fact, the information on the distribution is obtained directly from the dataset itself. This is the most fundamental idea of non-parametric statistics! Bootstrapping, for example, is one popular non-parametric test that relies on "random sampling with replacement". This test can be used to assess stability of your statistics and it is pretty straightforward. Another test called permutation test builds, rather than assumes, the sampling distribution by 'shuffling' the observed data. Unlike bootstrapping, we use only the existing data set and the shuffles are without replacement.

I took this from Matlab website. Suppose you have a set of LSAT scores and  GPA. You can easily plot the data and compute the correlation. It is found that there is a positive relationship between LSAT and GPA, r = +0.78. Now although it may seem large, it is unknown whether this observation is statistically significant or just happens by chance. Using the bootstrp function you can resample the two vectors as many times as you like and consider the variation in the resulting correlation coefficients. If the finding is true, you would expect that bootstrapping process yields more or less the same correlation value. Refer to the following histogram. Bootstrapping gives a distribution of correlation coefficients that is centered near r = +0.78. Using permutation strategy, we can also show how the correlation found is not by chance. We permute the LSAT-GPA pairs for 1000 times, each time is followed by computing the correlation. The correlation values of the permuted data sets should collapse since the relationship has been "broken" following shuffling.

Fig 4: Distribution of correlation coefficients with bootstrapping and the location of the original r = +0.78 (red line).

During the period when I first learnt statistics, we didn't care much about such assumptions and treated them as dogma. Now, I'm learning to accept the presence of non-parametric statistical tests used when the assumptions are not met, that is, to replace the usual t-test etc. The principle behind those tests is assigning ranking. High scores will be represented by large ranks, low scores by small ranks. The analysis is then performed on the ranks rather than the original raw data. In other words, non-parametric tests are useful when dealing with ordinal data. Example of ordinal dataset includes: the satisfaction rating from 1 - 5, the order of athletes based on speed, and the school ranking.

When we are interested in comparing two independent datasets, we need to use the Mann-Whitney test or the Wilcoxon's rank-sum test. This is the counterpart of the independent t-test. In R, you can do this by:  wilcox.test(x,y,____), some parameters inside the brackets have to be defined. In Matlab, use the following instead ranksum(x,y). When we deal with repeated measures or paired sample, then we have to use the Wilcoxon signed-rank test. For example: testing the performance before and after drug consumption. You can use the same command as in (1) in R but by defining "paired = TRUE". In Matlab, however, the command will be different now, signrank(x,y). If we are interested in testing more independent groups, we have to do non-parametric version of the regular ANOVA, e.g. Kruskal-Wallis test and Friedman's ANOVA for repeated-measure designs. Other non-parametric analysis such as the Spearman's correlation coefficient (the counterpart of Pearson's r) and Kendall tau coefficient.



Final note: ANOVA can be understood from the perspective of Linear Modeling and in itself is a very huge topic, covering multiple regression, mixed-effects or hierarchical modeling, and generalized linear models. The mathematical discussions on this theme will be beyond this blog post.....

[Reference #2]: Andy Field et al, "Discovering Statistics using R" (2012).

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.