Developing Leverage Statistics Manually in R

Question

I've asked this question on Stack Overflow but think it might be better here.

I currently have the below data frame df with the regression equation

regression_eq <- -0.09999975 * df$x + 1.999999 * df$y and am trying to determine the leverage points by hand. I have now reviewed many sources and am still stuck. I understand that typically you would use the lm function then develop the diagnostics plots with R but was hoping to do this by hand as I developed results separately.

        response        x             y       xx            xy           yy      xxx
1 -0.1999981 2.000000 -4.794927e-09 4.000000 -9.589855e-09 2.299133e-17 8.000000
2 -0.2796748 1.997601 -3.995733e-02 3.990411 -7.981882e-02 1.596588e-03 7.971252
3 -0.3590789 1.994407 -7.981885e-02 3.977661 -1.591913e-01 6.371049e-03 7.933076
4 -0.4381798 1.990421 -1.195688e-01 3.961775 -2.379922e-01 1.429669e-02 7.885600
5 -0.5169470 1.985645 -1.591913e-01 3.942786 -3.160975e-01 2.534188e-02 7.828973
6 -0.5953499 1.980083 -1.986710e-01 3.920729 -3.933850e-01 3.947016e-02 7.763370

If it is possible to convert this manually in R please let me know as I have found results that I am unsure of thus far and want to clarify if they are correct.

Thanks in advance.

Answer 1

From Wikipedia the leverage score of each element (row) of the independent variable is the diagonal of the hat matrix:

$\mathbf{H} = \mathbf{X(X^{T}X)^{-1}X^{T}}$

where $\mathbf{X}$ is the design matrix.

Here I reproduce your example. For comparison with the manual calculation, I fit the model with lm and use the R function hatvalues to extract the leverage:

df <- read.table(text="
  response        x             y
-0.1999981 2.000000 -4.794927e-09
-0.2796748 1.997601 -3.995733e-02
-0.3590789 1.994407 -7.981885e-02
-0.4381798 1.990421 -1.195688e-01
-0.5169470 1.985645 -1.591913e-01
-0.5953499 1.980083 -1.986710e-01
", header= TRUE)

fit <- lm(response ~ x + y, data= df)
leverage <- hatvalues(fit)
leverage
   1    2    3    4    5    6 
0.82 0.31 0.37 0.37 0.31 0.82 

This is instead the manual calculation returning the same results:

# Design matrix. Can be done also with model.matrix(~x + y)
X <- cbind(intercept= 1, x= df$x, y= df$y)
X
     intercept x        y
[1,]         1 2 -4.8e-09
[2,]         1 2 -4.0e-02
[3,]         1 2 -8.0e-02
[4,]         1 2 -1.2e-01
[5,]         1 2 -1.6e-01
[6,]         1 2 -2.0e-01

H <- X %*% solve(t(X) %*% X) %*% t(X)
diag(H)
[1] 0.82 0.31 0.37 0.37 0.31 0.82