View source: R/EMM_functions_cpp.R

EM3.linker.cpp | R Documentation |

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.

```
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
)
```

`y0` |
A |

`X0` |
A |

`ZETA` |
A list of variance (relationship) matrix (K; |

`Zs0` |
A list of design matrices (Z; |

`Ws0` |
A list of low rank matrices (W; |

`Gammas0` |
A list of matrices for weighting SNPs (Gamma; |

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

`eigen.G` |
A list with - $values
Eigen values - $vectors
Eigen vectors
The result of the eigen decompsition of |

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

`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'; |

`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, |

- $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`

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.

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

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.