RAINBOWR (version 0.1.35)

EM3.linker.cpp: Equation of mixed model for multi-kernel (fast, for limited cases)

Description

This function solves multi-kernel mixed model using fastlmm.snpset approach (Lippert et al., 2014). This function can be used only when the kernels other than genomic relationship matrix are linear kernels.

Usage

EM3.linker.cpp(
  y0,
  X0 = NULL,
  ZETA = NULL,
  Zs0 = NULL,
  Ws0,
  Gammas0 = lapply(Ws0, function(x) diag(ncol(x))),
  gammas.diag = TRUE,
  X.fix = TRUE,
  eigen.SGS = NULL,
  eigen.G = NULL,
  n.core = 1,
  tol = NULL,
  bounds = c(1e-06, 1e+06),
  optimizer = "nlminb",
  traceInside = 0,
  n.thres = 450,
  spectral.method = NULL,
  REML = TRUE,
  pred = TRUE,
  return.u.always = TRUE,
  return.u.each = TRUE,
  return.Hinv = TRUE
)

Value

$y.pred

The fitting values of y \(y = X\beta + Zu\)

$Vu

Estimator for \(\sigma^2_u\), all of the genetic variance

$Ve

Estimator for \(\sigma^2_e\)

$beta

BLUE(\(\beta\))

$u

BLUP(Sum of \(Zu\))

$u.each

BLUP(Each \(u\))

$weights

The proportion of each genetic variance (corresponding to each kernel of ZETA) to Vu

$LL

Maximized log-likelihood (full or restricted, depending on method)

$Vinv

