In variani/biglmmz: Low Rank Linear Mixed Models (LMMs) Powered by 'bigstatsr'

knitr::opts_chunk$set(
  comment = "#>",  
  fig.path = "man/figures/README-",
  out.width = "100%"
)

biglmmz

Low-rank linear mixed models (LMMs) powered by bigstatsr.

Installation

# install.packages("remotes")
remotes::install_github("variani/biglmmz")

Examples

We first check whether the low-rank linear mixed model (LMM) recovers the heritability of a quantitative trait on simulated data (1500 individuals, 200 genetic markers, heritability 80%).

library(biglmmz)

## load simulated genotypes
(G <- attach_example200())
(N <- nrow(G))
(M <- ncol(G))

G[1:5, 1:10]

## simulate a phenotype with heritability h2 = 0.8 
h2 <- 0.8
# generate effect sizes under the polygenic model
b <- rnorm(M, 0, sqrt(h2/M))
# pre-compute mean and sd values for each genotype
stats <- big_scale()(G) 
# compute the matrix-vector product, Z b
Zb <- big_prodVec(G, b, center = stats$center, scale = stats$scale)
y <- Zb + rnorm(N, 0, sqrt(1 - h2))

## fit low-rank linear mixed model
mod <- biglmmg(y, G = G)
# check the estimate of h2
mod$gamma

We next compute the effective sample size (ESS) multiplier for the LMM. See the pre-print.

Calling the ess function:

# varianace components = s2 * c(h2, 1 - h2)
h2 <- mod$gamma
s2 <- mod$s2

ess(G, h2 = h2, s2 = s2)

Calculating the ESS manually:

# EVD on K = Z'Z/M, where Z is a matrix of scaled genotypes G
K <- big_crossprodSelf(G, fun.scaling = big_scale_grm(M = M))[]
# EVD of K
lambdas <- eigen(K, only.values = TRUE)$values

# the multiplier
mult <- (1/s2) * (sum(1/(h2*lambdas + (1-h2))) + (N-M)/(1-h2)) / N

## print results
(res <- data.frame(N = N, M = M, h2_hat = h2, s2 = s2, mult = mult))