lmmvar-vignette

In this document, we provide a short tutorial on how to use the \(\texttt{lmmvar}\) (Fast and reliable confidence intervals for a variance component or proportion) software in R. If you encounter any errors or strange behavior, please report it as an issue at https://github.com/yqzhang5972/lmmvar.

Model

Suppose, for a design matrix \(X \in \mathbb{R}^{n\times p}\) and known positive, semi-definite \(K \in \mathbb{R}^{n\times n}\), \[ Y \sim \mathrm{N}_n(X\beta, \sigma_g^2 K + \sigma_e^2 I_n). \] We are interested in test-statistics and confidence intervals for \(h^2 = \sigma_g^2 / (\sigma_g^2+\sigma_e^2)\), in some settings called heritability, and \(\sigma^2 = \sigma_g^2+\sigma_e^2\).

Installation

# install.packages("devtools")
# library(devtools)
# devtools::install_github("yqzhang5972/lmmvar")
library(lmmvar)

Example

# Function to simulate Y = Z U + E, where K = ZZ'.
# \beta = 0 WLOG since using REML
simData <- function(sigma2_g, sigma2_e = 1, Z) {  
  n = dim(Z)[1]
  m = dim(Z)[2] 
  e = rnorm(n, mean = 0, sd = sqrt(sigma2_e))
  u = rnorm(m, mean = 0, sd = sqrt(sigma2_g))
  y = Z %*% u + e
  return(y)
}

# Settings
n <- 100
p <- 10
h2 <- 0.5
s2 <- 1 

# K with autoregressive structure
K <- 0.9^abs(outer(1:n, 1:n, FUN = "-"))
K_eig <- eigen(K, symmetric = TRUE)
V <- K_eig$vectors        # K's eigenvectors
lambdas <- K_eig$values   # K's eigenvalues

# Ensure ZZ' = K
Z <-  K_eig$vectors %*% diag(sqrt(lambdas))

# Simulate dataset
set.seed(0)
X <- matrix(rnorm(n*p, 0, 1), nrow = n)
y <- simData(sigma2_g = h2 * s2, sigma2_e = (1 - h2) * s2, Z)

Xnew <- crossprod(V, X)         # transform X with V'X
ynew <- crossprod(V, y)         # transform y with V'y

Compute the 1d and 2d test statistics, 1d confidence interval and 2d test statistics matrix

# score test statistic of h2 at true value
varRatioTest1d(h2, ynew, Xnew, lambdas)   
## [1] 0.4845003
# score test statistic of (h2,sigma^2_p) at true value
varRatioTest2d(h2, s2, ynew, Xnew, lambdas) 
## [1] 0.6602168
# 95% confidence interval of h2 treating sigma^2 as nuisance parameter
confInv(ynew, Xnew, lambdas) 
## [1] 0.2523193 0.7092285

# Matrix of test statistics at a grid of h2 from 0 to 1 and sigma^2 from 0.0001
# to 3.
cr <- confReg(range_p = c(1e-4,3), ynew, Xnew, lambdas) 
library(ggplot2)
library(dplyr)
library(latex2exp)

# Plot joint test-statistic for different values of \sigma, fixed value of h^2
data.frame(p = seq(1e-4,3, length.out = 200)[30:150], teststatistic = cr[100,30:150]) %>%
  ggplot(aes(x = p, y = teststatistic))  + theme_bw() + #geom_point(size=0.8) +
  geom_line() + xlab(TeX("$sigma^2$")) + ylab("Test-statistic") +
  ggtitle(TeX("$h^2=0.5$")) + theme(plot.title = element_text(hjust = 0.5)) +
  geom_hline(yintercept=qchisq(0.95, 1), color = "red")+
  theme(axis.title=element_text(size=16), 
        plot.title = element_text(size = 18), axis.text = element_text(size=12))