Hat Matrix of a Smoother

Compute the hat matrix or smoother matrix, of ‘any’ (linear) smoother, smoothing splines, by default.

Usage

hatMat(x, trace= FALSE,
       pred.sm = function(x, y, ...)
                 predict(smooth.spline(x, y, ...), x = x)$y,
       ...)

Arguments

x: numeric vector or matrix.
trace: logical indicating if the whole hat matrix, or only its trace, i.e. the sum of the diagonal values should be computed.
pred.sm: a function of at least two arguments (x,y) which returns fitted values, i.e. \(\hat{y}\), of length compatible to x (and y).
...: optionally further arguments to the smoother function pred.sm.

Value

The hat matrix \(H\) (if trace = FALSE as per default) or a number, \(tr(H)\), the trace of \(H\), i.e., \(\sum_i H_{ii}\).

Note that dim(H) == c(n, n) where n <- length(x) also in the case where some x values are duplicated (aka ties).

References

Hastie and Tibshirani (1990). Generalized Additive Models. Chapman & Hall.

Author

Martin Maechler maechler@stat.math.ethz.ch

Examples

require(stats) # for smooth.spline() or loess()

x1 <- c(1:4, 7:12)
H1 <- hatMat(x1, spar = 0.5) # default : smooth.spline()

matplot(x1, H1, type = "l", main = "columns of smoother hat matrix")


## Example 'pred.sm' arguments for hatMat() :
pspl <-  function(x,y,...) predict(smooth.spline(x,y, ...), x = x)$y
pksm <-  function(x,y,...) ksmooth(sort(x),y, "normal", x.points=x, ...)$y
## Rather than ksmooth():
if(require("lokern"))
  pksm2 <- function(x,y,...) glkerns(x,y, x.out=x, ...)$est
#> Loading required package: lokern




## Explaining 'trace = TRUE'
all.equal(sum(diag((hatMat(c(1:4, 7:12), df = 4)))),
                    hatMat(c(1:4, 7:12), df = 4, trace = TRUE), tol = 1e-12)
#> [1] TRUE

## ksmooth() :
Hk <- hatMat(x1, pr = pksm, bandwidth = 2)
cat(sprintf("df = %.2f\n", sum(diag(Hk))))
#> df = 6.06
image(Hk)

Matrix::printSpMatrix(as(round(Hk, 2), "sparseMatrix"))
#>                                                        
#>  [1,] 0.70 0.28 0.02 .    .    .    .    .    .    .   
#>  [2,] 0.22 0.55 0.22 0.01 .    .    .    .    .    .   
#>  [3,] 0.01 0.22 0.55 0.22 .    .    .    .    .    .   
#>  [4,] .    0.02 0.28 0.70 .    .    .    .    .    .   
#>  [5,] .    .    .    .    0.70 0.28 0.02 .    .    .   
#>  [6,] .    .    .    .    0.22 0.55 0.22 0.01 .    .   
#>  [7,] .    .    .    .    0.01 0.22 0.54 0.22 0.01 .   
#>  [8,] .    .    .    .    .    0.01 0.22 0.54 0.22 0.01
#>  [9,] .    .    .    .    .    .    0.01 0.22 0.55 0.22
#> [10,] .    .    .    .    .    .    .    0.02 0.28 0.70

##--->  see demo("hatmat-ex")  for more (and larger) examples