Multiple Linear Regression

Section 15 Multiple linear regression. Let us consider a model Yi = �1 Xi1 + . . . + �p Xip + χi where random noise variables χ1 , . . . , χn are i.i.d. N (0, π 2 ). We can write this in a matrix form Y = X� + χ, where Y and χ are n × 1 vectors, � is p × 1 vector and X is n × p matrix. We will denote the columns of matrix X by X1 , . . . Xp , i.e. X = (X1 , . . . , Xp ) and we will assume that these columns are linearly independent. If they are not linearly independent, we can not reconstruct parameters � from X and Y even if there is no noise χ. In simple linear regression this would correspond to all Xs being equal and we can not estimate a line from observations only at one point. So from now on we will assume that n > p and the rank of matrix X is equal to p. To estimate unknown parameters � and π we will use maximum likelihood estimators. Lemma 1. The MLE of � and π 2 are given by: �ˆ = (X T X)−1 X T Y and π̂ 2 = 1 1 |Y − X�ˆ|2 = |Y − X(X T X)−1 X T Y |2 . n n Proof. The p.d.f. of Yi is fi (x) = ∼ ⎟ 1 ⎠ 1 exp − 2 (x − �1 Xi1 − . . . − �p Xip )2 2π 2απ and, therefore, the likelihood function is n � n ⎠ ⎟ 1 � 1 ⎠n ∼ exp − 2 fi (Yi ) = (Yi − �1 Xi1 − . . . − �p Xip )2 2π i=1 2απ i=1 ⎟ 1 ⎟ 1 ⎠n ⎠ 2 ∼ exp − 2 |Y − X�| . = 2π 2απ ⎟ To maximize the likelihood function, first, we need to minimize |Y − X�|2 . If we rewrite the norm squared using scalar product: |Y − X�| 2 = (Y − p � i=1 � i Xi , Y − = (Y, Y ) − 2 p � p � � i Xi ) i=1 �i (Y, Xi ) + i=1 p � �i �j (Xi , Xj ). i,j=1 Then setting the derivatives in each �i equal to zero −2(Y, Xi ) + 2 we get (Y, Xi ) = p � i=1 p � �j (Xi , Xj ) = 0 j=1 �j (Xi , Xj ) for all i � p. In matrix notations this can be written as X T Y = X T X�. Matrix X T X is a p × p matrix. Is is invertible since by assumption X has rank p. So we can solve for � to get the MLE �ˆ = (X T X)−1 X T Y. It is now easy to minimize over π to get π̂ 2 = 1 1 |Y − X�ˆ|2 = |Y − X(X T X)−1 X T Y |2 . n n To do statistical inference we need to compute the joint distribution of these estimates. We will prove the following. Theorem. We have ⎟ ⎠ nπ̂ 2 �ˆ � N �, π 2 (X T X)−1 , � �2n−p π2 and estimates �ˆ and π̂ 2 are independent. Proof. First of all, let us rewrite the estimates in terms of random noise χ using Y = X� + χ. We have �ˆ = (X T X)−1 X T Y = (X T X)−1 X T (X� + χ) = (X T X)−1 (X T X)� + (X T X)−1 X T χ = � + (X T X)−1 X T χ and since Y − X(X T X)−1 X T Y = X� + χ − X(X T X)−1 X T (X� + χ) = X� + χ − X� − X(X T X)−1 X T χ = (I − X(X T X)−1 X T )χ we have π̂ 2 = 1 |(I − X(X T X)−1 X T )χ|2 . n Since �ˆ is a linear transformation of a normal vector χ it will also be normal with mean E�ˆ = E(� + (X T X)−1 X T χ) = � and covariance matrix E(�ˆ − �)(�ˆ − �)T = = = = E(X T X)−1 X T χχT X(X T X)−1 (X T X)−1 X T EχχT X(X T X)−1 (X T X)−1 X T (π 2 I)X(X T X)−1 π 2 (X T X)−1 (X T X)(X T X)−1 = π 2 (X T X)−1 . This proves that �ˆ � N (�, π 2 (X T X)−1 ). To prove that �ˆ and π̂ 2 are independent and to find the distribution of nπ̂ 2 /π 2 we will use the following trick. This trick can also be very useful computationally since it will relate all quantities of interest expressed in terms of n×p matrix X to quantities expressed in terms of a certain p × p matrix R which can be helpful when n is very large compared to p. We would like to manipulate the columns of matrix X to make them orthogonal to each other, which can be done by Gram-Schmidt orthogonalization. In other words, we want to represent matrix X as X = X0 R where X0 is n × p matrix with columns X01 , . . . , X0p that are orthogonal to each other and, moreover, form an orthonormal basis, and matrix R is p × p invertible (and upper triangular) matrix. In Matlab this can be done using economy size QR factorization [X0, R]=qr(X,0). The fact that columns of X0 are orthonormal implies that X0T X0 = I - a p × p identity matrix. Let us replace X by X0 R everywhere in the estimates. We have (X T X)−1 X T = (RT X0T X0 R)−1 RT X0T = (RT R)−1 RT X0T = R−1 (RT )−1 RT = R−1 X0T , X(X T X)−1 X T = X0 R(RT X02 X0 R)−1 RT X0T = X0 RR−1 (RT )−1 RT X0T = X0 X0T . As a result �ˆ − � = R−1 X0T χ and nπ̂ 2 = |(I − X0 X0T )χ|2 . (15.0.1) By construction p columns of X0 , which are also the rows of X0T , are orthonormal. Therefore, we can choose the last n − p rows of a n × n matrix ⎛ T � X0 A= ··· to make A an orthogonal matrix, we just need to choose them to complete, together with rows of X0T , the orthonormal basis in Rn . Let us define a vector � � ⎞ ⎞ g1 χ1 � g2 ⎜ ⎛ T � � χ2 ⎜ X0 � � ⎜ ⎜ g = Aχ, i.e. � .. ⎜ = � .. ⎜ . · · · � . ⎝ � . ⎝ gn χn Since χ is a vector of i.i.d. standard normal, we proved before that its orthogonal transfor­ mation g will also be a vector of independent N (0, π 2 ) random variables g1 , . . . , gn . First of all, since � � ⎞ ⎞ g1 χ1 ⎜ ⎜ � � ĝ := � ... ⎝ = X0T � ... ⎝ , gp χn we have � g1 � ⎜ �ˆ − � = R−1 X0T χ = R−1 � ... ⎝ = R−1 ĝ. gp (15.0.2) |(I − X0 X0T )χ|2 = gp2+1 + . . . + gn2 . (15.0.3) ⎞ Next, we will prove that First of all, orthogonal transformation preserves lengths, so |g |2 = |Aχ|2 = |χ|2 . On the other hand, let us write |χ|2 = χT χ and break χ into a sum of two terms χ = X0 X0T χ + (I − X0 X0T )χ. Then we get ⎟ ⎠⎟ ⎠ |g|2 = |χ|2 = χT χ = χT X0 X0T + χT (I − X0 X0T ) X0 X0T χ + (I − X0 X0T )χ . When we multiply all the terms out we will use that X0T X0 = I since the matrix X0T X0 consists of scalar products of columns of X0 which are orthonormal. This also implies that X0 X0T (I − X0 X0T ) = X0 X0T − X0 IX0T = 0. Using this we get |g |2 = |χ|2 = χT X0 X0T χ + χT (I − X0 X0T )(I − X0 X0T )χ = |X0T χ|2 + |(I − X0 X0T )χ|2 = |ĝ|2 + |(I − X0 X0T )χ|2 because ĝ = X0T χ so we finally proved that |(I − X0 X0T )χ|2 = |g|2 − |ĝ|2 = g12 + . . . + gn2 − g12 − . . . − gp2 = gp2+1 + . . . + gn2 which is (15.0.3). This proves that nπ̂ 2 /π 2 � �2n−p and it is also independent of �ˆ which depends only on g1 , . . . , gp by (15.0.2). Let us for convenience write down equation (15.0.2) as a separate result. Lemma 2. Given a decomposition X = X0 R with n × p matrix X0 with orthonormal columns and invertible (upper triangular) p × p matrix R we can represent ⎞ � g1 � ⎜ �ˆ − � = R−1 ĝ = R−1 � ... ⎝ gp for independent N (0, π 2 ) random variables g1 , . . . , gp . Confidence intervals and t-tests for linear combination of parameters �. Let us consider a linear combination c 1 �1 + . . . + c p �p = c T � where c = (c1 , . . . , cp )T . To construct confidence intervals and t-tests for this linear combinaˆ Clearly, it has a normal distribution with tion we need to write down a distribution of cT �. T ˆ T mean Ec � = c � and variance E(cT (�ˆ − �))2 = EcT (�ˆ − �)(�ˆ − �)T c = cT Cov(�ˆ)c = π 2 cT (X T X)−1 c. Therefore, cT (�ˆ − �) � � N(0, 1) π 2 cT (X T X)−1 c and using that nπ̂ 2 /π 2 � �2n−p we get � c (�ˆ − �) � π 2 cT (X T X)−1 c T � 1 nπ̂ 2 = cT (�ˆ − �) 2 n−p π � n−p nπ̂ 2 cT (X T X)−1 c � tn−p . To obtain the distribution of one parameter �ˆi we need to choose a vector c that has all zeros and 1 in the ith coordinate. Then we get � n−p ˆ (�i − �i ) � tn−p . 2 nπ̂ ((X T X)−1 )ii Here ((X T X)−1 )ii is the ith diagonal element of the matrix (X T X)−1 . This is a good time to mention how the quality of estimation of � depends on the choice of X. For example, we mentioned before that the columns of X should be linearly independent. What happens if some of them are nearly collinear? Then some eigenvalues of (X T X) will be ’small’ (in some sense) and some eigenvalues of (X T X)−1 will be ’large’. (Small and large here are relative terms because the size of the matrix also grows with n.) As a result, the confidence intervals for some parameters will get very large too which means that their estimates are not very accurate. To improve the quality of estimation we need to avoid using collinear predictors. We will see this in the example below. Joint confidence set for � and F -test. By Lemma 2, R(�ˆ − �) = ĝ and, therefore, g12 + . . . + gp2 = |ĝ|2 = ĝ T ĝ = (�ˆ − �)T RT R(�ˆ − �) = (�ˆ − �)T X T X(�ˆ − �). Since gi � N (0, π 2 ) this proves that (�ˆ − �)T X T X(�ˆ − �) � �2p . π2 Using that nπ̂ 2 /π 2 � �2n−p gives (�ˆ − �)T X T X(�ˆ − �) � nπ̂ 2 (n − p) ˆ = (� − �)T X T X(�ˆ − �) � Fp,n−p . 2 2 pπ (n − p)π npˆ π2 If we take c such that Fp,n−p (0, c� ) = � then (n − p) ˆ (� − �)T X T X(�ˆ − �) � c� npπ̂ 2 (15.0.4) defines a joint confidence set for all parameters � simultaneously with confidence level �. Suppose that we want to test a hypothesis about all parameters simultaneously, for example, H0 : � = �0 . Then we consider a statistic F = (n − p) ˆ (� − �0 )T X T X(�ˆ − �0 ), 2 npπ̂ (15.0.5) which under null hypothesis has Fp,n−p distribution, and define a decision rule by � H0 : F � c β= H1 : F > c, where a threshold c is determined by Fp,n−p (c, �) = � - a level of significance. Of course, this test is equivalent to checking if vector �0 belongs to a confidence set (15.0.4)! (We just need to remember that confidence level = 1 - level of significance.) Simultaneous confidence set and F -test for subsets of �. Let s = {i1 , . . . , ik } ≤ {1, . . . , p} be s subset of size k � p of indices {1, . . . , p} and let �s = (�i1 , . . . , �ik )T be a vector that consists of the corresponding subset of parameters �. Suppose that we would like to test the hypothesis H0 : �s = �s0 for some given vector �s0 , for example, �s0 = 0. Let �ˆs be a corresponding vector of estimates. Let ⎟ ⎠ −1 �s = (X T X)i,j i,j�s be a k × k submatrix of (X T X)−1 with row and column indices in the set s. By the above Theorem, the joint distribution of �ˆs is �ˆs � N(�s , π 2 �s ). 1/2 Let A = �s , i.e. A is a symmetric k × k matrix such that �s = AAT . As a result, a centered vector of estimates can be represented as �ˆs − �s = Ag, where g = (g1 , . . . , gk )T are independent N(0, π 2 ). Therefore, g = A−1 (�ˆs − �s ) and the rest is similar to the above argument. Namely, g12 + . . . + gk2 = |g|2 = g T g = (�ˆs − �s )T (A−1 )T A−1 (�ˆs − �s ) 2 2 ˆ = (�ˆs − �s )T (AAT )−1 (�ˆs − �s ) = (�ˆs − �s )T �−1 s (�s − �s ) � π �k . As before we get (n − p) ˆ ˆ (�s − �s )T �−1 s (�s − �s ) � Fk,n−p nkˆ π2 and we can now construct a simultaneous confidence set and F -tests. F = Remark. Matlab regression function ’regress’ assumes that a matrix X of explanatory variables will contain a first column of ones that corresponds to an ”intercept” parameter � 1 . The F -statistic output by ’regress’ corresponds to F -test about all other ”slope” parameters: H0 : �2 = . . . = �p = 0. In this case s = {2, 3, . . . , p}, k = p − 1 and F = (n − p) ˆT −1 ˆ � � �s � Fp−1,n−p . n(p − 1)π̂ 2 s s Example. Let us take a look at the ’cigarette’ dataset from previous lecture. We saw that tar, nicotine and carbon monoxide content are positively correlated and any pair is well described by a simple linear regression. Suppose that we would like to predict carbon monoxide as a linear function of both tar and nicotine content. We create a 25 × 3 matrix X: X=[ones(25,1),tar,nic]; We introduce a first column of ones to allow an intercept parameter �1 in our multiple linear regression model: COi = �1 + �2 Tari + �3 Nicotini + χi . If we perform a multiple linear regression: [b,bint,r,rint,stats] = regress(carb,X); We get the estimates of parameters and 95% confidence intervals for each parameter b = 3.0896 0.9625 -2.6463 bint = 1.3397 0.4717 -10.5004 4.8395 1.4533 5.2079 and, in order, R2 -statistic, F -statistic from (15.0.5), p-value for this statistic Fp,n−p (F, +�) = F3,25−3 (F, +�) and the estimate of variance π̂ 2 : stats = 0.9186 124.1102 0.000 1.9952. First of all, we see that high R2 means that linear model explain most of the variability in the data and small p-value means that we reject the hypothesis that all parameters are equal to zero. On the other hand, simple linear regression showed that carbon monoxide had a positive correlation with nicotine and now we got �ˆ3 = −2.6463. Also, notice that the confidence interval for �3 is very poor. The reason for this is that tar and nicotine are nearly collinear. Because of this the matrix ⎞ � 0.3568 0.0416 −0.9408 0.0281 −0.4387 ⎝ (X T X)−1 = � 0.0416 −0.9408 −0.4387 7.1886 has relatively large last diagonal diagonal value. We recall that Theorem gives that the variance of estimate �ˆ3 is 7.1886π 2 and we also see that the estimate of π 2 is π̂ 2 = 1.9952. As a result the confidence interval for �3 is rather poor. Of course, looking at linear combinations of tar and nicotine as new predictors does not make sense because they lose their meaning, but for the sake of illustrations let us see what would happen if our predictors were not nearly collinear but, in fact, orthonormal. Let us use economic QR decomposition [X0,R]=qr(X,0) a new matrix of predictor X0 with orthonormal columns that are some linear combinations of tar and nicotine. Then regressing carbon monoxide on these new predictors [b,bint,r,rint,stats] = regress(carb,X0); we would get b = -62.6400 -22.2324 0.9870 bint = -65.5694 -25.1618 -1.9424 -59.7106 -19.3030 3.9164 all confidence intervals of the same relatively better size. Example. The following data presents per capita income of 20 countries for 1960s. Also presented are the percentages of labor force employed in agriculture, industry and service for each country. (Data source: lib.stat.cmu.edu/DASL/Datafiles/oecdat.html) COUNTRY PCINC AGR IND SER CANADA SWEEDEN SWITZERLAND LUXEMBOURG U. KINGDOM DENMARK W. GERMANY FRANCE BELGUIM NORWAY ICELAND NETHERLANDS AUSTRIA IRELAND ITALY JAPAN GREECE SPAIN PORTUGAL TURKEY 1536 1644 1361 1242 1105 1049 1035 1013 1005 977 839 810 681 529 504 344 324 290 238 177 13 14 11 15 4 18 15 20 6 20 25 11 23 36 27 33 56 42 44 79 43 53 56 51 56 45 60 44 52 49 47 49 47 30 46 35 24 37 33 12 45 33 33 34 40 37 25 36 42 32 29 40 30 34 28 32 20 21 23 9 We can perform simple linear regression of income on each of the other explanatory variables or multiple linear regression on any pair of the explanatory variables. Fitting simple linear regression of income vs. percent of labor force in agriculture, industry and service: polytool(agr,income,1), etc., produces figure 15.1. Next, we perform statistical inference using ’regress’ function. Statistical analysis of linear regression fit of income vs. percent of labor force in agriculture: [b,bint,r,rint,stats]=regress(income,[ones(20,1),agr]) b = 1317.9 -18.9 stats = bint = 1094.7 -26.0 0.6315 30.8472 1541.1 -11.7 2.8e-005 74596. 2500 2000 1500 1500 1000 1000 500 500 0 0 −500 −500 −1000 10 20 30 40 50 60 70 80 40 45 15 20 25 30 35 40 45 50 55 60 2500 2000 1500 1000 500 0 −500 −1000 10 15 20 25 30 35 Figure 15.1: Linear regression of income on percent of labor force in agriculture, industry and service. For income vs. percent of labor force in industry [b,bint,r,rint,stats]=regress(income,[ones(20,1),ind]); b = -359.3115 27.4905 bint = -907.1807 15.3058 stats = 0.5552 22.4677 0.0002 188.5577 39.6751 90042 and for income vs. labor force in service [b,bint,r,rint,stats]=regress(income,[ones(20,1),serv]); b = -264.5199 35.3024 bint = -858.0257 16.8955 328.9858 53.7093 stats = 0.4742 16.2355 0.0008 106430. We see that in all three cases, the hypotheses that parameters of least-squares line are both zero can be rejected at conventional level of significance � = 0.05. Looking at the confidence intervals for the estimates of slopes we observe that the correlation of income with percent of labor force in agriculture is negative, and other two correlations are positive. We can also perform a multiple regression on any two explanatory variables. We can not perform multiple linear regression with all three explanatory variables because they add up to 100%, i.e. they are linearly dependent. If we create a predictor matrix X=[ones(20,1),agr,ind]; and perform multiple linear regression [b,bint,r,rint,stats]=regress(income,X); we get b = 1272.1 -18.4 0.8 stats = bint = -632.6 -39.1 -31.4 0.6316 14.5703 3176.9 2.3 32.9 0.0002 78972 Of course, one can find many shortcomings of this model. For example, having the entire population in agriculture results in prediction of 1272.1 − 1840 < 0 negative income per capita.