4.3 Paired vectors
4.3.2 Functional relations: linear regression
Instead of comparing just the means of the size and weight ratings, or comparing their distributions by means of boxplots, we can graph the individual data points in a scatter-plot, as shown in the left panel of Figure 4.8.
91
DRAFT
> plot(ratings$meanWeightRating, ratings$meanSizeRating, + xlab = "mean weight rating", ylab = "mean size rating")
What this panel shows is that the data points pattern into a nearly straight line. In other words, we observe an exceptionally clear linear functional relation between estimated sizeandweight. This functional relation can be visualized by means of a line drawn through the scatter of data points in such a way that the line is as close as possible to each of these data points. The question that arises here is how to obtain this regression line. In order to answer this question, we begin with recapitulating how a straight line is characterized.
Slope and intercept
Consider the two lines shown in Figure 4.9. For the dashed line, theINTERCEPTis2and the slope−2. For the dotted line, the intercept is−2and the slope1. It is easy to see that the intercept is theY-coordinate of the line where it crosses the vertical axis. The slope of the line specifies the direction of the line in terms of how far you have to move along the horizontal axis for a unit change in the vertical direction. For the dashed line, two units down corresponds to one unit to the right. Using∆yand∆xto denote the change iny corresponding to a change inx, we find that we have a slope of∆y/∆x=−2/1 =−2. For the dotted line, moving two units up corresponds with moving two units to the right, so the slope is∆y/∆x= 2/2 = 1.
The functionabline()adds parametrically specified lines to a plot. It takes two arguments, first the intercept, and then the slope. This is illustrated by the following code, which produces Figure 4.9.
> plot(c(-4, 4), c(-4, 4), xlab = "x", ylab = "y", type = "n")
# set up the plot region
> abline(2, -2, lty = 2) # add the lines
> abline(-2, 1, lty = 3)
> abline(h = 0) # and add the axes
> abline(v = 0)
> abline(h = -2, col = "grey") # and ancillary lines in grey
> abline(h = 2, col = "grey")
> abline(v = 1, col = "grey", lty = 2)
> abline(v = 2, col = "grey", lty = 2)
The right panel of Figure 4.8 shows a straight line that has been drawn through the data points in such a way that all the data points are as close to the line as possible. Its intercept is0.527and its slope is0.926.
> plot(ratings$meanWeightRating, ratings$meanSizeRating, + xlab = "mean weight rating", ylab = "mean size rating", + col = "darkgrey")
> abline(0.527, 0.926)
The question, of course, is how to determine this slope and intercept.
92
DRAFT
−4 −2 0 2 4
−4−2024
x
y
Figure 4.9: Straight lines are defined by intercept and slope.
Estimating slope and intercept
We estimate slope and intercept with the help of the function for linear modeling,lm().
This function needs to be told what variable is the dependent variable (the variable on theY-axis) and what variable is the predictor (the variable on theX-axis). We provide this information by means of a formula that we supply as the first argument tolm().
> ratings.lm = lm(meanSizeRating ˜ meanWeightRating, data = ratings) The formula specifies thatmeanSizeRatingis to be modelled as a function of, or de-pending on,meanWeightRating. The second argument tellsRto look for these two variables in the data frameratings. The output oflm()is aLINEAR MODELobject that we name after its input and the function that created it. By typingratings.lmto the prompt, we get to see the coefficients of the desiredLEAST SQUARES REGRESSION LINE. (The termLEAST SQUARESrefers to the way in which slope and intercept are estimated, namely, by minimizing the squared vertical differences between the data points and the line.)
> ratings.lm Call:
lm(formula = meanSizeRating ˜ meanWeightRating, data = ratings) Coefficients:
(Intercept) meanWeightRating
0.5270 0.9265
We can extract from the model object a vector with just the intercept and slope with the functioncoef(), which returns the model’sCOEFFICIENTS.
93
DRAFT
> coef(ratings.lm)
(Intercept) meanWeightRating 0.5269981 0.9264743
In order to add this regression line to our scatterplot, we simply type
> abline(ratings.lm)
asabline()is smart enough to extract slope and intercept from the linear model object by itself.
Correlation
You now know how to estimate the intercept and the slope of a regression line. There is much more to be learned from a linear model than just this. We illustrate this by looking in some more detail at scatterplots of paired standard normal random variables. Each panel of Figure 4.10 plots random samples of such paired vectors. The technical name for such paired distributions is aBIVARIATE STANDARD NORMAL DISTRIBUTION. The dashed line in these panels represents the lineY =X, the solid line the regression line. In the upper left panel, we have a scatter of points roughly in the form of a disc. Many points are far away from the regression line, which happens to have a negative slope. The upper right panel also shows a wide scatter, but here the regression line has a positive slope. The points in the lower left panel are somewhat more concentrated and closer to the regression line. The regression line itself is becoming more similar to the lineY =X. Finally, the lower right panel has a regression line that has crept even closer to the dashed line, and the data points are again much closer to the regression line.
The technical term for the degree to which the data points cluster around the regres-sion line isCORRELATION. This degree of correlation is quantified by means of aCORRE -LATION COEFFICIENT. The correlation coefficient of a given population is denoted byρ, and that of a sample from that population byr. The correlation coefficient is bounded by−1(a perfect negative correlation) and+1(a perfect positive correlation). When the correlation is−1or+1, all the data points lie exactly on the regression line, and in that case the regression line is equal to the lineY =−XandY =Xrespectively. This is a limiting case that never arises in practice.
The sample correlationrfor each of the four scatterplots in Figure 4.10 is listed above each panel, and varies from−0.07in the upper left to0.89in the lower right. You can regardras a measure of how useful it is to fit a straight line to the data. Ifris close to zero, the regression line does not help at all to predict where the data points will be for a given value of the predictor variable. This is easy to see by comparing the upper and lower panels. In the upper panels, the scatter along theY-axis is very large for almost all values ofX. For a large majority of observed data points, the predicted value (somewhere on the regression line) is going to be way off. This changes in the lower panels, where the regression line starts to become predictive. Another way of thinking aboutris that it tells us something about how much of the scatter we get a handle on. In fact, the appro-priate measure for evaluating how much of the scatter is accounted for, or explained, by
94
DRAFT
−3 −2 −1 0 1 2
−3−2−1012
x
y
r = −0.067
−2 −1 0 1 2
−2−1012
x
y
r = 0.2
−3 −2 −1 0 1 2
−2−1012
x
y
r = 0.61
−3 −2 −1 0 1 2
−3−2−1012
x
y
r = 0.888
Figure 4.10: Scatterplots for four paired standard normal random variables with different population correlations (ρ = 0.0,0.2,0.5,0.9). The correlations shown above each panel are the sample correlations.
95
DRAFT
the model is notritself, butr2(often denoted byR2). More precisely,R2quantifies the proportion of the variance in the data that is captured and explained by the regression model.
Let’s pause for a moment to think about what it means to explain variance. When we try to fit a model to a data set, the goal is to be able to predict what the value of the dependent variable is given the predictors. The better we succeed in predicting, the better the predictors succeed in explaining the variability in the dependent variable. When you are in bad luck, with lousy predictors, there is little variability that your model explains.
In that case, the values of the dependent variable jump around almost randomly. In this situation,R2will be close to zero. The better the model, the smaller the random variation, the variation that we do not yet understand, will be, and the closerR2will be to one.
Scatterplots like those shown in the panels of Figure 4.10 can be obtained with the help ofmvrnormplot.fnc(),
> mvrnormplot.fnc(r = 0.9)
a convenience functionmvrnormplot.fnc()defined in thelanguageRpackage. As we are dealing with random numbers, your output will be somewhat different each time you run this code, even for the samer. You should try outmvrnormplot.fnc() with different values ofrto acquire some intuitions about what correlations of different strengths look like.
Summarizing a linear model object
We return to our running example of mean size and weight ratings. Recall that we created a linear model object,ratings.lm, and extracted the coefficients of the regression line from this object. If we summarize the model withsummary(), we obtain much more detailed information, including information aboutR2.
> summary(ratings.lm) Call:
lm(formula = meanSizeRating ˜ meanWeightRating, data = ratings) Residuals:
Min 1Q Median 3Q Max
-0.096368 -0.020285 0.002058 0.024490 0.075310 Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.526998 0.010632 49.57 <2e-16 meanWeightRating 0.926474 0.003837 241.45 <2e-16 Residual standard error: 0.03574 on 79 degrees of freedom Multiple R-Squared: 0.9986, Adjusted R-squared: 0.9986 F-statistic: 5.83e+04 on 1 and 79 DF, p-value: < 2.2e-16 In what follows, we will walk through this summary line by line.
The first thing that the summary does is remind us of how the object was created. We then get a brief summary of the distribution of the residuals. We postpone to Chapter 6
96
DRAFT
the discussion of what the residuals are and why they are so important to be mentioned in the summary of the model.
Next, we see a table with the coefficients of the model: a coefficient for the intercept (0.527) and a coefficient for the slope (0.926). Each coefficient comes with three other numbers: its standard error, at-value, and ap-value. Thep-value tells us whether the coefficient is significantly different from zero. If the coefficient for a predictor is zero, there is no relation at all between the predictor and the dependent variable, in which case it is worthless as a predictor. In order to ascertain whether a coefficient is significantly different from zero, and hence potentially useful, a two-tailedt-test is carried out, using thet-value and the associated degrees of freedom (79, this number is listed further down in the summary). Thet-value itself is the value of the coefficient divided by its standard error. This standard error is a measure of how sure we are about the estimate of the coefficient. The smaller the standard error, the smaller the confidence interval around the estimate, the less likely that zero will be included in the acceptance region, and hence the smaller the probability that it might just as well be zero.
Sometimes, it is useful to be able to access the different parts of the summary. You can identify the components of the summary withnames(summary(ratings.lm)), and we can extract these components from the summary with the help of the$operator. For instance, we obtain the table of coefficients with$coef.
> summary(ratings.lm)$coef
Estimate Std. Error t value Pr(>|t|) (Intercept) 0.5269981 0.010632282 49.56585 2.833717e-61 meanWeightRating 0.9264743 0.003837106 241.45129 4.380725e-115
Because this table is a matrix, we can access thet-values or the estimates of the coefficients themselves:
> summary(ratings.lm)$coef[ ,3]
(Intercept) meanWeightRating 49.56585 241.45129
> summary(ratings.lm)$coef[ ,1]
(Intercept) meanWeightRating 0.5269981 0.9264743
Sincesummary(ratings.lm)$coefis not a data frame, we cannot reference columns by name with the $ operator, unfortunately. To do so, we first have to convert it explicitly into a data frame.
> data.frame(summary(ratings.lm)$coef)$Estimate [1] 0.5269981 0.9264743
Let’s return to the summary, and proceed to its last three lines. TheRESIDUAL STAN -DARD ERRORis a measure of how unsuccesful the model is, it gauges the variability in the dependent variable that we can’t handle through the predictor variables. The better a model is, the smaller its residual standard error will be. The next line states that the
97
DRAFT
multipleR-squared equals0.9986. ThisR-squared is the squared correlation coefficient, r2, which quantifies, on a scale from0to1, the proportion of the variance that the model explains. We get the value of the correlation coefficientrby taking the square root of 0.9986, which is0.9993. This is a bit cumbersome, but, fortunately, there are quicker ways of calculatingr. The functioncor()returns the correlation coefficient,
> cor(ratings$meanSizeRating, ratings$meanWeightRating) [1] 0.9993231
andcor.test()provides the correlation coefficient and also tests whether it is signifi-cantly different from zero. It also lists a95% confidence interval.
> cor.test(ratings$meanSizeRating, ratings$meanWeightRating) Pearson’s product-moment correlation
data: ratings$meanSizeRating and ratings$meanWeightRating t = 241.4513, df = 79, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval:
0.9989452 0.9995657 sample estimates:
cor 0.9993231
There is also a distribution free, non-parametric correlation test, which does not depend on the input vectors being approximally normally distributed, the Spearman correlation test, which is based on the ranks of the observations in the two vectors. It is carried out bycor.test()when you specify the optionmethod="spearman":
> cor.test(ratings$meanSizeRating, ratings$meanWeightRating, + method = "spearman")
Spearman’s rank correlation rho
data: ratings$meanSizeRating and ratings$meanWeightRating S = 118, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0 sample estimates:
rho 0.9986676
Warning message: p-values may be incorrect due to ties
We could have avoided the warning message by adding some jitter to the ratings, but given the very lowp-value, this is superfluous. The Spearman correlation coefficient is often referenced asrs.
Returning to the summary ofratings.lmand leaving the discussion of the adjusted R-squared to Chapter 6, we continue with the last line, which lists anF-value. This F value goes with an overall test of whether the linear model as a whole succeeds in explaining a significant portion of the variance. Given the smallp-value listed in the summary, there is no question about lack of statistical significance.
98
DRAFT
0 500 1000 1500
0200400600800
frequency of singular
frequency of plural
2 3 4 5 6 7
01234567
log(frequency of singular+1)
log(frequency of plural+1)
Figure 4.11: Scatterplots for singular and plural frequency with regression lines. Solid lines represent ordinary least squares regression on all data points, the dashed line rep-resents an ordinary least squares regression with 4 outliers excluded, and dotted lines represents robust regression lines obtained withlmsreg().
Problems and pitfalls of linear regression
Now that we have seen how to fit a linear model to a data set with paired vectors, we proceed to two more complex examples that illustrate some of the problems and pitfalls of linear regression. First consider the left panel of Figure 4.11, which plots the frequency of the plural against the frequency of the singular for the81nouns for animals and plants in theratingsdata frame. The problem that we are confronted with here is that there is a cluster of observations near the origin combined with a handful of atypical points with very high values. The presence of suchOUTLIERSmay mislead the algorithm that estimates the coefficients of the linear model. If we fit a linear model to these data points, we obtain the solid line. But if we exclude just the4words with singular frequencies greater than500, and then refit the model, we obtain the dashed line. The two lines tell a rather different story which suggests that these4words are atypical with respect to the lower-frequency words. There are various regression techniques that are more robust with respect to outliers than islm(). The dotted line illustrates thelmsreg()function, which, unfortunately, does not tell us whether the predictors are significant. From the graph we can tell that it considers rather different words to be outliers, namely, the words with high plural frequency but singular frequency less than500.
Before we move on to a better solution for this regression problem, let’s first review the code for the left panel of Figure 4.11.
> plot(ratings$FreqSingular, ratings$FreqPlural)
99
DRAFT
> abline(lm(FreqPlural ˜ FreqSingular, data = ratings), lty = 1)
> abline(lm(FreqPlural ˜ FreqSingular,
+ data = ratings[ratings$FreqSingular < 500, ]), lty = 2) In order to have access tolmsreg(), we must first load theMASSpackage.
> library(MASS)
> abline(lmsreg(FreqPlural ˜ FreqSingular, data = ratings), lty = 3) The problem illustrated in the left panel of Figure 4.11 is that word frequency distri-butions are severely skewed. There are many low-probability words and relatively few high-probability words. This skewness poses a technical problem tolm(). A few high-probability outliers become overly influential, and shift the slope and intercept to such an extent that it becomes suboptimal for the majority of data points. The technical solution is to apply a logarithmic transformation in order to remove at least a substantial amount of this skewness by bringing many straying outliers back into the fold. The right panel of Figure 4.11 visualizes these benefits of the logarithmic transforms. We now have a regres-sion line that captures the main trend in the data quite well. The robust regresregres-sion line has nearly the same slope, albeit a slightly higher intercept. It is influenced less by the four data points with exceptionally low plural frequencies given their singular frequencies, which have a small but apparently somewhat disproportionate effect onlm()’s estimate of the intercept. In Chapter 6, we will discuss in more detail how undue influence of po-tential outliers can be detected and what measures can be taken to protect one’s model against them.
The second example addresses the relation between the mean familiarity rating and mean size rating for our81nouns in theratingsdata set. The question of interest is whether it is possible to predict how heavy people think an object is from how frequently they think the name for that object is used in the language. We address this question with lm(),
> ratings.lm = lm(meanSizeRating ˜ meanFamiliarity, data = ratings) extract the table of coefficients from the summary, and round it to four decimal digits.
> round(summary(ratings.lm)$coef, 4)
Estimate Std. Error t value Pr(>|t|) (Intercept) 3.7104 0.4143 8.9549 0.0000 meanFamiliarity -0.2066 0.1032 -2.0014 0.0488
The summary presents a negative coefficient formeanFamiliaritythat is just signif-icant at the5% level. This suggests that objects that participants judge to have more familiar names in the language receive somewhat lower size ratings.
This conclusion is, however, unwarranted as there are lots of things wrong with this analysis. But this becomes apparent only by graphical inspection of the data and of the predictions of the model. Let’s make a scatterplot of the data, the first thing that we should have done anyway. The scatterplot smoother (lowess()) shown in the upper left
100
DRAFT
Figure 4.12: Scatterplots for mean rated size as a function of mean familiarity, with scat-terplot smoothers (left panels) and linear (upper right, lower left) and quadratic (lower right) fits. The upper panels show fits to all data points, the lower panels show fits to the words for plants (p) and animals (a) separately.
101
DRAFT
panel of Figure 4.12 suggests a negative correlation, but what is worrying is that there are no points close to the line in the center of the graph. The same holds for the regression line for the model that we just fitted to the data withlm(), as shown in the upper right panel.
If you look carefully at the scatterplots, you can see that there seem to be two separate strands of data points, one with higher size ratings, and one with lower size ratings. This
If you look carefully at the scatterplots, you can see that there seem to be two separate strands of data points, one with higher size ratings, and one with lower size ratings. This