# EM3.cpp: Equation of mixed model for multi-kernel (slow, general... In RAINBOWR: Genome-Wide Association Study with SNP-Set Methods

 EM3.cpp R Documentation

## Equation of mixed model for multi-kernel (slow, general version)

### Description

This function solves the following multi-kernel linear mixed effects model.

y = X \beta + \sum _{l=1} ^ {L} Z _ {l} u _ {l} + \epsilon

where Var[y] = \sum _{l=1} ^ {L} Z _ {l} K _ {l} Z _ {l}' \sigma _ {l} ^ 2 + I \sigma _ {e} ^ {2}.

### Usage

EM3.cpp(
y,
X0 = NULL,
ZETA,
eigen.G = NULL,
eigen.SGS = NULL,
tol = NULL,
n.core = NA,
optimizer = "nlminb",
traceInside = 0,
n.thres = 450,
REML = TRUE,
pred = TRUE,
return.u.always = TRUE,
return.u.each = TRUE,
return.Hinv = TRUE
)


### Arguments

 y 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 matrices and its design matrices of random effects. You can use more than one kernel matrix. For example, ZETA = list(A = list(Z = Z.A, K = K.A), D = list(Z = Z.D, K = K.D)) (A for additive, D for dominance) Please set names of lists "Z" and "K"! eigen.G A list with $valuesEigen values$vectorsEigen 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. eigen.SGS A list with $valuesEigen values$vectorsEigen 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. 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. n.core Setting n.core > 1 will enable parallel execution on a machine with multiple cores. 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. 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}.

### 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

### 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.

### Examples



### Import RAINBOWR
require(RAINBOWR)

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) & epistatic relationship matrix K.A <- calcGRM(genoMat = x) K.AA <- K.A * K.A ### additive x additive epistatic effects ### 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), AA = list(Z = Z, K = K.AA)) ### Solve multi-kernel linear mixed effects model (2 random efects) EM3.res <- EM3.cpp(y = pheno.mat, X0 = NULL, ZETA = ZETA) (Vu <- EM3.res$Vu)   ### estimated genetic variance
(Ve <- EM3.res$Ve) ### estimated residual variance (weights <- EM3.res$weights)   ### estimated proportion of two genetic variances

(beta <- EM3.res$beta) ### Here, this is an intercept. u.each <- EM3.res$u.each   ### estimated genotypic values (additive, additive x additive)
See(u.each)

### Perform genomic prediction with 10-fold cross validation (multi-kernel)
noNA <- !is.na(c(pheno.mat))   ### NA (missing) in the phenotype data

phenoNoNA <- pheno.mat[noNA, , drop = FALSE]   ### remove NA
ZETANoNA <- ZETA
ZETANoNA <- lapply(X = ZETANoNA, FUN = function (List) {
List$Z <- List$Z[noNA, ]

return(List)
})   ### remove NA

nFold <- 10    ### # of folds
nLine <- nrow(phenoNoNA)
idCV <- sample(1:nLine %% nFold)   ### assign random ids for cross-validation
idCV[idCV == 0] <- nFold

yPred <- rep(NA, nLine)

for (noCV in 1:nFold) {
print(paste0("Fold: ", noCV))
yTrain <- phenoNoNA
yTrain[idCV == noCV, ] <- NA   ### prepare test data

EM3.resCV <- EM3.cpp(y = yTrain, X0 = NULL, ZETA = ZETANoNA)   ### prediction
yTest <-  EM3.resCV\$y.pred     ### predicted values

yPred[idCV == noCV] <- yTest[idCV == noCV]
}

### Plot the results
plotRange <- range(phenoNoNA, yPred)
plot(x = phenoNoNA, y = yPred,xlim = plotRange, ylim = plotRange,
xlab = "Observed values", ylab = "Predicted values",
main = "Results of Genomic Prediction (multi-kernel)",
cex.lab = 1.5, cex.main = 1.5, cex.axis = 1.3)
abline(a = 0, b = 1, col = 2, lwd = 2, lty = 2)
R2 <- cor(x = phenoNoNA[, 1], y = yPred) ^ 2
text(x = plotRange - 10,
y = plotRange + 10,
paste0("R2 = ", round(R2, 3)),
cex = 1.5)



RAINBOWR documentation built on Sept. 12, 2023, 9:08 a.m.