coxmeg_gds | R Documentation |
coxmeg_gds
first estimates the variance component under a null model with only cov if tau is not given, and then analyzing each SNP in the gds file.
coxmeg_gds( gds, pheno, corr, type, cov = NULL, tau = NULL, snp.id = NULL, maf = 0.05, min_tau = 1e-04, max_tau = 5, eps = 1e-06, order = NULL, detap = NULL, opt = "bobyqa", score = FALSE, threshold = 0, solver = NULL, spd = TRUE, mc = 100, verbose = TRUE )
gds |
A GDS object created by |
pheno |
A data.frame containing the following four columns, family ID, individual ID, time and status. The status is a binary variable (1 for events/0 for censored). Missing values should be denoted as NA. |
corr |
A relatedness matrix or a List object of matrices if there are multiple relatedness matrices. They can be a matrix or a 'dgCMatrix' class in the Matrix package. The matrix (or the sum if there are multiple) must be symmetric positive definite or symmetric positive semidefinite. The order of subjects must be consistent with that in pheno. |
type |
A string indicating the sparsity structure of the relatedness matrix. Should be 'bd' (block diagonal), 'sparse', or 'dense'. See details. |
cov |
An optional data.frame of covariates. Same as |
tau |
An optional positive value or vector for the variance component(s). If tau is given, the function will skip estimating the variance component, and use the given tau to analyze the SNPs. |
snp.id |
An optional vector of snp.id (or variant.id). Only these SNPs will be analyzed. Default is |
maf |
An optional positive value. All SNPs with MAF<maf in the GDS file will not be analyzed. If |
min_tau |
An optional positive value indicating the lower bound in the optimization algorithm for the variance component tau. Default is 1e-4. |
max_tau |
An optional positive value indicating the upper bound in the optimization algorithm for the variance component tau. Default is 5. |
eps |
An optional positive value indicating the relative convergence tolerance in the optimization algorithm. Default is 1e-6. A smaller value (e.g., 1e-8) can be used for better precision of the p-values in the situation where most SNPs under investigation have a very low minor allele count (<5). |
order |
An optional integer value starting from 0. Only effective when |
detap |
An optional string indicating whether to use an approximation for log-determinant. Can be 'exact', 'diagonal', 'gkb', or 'slq'. Default is NULL, which lets the function select a method based on 'type' and other information. See details. |
opt |
An optional logical scalar for the Optimization algorithm for estimating the variance component(s). Can be one of the following values: 'bobyqa', 'Brent', 'NM', or 'L-BFGS-B' (only for >1 variance components). Default is 'bobyqa'. |
score |
An optional logical value indicating whether to perform a score test. Default is FALSE. |
threshold |
An optional non-negative value. If threshold>0, coxmeg_gds will reestimate HRs for those SNPs with a p-value<threshold by first estimating a variant-specific variance component. Default is 0. |
solver |
An optional binary value taking either 1 or 2. Default is 1. See details. |
spd |
An optional logical value indicating whether the relatedness matrix is symmetric positive definite. Default is TRUE. See details. |
mc |
An optional integer scalar specifying the number of Monte Carlo samples used for approximating the log-determinant when |
verbose |
An optional logical value indicating whether to print additional messages. Default is TRUE. |
beta: The estimated coefficient for each predictor in X.
HR: The estimated HR for each predictor in X.
sd_beta: The estimated standard error of beta.
p: The p-value of each SNP.
tau: The estimated variance component.
rank: The rank of the relatedness matrix.
nsam: Actual sample size.
corr
The subjects in corr
must be in the same order as in pheno
and cov
.
type
Specifying the type of the relatedness matrix (whether it is block-diagonal, general sparse, or dense). In the case of multiple relatedness matrices, it refers to the type of the sum of these matrices.
"bd" - used for a block-diagonal relatedness matrix, or a sparse matrix the inverse of which is also sparse.
"sparse" - used for a general sparse relatedness matrix the inverse of which is not sparse.
"dense" - used for a dense relatedness matrix.
spd
When spd=TRUE
, the relatedness matrix is treated as SPD. If the matrix is SPSD or not sure, use spd=FALSE
.
solver
Specifying which method is used to solve the linear system in the optimization algorithm.
"1" - Cholesky decompositon (RcppEigen:LDLT) is used to solve the linear system.
"2" - PCG is used to solve the linear system. When type='dense'
, it is recommended to set solver=2
to have better computational performance.
detap
Specifying the method to compute the log-determinant for estimating the variance component(s).
"exact" - the exact log-determinant is computed for estimating the variance component.
"diagonal" - uses diagonal approximation and is only effective for a sparse relatedness matrix.
"slq" - uses stochastic lanczos quadrature approximation. It uses the Lanczos algorithm to compute the weights and nodes. When type is 'bd' or 'sparse', it is often faster than 'gkb' and has the same accuracy. When type='dense', it is fater than 'gkb' by around half, but can be inaccurate if the relatedness matrix is (almost) singular.
"gkb" - uses stochastic lanczos quadrature approximation. It uses the Golub-Kahan bidiagonalization algorithm to compute the weights and nodes. It is robust against an (almost) singular relatedness matrix when type='dense', but it is generally slower than 'slq'.
Liang He, Stephanie Gogarten
library(Matrix) library(MASS) library(coxmeg) ## build a block-diagonal relatedness matrix n_f <- 600 mat_list <- list() size <- rep(5,n_f) offd <- 0.5 for(i in 1:n_f) { mat_list[[i]] <- matrix(offd,size[i],size[i]) diag(mat_list[[i]]) <- 1 } sigma <- as.matrix(bdiag(mat_list)) ## Example data files pheno.file = system.file("extdata", "ex_pheno.txt", package = "coxmeg") cov.file = system.file("extdata", "ex_cov.txt", package = "coxmeg") bed = system.file("extdata", "example_null.bed", package = "coxmeg") bed = substr(bed,1,nchar(bed)-4) ## Read phenotype and covariates pheno <- read.table(pheno.file, header=FALSE, as.is=TRUE, na.strings="-9") cov <- read.table(cov.file, header=FALSE, as.is=TRUE) ## SNPRelate GDS object gdsfile <- tempfile() SNPRelate::snpgdsBED2GDS(bed.fn=paste0(bed,".bed"), fam.fn=paste0(bed,".fam"), bim.fn=paste0(bed,".bim"), out.gdsfn=gdsfile, verbose=FALSE) gds <- SNPRelate::snpgdsOpen(gdsfile) ## Estimate variance component under a null model re <- coxmeg_gds(gds,pheno,sigma,type='bd',cov=cov,detap='diagonal',order=1) re SNPRelate::snpgdsClose(gds) unlink(gdsfile) ## SeqArray GDS object gdsfile <- tempfile() SeqArray::seqBED2GDS(bed.fn=paste0(bed,".bed"), fam.fn=paste0(bed,".fam"), bim.fn=paste0(bed,".bim"), out.gdsfn=gdsfile, verbose=FALSE) gds <- SeqArray::seqOpen(gdsfile) ## Estimate variance component under a null model re <- coxmeg_gds(gds,pheno,sigma,type='bd',cov=cov,detap='diagonal',order=1) re SeqArray::seqClose(gds) unlink(gdsfile)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.