Download Data Analysis Introduction To The Class, Lecture Slide - Engineering and more Slides Advanced Data Analysis in PDF only on Docsity!
Regression: Predicting and Relating Quantitative
Features
36-402, Advanced Data Analysis
11 January 2011
Reading: Chapter 1 in Faraway, especially up to p. 17. Optional Readings: chapter 1 of Berk; ch. 20 in Adler (through p. 386).
Contents
1 Statistics, Data Analysis, Regression 1
2 Guessing the Value of a Random Variable 3 2.1 Estimating the Expected Value................... 3
3 The Regression Function 4 3.1 Some Disclaimers........................... 4
4 Estimating the Regression Function 7 4.1 The Bias-Variance Tradeoff..................... 8 4.2 The Bias-Variance Trade-Off in Action............... 10 4.3 Ordinary Least Squares Linear Regression as Smoothing..... 10
5 Linear Smoothers 13 5.1 k-Nearest-Neighbor Regression................... 13 5.2 Kernel Smoothers........................... 15
1 Statistics, Data Analysis, Regression
Statistics is the science which uses mathematics to study and improve ways of drawing reliable inferences from incomplete, noisy, corrupt, irreproducible and otherwise imperfect data. The subject of most sciences is some aspect of the world around us, or within us. Psychology studies minds; geology studies the Earth’s composition and form; economics studies production, distribution and exchange; mycology studies mushrooms. Statistics does not study the world, but some of the ways
we try to understand the world — some of the intellectual tools of the other sciences. Its utility comes indirectly, through helping those other sciences. This utility is very great, because all the sciences have to deal with imperfect data. Data may be imperfect because we can only observe and record a small fraction of what is relevant; or because we can only observe indirect signs of what is truly relevant; or because, no matter how carefully we try, our data always contain an element of noise. Over the last two centuries, statistics has come to model all such imperfections by modeling them as random processes, and probability has become so central to statistics that we introduce random events deliberately (as in sample surveys).^1 Statistics, then, uses probability to model inference from data. We try to mathematically understand the properties of different procedures for drawing inferences: Under what conditions are they reliable? What sorts of errors do they make, and how often? What can they tell us when they work? What are signs that something has gone wrong? Like some other sciences, such as engineering, medicine and economics, statistics aims not just at understanding but also at improvement: we want to analyze data better, more reliably, with fewer and smaller errors, under broader conditions, faster, and with less mental effort. Sometimes some of these goals conflict — a fast, simple method might be very error-prone, or only reliable under a narrow range of circumstances. One of the things that people most often want to know about the world is how different variables are related to each other, and one of the central tools statistics has for learning about relationships is regression.^2 In 36-401, you learned how to do linear regression, learned about how it could be used in data analysis, and learned about its properties. In this class, we will build on that foundation, extending beyond simple linear regression in many directions, to answer many questions about how variables are related to each other. This is intimately related to prediction. Being able to make predictions isn’t the only reason we want to understand relations between variables, but prediction tests our knowledge of relations. (If we misunderstand, we might still be able to predict, but it’s hard to see how we could understand and not be able to predict.) So before we go beyond linear regression, we will first look at prediction, and how to predict one variable from nothing at all. Then we will look at predictive relationships between variables, and see how linear regression is just one member of a big family of smoothing methods, all of which are available to us. (^1) Two excellent, but very different, histories of statistics are Hacking (1990) and Porter (1986). (^2) The origin of the name is instructive. It comes from 19th century investigations into the relationship between the attributes of parents and their children. People who are taller (heavier, faster,... ) than average tend to have children who are also taller than average, but not quite as tall. Likewise, the children of unusually short parents also tend to be closer to the average, and similarly for other traits. This came to be called “regression towards the mean”, or even “regression towards mediocrity”, and the word stuck.
with a common expected value. Even if they are correlated, but the correlations decay fast enough, all that changes is the rate of convergence. So “sit, wait, and average” is a pretty reliable way of estimating the expectation value.
3 The Regression Function
Of course, it’s not very useful to predict just one number for a variable. Typ- ically, we have lots of variables in our data, and we believe they are related somehow. For example, suppose that we have data on two variables, X and Y , which might look like Figure 1. The feature Y is what we are trying to predict, a.k.a. the dependent variable or output or response, and X is the predic- tor or independent variable or covariate or input. Y might be something like the profitability of a customer and X their credit rating, or, if you want a less mercenary example, Y could be some measure of improvement in blood cholesterol and X the dose taken of a drug. Typically we won’t have just one input feature X but rather many of them, but that gets harder to draw and doesn’t change the points of principle. Figure 2 shows the same data as Figure 1, only with the sample mean added on. This clearly tells us something about the data, but also it seems like we should be able to do better — to reduce the average error — by using X, rather than by ignoring it. Let’s say that the we want our prediction to be a function of X, namely f (X). What should that function be, if we still use mean squared error? We can work this out by using the law of total expectation, i.e., the fact that E [U ] = E [E [U |V ]] for any random variables U and V.
MSE(f (X)) = E
[
(Y − f (X))^2
]
= E
[
E
[
(Y − f (X))^2 |X
]]
= E
[
Var [Y |X] + (E [Y − f (X)|X])^2
]
When we want to minimize this, the first term inside the expectation doesn’t depend on our prediction, and the second term looks just like our previous optimization only with all expectations conditional on X, so for our optimal function r(x) we get r(x) = E [Y |X = x] (14)
In other words, the (mean-squared) optimal conditional prediction is just the conditional expected value. The function r(x) is called the regression func- tion. This is what we would like to know when we want to predict Y.
3.1 Some Disclaimers
It’s important to be clear on what is and is not being assumed here. Talking about X as the “independent variable” and Y as the “dependent” one suggests
0.0 0.2 0.4 0.6 0.8 1.
x
y
plot(all.x,all.y,xlab="x",ylab="y") rug(all.x,side=1,col="grey") rug(all.y,side=2,col="grey")
Figure 1: Scatterplot of the example data. (These are made up.) The rug commands add horizontal and vertical ticks to the axes to mark the location of the data (in grey so they’re less strong than the main tick-marks). This isn’t necessary but is often helpful. The data are in the example.dat file.
a causal model, which we might write
Y ← r(X) + (15)
where the direction of the arrow, ←, indicates the flow from causes to effects, and is some noise variable. If the gods of inference are very, very kind, then would have a fixed distribution, independent of X, and we could without loss of generality take it to have mean zero. (“Without loss of generality” because if it has a non-zero mean, we can incorporate that into r(X) as an additive constant.) This is the kind of thing we saw with the factor model. However, no such assumption is required to get Eq. 14. It works when predicting effects from causes, or the other way around when predicting (or “retrodicting”) causes from effects, or indeed when there is no causal relationship whatsoever between X and Y. It is always true that
Y |X = r(X) + η(X) (16)
where η(X) is a noise variable with mean zero, but as the notation indicates the distribution of the noise generally depends on X. It’s also important to be clear that when we find the regression function is a constant, r(x) = r 0 for all x, that this does not mean that X and Y are independent. If they are independent, then the regression function is a constant, but turning this around is the logical fallacy of “affirming the consequent”.^3
4 Estimating the Regression Function
We want to find the regression function r(x) = E [Y |X = x], and what we’ve got is a big set of training examples, of pairs (x 1 , y 1 ), (x 2 , y 2 ),... (xn, yn). How should we proceed? If X takes on only a finite set of values, then a simple strategy is to use the conditional sample means:
r̂(x) =
{i : xi = x}
i:xi=x
yi (17)
By the same kind of law-of-large-numbers reasoning as before, we can be confi- dent that ̂r(x) → E [Y |X = x]. Unfortunately, this only works if X has only a finite set of values. If X is continuous, then in general the probability of our getting a sample at any particular value is zero, is the probability of getting multiple samples at exactly the same value of x. This is a basic issue with estimating any kind of function from data — the function will always be undersampled, and we need to fill
(^3) As in combining the fact that all human beings are featherless bipeds, and the observation that a cooked turkey is a featherless biped, to conclude that cooked turkeys are human beings. An econometrician stops there; an econometrician who wants to be famous writes a best-selling book about how this proves that Thanksgiving is really about cannibalism.
in between the values we see. We also need to somehow take into account the fact that each yi is a sample from the conditional distribution of Y |X = xi, and so is not generally equal to E [Y |X = xi]. So any kind of function estimation is going to involve interpolation, extrapolation, and smoothing. Different methods of estimating the regression function — different regres- sion methods, for short — involve different choices about how we interpolate, extrapolate and smooth. This involves our making a choice about how to ap- proximate r(x) by a limited class of functions which we know (or at least hope) we can estimate. There is no guarantee that our choice leads to a good ap- proximation in the case at hand, though it is sometimes possible to say that the approximation error will shrink as we get more and more data. This is an extremely important topic and deserves an extended discussion, coming next.
4.1 The Bias-Variance Tradeoff
Suppose that the true regression function is r(x), but we use the function ̂r to make our predictions. Let’s look at the mean squared error at X = x in a slightly different way than before, which will make it clearer what happens when we can’t use r to make predictions. We’ll begin by expanding (Y − ̂r(x))^2 , since the MSE at x is just the expectation of this.
(Y − ̂r(x))^2 (18) = (Y − r(x) + r(x) − ̂r(x))^2 = (Y − r(x))^2 + 2(Y − r(x))(r(x) − ̂r(x)) + (r(x) − ̂r(x))^2 (19)
We saw above (Eq. 16) that Y − r(x) = η, a random variable which has ex- pectation zero (and is uncorrelated with x). When we take the expectation of Eq. 19, nothing happens to the last term (since it doesn’t involve any random quantities); the middle term goes to zero (because E [Y − r(x)] = E [η] = 0), and the first term becomes the variance of η. This depends on x, in general, so let’s call it σ^2 x. We have
MSE(̂r(x)) = σ^2 x + ((r(x) − r̂(x))^2 (20)
The σ^2 x term doesn’t depend on our prediction function, just on how hard it is, intrinsically, to predict Y at X = x. The second term, though, is the extra error we get from not knowing r. (Unsurprisingly, ignorance of r cannot improve our predictions.) This is our first bias-variance decomposition: the total MSE at x is decomposed into a (squared) bias r(x) − ̂r(x), the amount by which our predictions are systematically off, and a variance σ^2 x, the unpredictable, “statistical” fluctuation around even the best prediction. All of the above assumes that ̂r is a single fixed function. In practice, of course, ̂r is something we estimate from earlier data. But if those data are ran- dom, the exact regression function we get is random too; let’s call this random function R̂ n, where the subscript reminds us of the finite amount of data we used to estimate it. What we have analyzed is really MSE( R̂ n(x)|R̂ n = ̂r), the
Again, consistency depends on how well the method matches the actual data- generating process, not just on the method, and again, there is a bias-variance trade-off. There can be multiple consistent methods for the same problem, and their biases and variances don’t have to go to zero at the same rates.
4.2 The Bias-Variance Trade-Off in Action
Let’s take an extreme example: we could decide to approximate r(x) by a con- stant r 0. The implicit smoothing here is very strong, but sometimes appropri- ate. For instance, it’s appropriate when r(x) really is a constant! Then trying to estimate any additional structure in the regression function is just so much wasted effort. Alternately, if r(x) is nearly constant, we may still be better off approximating it as one. For instance, suppose the true r(x) = r 0 + a sin (νx), where a 1 and ν 1 (Figure 3 shows an example). With limited data, we can actually get better predictions by estimating a constant regression function than one with the correct functional form.
4.3 Ordinary Least Squares Linear Regression as Smooth-
ing
Let’s revisit ordinary least-squares linear regression from this point of view. Let’s assume that the independent variable X is one-dimensional, and that both X and Y are centered (i.e. have mean zero) — neither of these assumptions is really necessary, but they reduce the book-keeping. We choose to approximate r(x) by α + βx, and ask for the best values a, b of those constants. These will be the ones which minimize the mean-squared error.
M SE(α, β) = E
[
(Y − α − βX)^2
]
= E
[
(Y − α − βX)^2 |X
]
= E
[
Var [Y |X] + (E [Y − α − βX|X])^2
]
= E [Var [Y |X]] + E
[
(E [Y − α − βX|X])^2
]
The first term doesn’t depend on α or β, so we can drop it for purposes of optimization. Taking derivatives, and then brining them inside the expectations,
∂M SE ∂α
= E [2(Y − α − βX)(−1)] (31) E [Y − a − bX] = 0 (32) a = E [Y ] − bE [X] = 0 (33)
using the fact that X and Y are centered; and,
∂M SE ∂β
= E [2(Y − α − βX)(−X)] (34)
0.0 0.2 0.4 0.6 0.8 1.
x
y
ugly.func = function(x) {1 + 0.01sin(100x)} r = runif(100); y = ugly.func(r) + rnorm(length(r),0,0.5) plot(r,y,xlab="x",ylab="y"); curve(ugly.func,add=TRUE) abline(h=mean(y),col="red") sine.fit = lm(y ~ 1+ sin(100r)) curve(sine.fit$coefficients[1]+sine.fit$coefficients[2]sin(100*x), col="blue",add=TRUE)
Figure 3: A rapidly-varying but nearly-constant regression function; Y = 1 + 0 .01 sin 100x + , with ∼ N (0, 0 .1). (The x values are uniformly distributed between 0 and 1.) Red: constant line at the sample mean. Blue: estimated function of the same form as the true regression function, i.e., r 0 + a sin 100x. If the data set is small enough, the constant actually generalizes better — the bias of using the wrong functional form is smaller than the additional variance from the extra degrees of freedom. Here, the root-mean-square (RMS) error of the constant on new data is 0.50, while that of the estimated sine function is 0 .51 — using the right function actually hurts us!
where s^2 X is the sample variance of X. In words, our prediction is a weighted average of the observed values yi of the dependent variable, where the weights are proportional to how far xi is from the center, relative to the variance, and proportional to the magnitude of x. If xi is on the same side of the center as x, it gets a positive weight, and if it’s on the opposite side it gets a negative weight. Figure 4 shows the data from Figure 1 with the least-squares regression line added. It will not escape your notice that this is very, very slightly different from the constant regression function; the coefficient on X is 6. 3 × 10 −^3. Visually, the problem is that there should be a positive slope in the left-hand half of the data, and a negative slope in the right, but the slopes are the densities are balanced so that the best single slope is zero.^7 Mathematically, the problem arises from the somewhat peculiar way in which least-squares linear regression smoothes the data. As I said, the weight of a data point depends on how far it is from the center of the data, not how far it is from the point at which we are trying to predict. This works when r(x) really is a straight line, but otherwise — e.g., here — it’s a recipe for trouble. However, it does suggest that if we could somehow just tweak the way we smooth the data, we could do better than linear regression.
5 Linear Smoothers
The sample mean and the linear regression line are both special cases of linear smoothers, which are estimates of the regression function with the following form: ̂ r(x) =
i
yi ŵ (xi, x) (47)
The sample mean is the special case where ŵ (xi, x) = 1/n, regardless of what xi and x are. Ordinary linear regression is the special case where ŵ(xi, x) = (xi/ns^2 X )x. Both of these, as remarked, ignore how far xi is from x.
5.1 k-Nearest-Neighbor Regression
At the other extreme, we could do nearest-neighbor regression:
ŵ (xi, x) =
1 xi nearest neighbor of x 0 otherwise
This is very sensitive to the distance between xi and x. If r(x) does not change too rapidly, and X is pretty thoroughly sampled, then the nearest neighbor of x among the xi is probably close to x, so that r(xi) is probably close to r(x).
(^7) The standard test of whether this coefficient is zero is about as far from rejecting the null hypothesis as you will ever see, p = 0.95. You should remember this the next time you look at regression output.
0.0 0.2 0.4 0.6 0.8 1.
x
y
fit.all = lm(all.y~all.x) abline(fit.all)
Figure 4: Data from Figure 1, with a horizontal line at the mean (dotted) and the ordinary least squares regression line (solid). If you zoom in online you will see that there really are two lines there. (The abline adds a line to the current plot with intercept a and slope b; it’s set up to take the appropriate coefficients from the output of lm.
0.0 0.2 0.4 0.6 0.8 1.
x
y
library(knnflex) all.dist = knn.dist(c(all.x,seq(from=0,to=1,length.out=100))) all.nn1.predict = knn.predict(1:110,111:210,all.y,all.dist,k=1) abline(h=mean(all.y),lty=2) lines(seq(from=0,to=1,length.out=100),all.nn1.predict,col="blue") all.nn3.predict = knn.predict(1:110,111:210,all.y,all.dist,k=3) lines(seq(from=0,to=1,length.out=100),all.nn3.predict,col="red") all.nn5.predict = knn.predict(1:110,111:210,all.y,all.dist,k=5) lines(seq(from=0,to=1,length.out=100),all.nn5.predict,col="green") all.nn20.predict = knn.predict(1:110,111:210,all.y,all.dist,k=20) lines(seq(from=0,to=1,length.out=100),all.nn20.predict,col="purple")
Figure 5: Data points from Figure 1 with horizontal dashed line at the mean and the k-nearest-neighbor regression curves for k = 1 (blue), k = 3 (red), k = 5 (green) and k = 20 (purple). Note how increasing k smoothes out the regression line, and pulls it back towards the mean. (k = 100 would give us back the dashed horizontal line.)
x^2 K(0, x)dx < ∞
These conditions together (especially the last one) imply that K(xi, x) → 0 as |xi − x| → ∞. Two examples of such functions are the density of the Unif(−h/ 2 , h/2) distribution, and the density of the standard Gaussian N (0,
h) distribution. Here h can be any positive number, and is called the bandwidth. The Nadaraya-Watson estimate of the regression function is
̂ r(x) =
i
yi K(xi, x) ∑ j K(xj^ , x)^
i.e., in terms of Eq. 47,
ŵ(xi, x) = K(xi, x) ∑ j K(xj^ , x)^
(Notice that here, as in k-NN regression, the sum of the weights is always 1. Why?)^10 What does this achieve? Well, K(xi, x) is large if xi is close to x, so this will place a lot of weight on the training data points close to the point where we are trying to predict. More distant training points will have smaller weights, falling off towards zero. If we try to predict at a point x which is very far from any of the training data points, the value of K(xi, x) will be small for all xi, but it will typically be much, much smaller for all the xi which are not the nearest neighbor of x, so ŵ (xi, x) ≈ 1 for the nearest neighbor and ≈ 0 for all the others.^11 That is, far from the training data, our predictions will tend towards nearest neighbors, rather than going off to ±∞, as linear regression’s predictions do. Whether this is good or bad of course depends on the true r(x) — and how often we have to predict what will happen very far from the training data. Figure 6 shows our running example data, together with kernel regression estimates formed by combining the uniform-density, or box, and Gaussian ker- nels with different bandwidths. The box kernel simply takes a region of width h around the point x and averages the training data points it finds there. The Gaussian kernel gives reasonably large weights to points within h of x, smaller ones to points within 2h, tiny ones to points within 3h, and so on, shrinking like e−(x−xi)
(^2) / 2 h
. As promised, the bandwidth h controls the degree of smooth- ing. As h → ∞, we revert to taking the global mean. As h → 0, we tend to
(^10) What do we do if K(xi, x) is zero for some xi? Nothing; they just get zero weight in the average. What do we do if all the K(xi, x) are zero? Different people adopt different conventions; popular ones are to return the global, unweighted mean of the yi, to do some sort of interpolation from regions where the weights are defined, and to throw up our hands and refuse to make any predictions (computationally, return NA). (^11) Take a Gaussian kernel in one dimension, for instance, so K(xi, x) ∝ e−(xi−x)^2 / 2 h^2. Say
xi is the nearest neighbor, and |xi − x| = L, with L h. So K(xi, x) ∝ e−L (^2) / 2 h 2 , a small number. But now for any other xj , K(xi, x) ∝ e−L^2 /^2 h^2 e−(xj^ −xi)L/^2 h^2 e−(xj^ −xi)^2 /^2 h^2 e−L^2 /^2 h^2. — This assumes that we’re using a kernel like the Gaussian, which never quite goes to zero, unlike the box kernel.
0.0 0.2 0.4 0.6 0.8 1.
x
y
plot(all.x,all.y,xlab="x",ylab="y") rug(all.x,side=x,col="grey") rug(all.y,side=y,col="grey") lines(ksmooth(all.x, all.y, "box", bandwidth=2),col="blue") lines(ksmooth(all.x, all.y, "box", bandwidth=1),col="red") lines(ksmooth(all.x, all.y, "box", bandwidth=0.1),col="green") lines(ksmooth(all.x, all.y, "normal", bandwidth=2),col="blue",lty=2) lines(ksmooth(all.x, all.y, "normal", bandwidth=1),col="red",lty=2) lines(ksmooth(all.x, all.y, "normal", bandwidth=0.1),col="green",lty=2)
Figure 6: Data from Figure 1 together with kernel regression lines. Solid colored lines are box-kernel estimates, dashed colored lines Gaussian-kernel estimates. Blue, h = 2; red, h = 1; green, h = 0.5; purple, h = 0.1 (per the definition of bandwidth in the ksmooth function). Note the abrupt jump around x = 0.75 in the box-kernel/h = 0.1 (solid purple) line — with a small bandwidth the box kernel is unable to interpolate smoothly across the break in the training data, while the Gaussian kernel can.
Exercises
These are for you to think through, not to hand in.
- Suppose we use the mean absolute error instead of the mean squared error: MAE(a) = E [|Y − a|] (52) Is this also minimized by taking a = E [Y ]? If not, what value ˜r minimizes the MAE? Should we use MSE or MAE to measure error?
- Derive Eqs. 41 and 42 by minimizing Eq. 40.
- What does it mean for Gaussian kernel regression to approach nearest- neighbor regression as h → 0? Why does it do so? Is this true for all kinds of kernel regression?
