Introduction to EMMREML

Introduction

The EMMREML package provides an improved mixed model solver that handles a variety of covariance structures for genetic data analysis. It is particularly useful in plant and animal breeding, quantitative genetics, and genome-wide association studies (GWAS).

This vignette introduces the main functionality of the package and shows how to use it through practical examples.

Key Features

The package provides several key functions:

emmreml: Solves mixed models with known covariance structure using the ‘EMMA’ algorithm
emmremlMultiKernel: Fits models with multiple random components with known covariance structures
emmremlMultivariate: Solves multivariate Gaussian mixed models with known covariance structures
gwas: Performs genome-wide association studies with mixed models

Additionally, the package includes functionality for:

Computing and manipulating relationship/kinship matrices
Visualizing results
Managing marker data

Installation

You can install the released version of EMMREML from CRAN:

#install.packages("EMMREML")

Or the development version from GitHub:

# install.packages("devtools")
#devtools::install_github("denizakdemir/EMMREML")

Basic Usage: Mixed Models with Known Covariance

The core functionality of EMMREML is fitting mixed models with known covariance structures. Let’s start with a simple example:

library(EMMREML)
#> Loading required package: Matrix

Example Data

First, let’s generate some example data for our demonstration:

# Generate example data
set.seed(123)
n <- 100 # number of samples

# Create fixed effects
X <- cbind(1, rnorm(n))
colnames(X) <- c("intercept", "x")

# Create a realistic kinship matrix
Z <- diag(n)
K_raw <- matrix(rnorm(n*20), n, 20)
K <- tcrossprod(K_raw)
diag(K) <- diag(K) + 1e-6  # ensure positive definiteness
K <- K / mean(diag(K))  # Scale K

# Simulate data from the model
u <- crossprod(chol(K), rnorm(n)) * sqrt(2)  # random effect with Var=2*K
e <- rnorm(n)  # residual with Var=1
y <- 10 + 2 * X[, 2] + u + e

# Create data frame
df <- data.frame(y = y, x = X[, 2], id = 1:n)

Fitting a Mixed Model

We can fit a mixed model using the emmreml function:

# Fit model
result <- emmreml(y = df$y, 
                  X = cbind(1, df$x),  # Intercept + x
                  Z = diag(nrow(df)),  # Identity for observations
                  K = K)

# Display results
cat("Variance components:\n")
#> Variance components:
cat("Genetic variance (Vu):", result$Vu, "\n")
#> Genetic variance (Vu): 0.7811974
cat("Residual variance (Ve):", result$Ve, "\n")
#> Residual variance (Ve): 0.8963003
cat("\nFixed effects:\n")
#> 
#> Fixed effects:
names(result$betahat) <- c("Intercept", "x")
print(result$betahat)
#>          [,1]
#> [1,] 9.848929
#> [2,] 2.128694
#> attr(,"names")
#> [1] "Intercept" "x"

Multi-Kernel Models

The emmremlMultiKernel function allows you to fit models with multiple random components. This is useful when you have multiple sources of genetic effects or environmental factors that you want to model separately.

# Create a second kinship matrix
set.seed(456)
K2_raw <- matrix(rnorm(n*15), n, 15)
K2 <- tcrossprod(K2_raw)
diag(K2) <- diag(K2) + 1e-6
K2 <- K2 / mean(diag(K2))

# Simulate new data with two random components
u1 <- crossprod(chol(K), rnorm(n)) * sqrt(1.5)  # first random effect
u2 <- crossprod(chol(K2), rnorm(n)) * sqrt(0.8)  # second random effect
e <- rnorm(n) * sqrt(0.7)  # residual
y_mk <- 10 + 2 * X[, 2] + u1 + u2 + e

# Create data frame
df_mk <- data.frame(y = y_mk, x = X[, 2], id = 1:n)

# List of kernel matrices
K_list <- list(K1 = K, K2 = K2)

# Fit multi-kernel model
mk_result <- emmremlMultiKernel(y = df_mk$y,
                               X = cbind(1, df_mk$x),
                               Z = diag(nrow(df_mk)),
                               K = K_list)

# Display results
cat("Variance components:\n")
print(mk_result$varcomp)
cat("\nFixed effects:\n")
names(mk_result$betahat) <- c("Intercept", "x")
print(mk_result$betahat)

Multivariate Models

EMMREML can also handle multivariate mixed models, where you analyze multiple response variables simultaneously using the emmremlMultivariate function:

# Generate multivariate example data
p <- 3  # number of traits

# Create covariance matrices
Vg <- matrix(c(2, 1, 0.5, 1, 1.5, 0.75, 0.5, 0.75, 1), p, p)
Ve <- matrix(c(1, 0.2, 0.1, 0.2, 0.8, 0.15, 0.1, 0.15, 0.5), p, p)

# Simulate data
B <- matrix(c(10, 2, 5, 1, 3, 0.5), 2, p)  # fixed effects for 3 traits
G <- matrix(rnorm(n*p), n, p) %*% chol(Vg)  # random effects
E <- matrix(rnorm(n*p), n, p) %*% chol(Ve)  # errors
Y <- X %*% B + G + E

# Create data frame
df_mv <- data.frame(y1 = Y[,1], y2 = Y[,2], y3 = Y[,3], x = X[,2], id = 1:n)

# Fit multivariate model
mv_result <- emmremlMultivariate(formula = cbind(y1, y2, y3) ~ x, 
                               random = ~ id, 
                               data = df_mv, 
                               K = K)

# Display results
cat("Genetic covariance matrix (Vg):\n")
print(mv_result$Vg)
cat("\nError covariance matrix (Ve):\n")
print(mv_result$Ve)
cat("\nFixed effects:\n")
print(mv_result$beta)

Genome-Wide Association Studies (GWAS)

The package includes functionality for performing GWAS with mixed models, which helps control for population structure:

# Generate GWAS example data
p <- 1000  # number of markers
markers <- matrix(sample(c(0, 1, 2), n*p, replace=TRUE), nrow=n, ncol=p)
colnames(markers) <- paste0("SNP", 1:p)
    
# Select some causal variants
causal_idx <- c(100, 300, 700)
betas <- rep(0, p)
betas[causal_idx] <- c(0.5, 0.4, 0.3)
    
# Create genetic values from markers
g <- markers %*% betas
    
# Create phenotype for GWAS
y_gwas <- 10 + g + u + e
df_gwas <- data.frame(y = y_gwas, x = X[, 2], id = 1:n)

# Run GWAS
gwas_result <- gwas(
  formula = y ~ x, 
  random = ~ id, 
  data = df_gwas,
  markers = markers,
  K = K,
  accelerated = TRUE
)

# Display top results
head(gwas_result$results[order(gwas_result$results$p_value),], 10)

Working with Kinship Matrices

The kinship module provides functions for processing marker data and calculating different types of kinship matrices:

# Calculate a VanRaden kinship matrix
K_vr <- calculate_kinship(markers, method = "vanraden1", coding = "012")

# Ensure positive definiteness
diag(K_vr) <- diag(K_vr) + 1e-4

# Use this kinship matrix in a mixed model
mod <- emmreml(y = df$y, 
              X = cbind(1, df$x),
              Z = diag(nrow(df)),
              K = K_vr)

Conclusion

The EMMREML package provides a comprehensive set of tools for fitting mixed models with known covariance structures. It is particularly useful for genetic data analysis, including GWAS and multivariate trait analysis. The package offers optimized implementations for improved performance, especially when dealing with identity kinship matrices.

For more detailed information, refer to the function documentation using the ? operator, e.g., ?emmreml.