The inverse of \(V = Vu \times ZKZ' + Ve \times I\)

$Hinv

The inverse of \(H = ZKZ' + \lambda I\)

Arguments

y0

A \(n \times 1\) vector. A vector of phenotypic values should be used. NA is allowed.

X0

A \(n \times p\) matrix. You should assign mean vector (rep(1, n)) and covariates. NA is not allowed.

ZETA

A list of variance (relationship) matrix (K; \(m \times m\)) and its design matrix (Z; \(n \times m\)) of random effects. You can use only one kernel matrix. For example, ZETA = list(A = list(Z = Z, K = K)) Please set names of list "Z" and "K"!

Zs0

A list of design matrices (Z; \(n \times m\) matrix) for Ws. For example, Zs0 = list(A.part = Z.A.part, D.part = Z.D.part)

Ws0

A list of low rank matrices (W; \(m \times k\) matrix). This forms linear kernel \(K = W \Gamma W'\). For example, Ws0 = list(A.part = W.A, D.part = W.D)

Gammas0

A list of matrices for weighting SNPs (Gamma; \(k \times k\) matrix). This forms linear kernel \(K = W \Gamma W'\). For example, if there is no weighting, Gammas0 = lapply(Ws0, function(x) diag(ncol(x)))

gammas.diag

If each Gamma is the diagonal matrix, please set this argument TRUE. The calculationtime can be saved.

X.fix

If you repeat this function and when X0 is fixed during iterations, please set this argument TRUE.

eigen.SGS

A list with

$values

Eigen values

$vectors

Eigen vectors

The result of the eigen decompsition of \(SGS\), where \(S = I - X(X'X)^{-1}X'\), \(G = ZKZ'\). You can use "spectralG.cpp" function in RAINBOWR. If this argument is NULL, the eigen decomposition will be performed in this function. We recommend you assign the result of the eigen decomposition beforehand for time saving.

eigen.G

A list with

$values

Eigen values

$vectors

Eigen vectors

The result of the eigen decompsition of \(G = ZKZ'\). You can use "spectralG.cpp" function in RAINBOWR. If this argument is NULL, the eigen decomposition will be performed in this function. We recommend you assign the result of the eigen decomposition beforehand for time saving.

n.core

Setting n.core > 1 will enable parallel execution on a machine with multiple cores.

tol

The tolerance for detecting linear dependencies in the columns of G = ZKZ'. Eigen vectors whose eigen values are less than "tol" argument will be omitted from results. If tol is NULL, top 'n' eigen values will be effective.

bounds

Lower and upper bounds for weights.

optimizer

The function used in the optimization process. We offer "optim", "optimx", and "nlminb" functions.

traceInside

Perform trace for the optimzation if traceInside >= 1, and this argument shows the frequency of reports.

n.thres

If \(n >= n.thres\), perform EMM1.cpp. Else perform EMM2.cpp.

spectral.method

The method of spectral decomposition. In this function, "eigen" : eigen decomposition and "cholesky" : cholesky and singular value decomposition are offered. If this argument is NULL, either method will be chosen accorsing to the dimension of Z and X.

REML

You can choose which method you will use, "REML" or "ML". If REML = TRUE, you will perform "REML", and if REML = FALSE, you will perform "ML".

pred

If TRUE, the fitting values of y is returned.

return.u.always

If TRUE, BLUP (`u`; \(u\)) will be returned.

return.u.each

If TRUE, the function also computes each BLUP corresponding to different kernels (when solving multi-kernel mixed-effects model). It takes additional time compared to the one with `return.u.each = FALSE`.

return.Hinv

If TRUE, \(H ^ {-1} = (Var[y] / \sum _{l=1} ^ {L} \sigma _ {l} ^ 2) ^ {-1}\) will be computed. It also returns \(V ^ {-1} = (Var[y]) ^ {-1}\).

References

Kang, H.M. et al. (2008) Efficient Control of Population Structure in Model Organism Association Mapping. Genetics. 178(3): 1709-1723.

Zhou, X. and Stephens, M. (2012) Genome-wide efficient mixed-model analysis for association studies. Nat Genet. 44(7): 821-824.

Lippert, C. et al. (2014) Greater power and computational efficiency for kernel-based association testing of sets of genetic variants. Bioinformatics. 30(22): 3206-3214.

Examples

# \dontshow{
  ### Import RAINBOWR
  require(RAINBOWR)
  
  ### Load example datasets
  data("Rice_Zhao_etal")
  Rice_geno_score <- Rice_Zhao_etal$genoScore
  Rice_geno_map <- Rice_Zhao_etal$genoMap
  Rice_pheno <- Rice_Zhao_etal$pheno

    
  ### Select one trait for example
  trait.name <- "Flowering.time.at.Arkansas"
  y <- as.matrix(Rice_pheno[1:30, trait.name, drop = FALSE])
  # use first 30 accessions
  
  ### Remove SNPs whose MAF <= 0.05
  x.0 <- t(Rice_geno_score)
  MAF.cut.res <- MAF.cut(x.0 = x.0, map.0 = Rice_geno_map)
  x <- MAF.cut.res$x
  map <- MAF.cut.res$map
  
  
  ### Estimate additive genomic relationship matrix (GRM)
  K.A <- calcGRM(genoMat = x)
  
  
  ### Modify data
  Z <- design.Z(pheno.labels = rownames(y),
                geno.names = rownames(K.A))  ### design matrix for random effects
  pheno.mat <- y[rownames(Z), , drop = FALSE]
  ZETA <- list(A = list(Z = Z, K = K.A))
  
  
  ### Including the first 20 SNPs
  W.A <- x[, 1:20]    ### marker genotype data of first 20 SNPs
  
  Zs0 <- list(A.part = Z)
  Ws0 <- list(A.part = W.A)       ### This will be regarded as linear kernel
  ### for the variance-covariance matrix of another random effects.
  
  
  ### Solve multi-kernel linear mixed effects model (2 random efects)
  EM3.linker.res <- EM3.linker.cpp(y0 = pheno.mat, X0 = NULL, ZETA = ZETA,
                                   Zs0 = Zs0, Ws0 = Ws0)
  Vu <- EM3.linker.res$Vu   ### estimated genetic variance
  Ve <- EM3.linker.res$Ve   ### estimated residual variance
  weights <- EM3.linker.res$weights   ### estimated proportion of two genetic variances
  herit <- Vu * weights / (Vu + Ve)   ### genomic heritability (all chromosomes, chromosome 12)
  
  beta <- EM3.linker.res$beta   ### Here, this is an intercept.
  u.each <- EM3.linker.res$u.each   ### estimated genotypic values (all chromosomes, chromosome 12)
# }


# \donttest{
  ### Import RAINBOWR
  require(RAINBOWR)
  
  ### Load example datasets
  data("Rice_Zhao_etal")
  Rice_geno_score <- Rice_Zhao_etal$genoScore
  Rice_geno_map <- Rice_Zhao_etal$genoMap
  Rice_pheno <- Rice_Zhao_etal$pheno
  
  ### View each dataset
  See(Rice_geno_score)
  See(Rice_geno_map)
  See(Rice_pheno)
  
  ### Select one trait for example
  trait.name <- "Flowering.time.at.Arkansas"
  y <- as.matrix(Rice_pheno[, trait.name, drop = FALSE])
  
  ### Remove SNPs whose MAF <= 0.05
  x.0 <- t(Rice_geno_score)
  MAF.cut.res <- MAF.cut(x.0 = x.0, map.0 = Rice_geno_map)
  x <- MAF.cut.res$x
  map <- MAF.cut.res$map
  
  
  ### Estimate additive genomic relationship matrix (GRM)
  K.A <- calcGRM(genoMat = x)
  
  
  ### Modify data
  Z <- design.Z(pheno.labels = rownames(y),
                geno.names = rownames(K.A))  ### design matrix for random effects
  pheno.mat <- y[rownames(Z), , drop = FALSE]
  ZETA <- list(A = list(Z = Z, K = K.A))
  
  
  ### Including the additional linear kernel for chromosome 12
  chrNo <- 12
  W.A <- x[, map$chr == chrNo]    ### marker genotype data of chromosome 12
  
  Zs0 <- list(A.part = Z)
  Ws0 <- list(A.part = W.A)       ### This will be regarded as linear kernel
  ### for the variance-covariance matrix of another random effects.
  
  
  ### Solve multi-kernel linear mixed effects model (2 random efects)
  EM3.linker.res <- EM3.linker.cpp(y0 = pheno.mat, X0 = NULL, ZETA = ZETA,
                                   Zs0 = Zs0, Ws0 = Ws0)
  (Vu <- EM3.linker.res$Vu)   ### estimated genetic variance
  (Ve <- EM3.linker.res$Ve)   ### estimated residual variance
  (weights <- EM3.linker.res$weights)   ### estimated proportion of two genetic variances
  (herit <- Vu * weights / (Vu + Ve))   ### genomic heritability (all chromosomes, chromosome 12)
  
  (beta <- EM3.linker.res$beta)   ### Here, this is an intercept.
  u.each <- EM3.linker.res$u.each   ### estimated genotypic values (all chromosomes, chromosome 12)
  See(u.each)
# }