## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.2, dpi = 96, out.width = "100%")

## ----example------------------------------------------------------------------
library(matsketch)
set.seed(11)
dat <- sim_genomic(n = 1500, p = 3000, h2 = 0.5, pops = 4)
G <- grm_matrix(dat$M)
fit <- reml_sketch(dat$y, G)
fit

## ----history------------------------------------------------------------------
fit$history
plot(fit)

## ----compare------------------------------------------------------------------
exact <- reml_exact(dat$y, as.matrix(G))
rbind(sketched = c(fit$sigma2, h2 = fit$h2, se_h2 = fit$se[["h2"]]),
      exact = c(exact$sigma2, h2 = exact$h2, se_h2 = exact$se[["h2"]]))

## ----spectrum-----------------------------------------------------------------
ev <- eigen(as.matrix(G), symmetric = TRUE, only.values = TRUE)$values
mu <- fit$sigma2[["residual"]] / fit$sigma2[["genetic"]]
keep <- ev > 1e-8
plot(which(keep), ev[keep], log = "y", pch = 19, cex = 0.4, col = "#0072B2",
     xlab = "index", ylab = "eigenvalue of G")
abline(v = 100.5, lty = 2, col = "grey50")
abline(h = mu, lty = 3, lwd = 2, col = "#D55E00")
legend("topright", c("preconditioner rank", "mu at the estimate"),
       lty = c(2, 3), lwd = c(1, 2), col = c("grey50", "#D55E00"),
       bty = "n")

## ----accuracy-----------------------------------------------------------------
acc <- read.csv(system.file("extdata", "reml-accuracy.csv",
                            package = "matsketch"))
z <- (acc$sketch_h2 - acc$exact_h2) / acc$exact_se
summary(z)

## ----accuracy-plot, fig.height = 3.6------------------------------------------
op <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 1, 1))
cols <- ifelse(acc$true_h2 < 0.5, "#0072B2", "#D55E00")
plot(acc$exact_h2, acc$sketch_h2, pch = 19, col = cols, asp = 1,
     xlab = "exact REML estimate", ylab = "sketched REML estimate")
abline(0, 1, lty = 2)
legend("topleft", c("true h2 = 0.3", "true h2 = 0.6"), pch = 19,
       col = c("#0072B2", "#D55E00"), bty = "n")
hist(z, breaks = 12, col = "grey80", border = "white", main = "",
     xlab = "(sketched - exact) / exact SE")
par(op)

## ----scaling------------------------------------------------------------------
sc <- read.csv(system.file("extdata", "reml-scaling.csv",
                           package = "matsketch"))
secs <- with(sc, tapply(seconds, list(n, method), sum))
secs <- secs[, c("form_G", "exact", "eigen", "sketch_dense", "sketch_lazy")]
round(secs, 1)

## ----scaling-plot-------------------------------------------------------------
tot <- cbind(
  `exact REML` = secs[, "form_G"] + secs[, "exact"],
  `eigendecomposition only` = secs[, "form_G"] + secs[, "eigen"],
  `sketch, G formed` = secs[, "form_G"] + secs[, "sketch_dense"],
  `sketch, grm_matrix()` = secs[, "sketch_lazy"]
)
n <- as.numeric(rownames(secs))
cols <- c("#999999", "#0072B2", "#E69F00", "#D55E00")
matplot(n, tot / 60, log = "xy", type = "b", pch = 19, lty = 1, lwd = 2,
        col = cols, xlab = "individuals", ylab = "minutes")
legend("topleft", colnames(tot), col = cols, lwd = 2, pch = 19, bty = "n")

## ----memory-------------------------------------------------------------------
mem <- with(sc[sc$method %in% c("exact", "sketch_dense", "sketch_lazy"), ],
            tapply(memory_gb, list(n, method), sum))
round(mem, 2)