R Under development (unstable) (2023-11-03 r85470 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> if (.Platform$OS.type != "windows") withAutoprint({
+ 
+     library(lme4)
+     source(system.file("testdata", "lme-tst-funs.R", package="lme4", mustWork=TRUE))# -> unn()
+ 
+     m <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
+     bruteForceHat <- function(object) {
+         with(getME(object, c("Lambdat", "Lambda", "Zt", "Z", "q", "X")), {
+             ## cp:= the cross product block matrix in (17) and (18):
+             W <- Diagonal(x = weights(object))
+             I <- Diagonal(q)
+             A.21 <- t(X) %*% W %*% Z %*% Lambda
+             cp <- rbind(cbind(Lambdat %*% Zt %*% W %*% Z %*% Lambda + I, t(A.21)),
+                         cbind(A.21, t(X) %*% W %*% X))
+             mm <- cbind(Z %*% Lambda, X)
+             ## a bit efficient: both cp and mm are typically quite sparse
+             ## mm %*% solve(as.matrix(cp)) %*% t(mm)
+             mm %*% solve(cp, t(mm), sparse=FALSE)
+         })
+     }
+ 
+ 
+     str(H <- bruteForceHat(m))
+ 
+     set.seed(7)
+     ii <- sample(nrow(sleepstudy), 500, replace=TRUE)
+     m2 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy[ii, ])
+ 
+     
+     stopifnot(all.equal(diag(H),
+                         unn(hatvalues(m)),  tol= 1e-14),
+               all.equal(diag(bruteForceHat(m2)),
+                         unn(hatvalues(m2)), tol= 1e-14)
+               )
+ }) ## skip on windows (for speed)
> 
> proc.time()
   user  system elapsed 
   0.10    0.07    0.17