{"word_count": 670, "character_count": 7745}
Last week, in our STT5100 (applied linear models) class, I've introduce the hat matrix, and the notion of leverage. In a classical regression model, [latex]\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}[/latex] (in a matrix form), the ordinary least square estimator of parameter [latex]\boldsymbol{\beta}[/latex] is [latex display="true"]\widehat{\boldsymbol{\beta}}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{y}[/latex]The prediction can then be written[latex display="true"]\widehat{\boldsymbol{y}}=\boldsymbol{X}\widehat{\boldsymbol{\beta}}=\underbrace{\color{blue}{\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top}}_{\color{blue}{\boldsymbol{H}}}\boldsymbol{y}[/latex]where [latex]\color{blue}{\boldsymbol{H}}[/latex] is called the hat matrix.The matrix is idempotent, i.e. [latex display="true"]\boldsymbol{H}^2={\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\textcolor{grey}{\boldsymbol{X}^\top{\boldsymbol{X}}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}}\boldsymbol{X}^\top}={\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top}=\boldsymbol{H}[/latex]so it can be interpreted as a projection matrix. Furthermore, since[latex]\boldsymbol{H}\boldsymbol{X}=\boldsymbol{X}[/latex] (just do the maths), the projection is on a subspace that contains all the linear combinations of columns of [latex]\boldsymbol{X}[/latex]. One can also observe that [latex]\mathbb{I}-\boldsymbol{H}[/latex] is also a projection matrix. And we can write[latex display="true"]\boldsymbol{y}=\underbrace{\boldsymbol{H}\boldsymbol{y}}_{\widehat{\boldsymbol{y}}}+\underbrace{(\mathbb{I}-\boldsymbol{H})\boldsymbol{y}}_{\widehat{\boldsymbol{\varepsilon}}}[/latex]where [latex]\widehat{\boldsymbol{y}}[/latex] is the orthogonal projection of [latex]\boldsymbol{y}[/latex] on the (linear) space of linear combinations of columns of [latex]\boldsymbol{X}[/latex], and [latex]\widehat{\boldsymbol{y}}\perp\widehat{\boldsymbol{\varepsilon}}[/latex], which gives the classical interpretation of residuals, being unpredictible (at least with a linear model using variables [latex]\boldsymbol{X}[/latex]).Let's move a bit faster now (we've seen many other properties last week), and consider elements on the diagonal of matrix [latex]\boldsymbol{H}[/latex]. Recall that we haveso entry [latex]\boldsymbol{H}_{i,i}[/latex] is a measure of the influence of entry [latex]\boldsymbol{y}_i[/latex] on its prediction latex]\widehat{\boldsymbol{y}}_i[/latex].We have seen that[latex display="true"]\sum_{i=1}^n\boldsymbol{H}_{i,i}=\text{trace}(\boldsymbol{H})=\text{trace}(\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top)[/latex]which can be written[latex display="true"]\sum_{i=1}^n\boldsymbol{H}_{i,i}=\text{trace}\boldsymbol{X}^\top(\boldsymbol{X}(\boldsymbol{X}^\top\boldsymbol{X})^{-1})=\text{trace}(\mathbb{I})=p[/latex]where classically [latex]p=k+1[/latex], where [latex]k[/latex] is the number of explanatory variables. Further, since [latex]\boldsymbol{H}[/latex] is idempotent, we can write (from [latex]\boldsymbol{H}=\boldsymbol{H}^2[/latex]) that[latex display="true"]\boldsymbol{H}_{i,i}=\boldsymbol{H}_{i,i}^2 + \sum_{j\neq i}\boldsymbol{H}_{i,j}\boldsymbol{H}_{j,i}=\boldsymbol{H}_{i,i}^2 + \sum_{j\neq i}\boldsymbol{H}_{i,j}^2[/latex]One the one hand, since the second term is positive [latex]\boldsymbol{H}_{i,i}\geq\boldsymbol{H}_{i,i}^2[/latex], i.e. [latex]1\geq\boldsymbol{H}_{i,i}[/latex]. And since both terms are positive, then [latex]\boldsymbol{H}_{i,i}\in[0,1][/latex]. And there was a question in the course on the sharpeness of the bounds.Using Anscombe's dataset, we've seen that it was possible to get a leverage of 1. Using something rather similardf = data.frame(x = c(rep(1,10),6), y = c(1:10,8))plot(df)we obtainmodel = lm(y~x,data=df)abline(model,col="red",lwd=2)H = lm.influence(model)$hatplot(1:11,H,type="h")The very last observation, the one one the right, is here extremely influencial : if we remove it, the model is completely different ! And here, we reach the upper bound, [latex]\boldsymbol{H}_{11,11}=1[/latex]. Observe that all other points are equally influencial, and because on the constraint on the trace of the matrix, [latex]\boldsymbol{H}_{i,i}=1/10[/latex] when [latex]i\in\{1,2,\cdots,10\}[/latex].Now, what about the lower bound ? In order to have some sort of "non-influencial" observations, consider the two following case. the case where one observation (below the first one) is such that [latex]\widehat{\boldsymbol{y}}_{i}=\boldsymbol{y}_{i}[/latex] (perfect prediction) the case where one observation (below the tenth one) is such that [latex]\boldsymbol{x}_{i}=\overline{\boldsymbol{x}}[/latex] and [latex]\boldsymbol{y}_{i}=\overline{\boldsymbol{y}}[/latex] (from the first order condition - or normal equation), the fitted regression line always go through point [latex](\overline{\boldsymbol{x}},\overline{\boldsymbol{y}})[/latex]Let us move two observations from our dataset,mean(c(4,rep(1,8),6))[1] 1.8df = data.frame(x = c(4,rep(1,8),6,1.8),y = c(predict(model,newdata=data.frame(x=4)),2:9,8,predict(model,newdata=data.frame(x=1.8))))We now haveIf we compute the leverages, we obtainmodel = lm(y~x,data=df)H = lm.influence(model)$hatplot(1:11,H,type="h")so, for the first observation, its leverage actually increased (the blue part), and for the tenth one, we have the lowest influence, but it is not zero. Is it possible to reach zero ?Here, observe that for the tenth observation, [latex]\boldsymbol{H}_{i,i}=1/n[/latex]. And actually, that's the best we can do... We can prove that, in the case of a simple regression (as above)[latex display="true"]\boldsymbol{H}_{i,i}=\frac{1}{n}+\frac{(x_i-\overline{x})^2}{n\text{Var}(x)}[/latex]which is minimum when [latex]x_i=\overline{x}[/latex], and then [latex]\boldsymbol{H}_{i,i}=1/n[/latex], otherwise [latex]\boldsymbol{H}_{i,i}&#62;1/n[/latex]. And this property is also valid in a multiple regression (as soon as an intercept is included in the regression - which should always be the case). To prove that result, let [latex]\tilde{\boldsymbol{X}}[/latex] denote the matrix of centered variables [latex]\boldsymbol{X}[/latex], then we can prove that [latex display="true"]\boldsymbol{H}_{i,i}=\frac{1}{n}+\big[\tilde{\boldsymbol{X}}(\tilde{\boldsymbol{X}}^\top\tilde{\boldsymbol{X}})^{-1}\tilde{\boldsymbol{X}}^\top\big]_{i,i}[/latex](which is basically a matrix version of the previous equation).I can maybe add another comment on Anscombe's data. We've seen that on the right that we did reach 1. But I did not prove it. One way to prove it is actually to focus on the remaining [latex]n-1[/latex] points, on the left. Those have all the same [latex]x[/latex] values. We can prove that if [latex]\boldsymbol{X}_{i_1}=\boldsymbol{X}_{i_2}[/latex], then [latex display="true"]\boldsymbol{H}_{i_1,i_2}=\boldsymbol{X}_{i_1}^\top(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}_{i_2}=\boldsymbol{H}_{i_1,i_1}[/latex]hence, using the relationship obtained since the hat matrix is idempotent[latex display="true"]\boldsymbol{H}_{i_1,i_1}=2\boldsymbol{H}_{i_1,i_1}^2+\sum_{j\notin\{i_1,i_2\}}\boldsymbol{H}_{i_1,j}^2[/latex]thus, we now have[latex display="true"]\boldsymbol{H}_{i_1,i_1}\big(1-2\boldsymbol{H}_{i_1,i_1}\big)&#62;0[/latex]i.e. [latex]\boldsymbol{H}_{i_1,i_1}\in[0,1/2][/latex], where the upper bound becomes [latex]1/(n-1)[/latex] "duplicates". So for [latex]n-1[/latex] [latex]\boldsymbol{H}_{i,i}[/latex]'s, we have values below [latex]1/(n-1)[/latex], the last one should be below [latex]1[/latex] and the sum has to be [latex]k=2[/latex] . So we have the value of the [latex]n[/latex] [latex]\boldsymbol{H}_{i,i}[/latex]'s.&#160;