Linear Model I: Simple Regression

9.07 Introduction to Statistics for Brain and Cognitive Sciences Emery N. Brown Lecture 13 Linear Model I: Simple Regression I. Objectives Understand the simple linear regression model and its assumptions. Understand the relation between fitting the simple linear regression model by maximum likelihood and the method-of-least squares. Understand the properties of the parameter estimates, how we compute confidence intervals for them and how we test hypotheses about them. Understand how we assess model goodness-of-fit. Understand the Pythagorean formulation of regression analysis. The linear model is the most used tool in modern statistical analysis. It also provides an important conceptual framework for constructing more complex models. For example, in neuroscience these methods are at the heart of most functional neuroimaging data analysis paradigms. In the next three lectures we will present the linear model by studying simple regression, multiple regression and analysis of variance (ANOVA). II. Simple Linear Regression Model A. Motivation for the Simple Linear Regression Model To motivate the simple linear regression problem, we consider the following examples from the neuroscience literature. Figure 14.1. Relation between Conduction Velocity and Axon Diameter. Replotted from Hursh (1939). page 2: 9.07 Lecture 14: Simple Linear Regression Example 14.1. The Relation between Axon Conduction Velocity and Axon Diameter. Hursh (1939) presented data from an adult cat relating neuron conduction velocities to the diameters of axons. Hursh reported 67 velocity and diameter measurement pairs. He measured maximal velocity among fibers in several nerve bundles and then measured the diameter of the largest fiber in the bundle. The data are replotted in Figure 9.1. Diameter is reported in microns and velocity is reported in meters per second. We see that there is a strong positive relation in that as the axon diameter increases, so does the conduction velocity. Possible questions we may wish to answer include: How strong is the relation between conduction velocity and axon diameter? Can we quantify how conduction velocity changes as a function of diameter? Figure 14.2. Time-Series Plot of first 500 observations of the MEG sensor background noise measurements. Example 3.2 (continued). In this example, we have until now only considered the distributional properties of the MEG measurements. In so doing, we have treated these observations as if they were independent. If we plot the first 500 observations of the time-series (~ 833 msec), we see that the measurements do not appear to be independent (Figure 14.2). Indeed, there seems to be almost an oscillatory pattern in the time-series. Furthermore, if we plot xt versus xt−1 , we see that there is a strong linear relation between adjacent measurements (Fig. 14.3). When xt −1 is large, xt is also large and when xt −1 is small, xt also tends to be small. In fact, the correlation coefficient between xt and xt−1 , which we define below, is 0.65. Is this relation really there and if so, could this represent a systematic distortion in the local magnetic field? page 3: 9.07 Lecture 14: Simple Linear Regression Figure 14.3. Plot of MEG background noise values xt versus xt −1. Courtesy of Science. Used with permission. Figure 14.4. Plot of Change in Neural Activity during Learning Experiments versus Learning Trial (Wirth et al. 2003). Example 14.2. Defining Neural Correlates of Behavioral Learning. Wirth et al. (2003) studied the relation between a monkey’s performance learning a location-scene association task and changes in neural activity in the animal’s hippocampus (Figure 14.4). They were interested in the question of how the timing of changes in the animal’s neural activity during the task related to the trial in the experiment when the animal learned the task. Change in neural activity was defined as a significant increase in neural activity above baseline. Learning was defined by performance better than chance with a high degree of certainty. The 45 degree line in Figure 14.4 is the line which would mean that the neural change and the behavioral change occurred page 4: 9.07 Lecture 14: Simple Linear Regression on the same trial. Is this the most accurate description of the relation between neural change and learning for these data? B. Simple Linear Regression Model Assumptions To formulate a statistical framework to study these problems, we assume we have data consisting of pairs of observations that we denote as ( x1 , y1 ), ..., ( xn , yn ). For example, in Example 14.1, the variable yi is axon velocity and the variable xi is axon diameter. Let us assume that there is a linear relation between xi and yi and write it as yi = α + β xi . (14.1) To make the relation in Eq. 14.1 into a statistical model, we will make the following assumptions: i) E [ yi | xi ] = α + β xi for i = 1,. ... , n. ii) The xi 's are fixed non-random measurements termed covariates, regressors or carriers. iii) The yi 's are independent Gaussian random variables with mean α + β xi and variance σ 2 . Equation 14.1 and these three assumptions are often summarized as yi = α + β xi + ε i (14.2) where the ε i 's are independent Gaussian random variables with mean 0 and variance σ 2 . Equation 14.2 defines a simple linear regression model. It is simple because only one variable xi is being used to describe or predict yi and it is linear because the relation between yi and xi is assumed to be linear. C. Model Parameter Estimation Our objective is to estimate the parameters α , β and σ 2 . Because yi is assumed to have a Gaussian distribution conditional on xi , a logical approach is to use maximum likelihood estimation. For these data the joint probability density (likelihood) is ( ) ( L α , β , σ 2 | y, x = f y | α , β , σ 2 , x ) n = ∏ f ( y | α + β x ,σ ) i 2 i (14.3) i=1 n ⎧⎪ 1 ⎛ 1 ⎞2 =⎜ exp ⎨− 2 ⎟ ⎝ 2πσ ⎠ ⎪⎩ 2 n ∑ ( yi − α − β xi )2 ⎫⎪ i=1 σ2 ⎬ ⎪⎭ where y = ( y1, ..., yn ) and x = ( x1 ,..., xn ) . The log likelihood is ( ) ( ) ∑( n 1 logf y | x, α , β , σ 2 = − log 2πσ 2 − 2 2 n i=1 yi − α − β xi ) σ2 2 . (14.4) page 5: 9.07 Lecture 14: Simple Linear Regression To compute the maximum likelihood estimates of α , β and σ 2 we differentiate Eq. 14.4 with respect to α and β and σ 2 ( ∂ log f y | x, α , β , σ 2 ( ∂α )= n ∑ ( yi − α − β xi ) σ2 i=1 ∂ log f y | x, α , β , σ 2 ∂β )= n ∑ ( yi − α − β xi ) xi σ2 i=1 ( ∂logf y | x, α , β , σ 2 ∂σ 2 )=− n 2σ 2 1 2 + n ∑ ( yi − α − β xi )2 . i=1 (σ 2 )2 (14.5) (14.6) (14.7) Setting the derivatives equal to zero in Eqs. 14.5 and 14.6 yields the normal equations αn + β n n ∑x = ∑ y i i=1 α n n ∑x + β∑ i i=1 (14.8) i i=1 xi2 n = i=1 ∑x y . i i (14.9) i=1 In matrix form Eqs. 14.8 and 14.9 become ⎡ ⎢ n ⎢ ⎢ n ⎢ x i ⎢ ⎣ i=1 ∑ ⎤ ⎡ n ⎤ xi ⎥ yi ⎥ ⎢ ⎥ ⎡α ⎤ ⎢ i=1 ⎥ i=1 . =⎢ ⎥ ⎥ ⎢ ⎥ n n β ⎢ ⎥ 2⎥⎣ ⎦ xi ⎥ ⎢ xi yi ⎥ i=1 ⎦ ⎣ i=1 ⎦ n ∑ ∑ ∑ ∑ (14.10) Solving for α and β yields ⎡ ⎢ n ⎡α ⎤ ⎢ ⎢β ⎥ = ⎢ n ⎣ ⎦ ⎢ ⎢ xi ⎣ i=1 ∑ ⎤ xi ⎥ ⎥ i=1 ⎥ n xi2 ⎥⎥ i=1 ⎦ n ∑ ∑ −1 ⎡ n ⎤ yi ⎥ ⎢ ⎢ i=1 ⎥ ⎢ n ⎥. ⎢ xy⎥ i i⎥ ⎢ ⎣ i=1 ⎦ ∑ ∑ The solutions for β̂ and α̂ , which are the maximum likelihood estimates of β and α are respectively (14.11) page 6: 9.07 Lecture 14: Simple Linear Regression n βˆ = ∑ ⎛ ⎜ ⎜ xi yi − ⎝ ⎞⎛ n i ∑ i n i=1 n ⎛ ⎜ ⎜ 2 x −⎝ i=1 i i=1 i=1 n ⎞ n ∑ x ⎟⎠⎟ ⎜⎝⎜ ∑ y ⎟⎠⎟ ⎞ xi ⎟ ⎟ i=1 ⎠ n n ∑ i = 2 ∑ ( x − x )( y − y ) i i=1 n ∑(x − x ) (14.12) 2 i i=1 α̂ = y − β̂ x. (14.13) ŷi = αˆ + βˆ xi = y − βˆ x + βˆ xi = y + β̂ ( xi − x ) , (14.14) We may write the estimate of yi as for i = 1,..., n. If we return to Eq. 14.7 we can compute the maximum likelihood estimate of σ 2 by first substituting in β̂ and α̂ for β and α . Setting the left hand side of Eq. 14.7 to zero we obtain the maximum likelihood estimate of σ 2 σ̂ 2 = n−1 n ∑ ( y − α̂ − βˆ x ) i i 2 (14.15) i=1 which is what we would have predicted based on our analyses in Lecture 9. Remark 14.1. Estimating α and β by maximum likelihood under an assumption of Gaussian errors is equivalent to the method of least-squares. In least-squares estimation, we find the values of α and β which minimize the sum of the squared deviations of the yi 's from the estimated regression line. The method of least-squares is a form of method-of-moments estimation as we will show below. Remark 14.2. Other metrics could be used to estimate α and β , such as minimizing the sum of the absolute deviations. This would be defined as n−1 n ∑| y − α − β x |. i i i=1 Remark 14.3. The estimate α̂ and Eq. 14.14 show that every regression line goes through the point ( x , y ) . Remark 14.4. The residuals, yi − ŷi , are the components in the data which the model does not explain. They are estimates of the ε i ’s and we often write εˆi = yi − ŷi . We note that ŷi = y + β̂ ( xi − x ) , (14.16) and summing we have n n n ∑( y − ŷ ) = ∑( y − y ) − β̂ ∑( x − x ) = 0. i i=1 i i i=1 i i=1 (14.17) page 7: 9.07 Lecture 14: Simple Linear Regression Remark 14.5. There is an important Pythagorean relation between the sum of squared deviation in the data about its mean, the sum of squared deviations of the regression estimates about the mean of the data, and the sum of squared residuals. We derive it now. By Remark 14.4, the i th residual is yi − yˆi = ( yi − y ) − ( yˆi − y ). (14.18) Squaring and summing both sides of Eq. 14.18, we obtain n n ∑ ( y − ŷ ) = ∑{( y − y ) − ( ŷ − y )} i i 2 i=1 n = ∑ 2 i (14.19) i i=1 ( yi − y ) + 2 i=1 n ∑ n ∑ ( y − y )( ŷ − y ). ( ŷi − y ) − 2 2 i=1 i (14.20) i i=1 Now if we analyze the last term in Eq. 14.20 using Eq. 14.16 we find that n n ∑( y − y ) β̂ ( x − x ) = −2β̂ ∑ ( y − y )( x − x ) −2 i i i i=1 (14.21) i i=1 = −2β̂ 2 n ∑( x − x ) 2 i i=1 n ∑ ( ŷ − y ) , = −2 2 i i=1 and Eq. 14.20 becomes n n n i=1 i=1 i=1 2 2 2 ∑ ( yi − ŷi ) = ∑( yi − y ) − ∑ ( ŷi − y ) , (14.22) which yields the desired Pythagorean Relation n n n i=1 i=1 i=1 2 2 2 ∑( yi − y ) = ∑( ŷi − y ) + ∑( yi − ŷi ) . (14.23) It states that Total sum of squares (TSS) = Explained sum of squares (ESS) + Residual sum of squares (RSS) We will make extensive use of this relation in our analyses of our regression fits. page 8: 9.07 Lecture 14: Simple Linear Regression D. The Relation between Correlation and Regression Let us define the sample correlation coefficient between x and y as n ∑ ( x − x )( y − y ) i i i=1 rxy = ⎡ 2 ⎢ ( xi − x ) ⎣⎢ i=1 n ∑ n 1 ⎤2 . (14.24) ∑ ( y − y ) ⎥⎥⎦ 2 i i=1 It is the method-of-moments estimate of the theoretical correlation coefficient between x and y which is defined as Cov( X , Y ) ρ xy = 1 [Var ( X )Var (Y )] 2 , (14.25) where Cov( X , Y ) = E[( X − E( X ))(Y − E(Y ))] and where the expectation is taken with respect to f ( X , Y ), the joint distribution of X and Y. If X and Y have a bivariate Gaussian distribution, then it is easy to show that rxy is also the maximum likelihood estimate of ρ xy . It is also straight forward to show that we always have −1 < rxy < 1. Recall we showed in Lecture 5 −1 < ρxy < 1. As we suggested in Lecture 4, when we use the term correlation we will be using it in the technical sense defined by either Eqs. 14.24 or 14.25. To work out the relation between the sample correlation coefficient and the slope parameter estimate, we recall that the slope parameter estimate (Eq. 14.12) can be expressed as n ∑ ( x − x )( y − y ) i βˆ = i i=1 n ∑(x − x ) . (14.26) 2 i i=1 Hence, n ∑ ( y − y )( x − x ) i βˆ = i i=1 n ∑(x − x ) 2 i = σˆ xy σˆ x2 i=1 1 2 ⎧ n 2⎫ ⎪ ( yi − y ) ⎪ σˆ y ⎪ ⎪ = ⎨ i=1 ⎬ rxy = rxy n σˆ x ⎪ ( xi − x )2 ⎪⎪ ⎪ ⎩ i=1 ⎭ ∑ ∑ (14.27) page 9: 9.07 Lecture 14: Simple Linear Regression where σ̂ xy = n−1 n ∑( y − y )( x − x ). We see that the regression coefficient is a scaled version of the i i i=1 sample correlation coefficient or simply the sample covariance of x and y divided by the sample variance of x. In this way, we see that the least-squares estimates are method-ofmoments estimates. E. Distributions of the Parameter Estimates The parameter estimates α̂ and β̂ are functions of the data. Therefore, they are random variables and have distributions. It is straight forward to show that these distributions are Gaussian. Therefore, it suffices to specify their means and variances because this is the minimal description needed to define Gaussian random variables. It is straight forward to show that E(β̂ ) = β and E(α̂ ) = α . The variances of the regression parameter estimates are Var ( β̂ ) = σ2 n ∑(x − x ) (14.28) 2 i i=1 ⎡ ⎢ 2 ⎢1 Var (α̂ ) = σ ⎢ + n ⎢ ⎢ ⎣ ⎤ ⎥ x ⎥ ⎥. n 2⎥ ( xi − x ) ⎥ i=1 ⎦ 2 ∑ (14.29) If we estimate σ 2 by ŝ2 = n(n − 2)−1σ̂ 2 (Eq. 14.15) then following the logic we used in Lecture 8, the 100%(1− δ ) confidence intervals for the parameters based on the t − distribution with n − 2 degrees of freedom are β̂ ± tn−2,1−δ / 2 sˆ (14.30) 1 ⎫2 ⎧⎪ n 2⎪ ⎨ ( xi − x ) ⎬ ⎪⎭ ⎪⎩ i=1 ∑ ⎧ ⎪ ⎪1 α̂ ± tn−2,1−δ / 2 ⎨ + ⎪n ⎪ ⎩ 1 ⎫2 ⎪ x2 ⎪ ⎬ ŝ n 2⎪ ( xi − x ) ⎪ i=1 ⎭ ∑ (14.31) There are n − 2 degrees of freedom for this t-distribution because we estimated two parameters in this analysis. Remark 14.6. We can also invert the above confidence intervals as we discussed in Lecture 12 to test hypotheses about the regression coefficients by constructing a t − test. page 10: 9.07 Lecture 14: Simple Linear Regression n ∑( x − x ) Remark 14.7. The variance of β̂ decreases as 2 i increases. Hence, if we can design i=1 an experiment in which we can choose the value of the xk 's, we will spread them out as far as possible across the relevant range to decrease the variance of the estimated slope. To construct a confidence interval for a predicted yk value at a given value, xk , we note that the Var ( yk ) = σ2 n + ( xk − x )2 σ 2 n ∑( x − x ) (14.32) 2 i i=1 and hence, the 100%(1− δ ) confidence interval is defined by ⎡ ⎢ ⎢1 ŷk ± tn−2,1−δ / 2 ⎢ + n ⎢ ⎢ ⎣ 1 ⎤2 2 ⎥ ( xk − x ) ⎥ ⎥ ŝ. n 2⎥ ( xi − x ) ⎥ i=1 ⎦ (14.33) ∑ Example 14.1. (continued). We fit the model in Eq. 14.1 to the axon diameter and velocity data in Fig. 14.1. The estimated regression line, along with the local 95% confidence intervals (Eq. 14.33), is shown in Fig. 14.5. The estimated regression line is consistent with the strong linear relation seen in the original data in Fig. 14.1. Most of the data lie close to the estimated regression line. Parameter Estimate Standard Error t-statistic α̂ = 1.44 0.14 -2.61 43.5 β̂ = -3.77 6.07 Table 14.1. Parameter Estimate Summary The parameter estimates and their corresponding standard errors and t − statistics are shown in Table 14.1. We reject the null hypothesis that there is no linear relation since the t-statistic for β̂ is 43.5 and has a p − value<<0.05. As there are 67 observations this t − statistic is essentially a z − statistic which has a value of 43.5!! Hence, this result is highly significant. More importantly, the 95% confidence interval for β̂ is approximately [5.79, 6.35]. The estimate of β̂ has units of meters per second per micron. This means that for every one micron increase in diameter, there is approximately a 6.07 meter per second increase in velocity. We also reject the null hypothesis that there is a zero intercept for the regression line since the t − statistic for α̂ is -2.61 and has a p-value<0.05. More importantly, the 95% confidence interval for α is approximately [-5.21, -2.33]. We conclude that the intercept is most likely negative given the limits of the confidence intervals. The intercept has units of meters per second. If this were a physical model, it would mean that at a zero diameter the velocity would be negative. This does page 11: 9.07 Lecture 14: Simple Linear Regression not make physical sense and shows the importance of not using a regression model beyond the range of the data where it is estimated. The smallest diameter in our sample is approximately 2 microns and the largest is approximately 20 microns. Figure 14.5. Fit and Confidence Intervals for a Simple Linear Regression Model of the Axon Conduction Velocity as a function of Axon Diameter. F. Model Goodness-of-Fit A crucial—if not the most crucial—step in a statistical analysis is measuring goodness-of fit. That is, how well does the model agree with the data. We previously used Q-Q plots to assess the degree to which elementary probability densities described sets of data. We discuss two statistical measures and three graphical measures of goodness-of-fit. 1. F-test Given the null hypothesis H 0 : β = 0, i.e. that there is no linear relation in x and y, we can test this hypothesis explicitly using an F − test. The F − statistic with p − 1 and n − p degrees of freedom is Fp −1, n − p = −1 ( p −1) ESS (n − p) −1 RSS ( p −1) −1 = (n − p) n ∑( y − y ) −1 2 i i=1 n ∑ ( y − yˆ ) , (14.34) 2 i i=1 where ESS is the explained sum of squares, RSS is the residual sum, p is the number of parameters in the model and n is the number of data points. We reject this null hypothesis for large values of the F − statistic. This suggests that if the amount of the variance in the data that the regression explains is large relative to the amount which is unexplained, then we reject the null hypothesis of no linear relation. page 12: 9.07 Lecture 14: Simple Linear Regression 2. R-Squared Another measure of goodness-of-fit is the Square of the Multiple Correlation Coefficient or R 2 . It is defined as N R2 = ESS = TSS ∑ ( ŷ − y ) 2 ∑( y − y ) 2 i i =1 N . (14.35) i i=1 The R 2 measures the fraction of variance in the data explained by the regression equation. By the Pythagorean relation in Eq. 14.23, we see that 0 < R2 ≤ 1. For the simple linear regression model R 2 = rxy2 . The R 2 and the F − statistic are related as R2 = F ( p −1) / ( n − p ) 1+ F ( p − 1) / ( n − p ) (14.36) As the F − statistic increases, the R 2 increases. Hence, the greater the magnitude of the F − statistic, the greater the R 2 . Example 14.1 (continued). We have from the Pythagorean relation in Eq. 14.23 ESS 65,041 RSS 2,233 TSS 67,274 and that F1,65 = 1,895 and R 2 = 0.967. Therefore, we reject the null hypothesis of no linear relation in the data with a p-value << 0.05. Furthermore, the R 2 suggests that the regressor, axon diameter can explain 97% of the variance in the axon conduction velocity. Remark 8. The square of the t − statistic with n − p degrees of freedom is the F − statistic with 1 and n − p degrees of freedom. To verify this we find in Table 1 that the ( t − statistic)2 = (43.5)2 = 1,893 which is the value of F1,65 . 3. Graphical Measures of Goodness-of-Fit a. Plot the raw data. Is the relation linear? (Fig. 14.1). b. Plot the residuals versus the covariate x. Is there lack of fit of the model? (Fig. 14.6) c. Plot the residuals versus the predicted values ŷ. Is there a relation? (Fig. 14.7) page 13: 9.07 Lecture 14: Simple Linear Regression Figure 14.6. Plot of Residuals versus Axon Diameter for the Simple Linear Regression Model of the Axon Diameter and Conduction Velocity. The residuals have no real discernible structure except that they seem to grow with the diameter of the axon suggesting that the assumption of homoscedasticity, i.e., that all of the observations have the same variance may be incorrect. When observations have different variances they are said to be heteroscedastic. We might model this by making the variance proportional to the diameter. Figure 14.7. Plot of the Residuals against the Predicted Velocity. The plot of the residuals versus the predicted velocity also shows the heteroscedasticity in the observations. This plot is performed because by the Pythagorean relation, the residuals and the predicted velocity estimates are orthogonal. Hence, if there is any relation between the two it suggests a way in which the model may be misspecified. Inference. Based on this regression analysis axon conduction velocity is strongly associated with (predicted by) axon diameter. This analysis also reveals that a one micron increase in page 14: 9.07 Lecture 14: Simple Linear Regression diameter is associated with a 6 meter per second change in velocity for axon diameters between 2 and 20 microns. Remark 14.9. Correlation does not mean causation. G. The Geometry of Regression Analysis (Method of Least-Squares) As we have shown, regression analysis has an intuitively appealing geometric interpretation defined primarily be Eq. 14.23. In particular, we can understand this geometry by considering two cases: i) Small Residual Error If the residual error is small it follows from the Pythogorean relation that the geometry of the analysis must have a small RSS relative to ESS. ii) Large Residual Error Similarly, if the residual error is large, it follows from the Pythogorean relation that the geometry of the analysis must have a large RSS relative to ESS. III. Summary The simple linear regression model can be used to study how well one variable may be described as a linear function of another variable. The Pythagorean relation is a useful construct for understanding how the analysis framework for the simple linear regression model is constructed. This Pythagorean relation will also be central to our analyses of multiple linear regression and ANOVA models. Acknowledgments I am grateful to Uri Eden for making the figures, to Julie Scott for technical assistance and to Jim Mutch for careful proofreading and comments. The data from the Hursh (1993) were provided by Rob Kass at Carnegie Mellon. Text References DeGroot MH, Schervish MJ. Probability and Statistics, 3rd edition. Boston, MA: Addison Wesley, 2002. Draper NR, Smith H. Applied Regression Analysis, 2nd ed. New York: Wiley, 1981. Rice JA. Mathematical Statistics and Data Analysis, 3rd edition. Boston, MA, 2007. Literature References Hursh JB. The properties of growing nerve fibers. American Journal of Physiology, 1939, 127: 131-39. Wirth S, Yanike M. Frank LM, Smith AC, Brown EN, Suzuki WA. Single neurons in the monkey hippocampus and learning of new associations. Science, 2003, 300: 1578-81.