Model Selection and Shrinkage A Project in Linear Modeling

Model Selection and Shrinkage: A Project in Linear Modeling 1 Introduction Linear Modeling, in a broad sense, involves using information that we have to predict outcomes in the future. As one can gather, this is concept is widely applicable to disciplines across the board. Big data is an up and coming catch phrase in the media these days. With the rise of technology and the Internet, we have more data at our fingertips than ever before. The question then becomes: How can we most effectively use it to our advantage? How can we filter it so that it is managable and still useful to us? How can we profit from or make use of its patterns? Data is cheap and easy to find. Effective tools for analysis and individuals who can use these tools, on the other hand, are neither cheap nor easy to find. This is where linear modeling comes in. Statisticians are working to develop algorithms that will efficiently select, compile, and sort through these data sets to most efficiently flesh out pertinent information. Using model selection, statisticians are able to parse through data sets, emphasizing the most pertinent information and diminishing or ignoring information which is repetitive or unrelated. In this way, our data sets are both interpretable and more computationally feasible so that prediction can move with the click of a button. As will be discussed in this paper, the most popular linear modeling technique is the Ordinary Least Squares technique, which minimizes the sum of the squared errors of the model. This technique, however, is unreliable when there are extreme outliers or when there is collinearity between variables in the model. Statisticians have been developing new modeling techniques to remedy these issues. This paper will explore the machinery behind these techniques, study their asymptotic behaviors, and run simulations to analyze their performance in various modeling scenarios. 2 2.1 Linear Models Some Notation As will be discussed, there exist many different modeling techniques at the fingertips of statisticians, each with a different approach to the problem at hand and each with their own pros and cons. The goal of the statistican is to use the one which is best. An illustrative example will help us see what this could mean. Suppose we want to model systolic blood pressure against age. In other words, this model will predict an individual’s blood pressure given their age. We collect data from existing patients and plot it, as in Figure 5. The line through the data shows one example of a predictive linear fit of the data. Input an age and the curve will output an expected blood pressure. 1 Figure 1: Blood Pressure vs. Age One obvious criterion for the performance of a linear model is to sum up the squares of the residuals, shown in the figure to be the distance between the predictive curve and each individual data point. We will discuss this in more detail later. In the above example we collected data on one predictor (age) and one outcome (blood pressure). In general , we have p predictors and list their observations in the obervational matrix as below.     X=    1 x11 x21 1 x12 x22 . . . . . . . . . 1 x1n x2n . . . . . . . . . . . . . xp1 . xp2 . . . . . . . xpn         The matrix X is n × P where P = p + 1. The first column of ones accounts for an intecept term β0 . We assume that the outcomes yj are distributed as yj = β0 + β1 x1j + ... + βp xpj + j = XjT β + j (1) 1 ≤ j ≤ n. The vector β is the P × 1 vector β = (β0 , β1 , ..., βp )T and the  terms are normally distributed with a mean of 0 and a standard deviation of σ. We write:  ∼ N (0, σ). 2 All modeling techniques discussed attempt to estimate the value of β through some estimator β̂. The hat notation indicates that the quantity is an estimate of β. Given a set of n observations, we apply a modeling technique which generates an estimate β̂ of the true β. Because β̂ depends on a set of observations, its value will change as the set of observations changes. In this sense, the estimator β̂ is a random variable. Often times we are concerned with the bias of such an estimator which captures the difference between the expected value of an estimator and the true value of the parameter which it estimates. In general, we say that an estimator, X̂ is unbiased iff E(X̂) = X. This leads right in to the variance of an estimator, defined as the expected value of the squared sampling deviations, that is, var(X̂) = E[(X̂ − E(X̂))2 ] = σ 2 (X̂) This tells us the average squared distance an estimate X̂ is away from its expected value, capturing a measure of spread of the values of X̂. We also examine the covariance matrix of the predictors. The covariance between two variables X1 and X2 is measured by σ(X1 , X2 ) = E{[X1 − E(X1 )][X2 − E(X2 )]T }. Covariance measures the percentage of variation in one variable that is explained by variation in the other. Thus we define the covariance matrix of the set of variables by    X1 − E(X1 )  X2 − E(X2 )            . 2   X1 − E(X1 ) X2 − E(X2 ) ... Xp − E(Xp )  σ (X) = E     .       . Xp − E(Xp )  σ 2 (X1 ) σ(X1 , X2 ) ... σ(X1 , Xp )  σ(X2 , X1 ) σ 2 (X2 ) ... σ(X2 , Xp )   . . . σ 2 (X) =   . . .   . . . 2 σ(Xp , X1 ) σ(Xp , X2 ) ... σ (Xp )     .    Noting σ(X1 , X2 ) = σ(X2 , X1 ) we see that covariance matrices are symmetric. These matricies are helpful at identifying highly correlated variables. Note that the correlation of X1 and X2 is p 2 defined as corr(X1 , X2 ) = σ(X1 , X2 )/ σ (X1 )σ 2 (X2 ). Correlation matricies are defined similarily and serve to standardize covariance matricies by making the diagonal entries equal to 1. The correlation matrix captures pairwise linear correlations between the predictor variables. For a more intuitve sense of this, we can imagine that in predicting a level of job satisfaction, years of 3 experience and annual income would be highly correlated predictors: when one changes, the other tends to change as well. This idea of correlation between predictors will play a large role in the performance of modeling techniques. 2.2 Ordinary Least Squares The Ordinary Least Squares (OLS) technique is the oldest and most popular method by which we measure other modeling techniques. Many of the techniques discussed in the paper involve some type of modification to the least squares approach. We will study it now to establish its machinery and its performance so that we have a basis for comparsion. The OLS approach attempts to minimize the squares of the error terms i = yi − XiT β as such: β̂ OLS = min β n X (i )2 i=1 where yi is the predicted outcome for the ith set of predictor observations xi . Minimizing this quantity over all possible values of β will produce an optimal P × 1 vector β̂ OLS which satisfies β̂ OLS = (X T X)−1 X T y and the OLS model then becomes ŷ = X β̂ OLS where again, a hat is placed over the y to indicate that it is only an estimate of our true observations y. We can see that our estimates are dimensionally correct by noting that the matrix X is n × P and the vector ŷ is n × 1. Each element of β̂ represents the estimated coefficient of a corresponding element of the true β. To foreshadow a future discussion, we note that if the matrix X T X is not invertible then problems arise in our β̂ OLS estimates. Theorem 1. Suppose yj is distributed as in (1). Then, β̂ OLS is an unbiased estimator of β. Proof. By definition, β̂ OLS = (X T X)−1 X T (Xβ + ). Thus, E[β̂ OLS ] = E[(X T X)−1 X T (Xβ + )] = E[β] + E[(X T X)−1 X T ] = β + (X T X)−1 X T E[] = β + 0 = β. (2) By noting that E is a linear operator with E[] = 0 and that xi is not a random variable by assumption. This is reassuring. It tells us that, on average, the β̂ coefficients will correctly estimate the β coefficents in the true model. On the other hand, the variance of β̂ can be large. Since variance captures the spread of an estimate, we can say that high variance in the β̂ coefficents results 4 Figure 2: Sampling Variations in β̂ against correlation between predictors in imprecise estimates of β. So although OLS produces unbiased esimates it can also have high variance in the estimates, causing them to be unreliable. When does the OLS estimate have high variance? As we will see, it will be when the predictors are highly correlated. Earlier we noted that the β̂ OLS are highly sensitive to the invertibility of X T X. We review the following theorem from Lay [5] to find where the OLS technique is unreliable. Theorem 2. The matrix AT A is invertible iff the columns of A are linearly independent. [5] If we have highly correlated predictors then there exists some dependancy between the colums of X. In this sense, there is a spectrum of the invertibility of a matrix. This causes interesting dynamics in the β̂ values. To illustrate this we ran some simulations in R to generate histograms of estimator coefficients in two models. Both models had two predictors (p=2). The first model had highly correlated predictors corr(X1 , X2 ) = .99, while the second had lower correlation between predictors corr(X1 , X2 ) = .2. The histograms in Figure 2 above reveal a widely spread distribution of estimates in the correlated case and a more narrowly spread distribution in the other case. The ranges along the x-axis are constant to illustrate the idea. Thus, when the predictors in a linear model are highly correlated, the OLS approach can give imprecise estimates. The OLS estimates are good in the sense that they are unbiased estimators. Their pitfall, however, is that their variance is highly dependant on the correlations between predictors. The other modeling techniques described below find unique ways to address this problem. 5 Figure 3: Biased estimators in Ridge Regression [4] 2.3 Performance Criterion: MSE Good estimators should have low bias and low variance. To capture both of these, we look at the Mean Square Error (MSE) for an estimater X̂ or a parameter X , defined as E[(X̂ − X)2 ]. Theorem 3. The Mean Square Error of X̂ is equal to the sum of the variance of X̂ and the square of the bias of X̂. Mathematically, E[(X̂ − X)2 ] = (variance of X̂) + (bias of X̂)2 . Proof. E[(X̂ − X)2 ] = E[{(X̂ − E(X̂)) + (E(X̂) − X)}2 ] = E[(X̂ − E(X̂))2 ] + E[(E(X̂) − X)2 ] + 2E[(X̂ − E(X̂))(E(X̂) − X))] = (variance of X̂) + (bias of X̂)2 + 2E[X̂E(X̂) − E[X̂]2 − X̂X + E[X̂]X) = (variance of X̂) + (bias of X̂)2 + 2(E[X̂]2 − E[X̂]2 − XE[X̂] + XE[X̂]) = (variance of X̂) + (bias of X̂)2 (3) So, we often use MSE as the performance criterion when doing model selection. We saw that the bias of the OLS estimates is zero, but their variance can be high. Some modeling techniques introduce bias into their estimates in a way that lowers their variance and minimizes the MSE. See Figure 3 for an illustration of this idea. 6 2.4 Model Shrinkage Including more predictors in our model will almost always improve the fit of the model. This kitchen sink approach, however, has huge pitfalls including but not limited to multiple collinear relationships between predictors. This causes the optimal β̂ coefficients to be highly sensitive to small changes in X and largely unstable due to the matrix X T X being nearly singular (barely invertible). More recently we have developed more complex measures for predictive power. Ideally, we want our model to emphasize the stronger predictors of the outcome variable and diminish or eliminate the other predictors. In some way, shape, or form, this is exactly what the techniques below attempt to do. 2.5 Subset Selection Why do we want to screen out data which may be useful in the model? For one, large data sets are time consuming and computationally inefficient. Having a model with fewer independent variables is also beneficial because it is easier to analyze and interpret. Furthermore, highly correlated variables often add little predictive power, but do negatively affect the sampling variation of the β coefficients. For the reasons discussed, it is desirable to screen out independent variables. We first look for the ones which are either not fundamental to the problem, are subject to large errors, or effectively duplicate other variables. This process is called subset selection. Its goal, then, is to find the ‘best’ subset of independent variables and while the criteria to measure this may be objective and mechanical, the process of screening ought not be such. With a set of P predictors, there are 2P −1 possible subsets of independent variables. We run OLS analysis on each of these subsets to find their respective MSE values and pick the subset of predictors which minimizes the MSE. Subset selection often generates models that perform very well, but the process takes a lot of computing time and energy, especially as P gets large. So, while subset selection often outperforms the other techniques, it is not efficient for larger models. Forward Stepwise Regression is an offshoot of subset selection that begins with all of the predictors and strategically removes one at a time from the model until the remaining model has optimized its MSE. This is a more computationally feasible approach, even with large P. As Kutner [4] notes, subjective judgment should be used to ensure that important predictors are not being thrown away. When considering a potential predictor we look to see whether its β coefficient is statistically significant. That is, we seek predictors with high t-scores where t∗i = β̂i /s(β̂i ) recalling that s(β̂i ) is the sampling variation of β̂i . So, predictors with larger values tend to have high statistical significance. A pivotal predictor may be ignored because of a narrow range of values, resulting in a statistically insignificant t-score. This is where judgement need be employed. 7 2.6 Ridge Regression Ridge Regression introduces a penalty contraint on the size of the model. Specifically, Ridge Regression controls how large the sum of the squares of the estimated coefficients. That is, ridge finds β̂ R that minimizes n X (yj − XiT β)2 subject to j=1 p X (βj )2 ≤ t. j=1 The parameter t ≥ 0 is a tuning parameter which controls the amount of shrinkage applied to the model. It can be shown that the ridge regression estimates satisfy β̂iR = 1 β̂ OLS 1+γ i where γ depends on the choice of t as above. 2.7 Lasso and its Extentions Lasso controls the size of the sum of the absolute value of the estimated coefficients. That is, it finds β̂ L which minimizes n X (yi − XiT β̂)2 subject to i=1 p X |β̂j | ≤ t. j=1 Again t ≥ 0 is a tuning parameter which controls the amount of shrinkage applied to the model. Solving for the lasso coefficients gives: β̂iL = sgn(β̂iOLS )(|β̂iOLS | − γ)+ where γ > 0 depends on the choice of t above. The + notation means that β̂iL = 0 for |β̂iOLS | < γ We can see the differences between Lasso and Ridge in figure 4 above. Both produce estimates which scale down the OLS coefficients. But while ridge regression scales the OLS coefficients continuously, Lasso sets smaller OLS coefficients to zero and thereby helps shrink and interpret the model. Therefore we note that Lasso is a nice combination of subset selection and Ridge regression: it shrinks the model both by eliminating some predictors (similar to subset selection) and by scaling the coefficients on the remaining predictors (similar to ridge). Adaptive and Group Lasso Adaptive Lasso is an extension of Lasso which helps to account for differences in the magnitudes of estimated β̂ coefficients. Adaptive Lasso minimizes 8 Figure 4: Model Shrinkage across techniques of Linear Modeling. The dotted line represents the OLS coefficients. [9] p n X X |βjaL | ≤t (yi − XiT β aL )2 subject to OLS | | β̂ j i=1 j=1 which increases the penalty for not eliminating predictors with smaller OLS coefficients and reduces the penalty for not scaling predictors with larger OLS coefficients. In this way, adaptive lasso has more of the advantages that subset selection has. 2.8 Choosing Optimal Contraints AIC For Lasso, choose the optimal tuning parameter (t constraint) we let t change across a range of values and choose the value that minimizes the AIC criteria according to: 1 AIC = ln( ||Y − X β̂||2 ) + 2df /n. n [10] where df is equal to the number of non-zero coefficents in the model. We see that the first term measures the accuracy (bias) of the model while the second term captures the size (variance) of the model. Under this procedure, larger models are therefore penalized. In this way, minimizing the 9 Figure 5: K fold cross validation, K=5 [2] AIC statistic will drive more coefficients to zero in the selected model and will balance bias and variance. For ridge regression, the process is more involved. AIC minimization is inappropriate for Ridge Regression. As discussed, the AIC procedure penalizes models that include more nonzero coefficients. Because Lasso sets some coefficients to zero, as t varies, so too will the number of nonzero estimated coefficients (df ). Since AIC chooses the value of t that minimizes the AIC criterion, it will find a balance between subset selection and ridge regression. Ridge Regression, on the other hand, uniformly scales the OLS coefficients but sets none equal to zero. This means that df is constant across all values of t. Thus, the AIC will be minimized for the smallest value in the range of t. In this case the Ridge solutions will be driven towards the OLS solutions because these minimize the bias of the model. K-Fold Cross Validation Instead, we use a technique called cross validation to optimize our t for ridge regression. Cross validation partitions the sample data into K equally sized complementary subsets. We remove one subset and fit the model using the data remaining in the K − 1 subsets. Then we calculate the prediction error by using the fitted model to predict the data from the removed subset. For an illustation, Figure 5 above shows data divided into 5 groups. For the ith part (i = 1, 2, 3, 4, 5), we fit the model to the other 4 parts of the data and calculate the prediction error of the model when predicting the ith part of the data. We repeat the process across all values of i and combine the values of prediction errors. Then, we let the t constraint vary, repeating the above process, and choose the value of t which minimizes the total sum of the squared residual errors. This becomes our optimal contraint. 3 3.1 Asymptotics Motivation As shown in Figure 3, some models introduce bias into their estimators to lower variance and improve overall prediction accuracy. We know that the OLS estimates whose solutions minimize the sum of the squares of the residual errors produce unbiased estimators. By contraining the size 10 Figure 6: Sampling Variations in β̂ OLS and β̂ L of the model (as Lasso and Ridge Regression do), we introduce some bias into the estimates. To illustrate this, we plotted the histograms of some OLS and Lasso estimates of β4 = 2 of 100 simulations, each with 20 observations from a population y = Xβ + 3. We had β = (2, 2, 2, 2, 0, 0, 0, 0)T and  is the standard normal. There was a pairwise correlation between predictors xi and xj of ρ|i−j| with ρ = .9. The results are shown in Figure 6. Note that the OLS estimates are widely spread and roughly centered around 2. The Lasso estimates, on the other hand, are more narrowly spread but their center is closer to 0. This is a perfect recreation of Figure 3. A well known result is that Lasso estimates are consistent. Essentially, they are biased with finite sample size (n < ∞) but their bias approaches zero as n approaches infinity. En route to proving this, we will first show that the OLS estimates are asymptotically normal. 3.2 Some Notation First, we will extend our notion of a sequence of real numbers to a sequence of random variables. The motivation for doing so should become clear to the reader before proceeding. Much of statistics is focused on estimating parameters (such as the mean, variance, etc.) of a population of an unknown distribution. We do so by collecting random observations X1 , X2 , ..., Xn (called random variables) from the population distribution and then use these to generate an estimate θ̂n = θ̂(X1 , ..., Xn ) of a population parameter θ. It should be noted that the estimate θ̂n is itself a random variable due to its dependence on the set of random variables. 11 Specifically, we attempt to find the analogous concept of a limit for a sequence of random variables. For example, let Xn be the number of heads observed after flipping a fair coin n times. One may expect that Xn /n
1/2 as n
∞. Upon closer inspection however, we find some glaring pitfalls with this notation: it is possible, albeit unlikely, that Xn /n = 1 even for large n. We can, however say that Xn /n is likely to be close to 1/2. This motivates our definition of convergence in probablity and in distribution. Definition 1. A sequence of random variables Xn converges to a constant c in probability, denoted p Xn → c if and only if for each  > 0 P [|Xn − c| < ]
1 as n
∞. Equivalently, P (|Xn − c| ≥ )
0. As Lehmann [6] notes, this condition may be hard to verify. We use the following lemma to find an equivalent condition for convergence in probability. Lemma 1. Chebyshev Inequality. For any random variable Xn and any constants  > 0 and c, E[(Xn − c)2 ] ≥ 2 P [|Xn − c| ≥ ] The proof is omitted but we can now find an alternative condition for convergence in probablity. p Theorem 4. A sufficient condition for Xn → c is that E[(Xn − c)2 ]
0. Proof. Suppose E[(Xn − c)2 ]
0. Fix  > 0. By the previous lemma we find P [|Xn − c| ≥ ] ≤ 1/2 E[(Xn − c)2 ] Letting n
∞ and noting the choice of  is arbitrary gives the desired result. Since it is often much easier to show that E(Xn − c)2 converges to zero, this theorem is very useful. We will discuss one last related topic (using suggestive notation for our parameters) before moving on. Definition 2. A sequence of estimators X̂n of a parameter X is consistent if p X̂n → X. Theorem 5. A sufficient condition for X̂n to be a consistent estimator of X is that both the bias and the variance of X̂n tend to zero as n
∞. Proof. Theorem 3 along with Theorem 4 give the desired result. 12 Theorem 1 shows that Mean Square Error, MSE = E[(X̂n − X)2 ] can be decomposed into terms having to with both the bias and the variance of Xn . Theorems 2 and 3 show that the limiting behavior of the MSE of Xn determines the consistency of Xn . We introduce a few concepts to expand our asymptotic analysis of sequences before visiting the limiting distributions of the Lasso estimators. Definition 3. Two sequences {an } and {bn } are asymptotically equivalent, denoted an ∼ bn , if an /bn
1 as n
∞. In the case where bn approaches a finite limit b 6= 0, we can say that an approaches the same limit. In the case that the limit b is 0 or ±∞, we must use more caution: {n} and {n2 } both approach ∞. They are not, however, equivalent since their ratio does not tend to 1. We thus need to look at order relations and rates of convergence. Definition 4. A sequence {an } has order smaller than the sequence {bn }, denoted an = o(bn ) as n
∞ if an /bn
0. Lemma 2. A sequence {an } satisfies an = o(1) if and only if an
0. Definition 5. Two sequences {an } and {bn } are said to be of the same order, denoted an  bn if an /bn 6= 0 is finite. That is, there exists positive numbers m, M, and N such that m < |an /bn | < M, n > N. Combining these gives rise to big O notation: Definition 6. A sequence {an } has order smaller than or equal to the sequence {bn }, denoted an = O(bn ) as n
∞ if an /bn is bounded. 3.3 The Central Limit Theorem and the Delta Method The Central Limit Theorem is one of the most widely cited theorem in statistics. It involves the asymptotic distribution of sample means and will be central to the analysis of behaviors of linear models in the remainder of this paper. Theorem 6. Central Limit Theorem [6]. Let Xi , i = 1, 2, ..., n be independently and identically distributed from an unknown distribution with E(Xi ) = ζ and V ar(Xi ) = σ 2 < ∞. Then, √ d n(X̄ − ζ)/σ → N (0, 1) or equivalently, √ d n(X̄ − ζ) → N (0, σ 2 ). This theorem is incredibly powerful because no matter what population distribution the Xi come from, the distribution of their sample means converges in distribution to the normal distribution. Some other useful tools are Taylor’s Theorem and the Delta Method. 13 Theorem 7. Taylor’s Theorem [6]. Suppose that f (x) has r derivatives at a the point a. Then, f (a + ∆) = f (a) + ∆f 0 (a) + ... + ∆r (r) f (a) + o(∆r ). r! Theorem 8. Delta Method [6]. If √ d n(Tn − θ) → N (0, τ 2 ) then, √ d n[f (Tn ) − f (θ)] → N (0, τ 2 [f 0 (θ)]2 ). Proof. By Taylor’s Theorem, with a = θ and ∆ = Tn − θ, f (Tn ) = f (θ) + (Tn − θ)f 0 (θ) + op (Tn − θ) and hence √ n[f (Tn ) − f (θ)] = √ √ n(Tn − θ)f 0 (θ) + op [ n(Tn − θ)]. The first term on the right hand side tends in distribution to N (0, τ 2 [f 0 (θ)]2 ). On the other hand, √ we know that n(Tn − θ) is bounded in probability and hence, the remainder tends to zero in probability. The result follows from previous results. 3.4 Asympotic Normality of the LS Estimators The paper by Knight and Fu [3] notes the apparently well known fact that the Ordinary Least Squares estimator β̂n satisfies √ d n(β̂n − β) → N (0, σ 2 C −1 ) under the conditions ∞ Cn = 1X (xi xi T )
C n i=1 where C is a nonneagtive definite matrix and 1 max xi T xi
0. n 1≤i≤n We will prove this here for the one dimensional (p = 1) case and note that the multivariate case is simlply a generalization shown using Slutzky’s Theorem. We have shown that β̂ OLS = (X T X)−1 X T y. In summation notation this translates to β̂ = n 1 X n xi xi T n −1  1 X i=1 14 n i=1 xi yi  −1   x̄ ȳ   n n β̂0   X  X xi  x y = = i i    x̄ β̂1 n n  1 i=1 i=1 (for the simple case where p=1) which, after inverting the matrix on the left and solving for β̂1 can be shown to give n X (xi − x̄)(yi − ȳ) β̂1 = i=1 n X = (xi − x̄)2 n X dni yi dni = i=1 i=1 (xi − x̄) n X (xi − x̄)2 i=1 We invoke the following theorem from Lehmann [6] to arrive at our desired result. Theorem 9 [6]. Let Y1 , ..., Yn be independently and identically distributed with E(Yi ) = 0 and V ar(Yi ) = σ 2 > 0. Then n X dni Yi d i=1 v → N (0, 1). uX u n σt dni 2 i=1 If we normalize our yi to be Yi = yi − ȳ then, since our β̂ is an unbiased estimator of β, the conditions of the above theorem are met and we find that: n X i=1 v uX u n t (xi − x̄)2  dni (yi − ȳ) i=1 n X  σ dni Yi i=1 d = v → N (0, 1). uX u n σt dni 2 i=1 Which is the univariate case of the result given in Knight [3]. 3.5 Asymptotic Consistency of Lasso Estimates and its Limiting Distribution Here, we again assume that the matrix C, defined above as ∞ Cn = 1X (xi xi T )
C n i=1 15 is nonsingular. We consider the asymptotic behavior of Lasso’s objective function by first defining a random variable Zn (β̂): p n 1X λn X T 2 Zn (β̂) = (yi − Xi β̂) + |β̂j | n n i=1 j=1 which, as we have seen, is minimized at β̂ = β̂nL . The following theorem tells us that β̂nL is a consistent estimator of β provided that λn = o(n) p Theorem 10 [3]. If C is nonsingular and λn /n
λ0 ≥ 0 then β̂nL → argmin(Z) where T Z(β̂) = (β̂ − β) C(β̂ − β) + λ0 p X |β̂j | j=1 Thus, if λn = o(n) then λ0 = 0 and argmin(Z) = β so that β̂nL is consistent. Proof. We will show that Zn (β̂), defined above, converges in probability to Z(β̂) + σ 2 . The result will follow by applying established previous results from Pollard [7]. To show convergence of Zn we translate it into matrix notation. Zn (β̂) = p n 1X λn X (yi − XiT β̂)2 + |β̂j | n n i=1 j=1 p λn X 1 T |β̂j | = (Y − X β̂) (Y − X β̂) + n n j=1 p 1 λn X T = (( + Xβ) − X β̂) (( + Xβ) − X β̂) + |β̂j | n n j=1 = 1 λn ( + X(β − β̂))T ( + X(β − β̂)) + n n p X |β̂j | j=1 p 1 T λn X T T = [  + 2X(β − β̂) + (β − β̂) X X(β − β̂)] + |β̂j | n n j=1 Now, we let n
∞ and note the following facts: σ 2 Inxn = var() = E(T ) + E()T E() = E(T ). 1 t p   → E[T ] = σ 2 Inxn (by the law of large numbers) n 16 1 T X X
C n 1 1 2 T X(β − β̂)
2E[ n n n n i=1 i=1 1X 1X i ] = E[i ] = 0 n n E[ n X i ]X T (β − β̂) = 0 (by the law of large numbers) i=1 λn
λo = 0 n so that p 2 T Zn (β̂) → σ + (β − β̂) C(β − β̂) + λ0 p X |β̂j | = Z(β̂) + σ 2 j=1 The pointwise convergence of Zn (β̂) to Z(β̂) + σ 2 allows us to conclude that p sup |Zn (β̂) − Z(β̂) − σ 2 | → 0 β̂∈K for any compact set K by Pollard [7] and that β̂nL = Op (1). It follows that p argmin(Zn ) → argmin(Z). [3] This tells us that the lasso estimates are consistent. To find their limiting distribution, we invoke the following theorem. √ Theorem 11 [3]. If λn / n
λ0 ≥ 0 and C is nonsingular then √ d n(β̂nL − β) → argmin(V ) where T T V (u) = u Cu − 2u W + λ0 p X j=1 and W ∼ N (0, σ 2 C). [3] 17 uj sgn(β̂j )     p p X X √ √ 2 Proof. Define Vn (u) = || − Xu/ n|| + λn |βj + uj / n| − ||||2 + λn |βj | and note √ that Vn (u) = Qn (β + u/ n) − Qn (β) where j=1 j=1 2 Qn (z) = ||Y − Xz|| + λn p X |zj | j=1 √ We find that argmin Vn = n(β̂n − β) by noting that √ Vn ( n(β̂n − β)) = Qn (β + (β̂n − β)) − Qn (β) = Qn (β̂n ) + C and that the objective function Qn is analogous to the Lasso objective function, minimized at β̂n . To show convergence of Vn , note that √ p √ X √ √ Qn (β + u/ n) − Qn (β) = |(Y − Xβ) − Xu/ n| |(Y − Xβ) − Xu/ n| + λn |βj + uj / n| − Qn (β) T j=1 √ = ||||2 − 2uT / nX T + uT X T Xu/n + λn p X j=1 = uT XT X n XT  u − 2uT √ n + λn p X p X √ |βj + uj / n| − (||||2 + λn |βj |) j=1 √ [|βj + uj / n| − |βj |] j=1 Then, letting n
∞ we apply the following facts: 1 T X X
C n XT  d √ →W n λ √n
λo n d So that Vn (u) → V (u), as defined above. 4 4.1 Simulations in R The Regression Scenario For the following simulations in R, we are in the usual regression situation. We have collected n observational data points from a set of p predictors ~x = (x1 , x2 , ..., xp )T along with from a 18 response/outcome variable y. In this case we will model the data without the intercept term so we are left with the augmented data matrix:     X| y =     x11 x21 x12 x22 . . . . . . x1n x2n . . . . . . . . . . . . . xp1 y1 . xp2 y2 . . . . . . . . . . xpn yn         In the following simulations, we compare the Ordinary Least Squares Estimates with the estimates of Lasso, Ridge Regression, Adaptive Lasso, and Forward Stepwise Regression. To optimize the tuning parameters we used five fold cross validation in ridge regression and AIC minimization in Lasso and its extensions. We also used cross validation in Lasso for comparative purposes. The steps procedure in R was used to perform forward stepwise regression. These examples are consistent with those published in Tibshirani’s paper Regression Shrinkage and Selection via the Lasso [9] in hopes of recreating his results. 4.2 Example 1 This example simulated 25 data sets, each with 20 observations from the population modeled as: y = Xβ + σ where β = (3, 1.5, 0, 0, 2, 0, 0, 0)T ,  is the standard normal, and σ = 3. There was a pairwise correlation between predictors xi and xj of ρ|i−j| with ρ = .5 Table 1 shows the median Mean Square Errors from each of the techniques. We use the median because it is a robust to outliers. We see that Lasso and Ridge Regression perform best. Table 2 shows the proportion of models from Lasso and forward stepwise regressions which include each β coefficient. We see that Ridge and Lasso perform comparably well and both outperform stepwise and OLS by a good measure. Table 2 shows that Lasso retains the correct coefficients more often than does forward stepwise. On average forward stepwise has closer to the correct number of zero coefficients (5) but suffers from high variablity in the nonzeo coefficients. Figure 7 below shows the inter quartile range of β estimates of the four linear models over the 25 data sets. We see that Lasso and forward stepwise are able to drive the zero β coefficents to zero very successfully. Table 1: Results for Example 1 Method Median MSE Avg. no. OLS 5.27 Ridge Regression 1.43 Lasso (AIC) 1.77 Lasso (cross-validation) 1.40 Adaptive Lasso 2.18 Forward Stepwise 2.45 19 of 0 coefficients 0 0 3.26 3.44 3.5 3.62 20 Figure 7: Estimates of the 8 coefficients in example 1. 1=OLS 2=Ridge 3=Lasso 4=Fwd Stepwise Table 2: Proportion of Models including β coefficients Coefficient Lasso Forward Stepwise β1 =3 1.0 .98 β2 =1.5 .82 .78 β3 =0 .22 .34 β4 =0 .22 .26 β5 =2 .98 .88 β6 =0 .16 .30 β7 =0 .22 .26 β8 =0 .24 .18 4.3 Example 2: Highly Correlated Predictors This example is similar to example 1 except we set ρ = .9 to increase the correlation between predictors. Table 3 shows the results and we see that ridge regression, Lasso and its extensions are able to handle the correlation issue by significantly outperforming the OLS estimates. Again we see that forward stepwise, on average, chooses the correct number of zero coefficients. Here, the variance problem in the nonzero coefficients in magnified, as its median mean square error is well above that of lasso and ridge regression. Table 3: Results for Example 2 Method Median MSE Avg. no. OLS 4.35 Ridge .27 Lasso (AIC) .64 Lasso (cross-validation) .75 Adaptive Lasso 1.92 Forward Stepwise 2.99 4.4 of 0 coefficients 0 0 5.76 5.2 4.92 5.08 Example 3: Sparse Signals Here we let β = (5, 0, 0, 0, 0, 0, 0, 0)T and σ = 2. As before, there are eight coefficients to estimate but only the first one is non-zero. In this sense the signals of β are sparse. This set up should be well suited for Lasso and its extensions. The results in Table 4 confirm this: Lasso’s median mean square error is unprecedentedly low. 4.5 Example 4: Smaller Signals Here we had the same set-up as in example 1 except βj = .85 for all 1 ≤ j ≤ 8. That is, the signals of β are equally weighted but relativiely smaller in magnitude. The results in Table 5 reveals that Ridge Regression outperforms the others. Lasso comes in second. OLS actually outperforms 21 Table 4: Results for Example 3 Method Median MSE Avg. no. OLS 4.71 Ridge 1.07 Lasso (AIC) .26 Lasso (cross-validation) .018 Adaptive Lasso 1.18 Forward Stepwise 2.41 of 0 coefficients 0 0 5.44 5.68 5.12 5.12 forward stepwise selection in this case. This could be due to the fact that model selection is useless in this case. Since all β elements are equal, it is impossible to forward stepwise to choose the strongest signals. Table 5: Results for Example 4 Method Median MSE Avg. no. of 0 coefficients OLS 6.6 0 Ridge 1.7 0 Lasso (AIC) 3.7 2.4 Lasso (cross-validation) 3.4 2.16 Adaptive Lasso 5.89 3.64 Forward Stepwise 6.93 4 4.6 Example 5: Larger Models Here, we are set up to analyze the performance of the estimators in larger models to better reflect a real world application. We simulated 25 data sets, each with 100 observations and 40 predictors. The coefficient vector was β = (0, 0, ..., 0, 2, 2, ..., 2, 0, 0, ..., 0, 2, 2, ..., 2)T and σ = 15. As in example 1, we had ρ = .5. Table 6 shows the performance results. We see that ridge outperforms the others but is doing no model selection. Lasso under cross validation has a reasonably low MSE and does a great job at model selection. Table 6: Results for Example 5 Method Median MSE Avg. no. OLS 164.99 Ridge 10.81 Lasso (AIC) 111.83 Lasso (cross-validation) 43.45 Adaptive Lasso 95.91 Forward Stepwise 305.22 22 of 0 coefficients 0 0 33.08 21.2 20.72 27.4 5 Conclusion We began our discussion by studying the Ordinary Least Squares technique. After some sensitivity analysis where we increased the correlation between predictors, we found that the OLS estimates are unreliable estimators when there is collinearity between predictors. We then touched on alternative modeling techniques that attempt to improve upon OLS and remedy its shortcomings. By constraining the size of the model, we introduce bias into our estimates in an attempt to lower the overall mean square error of the model. By applying the Central Limit Theorem and Taylor’s Theorem, we found that the distribution of OLS estimates are centered around the true parameter value and approach normality as n approaches infitity. After some legwork, we also found that the Lasso estimates are asymptotically consistent. That is, as n approaches infinity, the bias in the Lasso model approaches zero. Further, the distribution of the Lasso estimates is asymptotically normal. This is reassuring because with a lagre sample size, we are confident that the bias of the Lasso estimates is negligible. To test performance of the various models with finite sample size we ran simulations in R. For the most part, our mathematical intuitions were correct. In small models with sparse signals, Lasso outperformed all other techniques by a significant margin. Ridge was next followed by forward stepwise. In models with more signals and larger signals, Here Lasso and Ridge performed well (especially with highly correlated predictors) while forward stepwise and OLS did poorly. Finally, in models with a large number of smaller signals (no zero coefficients), Ridge Regression does best by far. Lasso still outperforms both OLS and forward stepwise. Overall, when the true model is unknown, we can conclude that Lasso will do the best job at balancing model selection with prediction accuracy. By using improved linear modeling techniques such as subset selection, ridge regression and Lasso, we are able to efficiently parse though large data sets (full of redundant and/or useless information) and accurately predict outcomes. These tools will continue to be improved upon and applied to larger and more complicated models. 23 References [1] Gleser et. al. The Limiting Distribution of Least Squares in Errors-in-Variables Regression Model , The Annals of Statistics, Mar. 1987. [2] Hastie, et. al. The Elements of Statistical Learning, Springer, 2009. [3] Knight, K. and Fu, W. Asympotics for Lasso-Type Estimators, University of Toronto, 2000. [4] Kutner, et. al. Applied Linear Statistical Models, 5th ed., McGraw-Hill, 2004. [5] Lay, David C. Linear Algebra and its Applications Pearson, 2006. [6] Lehmann, E.L. Elements of Large-Sample Theory, Springer, 1999. [7] Pollard, David Asymptotics for Least Absolute Deviation Regression Estimators, Yale University, 1991. [8] Srivastava, M. S. On Fixed-width Confidence Bounds for Regression Parameters and Mean Vector, Princeton University, 1967. [9] Tibshirani, Robert. Regression Shrinkage and Selection via the Lasso, University of Toronto, 1994. [10] Wang, H and Leng, C. A note on adaptive group lasso, Peking University, 2008. 24