Section 14.3 Linear Regression
Data that describes the joint distribution of two numeric variables can be represented with a scatter plot. The \(y\)-axis in this plot corresponds to the response and the \(x\)-axis corresponds to the explanatory variable. The regression line describes the linear trend of the response as a function of the explanatory variable. This line is characterized by a linear equation with an intercept and a slope that are computed from the data.
In the first subsection we present the computation of the regression linear equation from the data. The second subsection discusses regression as a statistical model. Statistical inference can be carried out on the basis of this model. In the context of the statistical model, one may consider the intercept and the slope of the regression model that is fitted to the data as point estimates of the model’s parameter. Based on these estimates, one may test hypotheses regarding the regression model and construct confidence intervals for parameters.
Subsection 14.3.1 Fitting the Regression Line
The
R function that fits the regression line to data is called “lm”, an acronym for Linear Model. The input to the function is a formula, with the response variable to the left of the tilde character and the explanatory variable to the right of it. The output of the function is the fitted linear regression model.
Let us apply the linear regression function to the data on the weight and the length of fish. The output of the function is saved by us in a object called “
fit”. Subsequently, the content of the object “fit” is displayed:
fit <- lm(y~x)
fit
## ## Call: ## lm(formula = y ~ x) ## ## Coefficients: ## (Intercept) x ## 4.616 1.427
When displayed, the output of the function “
lm” shows the formula that was used by the function and provides the coefficients of the regression linear equation. Observe that the intercept of the line is equal to 4.616. The slope of the line, the coefficient that multiplies “x” in linear equation, is equal to 1.427.
One may add the regression line to the scatter plot with the aid of the function “
abline”:
plot(y~x)
abline(fit)
The first expression produces the scatter plot of the data on fish. The second expression adds the regression line to the scatter plot. When the input to the graphical function “
abline” is the output of the function “lm” that fits the regression line, then the result is the addition of the regression line to the existing plot. The line that is added is the line characterized by the coefficients that are computed by the function “lm”. The coefficients in the current setting are 4.616 for the intercept and 1.427 for the slope.
The scatter plot and the added regression line are displayed in Figure 14.3.1. Observe that line passes through the points, balancing between the points that are above the line and the points that are below. The line captures the linear trend in the data.
Examine the line in Figure 14.3.1. When \(x=1\) then the \(y\) value of the line is slightly above 6. When the value of \(x\) is equal to 2, a change of one unit, then value of \(y\) is below 8, and is approximately equal to 7.5. This observation is consistent with the fact that the slop of the line is 1.427. The value of \(x\) is decreased by 1 when changing from \(x=1\) to \(x=0\text{.}\) Consequently, the value of \(y\) when \(x=0\) should decrease by 1.427 in comparison to its value when \(x=1\text{.}\) The value at \(x=1\) is approximately 6. Therefore, the value at \(x=0\) should be approximately 4.6. Indeed, we do get that the intercept is equal to 4.616.
The coefficients of the regression line are computed from the data and are hence statistics. Specifically, the slope of the regression line is computed as the ratio between the covariance of the response and the explanatory variable, divided by the variance of the explanatory variable. The intercept of the regression line is computed using the sample averages of both variables and the computed slope.
Start with the slope. The main ingredient in the formula for the slope, the numerator in the ratio, is the covariance between the two variables. The covariance measures the joint variability of two variables. Recall that the formula for the sample variance of the variable \(x\) is equal to:
\begin{equation*}
s^2 = \frac{\mbox{Sum of the squares of the deviations}}{\mbox{Number of values in the sample}-1} = \frac{\sum_{i=1}^n (x_i - \bar x)^2}{n-1}\;.
\end{equation*}
The formula of the sample covariance between \(x\) and \(y\) replaces the square of the deviations by the product of deviations. The product is between an \(y\) deviation and the parallel \(x\) deviation:
\begin{equation*}
\mbox{covariance} = \frac{\mbox{Sum of products of the deviations}}{\mbox{Number of values in the sample}-1} = \frac{\sum_{i=1}^n (y_i-\bar y)(x_i - \bar x)}{n-1}\;.
\end{equation*}
The function “
cov” computes the sample covariance between two numeric variables. The two variables enter as arguments to the function and the sample covariance is the output. Let us demonstrate the computation by first applying the given function to the data on fish and then repeating the computations without the aid of the function:
cov(y,x)
sum((y-mean(y))*(x-mean(x)))/9
## [1] 2.386111 ## [1] 2.386111
In both cases we obtained the same result. Notice that the sum of products of deviations in the second expression was divided by 9, which is the number of observations, minus 1.
The slope of the regression line is the ratio between the covariance and the variance of the explanatory variable.
The regression line passes through the point \((\bar x, \bar y)\text{,}\) a point that is determined by the means of the both the explanatory variable and the response. It follows that the intercept should obey the equation:
\begin{equation*}
\bar y = a + b\cdot \bar x \quad\Longrightarrow\quad a = \bar y - b\cdot \bar x\;,
\end{equation*}
The left-hand-side equation corresponds to the statement that the value of the regression line at the average \(\bar x\) is equal to the average of the response \(\bar y\text{.}\) The right-hand-side equation is the solution to the left-hand-side equation.
One may compute the coefficients of the regression model manually by computing first the slope as a ratio between the covariance and the variance of explanatory variable. The intercept can then be obtained by the equation that uses the computed slope and the averages of both variables:
b <- cov(x,y)/var(x)
a <- mean(y) - b*mean(x)
a
b
## [1] 4.616539 ## [1] 1.427396
Applying the manual method we obtain, after rounding up, the same coefficients that were produced by the application of the function “
lm” to the data.
As an exercise, let us fit the regression model to the data on the relation between the response “
horsepower” and the explanatory variable “engine.size”. Apply the function “lm” to the data and present the results:
fit.power <- lm(horsepower ~ engine.size, data=cars)
fit.power
## ## Call: ## lm(formula = horsepower ~ engine.size, data = cars) ## ## Coefficients: ## (Intercept) engine.size ## 6.641 0.770
The fitted regression model is stored in an object called “
fit.power”. The intercept in the current setting is equal to 6.6414 and the slope is equal to 0.7695.
Observe that one may refer to variables that belong to a data frame, provided that the name of the data frame is entered as the value of the argument “
data” in the function “lm”. Here we refer to variables that belong to the data frame “cars”.
Next we plot the scatter plot of the data and add the regression line:
plot(horsepower ~ engine.size, data=cars)
abline(fit.power)
The output of the plotting functions is presented in Figure 14.3.2. Again, the regression line describes the general linear trend in the data. Overall, with the increase in engine size one observes increase in the power of the engine.
Subsection 14.3.2 Inference
Up to this point we have been considering the regression model in the context of descriptive statistics. The aim in fitting the regression line to the data was to characterize the linear trend observed in the data. Our next goal is to deal with regression in the context of inferential statistics. The goal here is to produce statements on characteristics of an entire population on the basis of the data contained in the sample.
The foundation for statistical inference in a given setting is a statistical model that produces the sampling distribution in that setting. The sampling distribution is the frame of reference for the analysis. In this context, the observed sample is a single realization of the sampling distribution, one realization among infinitely many potential realizations that never take place. The setting of regression involves a response and an explanatory variable. We provide a description of the statistical model for this setting.
The relation between the response and the explanatory variable is such that the value of the later affects the distribution of the former. Still, the value of the response is not uniquely defined by the value of the explanatory variable. This principle also hold for the regression model of the relation between the response \(Y\) and the explanatory variable \(X\text{.}\) According to the model of linear regression the value of the expectation of the response for observation \(i\text{,}\) \(\Expec(Y_i)\text{,}\) is a linear function of the value of the explanatory variable for the same observation. Hence, there exist and intercept \(a\) and a slope \(b\text{,}\) common for all observations, such that if \(X_i = x_i\) then
\begin{equation*}
\Expec(Y_i) = a + b \cdot x_i\;.
\end{equation*}
The regression line can thus be interpreted as the average trend of the response in the population. This average trend is a linear function of the explanatory variable.
The intercept \(a\) and the slope \(b\) of the statistical model are parameters of the sampling distribution. One may test hypotheses and construct confidence intervals for these parameters based on the observed data and in relation to the sampling distribution.
Consider testing hypothesis. A natural null hypothesis to consider is the hypothesis that the slope is equal to zero. This hypothesis corresponds to statement that the expected value of the response is constant for all values of the explanatory variable. In other words, the hypothesis is that the explanatory variable does not affect the distribution of the response. One may formulate this null hypothesis as \(H_0:b = 0\) and test it against the alternative \(H_1: b \not= 0\) that states that the explanatory variable does affect the distribution of the response.
1
According to the model of linear regression, the only effect of the explanatory variable on the distribution of the response is via the expectation. If such an effect, according to the null hypothesis, is also excluded then the so called explanatory variable is not effecting at all the distribution of the response.
A test of the given hypotheses can be carried out by the application of the function “
summary” to the output of the function “lm”. Recall that the function “lm” was used in order to fit the linear regression to the data. In particular, this function was applied to the data on the relation between the size of the engine and the power that the engine produces. The function fitted a regression line that describes the linear trend of the data. The output of the function was saved in an object by the name “fit.power”. We apply the function “summary” to this object:
summary(fit.power)
## ## Call: ## lm(formula = horsepower ~ engine.size, data = cars) ## ## Residuals: ## Min 1Q Median 3Q Max ## -59.643 -12.282 -5.515 10.251 125.153 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 6.64138 5.23318 1.269 0.206 ## engine.size 0.76949 0.03919 19.637 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 23.31 on 201 degrees of freedom ## (2 observations deleted due to missingness) ## Multiple R-squared: 0.6574, Adjusted R-squared: 0.6556 ## F-statistic: 385.6 on 1 and 201 DF, p-value: < 2.2e-16
The output produced by the application of the function “
summary” is long and detailed. We will discuss this output in the next section. Here we concentrate on the table that goes under the title “Coefficients:”. The said table is made of 2 rows and 4 columns. It contains information for testing, for each of the coefficients, the null hypothesis that the value of the given coefficient is equal to zero. In particular, the second row may be used is order to test this hypothesis for the slope of the regression line, the coefficient that multiplies the explanatory variable.
Consider the second row. The first value on this row is 0.76949, which is equal (after rounding up) to the slope of the line that was fitted to the data in the previous subsection. However, in the context of statistical inference this value is the estimate of the slope of the population regression coefficient, the realization of the estimator of the slope.
2
The estimator of the slope is obtained via the application of the formula for the computation of the slope to the sample: \(\frac{1}{n-1}\sum_{i=1}^n (Y_i-\bar Y)(X_i - \bar X)/\frac{1}{n-1}\sum_{i=1}^n (X_i - \bar X)^2\text{.}\)
The second value is 0.03919. This is an estimate of the standard deviation of the estimator of the slope. The third value is the test statistic. This statistic is the ratio between the deviation of the sample estimate of the parameter (0.76949) from the value of the parameter under the null hypothesis (0), divided by the estimated standard deviation (0.03919): \((0.76949 - 0)/0.03919 = 0.76949/0.03919 = 19.63486\text{,}\) which is essentially the value given in the report.
3
Our computation involves rounding up errors, hence the small discrepancy between the value we computed and the value in the report.
The last value is the computed \(p\)-value for the test. It can be shown that the sampling distribution of the given test statistic, under the null distribution which assumes no slope, is asymptotically the standard Normal distribution. If the distribution of the response itself is Normal then the distribution of the statistic is the \(t\)-distribution on \(n-2\) degrees of freedom. In the current situation this corresponds to 201 degrees of freedom. The computed \(p\)-value is extremely small, practically eliminating the possibility that the slope is equal to zero.
4
Notice that the “
horsepower” measurement is missing for two observation. These observations are deleted for the analysis, leaving a total of \(n=203\) observations. The number of degrees of freedom is \(n-2 = 203-2=201\text{.}\)
The first row presents information regarding the intercept. The estimated intercept is 6.64138 with an estimated standard deviation of 5.23318. The value of the test statistic is 1.269 and the \(p\)-value for testing the null hypothesis that the intercept is equal to zero against the two sided alternative is 0.206. In this case the null hypothesis is not rejected since the \(p\)-value is larger than 0.05.
The report contains an inference for the intercept. However, one is advised to take this inference in the current case with a grain of salt. Indeed, the intercept is the expected value of the response when the explanatory variable is equal to zero. Here the explanatory variable is the size of the engine and the response is the power of that engine. The power of an engine of size zero is a quantity that has no physical meaning! In general, unless the intercept is in the range of observations (i.e. the value 0 is in the range of the observed explanatory variable) one should treat the inference on the intercept cautiously. Such inference requires extrapolation and is sensitive to the miss-specification of the regression model.
Apart from testing hypotheses one may also construct confidence intervals for the parameters. A crude confidence interval may be obtained by taking 1.96 standard deviations on each side of the estimate of the parameter. Hence, a confidence interval for the slope is approximately equal to \(0.76949 \pm 1.96\times 0.03919 = [0.6926776, 0.8463024]\text{.}\) In a similar way one may obtain a confidence interval for the slope: \(6.64138 \pm 1.96\times 5.23318 = [-3.615653, 16.89841]\text{.}\)
5
The warning message that was made in the context of testing hypotheses on the intercept should be applied also to the construction of confidence intervals. If the value 0 is not in the range of the explanatory variable then one should be careful when interpreting a confidence interval for the intercept.
Alternatively, one may compute confidence intervals for the parameters of the linear regression model using the function “
confint”. The input to this function is the fitted model and the output is a confidence interval for each of the parameters:
confint(fit.power)
## 2.5 % 97.5 % ## (Intercept) -3.6775989 16.9603564 ## engine.size 0.6922181 0.8467537
Observe the similarity between the confidence intervals that are computed by the function and the crude confidence intervals that were produced by us. The small discrepancies that do exist between the intervals result from the fact that the function “
confint” uses the \(t\)-distribution whereas we used the Normal approximation.
Subsection 14.3.3 R-squared and the Variance of Residuals
In this section we discuss the residuals between the values of the response and their estimated expected value according to the regression model. These residuals are the regression model equivalence of the deviations between the observations and the sample average. We use these residuals in order compute the variability that is not accounted for by the regression model. Indeed, the ratio between the total variability of the residuals and the total variability of the deviations from the average serves as a measure of the variability that is not explained by the explanatory variable. R-squared, which is equal to 1 minus this ratio, is interpreted as the fraction of the variability of the response that is explained by the regression model.
We start with the definition of residuals. Let us return to the artificial example that compared length of fish to their weight. The data for this example was given in Table 14.2.1 and was saved in the objects “
x” and “y”. The regression model was fitted to this data by the application of the function “lm” to the formula “y~x” and the fitted model was saved in an object called “fit”. Let us apply the function “summary” to the fitted model:
summary(fit)
The given report contains a table with estimates of the regression coefficients and information for conducting hypothesis testing. The report contains other information that is associated mainly with the notion of the residuals from regression line. Our current goal is to understand what is that other information.
The residual from regression for each observation is the difference between the value of the response for the observation and the estimated expectation of the response under the regression model. An observation is a pair \((x_i,y_i)\text{,}\) with \(y_i\) being the value of the response. The expectation of the response according to the regression model is \(a + b \cdot x_i\text{,}\) where \(a\) and \(b\) are the coefficients of the model. The estimated expectation is obtained by using, in the formula for the expectation, the coefficients that are estimated from the data. The residual is the difference between \(y_i\) and \(a + b \cdot x_i\text{.}\)
6
The estimated expectation of the response is also called the predicted response.
Consider an example. The first observation on the fish is \((4.5, 9.5)\text{,}\) where \(x_1 = 4.5\) and \(y_1 = 9.5\text{.}\) The estimated intercept is 4.6165 and the estimated slope is 1.4274. The estimated expectation of the response for the first variable is equal to
\begin{equation*}
4.6165 + 1.4274 \cdot x_1 = 4.6165 + 1.4274 \cdot 4.5 = 11.0398\;.
\end{equation*}
The residual is the difference between the observes response and this value:
\begin{equation*}
y_1 - (4.6165 + 1.4274 \cdot x_1) = 9.5 - 11.0398 = -1.5398\;.
\end{equation*}
The residuals for the other observations are computed in the same manner. The values of the intercept and the slope are kept the same but the values of the explanatory variable and the response are changed.
Consult the upper plot in Figure 14.3.3. This is a scatter plot of the data, together with the regression line in black and the line of the average in red. A vertical arrow extends from each data point to the regression line. The point where each arrow hits the regression line is associated with the estimated value of the expectation for that point. The residual is the difference between the value of the response at the origin of the arrow and the value of the response at the tip of its head. Notice that there are as many residuals as there are observations.
The function “
residuals” computes the residuals. The input to the function is the fitted regression model and the output is the sequence of residuals. When we apply the function to the object “fit”, which contains the fitted regression model for the fish data, we get the residuals:
residuals(fit)
Indeed, 10 residuals are produced, one for each observation. In particular, the residual for the first observation is -1.5397075, which is essentially the value that we obtained.
7
The discrepancy between the value that we computed and the value computed by the function results from rounding up errors. We used the vales of the coefficients that appear in the report. These values are rounded up. The function “
residuals” uses the coefficients without rounding.
Return to the report produced by the application of the function “
summary” to the fitted regression model. The first component in the report is the formula that identifies the response and the explanatory variable. The second component, the component that comes under the title “Residuals:”, gives a summary of the distribution of the residuals. This summary includes the smallest and the largest values in the sequence of residuals, as well as the first and third quartiles and the median. The average is not reported since the average of the residuals from the regression line is always equal to 0.
The table that contains information on the coefficients was discussed in the previous section. Let us consider the last 3 lines of the report.
The first of the three lines contains the estimated value of the standard deviation of the response from the regression model. If the expectations of the measurements of the response are located on the regression line then the variability of the response corresponds to the variability about this line. The resulting variance is estimated by the sum of squares of the residuals from the regression line, divided by the number of observations minus 2. A division by the number of observation minus 2 produces an unbiased estimator of the variance of the response about the regression model. Taking the square root of the estimated variance produces an estimate of the standard deviation:
sqrt(sum(residuals(fit)^2)/8)
The last computation is a manual computation of the estimated standard deviation. It involves squaring the residuals and summing the squares. This sum is divided by the number of observations minus 2 (\(10-2=8\)). Taking the square root produces estimate. The value that we get for the estimated standard deviation is 2.790787, which coincides with the value that appears in the first of the last 3 lines of the report.
The second of these lines reports the R-squared of the linear fit. In order to explain the meaning of R-squared let us consider Figure 14.3.3 once again. The two plots in the figure present the scatter plot of the data together with the regression line and the line of the average. Vertical black arrows that represent the residuals from the regression are added to the upper plot. The lower plot contains vertical red arrows that extend from the data points to the line of the average. These arrows represent the deviations of the response from the average.
Consider two forms of variation. One form is the variation of the response from its average value. This variation is summarized by the sample variance, the sum of the squared lengths of the red arrows divided by the number of observations minus 1. The other form of variation is the variation of the response from the fitted regression line. This variation is summarized by the sample variation of the residuals, the sum of squared lengths of the black arrows divided by the number of observations minus 1. The ratio between these two quantities gives the relative variability of the response that remains after fitting the regression line to the data.
The line of the average is a straight line. The deviations of the observations from this straight line can be thought of as residuals from that line. The variability of these residuals, the sum of squares of the deviations from the average divided by the number of observations minus 1, is equal to the sample variance.
The regression line is the unique straight line that minimizes the variability of its residuals. Consequently, the variability of the residuals from the regression, the sum of squares of the residuals from the regression divided by the number of observations minus 1, is the smallest residual variability produced by any straight line. It follows that the sample variance of the regression residuals is less than the sample variance of the response. Therefore, the ratio between the variance of the residuals and the variance of the response is less than 1.
R-squared is the difference between 1 and the ratio of the variances. Its value is between 0 and 1 and it represents the fraction of the variability of the response that is explained by the regression line. The closer the points are to the regression line the larger the value of R-squared becomes. On the other hand, the less there is a linear trend in the data the closer to 0 is the value of R-squared. In the extreme case of R-squared equal to 1 all the data point are positioned exactly on a single straight line. In the other extreme, a value of 0 for R-squared implies no linear trend in the data.
Let us compute manually the difference between 1 and the ratio between the variance of the residuals and the variance of the response:
1-var(residuals(fit))/var(y)
Observe that the computed value of R-squared is the same as the value “
Multiple R-squared: 0.4281” that is given in the report.
The report provides another value of R-squared, titled Adjusted R-squared. The difference between the adjusted and unadjusted quantities is that in the former the sample variance of the residuals from the regression is replaced by an unbiased estimate of the variability of the response about the regression line. The sum of squares in the unbiased estimator is divided by the number of observations minus 2. Indeed, when we re-compute the ratio using the unbiased estimate, the sum of squared residuals divided by \(10 - 2 = 8\text{,}\) we get:
1-(sum(residuals(fit)^2)/8)/var(y)
The value of this adjusted quantity is equal to the value “
Adjusted R-squared: 0.3566” in the report.
Which value of R-squared to use is a matter of personal taste. In any case, for a larger number of observations the difference between the two values becomes negligible.
The last line in the report produces an overall goodness of fit test for the regression model. In the current application of linear regression this test reduces to a test of the slope being equal to zero, the same test that is reported in the second row of the table of coefficients. The \(F\) statistic is simply the square of the \(t\) value that is given in the second row of the table. The sampling distribution of this statistic under the null hypothesis is the \(F\)-distribution on 1 and \(n-2\) degrees of freedom, which is the sampling distribution of the square of the test statistic for the slope. The computed \(p\)-value, “
8
In more complex applications of linear regression, applications that are not considered in this book, the test in the last line of the report and the tests of coefficients do not coincide.
p-value: 0.04003” is the identical (after rounding up) to the \(p\)-value given in the second line of the table.
Return to the R-squared coefficient. This coefficient is a convenient measure of the goodness of fit of the regression model to the data. Let us demonstrate this point with the aid of the “
cars” data. In Subsection Subsection 14.3.1 we fitted a regression model to the power of the engine as a response and the size of the engine as an explanatory variable. The fitted model was saved in the object called “fit.power”. A report of this fit, the output of the expression “summary(fit.power)” was also presented. The null hypothesis of zero slope was clearly rejected. The value of R-squared for this fit was 0.6574. Consequently, about 2/3 of the variability in the power of the engine is explained by the size of the engine.
Consider trying to fit a different regression model for the power of the engine as a response. The variable “
length” describes the length of the car (in inches). How well would the length explain the power of the car? We may examine this question using linear regression:
summary(lm(horsepower ~ length, data=cars))
We used one expression to fit the regression model to the data and to summarize the outcome of the fit.
A scatter plot of the two variables together with the regression line may be produced using the code:
plot(horsepower ~ length, data=cars)
abline(lm(horsepower ~ length, data=cars))
From the examination of the figure we may see that indeed there is a linear trend in the relation between the length and the power of the car. Longer cars tend to have more power. Testing the null hypothesis that the slope is equal to zero produces a very small \(p\)-value and leads to the rejection of the null hypothesis.
The length of the car and the size of the engine are both statistically significant in their relation to the response. However, which of the two explanatory variables produces a better fit?
An answer to this question may be provided by the examination of values of R-squared, the ratio of the variance of the response explained by each of the explanatory variable. The R-squared for the size of the engine as an explanatory variable is 0.6574, which is approximately equal to 2/3. The value of R-squared for the length of the car as an explanatory variable is 0.308, less than 1/3. It follows that the size of the engine explains twice as much of the variability of the power of the engine than the size of car and is a better explanatory variable.
