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

Description Usage Arguments Value References Examples

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

y = X β + ∑ _{l=1} ^ {L} Z _ {l} u _ {l} + ε

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

EM3.cpp(
  y,
  X0 = NULL,
  ZETA,
  eigen.G = NULL,
  eigen.SGS = NULL,
  tol = NULL,
  optimizer = "nlminb",
  traceInside = 0,
  n.thres = 450,
  REML = TRUE,
  pred = TRUE
)

`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 $values Eigen values $vectors Eigen vectors The result of the eigen decompsition of G = ZKZ'. You can use "spectralG.cpp" function in RAINBOW. 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 $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 RAINBOW. 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.
`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.

$y.pred: The fitting values of y y = Xβ + Zu
$Vu: Estimator for σ^2_u, all of the genetic variance
$Ve: Estimator for σ^2_e
$beta: BLUE(β)
$u: BLUP(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' + λ I

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.

  ### Import RAINBOW
  require(RAINBOW)
  
  ### 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) & 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, X = 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
  (herit <- Vu * weights / (Vu + Ve))   ### genomic heritability (additive, additive x additive)
  
  (beta <- EM3.res$beta)   ### Here, this is an intercept.
  u <- EM3.res$u   ### estimated genotypic values (additive, additive x additive)
  See(u)
  
  
  ### 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, X = 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[2] - 10,
       y = plotRange[1] + 10,
       paste0("R2 = ", round(R2, 3)), 
       cex = 1.5)